From fd27d140fbd2a0c1365454ba18059d346eca575b Mon Sep 17 00:00:00 2001 From: Jesse Nusbaumer Date: Tue, 24 Sep 2024 10:45:28 -0600 Subject: [PATCH 01/14] Change Liceense to Apache 2.0 Replace original, locally-developed license with Apache 2.0 license. --- LICENSE | 201 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 201 insertions(+) create mode 100644 LICENSE diff --git a/LICENSE b/LICENSE new file mode 100644 index 00000000..261eeb9e --- /dev/null +++ b/LICENSE @@ -0,0 +1,201 @@ + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright [yyyy] [name of copyright owner] + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. From 592fe830a0bf248456196555312320014c5f3b5f Mon Sep 17 00:00:00 2001 From: Jesse Nusbaumer Date: Tue, 24 Sep 2024 10:49:21 -0600 Subject: [PATCH 02/14] Remove old license file. --- LICENSE.txt | 34 ---------------------------------- 1 file changed, 34 deletions(-) delete mode 100644 LICENSE.txt diff --git a/LICENSE.txt b/LICENSE.txt deleted file mode 100644 index 8c3b9387..00000000 --- a/LICENSE.txt +++ /dev/null @@ -1,34 +0,0 @@ -Copyright (c) 2017, University Corporation for Atmospheric Research (UCAR) -All rights reserved. - -Developed by: - University Corporation for Atmospheric Research - National Center for Atmospheric Research - https://www.cesm.ucar.edu/working_groups/Atmosphere - -Permission is hereby granted, free of charge, to any person obtaining -a copy of this software and associated documentation files (the "Software"), -to deal with the Software without restriction, including without limitation -the rights to use, copy, modify, merge, publish, distribute, sublicense, -and/or sell copies of the Software, and to permit persons to whom -the Software is furnished to do so, subject to the following conditions: - - - Redistributions of source code must retain the above copyright notice, - this list of conditions and the following disclaimers. - - Redistributions in binary form must reproduce the above copyright notice, - this list of conditions and the following disclaimers in the documentation - and/or other materials provided with the distribution. - - Neither the names of [Name of Development Group, UCAR], - nor the names of its contributors may be used to endorse or promote - products derived from this Software without specific prior written permission. - -THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" -AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE -LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -POSSIBILITY OF SUCH DAMAGE. From 03f2251fc3f16cc699bed45c04fca31e6918dd0b Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Sun, 29 Sep 2024 20:35:44 -0600 Subject: [PATCH 03/14] Initial round of needed files for vertical diffusion solve. --- to-be-ccpp-ized/utilities/coords_1d.F90 | 152 +++ .../utilities/linear_1d_operators.F90 | 1181 +++++++++++++++++ .../utilities/shr_infnan_mod.F90.in | 406 ++++++ to-be-ccpp-ized/utilities/shr_kind_mod.F90 | 21 + to-be-ccpp-ized/utilities/shr_log_mod.F90 | 121 ++ .../utilities/shr_strconvert_mod.F90 | 167 +++ to-be-ccpp-ized/utilities/shr_sys_mod.F90 | 332 +++++ .../vertical_diffusion/vertical_diffusion.F90 | 136 ++ 8 files changed, 2516 insertions(+) create mode 100644 to-be-ccpp-ized/utilities/coords_1d.F90 create mode 100644 to-be-ccpp-ized/utilities/linear_1d_operators.F90 create mode 100644 to-be-ccpp-ized/utilities/shr_infnan_mod.F90.in create mode 100644 to-be-ccpp-ized/utilities/shr_kind_mod.F90 create mode 100644 to-be-ccpp-ized/utilities/shr_log_mod.F90 create mode 100644 to-be-ccpp-ized/utilities/shr_strconvert_mod.F90 create mode 100644 to-be-ccpp-ized/utilities/shr_sys_mod.F90 create mode 100644 to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion.F90 diff --git a/to-be-ccpp-ized/utilities/coords_1d.F90 b/to-be-ccpp-ized/utilities/coords_1d.F90 new file mode 100644 index 00000000..21475de0 --- /dev/null +++ b/to-be-ccpp-ized/utilities/coords_1d.F90 @@ -0,0 +1,152 @@ +module coords_1d + + ! This module defines the Coords1D type, which is intended to to cache + ! commonly used information derived from a collection of sets of 1-D + ! coordinates. + + use shr_kind_mod, only: r8 => shr_kind_r8 + + implicit none + private + save + + public :: Coords1D + + type :: Coords1D + ! Number of sets of coordinates in the object. + integer :: n = 0 + ! Number of coordinates in each set. + integer :: d = 0 + + ! All fields below will be allocated with first dimension "n". + ! The second dimension is d+1 for ifc, d for mid, del, and rdel, and + ! d-1 for dst and rdst. + + ! Cell interface coordinates. + real(r8), allocatable :: ifc(:,:) + ! Coordinates at cell mid-points. + real(r8), allocatable :: mid(:,:) + ! Width of cells. + real(r8), allocatable :: del(:,:) + ! Distance between cell midpoints. + real(r8), allocatable :: dst(:,:) + ! Reciprocals: 1/del and 1/dst. + real(r8), allocatable :: rdel(:,:) + real(r8), allocatable :: rdst(:,:) + contains + procedure :: section + procedure :: finalize + end type Coords1D + + interface Coords1D + module procedure new_Coords1D_from_fields + module procedure new_Coords1D_from_int + end interface + + contains + + ! Constructor to create an object from existing data. + function new_Coords1D_from_fields(ifc, mid, del, dst, & + rdel, rdst) result(coords) + real(r8), USE_CONTIGUOUS intent(in) :: ifc(:,:) + real(r8), USE_CONTIGUOUS intent(in) :: mid(:,:) + real(r8), USE_CONTIGUOUS intent(in) :: del(:,:) + real(r8), USE_CONTIGUOUS intent(in) :: dst(:,:) + real(r8), USE_CONTIGUOUS intent(in) :: rdel(:,:) + real(r8), USE_CONTIGUOUS intent(in) :: rdst(:,:) + type(Coords1D) :: coords + + coords = allocate_coords(size(ifc, 1), size(ifc, 2) - 1) + + coords%ifc = ifc + coords%mid = mid + coords%del = del + coords%dst = dst + coords%rdel = rdel + coords%rdst = rdst + + end function new_Coords1D_from_fields + + ! Constructor if you only have interface coordinates; derives all the other + ! fields. + function new_Coords1D_from_int(ifc) result(coords) + real(r8), USE_CONTIGUOUS intent(in) :: ifc(:,:) + type(Coords1D) :: coords + + coords = allocate_coords(size(ifc, 1), size(ifc, 2) - 1) + + coords%ifc = ifc + coords%mid = 0.5_r8 * (ifc(:,:coords%d)+ifc(:,2:)) + coords%del = coords%ifc(:,2:) - coords%ifc(:,:coords%d) + coords%dst = coords%mid(:,2:) - coords%mid(:,:coords%d-1) + coords%rdel = 1._r8/coords%del + coords%rdst = 1._r8/coords%dst + + end function new_Coords1D_from_int + + ! Create a new Coords1D object that is a subsection of some other object, + ! e.g. if you want only the first m coordinates, use d_bnds=[1, m]. + ! + ! Originally this used pointers, but it was found to actually be cheaper + ! in practice just to make a copy, especially since pointers can impede + ! optimization. + function section(self, n_bnds, d_bnds) + class(Coords1D), intent(in) :: self + integer, intent(in) :: n_bnds(2), d_bnds(2) + type(Coords1D) :: section + + section = allocate_coords(n_bnds(2)-n_bnds(1)+1, d_bnds(2)-d_bnds(1)+1) + + section%ifc = self%ifc(n_bnds(1):n_bnds(2),d_bnds(1):d_bnds(2)+1) + section%mid = self%mid(n_bnds(1):n_bnds(2),d_bnds(1):d_bnds(2)) + section%del = self%del(n_bnds(1):n_bnds(2),d_bnds(1):d_bnds(2)) + section%dst = self%dst(n_bnds(1):n_bnds(2),d_bnds(1):d_bnds(2)-1) + section%rdel = self%rdel(n_bnds(1):n_bnds(2),d_bnds(1):d_bnds(2)) + section%rdst = self%rdst(n_bnds(1):n_bnds(2),d_bnds(1):d_bnds(2)-1) + + end function section + + ! Quick utility to get allocate each array with the correct size. + function allocate_coords(n, d) result(coords) + integer, intent(in) :: n, d + type(Coords1D) :: coords + + coords%n = n + coords%d = d + + allocate(coords%ifc(coords%n,coords%d+1)) + allocate(coords%mid(coords%n,coords%d)) + allocate(coords%del(coords%n,coords%d)) + allocate(coords%dst(coords%n,coords%d-1)) + allocate(coords%rdel(coords%n,coords%d)) + allocate(coords%rdst(coords%n,coords%d-1)) + + end function allocate_coords + + ! Deallocate and reset to initial state. + subroutine finalize(self) + class(Coords1D), intent(inout) :: self + + self%n = 0 + self%d = 0 + + call guarded_deallocate(self%ifc) + call guarded_deallocate(self%mid) + call guarded_deallocate(self%del) + call guarded_deallocate(self%dst) + call guarded_deallocate(self%rdel) + call guarded_deallocate(self%rdst) + + contains + + subroutine guarded_deallocate(array) + real(r8), allocatable :: array(:,:) + + if (allocated(array)) deallocate(array) + + end subroutine guarded_deallocate + + end subroutine finalize + +end module coords_1d + \ No newline at end of file diff --git a/to-be-ccpp-ized/utilities/linear_1d_operators.F90 b/to-be-ccpp-ized/utilities/linear_1d_operators.F90 new file mode 100644 index 00000000..901200f4 --- /dev/null +++ b/to-be-ccpp-ized/utilities/linear_1d_operators.F90 @@ -0,0 +1,1181 @@ +module linear_1d_operators + + ! This module provides the type "TriDiagOp" to represent operators on a 1D + ! grid as tridiagonal matrices, and related types to represent boundary + ! conditions. + ! + ! The focus is on solving diffusion equations with a finite volume method + ! in one dimension, but other utility operators are provided, e.g. a second + ! order approximation to the first derivative. + ! + ! In order to allow vectorization to occur, as well as to avoid unnecessary + ! copying/reshaping of data in CAM, TriDiagOp actually represents a + ! collection of independent operators that can be applied to collections of + ! independent data; the innermost index is over independent systems (e.g. + ! CAM columns). + ! + ! A simple example: + ! ! First derivative operator + ! op = first_derivative(coords) + ! ! Convert data to its derivative (extrapolate at boundaries). + ! call op%apply(data) + ! + ! With explicit boundary conditions: + ! op = first_derivative(coords, & + ! l_bndry=BoundaryFixedFlux(), & + ! r_bndry=BoundaryFixedLayer(layer_distance)) + ! call op%apply(data, & + ! l_cond=BoundaryFlux(flux, dt, thickness), & + ! r_cond=BoundaryData(boundary)) + ! + ! Implicit solution example: + ! ! Construct diffusion matrix. + ! op = diffusion_operator(coords, d) + ! call op%lmult_as_diag(-dt) + ! call op%add_to_diag(1._r8) + ! ! Decompose in order to invert the operation. + ! decomp = TriDiagDecomp(op) + ! ! Diffuse data for one time step (fixed flux boundaries). + ! call decomp%left_div(data) + + use shr_kind_mod, only: r8 => shr_kind_r8 + use shr_log_mod, only: errMsg => shr_log_errMsg + use shr_sys_mod, only: shr_sys_abort + use coords_1d, only: Coords1D + + implicit none + private + save + + ! Main type. + public :: TriDiagOp + public :: operator(+) + public :: operator(-) + + ! Decomposition used for inversion (left division). + public :: TriDiagDecomp + + ! Multiplies by 0. + public :: zero_operator + + ! Construct identity. + public :: identity_operator + + ! Produce a TriDiagOp that is simply a diagonal matrix. + public :: diagonal_operator + + ! For solving the diffusion-advection equation with implicit Euler. + public :: diffusion_operator + public :: advection_operator + + ! Derivatives accurate to second order on a non-uniform grid. + public :: first_derivative + public :: second_derivative + + ! Boundary condition types. + public :: BoundaryType + public :: BoundaryZero + public :: BoundaryFirstOrder + public :: BoundaryExtrapolate + public :: BoundaryFixedLayer + public :: BoundaryFixedFlux + + ! Boundary data types. + public :: BoundaryCond + public :: BoundaryNoData + public :: BoundaryData + public :: BoundaryFlux + + ! TriDiagOp represents operators that can work between nearest neighbors, + ! with some extra logic at the boundaries. The implementation is a + ! tridiagonal matrix plus boundary info. + type :: TriDiagOp + private + ! The number of independent systems. + integer, public :: nsys + ! The size of the matrix (number of grid cells). + integer, public :: ncel + ! Super-, sub-, and regular diagonals. + real(r8), allocatable :: spr(:,:) + real(r8), allocatable :: sub(:,:) + real(r8), allocatable :: diag(:,:) + ! Buffers to hold boundary data; Details depend on the type of boundary + ! being used. + real(r8), allocatable :: left_bound(:) + real(r8), allocatable :: right_bound(:) + contains + ! Applies the operator to a set of data. + procedure :: apply => apply_tridiag + ! Given the off-diagonal elements, fills in the diagonal so that the + ! operator will have the constant function as an eigenvector with + ! eigenvalue 0. This is used internally as a utility for construction of + ! derivative operators. + procedure :: deriv_diag => make_tridiag_deriv_diag + ! Add/substract another tridiagonal from this one in-place (without + ! creating a temporary object). + procedure :: add => add_in_place_tridiag_ops + procedure :: subtract => subtract_in_place_tridiag_ops + ! Add input vector or scalar to the diagonal. + procedure :: scalar_add_tridiag + procedure :: diagonal_add_tridiag + generic :: add_to_diag => scalar_add_tridiag, diagonal_add_tridiag + ! Treat input vector (or scalar) as if it was the diagonal of an + ! operator, and multiply this operator on the left by that value. + procedure :: scalar_lmult_tridiag + procedure :: diagonal_lmult_tridiag + generic :: lmult_as_diag => & + scalar_lmult_tridiag, diagonal_lmult_tridiag + ! Deallocate and reset. + procedure :: finalize => tridiag_finalize + end type TriDiagOp + + interface operator(+) + module procedure add_tridiag_ops + end interface operator(+) + + interface operator(-) + module procedure subtract_tridiag_ops + end interface operator(-) + + interface TriDiagOp + module procedure new_TriDiagOp + end interface TriDiagOp + + ! + ! Boundary condition types for the operators. + ! + ! Note that BoundaryFixedLayer and BoundaryFixedFlux are the only options + ! supported for backwards operation (i.e. decomp%left_div). The others are + ! meant for direct application only (e.g. to find a derivative). + ! + ! BoundaryZero means that the operator fixes boundaries to 0. + ! BoundaryFirstOrder means a one-sided approximation for the first + ! derivative. + ! BoundaryExtrapolate means that a second order approximation will be used, + ! even at the boundaries. Boundary points do this by using their next- + ! nearest neighbor to extrapolate. + ! BoundaryFixedLayer means that there's an extra layer outside of the given + ! grid, which must be specified when applying/inverting the operator. + ! BoundaryFixedFlux is intended to provide a fixed-flux condition for + ! typical advection/diffusion operators. It tweaks the edge condition + ! to work on an input current rather than a value. + ! + ! The different types were originally implemented through polymorphism, but + ! PGI required this to be done via enum instead. + integer, parameter :: zero_bndry = 0 + integer, parameter :: first_order_bndry = 1 + integer, parameter :: extrapolate_bndry = 2 + integer, parameter :: fixed_layer_bndry = 3 + integer, parameter :: fixed_flux_bndry = 4 + + type :: BoundaryType + private + integer :: bndry_type = fixed_flux_bndry + real(r8), allocatable :: edge_width(:) + contains + procedure :: make_left + procedure :: make_right + procedure :: finalize => boundary_type_finalize + end type BoundaryType + + abstract interface + subroutine deriv_seed(del_minus, del_plus, sub, spr) + import :: r8 + real(r8), USE_CONTIGUOUS intent(in) :: del_minus(:) + real(r8), USE_CONTIGUOUS intent(in) :: del_plus(:) + real(r8), USE_CONTIGUOUS intent(out) :: sub(:) + real(r8), USE_CONTIGUOUS intent(out) :: spr(:) + end subroutine deriv_seed + end interface + + interface BoundaryZero + module procedure new_BoundaryZero + end interface BoundaryZero + + interface BoundaryFirstOrder + module procedure new_BoundaryFirstOrder + end interface BoundaryFirstOrder + + interface BoundaryExtrapolate + module procedure new_BoundaryExtrapolate + end interface BoundaryExtrapolate + + interface BoundaryFixedLayer + module procedure new_BoundaryFixedLayer + end interface BoundaryFixedLayer + + interface BoundaryFixedFlux + module procedure new_BoundaryFixedFlux + end interface BoundaryFixedFlux + + ! + ! Data for boundary conditions themselves. + ! + ! "No data" conditions perform extrapolation, if BoundaryExtrapolate was + ! the boundary type used to construct the operator. + ! + ! "Data" conditions contain extra data, which effectively extends the + ! system with an extra cell. + ! + ! "Flux" conditions contain prescribed fluxes. + ! + ! The condition you can use depends on the boundary type from above that + ! was used in the operator's construction. For BoundaryFixedLayer use + ! BoundaryData. For BoundaryFixedFlux use BoundaryFlux. For everything + ! else, use BoundaryNoData. + + ! The switches using this enumeration used to be unnecessary due to use of + ! polymorphism, but this had to be backed off due to insufficient PGI + ! support for type extension. + integer, parameter :: no_data_cond = 0 + integer, parameter :: data_cond = 1 + integer, parameter :: flux_cond = 2 + + type :: BoundaryCond + private + integer :: cond_type = no_data_cond + real(r8), allocatable :: edge_data(:) + contains + procedure :: apply_left + procedure :: apply_right + procedure :: finalize => boundary_cond_finalize + end type BoundaryCond + + ! Constructors for different types of BoundaryCond. + interface BoundaryNoData + module procedure new_BoundaryNoData + end interface BoundaryNoData + + interface BoundaryData + module procedure new_BoundaryData + end interface BoundaryData + + interface BoundaryFlux + module procedure new_BoundaryFlux + end interface BoundaryFlux + + ! Opaque type to hold a tridiagonal matrix decomposition. + ! + ! Method used is similar to Richtmyer and Morton (1967,pp 198-201), but + ! the order of iteration is reversed, leading to A and C being swapped, and + ! some differences in the indexing. + type :: TriDiagDecomp + private + integer :: nsys = 0 + integer :: ncel = 0 + ! These correspond to A_k, E_k, and 1 / (B_k - A_k * E_{k+1}) + real(r8), allocatable :: ca(:,:) + real(r8), allocatable :: ze(:,:) + real(r8), allocatable :: dnom(:,:) + contains + procedure :: left_div => decomp_left_div + procedure :: finalize => decomp_finalize + end type TriDiagDecomp + + interface TriDiagDecomp + module procedure new_TriDiagDecomp + end interface TriDiagDecomp + + contains + + ! Operator that sets to 0. + function zero_operator(nsys, ncel) result(op) + ! Sizes for operator. + integer, intent(in) :: nsys, ncel + + type(TriDiagOp) :: op + + op = TriDiagOp(nsys, ncel) + + op%spr = 0._r8 + op%sub = 0._r8 + op%diag = 0._r8 + op%left_bound = 0._r8 + op%right_bound = 0._r8 + + end function zero_operator + + ! Operator that does nothing. + function identity_operator(nsys, ncel) result(op) + ! Sizes for operator. + integer, intent(in) :: nsys, ncel + + type(TriDiagOp) :: op + + op = TriDiagOp(nsys, ncel) + + op%spr = 0._r8 + op%sub = 0._r8 + op%diag = 1._r8 + op%left_bound = 0._r8 + op%right_bound = 0._r8 + + end function identity_operator + + ! Create an operator that just does an element-wise product by some data. + function diagonal_operator(diag) result(op) + ! Data to multiply by. + real(r8), USE_CONTIGUOUS intent(in) :: diag(:,:) + + type(TriDiagOp) :: op + + op = TriDiagOp(size(diag, 1), size(diag, 2)) + + op%spr = 0._r8 + op%sub = 0._r8 + op%diag = diag + op%left_bound = 0._r8 + op%right_bound = 0._r8 + + end function diagonal_operator + + ! Diffusion matrix operator constructor. Given grid coordinates, a set of + ! diffusion coefficients, and boundaries, creates a matrix corresponding + ! to a finite volume representation of the operation: + ! + ! d/dx (d_coef * d/dx) + ! + ! This differs from what you would get from combining the first and second + ! derivative operations, which would be more appropriate for a finite + ! difference scheme that does not use grid cell averages. + function diffusion_operator(coords, d_coef, l_bndry, r_bndry) & + result(op) + ! Grid cell locations. + type(Coords1D), intent(in) :: coords + ! Diffusion coefficient defined on interfaces. + real(r8), USE_CONTIGUOUS intent(in) :: d_coef(:,:) + ! Objects representing the kind of boundary on each side. + class(BoundaryType), target, intent(in), optional :: l_bndry, r_bndry + ! Output operator. + type(TriDiagOp) :: op + + ! Selectors to implement default boundary. + class(BoundaryType), pointer :: l_bndry_loc, r_bndry_loc + ! Fixed flux is default, no allocation/deallocation needed. + type(BoundaryType), target :: bndry_default + + ! Level index. + integer :: k + + if (present(l_bndry)) then + l_bndry_loc => l_bndry + else + l_bndry_loc => bndry_default + end if + + if (present(r_bndry)) then + r_bndry_loc => r_bndry + else + r_bndry_loc => bndry_default + end if + + ! Allocate the operator. + op = TriDiagOp(coords%n, coords%d) + + ! d_coef over the distance to the next cell gives you the matrix term for + ! flux of material between cells. Dividing by cell thickness translates + ! this to a tendency on the concentration. Hence the basic pattern is + ! d_coef*rdst*rdel. + ! + ! Boundary conditions for a fixed layer simply extend this by calculating + ! the distance to the midpoint of the extra edge layer. + + select case (l_bndry_loc%bndry_type) + case (fixed_layer_bndry) + op%left_bound = 2._r8*d_coef(:,1)*coords%rdel(:,1) / & + (l_bndry_loc%edge_width+coords%del(:,1)) + case default + op%left_bound = 0._r8 + end select + + do k = 1, coords%d-1 + op%spr(:,k) = d_coef(:,k+1)*coords%rdst(:,k)*coords%rdel(:,k) + op%sub(:,k) = d_coef(:,k+1)*coords%rdst(:,k)*coords%rdel(:,k+1) + end do + + select case (r_bndry_loc%bndry_type) + case (fixed_layer_bndry) + op%right_bound = 2._r8*d_coef(:,coords%d+1)*coords%rdel(:,coords%d) / & + (r_bndry_loc%edge_width+coords%del(:,coords%d)) + case default + op%right_bound = 0._r8 + end select + + ! Above, we found all off-diagonals. Now get the diagonal. + call op%deriv_diag() + + end function diffusion_operator + + ! Advection matrix operator constructor. Similar to diffusion_operator, it + ! constructs an operator A corresponding to: + ! + ! A y = d/dx (-v_coef * y) + ! + ! Again, this is targeted at representing this operator acting on grid-cell + ! averages in a finite volume scheme, rather than a literal representation. + function advection_operator(coords, v_coef, l_bndry, r_bndry) & + result(op) + ! Grid cell locations. + type(Coords1D), intent(in) :: coords + ! Advection coefficient (effective velocity). + real(r8), USE_CONTIGUOUS intent(in) :: v_coef(:,:) + ! Objects representing the kind of boundary on each side. + class(BoundaryType), target, intent(in), optional :: l_bndry, r_bndry + ! Output operator. + type(TriDiagOp) :: op + + ! Selectors to implement default boundary. + class(BoundaryType), pointer :: l_bndry_loc, r_bndry_loc + ! Fixed flux is default, no allocation/deallocation needed. + type(BoundaryType), target :: bndry_default + + ! Negative derivative of v. + real(r8) :: v_deriv(coords%n,coords%d) + + if (present(l_bndry)) then + l_bndry_loc => l_bndry + else + l_bndry_loc => bndry_default + end if + + if (present(r_bndry)) then + r_bndry_loc => r_bndry + else + r_bndry_loc => bndry_default + end if + + ! Allocate the operator. + op = TriDiagOp(coords%n, coords%d) + + ! Construct the operator in two stages using the product rule. First + ! create (-v * d/dx), then -dv/dx, and add the two. + ! + ! For the first part, we want to interpolate to interfaces (weighted + ! average involving del/2*dst), multiply by -v to get flux, then divide + ! by cell thickness, which gives a concentration tendency: + ! + ! (del/(2*dst))*(-v_coef)/del + ! + ! Simplifying gives -v_coef*rdst*0.5, as seen below. + + select case (l_bndry_loc%bndry_type) + case (fixed_layer_bndry) + op%left_bound = v_coef(:,1) / & + (l_bndry_loc%edge_width+coords%del(:,1)) + case default + op%left_bound = 0._r8 + end select + + op%sub = v_coef(:,2:coords%d)*coords%rdst*0.5_r8 + op%spr = -op%sub + + select case (r_bndry_loc%bndry_type) + case (fixed_layer_bndry) + op%right_bound = v_coef(:,coords%d+1) / & + (r_bndry_loc%edge_width+coords%del(:,coords%d)) + case default + op%right_bound = 0._r8 + end select + + ! Above, we found all off-diagonals. Now get the diagonal. This must be + ! done at this specific point, since the other half of the operator is + ! not "derivative-like" in the sense of yielding 0 for a constant input. + call op%deriv_diag() + + ! The second half of the operator simply involves taking a first-order + ! derivative of v. Since v is on the interfaces, just use: + ! (v(k+1) - v(k))*rdel(k) + v_deriv(:,1) = v_coef(:,2)*coords%rdel(:,1) + + select case (l_bndry_loc%bndry_type) + case (fixed_layer_bndry) + v_deriv(:,1) = v_deriv(:,1) - v_coef(:,1)*coords%rdel(:,1) + end select + + v_deriv(:,2:coords%d-1) = (v_coef(:,3:coords%d) - & + v_coef(:,2:coords%d-1))*coords%rdel(:,2:coords%d-1) + + v_deriv(:,coords%d) = -v_coef(:,coords%d)*coords%rdel(:,coords%d) + + select case (r_bndry_loc%bndry_type) + case (fixed_layer_bndry) + v_deriv(:,coords%d) = v_deriv(:,coords%d) & + + v_coef(:,coords%d+1)*coords%del(:,coords%d) + end select + + ! Combine the two pieces. + op%diag = op%diag - v_deriv + + end function advection_operator + + ! Second order approximation to the first and second derivatives on a non- + ! uniform grid. + ! + ! Both operators are constructed with the same method, except for a "seed" + ! function that takes local distances between points to create the + ! off-diagonal terms. + function first_derivative(grid_spacing, l_bndry, r_bndry) result(op) + ! Distances between points. + real(r8), USE_CONTIGUOUS intent(in) :: grid_spacing(:,:) + ! Boundary conditions. + class(BoundaryType), intent(in), optional :: l_bndry, r_bndry + ! Output operator. + type(TriDiagOp) :: op + + op = deriv_op_from_seed(grid_spacing, first_derivative_seed, & + l_bndry, r_bndry) + + end function first_derivative + + subroutine first_derivative_seed(del_minus, del_plus, sub, spr) + ! Distances to next and previous point. + real(r8), USE_CONTIGUOUS intent(in) :: del_minus(:) + real(r8), USE_CONTIGUOUS intent(in) :: del_plus(:) + ! Off-diagonal matrix terms. + real(r8), USE_CONTIGUOUS intent(out) :: sub(:) + real(r8), USE_CONTIGUOUS intent(out) :: spr(:) + + real(r8) :: del_sum(size(del_plus)) + + del_sum = del_plus + del_minus + + sub = - del_plus / (del_minus*del_sum) + spr = del_minus / (del_plus*del_sum) + + end subroutine first_derivative_seed + + function second_derivative(grid_spacing, l_bndry, r_bndry) result(op) + ! Distances between points. + real(r8), USE_CONTIGUOUS intent(in) :: grid_spacing(:,:) + ! Boundary conditions. + class(BoundaryType), intent(in), optional :: l_bndry, r_bndry + ! Output operator. + type(TriDiagOp) :: op + + op = deriv_op_from_seed(grid_spacing, second_derivative_seed, & + l_bndry, r_bndry) + + end function second_derivative + + subroutine second_derivative_seed(del_minus, del_plus, sub, spr) + ! Distances to next and previous point. + real(r8), USE_CONTIGUOUS intent(in) :: del_minus(:) + real(r8), USE_CONTIGUOUS intent(in) :: del_plus(:) + ! Off-diagonal matrix terms. + real(r8), USE_CONTIGUOUS intent(out) :: sub(:) + real(r8), USE_CONTIGUOUS intent(out) :: spr(:) + + real(r8) :: del_sum(size(del_plus)) + + del_sum = del_plus + del_minus + + sub = 2._r8 / (del_minus*del_sum) + spr = 2._r8 / (del_plus*del_sum) + + end subroutine second_derivative_seed + + ! Brains behind the first/second derivative functions. + function deriv_op_from_seed(grid_spacing, seed, l_bndry, r_bndry) result(op) + ! Distances between points. + real(r8), USE_CONTIGUOUS intent(in) :: grid_spacing(:,:) + ! Function to locally construct matrix elements. + procedure(deriv_seed) :: seed + ! Boundary conditions. + class(BoundaryType), target, intent(in), optional :: l_bndry, r_bndry + ! Output operator. + type(TriDiagOp) :: op + + ! Selectors to implement default boundary. + class(BoundaryType), pointer :: l_bndry_loc, r_bndry_loc + ! Fixed flux is default, no allocation/deallocation needed. + type(BoundaryType), target :: bndry_default + + integer :: k + + if (present(l_bndry)) then + l_bndry_loc => l_bndry + else + l_bndry_loc => bndry_default + end if + + if (present(r_bndry)) then + r_bndry_loc => r_bndry + else + r_bndry_loc => bndry_default + end if + + ! Number of grid points is one greater than the spacing. + op = TriDiagOp(size(grid_spacing, 1), size(grid_spacing, 2) + 1) + + ! Left boundary condition. + call l_bndry_loc%make_left(grid_spacing, seed, & + op%left_bound, op%spr(:,1)) + + do k = 2, op%ncel-1 + call seed(grid_spacing(:,k-1), grid_spacing(:,k), & + op%sub(:,k-1), op%spr(:,k)) + end do + + ! Right boundary condition. + call r_bndry_loc%make_right(grid_spacing, seed, & + op%sub(:,op%ncel-1), op%right_bound) + + ! Above, we found all off-diagonals. Now get the diagonal. + call op%deriv_diag() + + end function deriv_op_from_seed + + ! Boundary constructors. Most simply set an internal flag, but + ! BoundaryFixedLayer accepts an argument representing the distance to the + ! location where the extra layer is defined. + + function new_BoundaryZero() result(new_bndry) + type(BoundaryType) :: new_bndry + + new_bndry%bndry_type = zero_bndry + + end function new_BoundaryZero + + function new_BoundaryFirstOrder() result(new_bndry) + type(BoundaryType) :: new_bndry + + new_bndry%bndry_type = first_order_bndry + + end function new_BoundaryFirstOrder + + function new_BoundaryExtrapolate() result(new_bndry) + type(BoundaryType) :: new_bndry + + new_bndry%bndry_type = extrapolate_bndry + + end function new_BoundaryExtrapolate + + function new_BoundaryFixedLayer(width) result(new_bndry) + real(r8), USE_CONTIGUOUS intent(in) :: width(:) + type(BoundaryType) :: new_bndry + + new_bndry%bndry_type = fixed_layer_bndry + new_bndry%edge_width = width + + end function new_BoundaryFixedLayer + + function new_BoundaryFixedFlux() result(new_bndry) + type(BoundaryType) :: new_bndry + + new_bndry%bndry_type = fixed_flux_bndry + + end function new_BoundaryFixedFlux + + ! The make_left and make_right methods implement the boundary conditions + ! using an input seed. + + subroutine make_left(self, grid_spacing, seed, term1, term2) + class(BoundaryType), intent(in) :: self + real(r8), USE_CONTIGUOUS intent(in) :: grid_spacing(:,:) + procedure(deriv_seed) :: seed + real(r8), USE_CONTIGUOUS intent(out) :: term1(:) + real(r8), USE_CONTIGUOUS intent(out) :: term2(:) + + real(r8) :: del_plus(size(term1)), del_minus(size(term1)) + + select case (self%bndry_type) + case (zero_bndry) + term1 = 0._r8 + term2 = 0._r8 + case (first_order_bndry) + ! To calculate to first order, just use a really huge del_minus (i.e. + ! pretend that there's a point so far away it doesn't matter). + del_plus = grid_spacing(:,1) + del_minus = del_plus * 4._r8 / epsilon(1._r8) + call seed(del_minus, del_plus, term1, term2) + case (extrapolate_bndry) + ! To extrapolate from the boundary, use distance from the nearest + ! neighbor (as usual) and the second nearest neighbor (with a negative + ! sign, since we are using two points on the same side). + del_plus = grid_spacing(:,1) + del_minus = - (grid_spacing(:,1) + grid_spacing(:,2)) + call seed(del_minus, del_plus, term1, term2) + case (fixed_layer_bndry) + ! Use edge value to extend the grid. + del_plus = grid_spacing(:,1) + del_minus = self%edge_width + call seed(del_minus, del_plus, term1, term2) + case (fixed_flux_bndry) + ! Treat grid as uniform, but then zero out the contribution from data + ! on one side (since it will be prescribed). + del_plus = grid_spacing(:,1) + del_minus = del_plus + call seed(del_minus, del_plus, term1, term2) + term1 = 0._r8 + case default + call shr_sys_abort("Invalid boundary type at "// & + errMsg(__FILE__, __LINE__)) + end select + + end subroutine make_left + + subroutine make_right(self, grid_spacing, seed, term1, term2) + class(BoundaryType), intent(in) :: self + real(r8), USE_CONTIGUOUS intent(in) :: grid_spacing(:,:) + procedure(deriv_seed) :: seed + real(r8), USE_CONTIGUOUS intent(out) :: term1(:) + real(r8), USE_CONTIGUOUS intent(out) :: term2(:) + + real(r8) :: del_plus(size(term1)), del_minus(size(term1)) + + select case (self%bndry_type) + case (zero_bndry) + term1 = 0._r8 + term2 = 0._r8 + case (first_order_bndry) + ! Use huge del_plus, analogous to how left boundary works. + del_minus = grid_spacing(:,size(grid_spacing, 2)) + del_plus = del_minus * 4._r8 / epsilon(1._r8) + call seed(del_minus, del_plus, term1, term2) + case (extrapolate_bndry) + ! Same strategy as left boundary, but reversed. + del_plus = - (grid_spacing(:,size(grid_spacing, 2) - 1) + & + grid_spacing(:,size(grid_spacing, 2))) + del_minus = grid_spacing(:,size(grid_spacing, 2)) + call seed(del_minus, del_plus, term1, term2) + case (fixed_layer_bndry) + ! Use edge value to extend the grid. + del_plus = self%edge_width + del_minus = grid_spacing(:,size(grid_spacing, 2)) + call seed(del_minus, del_plus, term1, term2) + case (fixed_flux_bndry) + ! Uniform grid, but with edge zeroed. + del_plus = grid_spacing(:,size(grid_spacing, 2)) + del_minus = del_plus + call seed(del_minus, del_plus, term1, term2) + term2 = 0._r8 + case default + call shr_sys_abort("Invalid boundary type at "// & + errMsg(__FILE__, __LINE__)) + end select + + end subroutine make_right + + subroutine boundary_type_finalize(self) + class(BoundaryType), intent(inout) :: self + + self%bndry_type = fixed_flux_bndry + if (allocated(self%edge_width)) deallocate(self%edge_width) + + end subroutine boundary_type_finalize + + ! Constructor for TriDiagOp; this just sets the size and allocates + ! arrays. + type(TriDiagOp) function new_TriDiagOp(nsys, ncel) + + integer, intent(in) :: nsys, ncel + + new_TriDiagOp%nsys = nsys + new_TriDiagOp%ncel = ncel + + allocate(new_TriDiagOp%spr(nsys,ncel-1), & + new_TriDiagOp%sub(nsys,ncel-1), & + new_TriDiagOp%diag(nsys,ncel), & + new_TriDiagOp%left_bound(nsys), & + new_TriDiagOp%right_bound(nsys)) + + end function new_TriDiagOp + + ! Deallocator for TriDiagOp. + subroutine tridiag_finalize(self) + class(TriDiagOp), intent(inout) :: self + + self%nsys = 0 + self%ncel = 0 + + if (allocated(self%spr)) deallocate(self%spr) + if (allocated(self%sub)) deallocate(self%sub) + if (allocated(self%diag)) deallocate(self%diag) + if (allocated(self%left_bound)) deallocate(self%left_bound) + if (allocated(self%right_bound)) deallocate(self%right_bound) + + end subroutine tridiag_finalize + + ! Boundary condition constructors. + + function new_BoundaryNoData() result(new_cond) + type(BoundaryCond) :: new_cond + + new_cond%cond_type = no_data_cond + ! No edge data, so leave it unallocated. + + end function new_BoundaryNoData + + function new_BoundaryData(data) result(new_cond) + real(r8), USE_CONTIGUOUS intent(in) :: data(:) + type(BoundaryCond) :: new_cond + + new_cond%cond_type = data_cond + new_cond%edge_data = data + + end function new_BoundaryData + + function new_BoundaryFlux(flux, dt, spacing) result(new_cond) + real(r8), USE_CONTIGUOUS intent(in) :: flux(:) + real(r8), intent(in) :: dt + real(r8), USE_CONTIGUOUS intent(in) :: spacing(:) + type(BoundaryCond) :: new_cond + + new_cond%cond_type = flux_cond + new_cond%edge_data = flux*dt/spacing + + end function new_BoundaryFlux + + ! Application of input data. + ! + ! When no data is input, assume that any bound term is applied to the + ! third element in from the edge for extrapolation. Boundary conditions + ! that don't need any edge data at all can then simply set the boundary + ! terms to 0. + + function apply_left(self, bound_term, array) result(delta_edge) + class(BoundaryCond), intent(in) :: self + real(r8), USE_CONTIGUOUS intent(in) :: bound_term(:) + real(r8), USE_CONTIGUOUS intent(in) :: array(:,:) + real(r8) :: delta_edge(size(array, 1)) + + select case (self%cond_type) + case (no_data_cond) + delta_edge = bound_term*array(:,3) + case (data_cond) + delta_edge = bound_term*self%edge_data + case (flux_cond) + delta_edge = self%edge_data + case default + call shr_sys_abort("Invalid boundary condition at "// & + errMsg(__FILE__, __LINE__)) + end select + + end function apply_left + + function apply_right(self, bound_term, array) result(delta_edge) + class(BoundaryCond), intent(in) :: self + real(r8), USE_CONTIGUOUS intent(in) :: bound_term(:) + real(r8), USE_CONTIGUOUS intent(in) :: array(:,:) + real(r8) :: delta_edge(size(array, 1)) + + select case (self%cond_type) + case (no_data_cond) + delta_edge = bound_term*array(:,size(array, 2)-2) + case (data_cond) + delta_edge = bound_term*self%edge_data + case (flux_cond) + delta_edge = self%edge_data + case default + call shr_sys_abort("Invalid boundary condition at "// & + errMsg(__FILE__, __LINE__)) + end select + + end function apply_right + + subroutine boundary_cond_finalize(self) + class(BoundaryCond), intent(inout) :: self + + self%cond_type = no_data_cond + if (allocated(self%edge_data)) deallocate(self%edge_data) + + end subroutine boundary_cond_finalize + + ! Apply an operator and return the new data. + function apply_tridiag(self, array, l_cond, r_cond) result(output) + ! Operator to apply. + class(TriDiagOp), intent(in) :: self + ! Data to act on. + real(r8), USE_CONTIGUOUS intent(in) :: array(:,:) + ! Objects representing boundary conditions. + class(BoundaryCond), target, intent(in), optional :: l_cond, r_cond + ! Function result. + real(r8) :: output(size(array, 1), size(array, 2)) + + ! Local objects to implement default. + class(BoundaryCond), pointer :: l_cond_loc, r_cond_loc + ! Default state is no data, no allocation/deallocation needed. + type(BoundaryCond), target :: cond_default + + ! Level index. + integer :: k + + if (present(l_cond)) then + l_cond_loc => l_cond + else + l_cond_loc => cond_default + end if + + if (present(r_cond)) then + r_cond_loc => r_cond + else + r_cond_loc => cond_default + end if + + ! Left boundary. + output(:,1) = self%diag(:,1)*array(:,1) + & + self%spr(:,1)*array(:,2) + & + l_cond_loc%apply_left(self%left_bound, array) + + do k = 2, self%ncel-1 + output(:,k) = & + self%sub(:,k-1)*array(:,k-1) + & + self%diag(:,k)*array(:,k ) + & + self%spr(:,k)*array(:,k+1) + end do + + ! Right boundary. + output(:,self%ncel) = & + self%sub(:,self%ncel-1)*array(:,self%ncel-1) + & + self%diag(:,self%ncel)*array(:,self%ncel) + & + r_cond_loc%apply_right(self%right_bound, array) + + end function apply_tridiag + + ! Fill in the diagonal for a TriDiagOp for a derivative operator, where + ! the off diagonal elements are already filled in. + subroutine make_tridiag_deriv_diag(self) + + class(TriDiagOp), intent(inout) :: self + + ! If a derivative operator operates on a constant function, it must + ! return 0 everywhere. To force this, make sure that all rows add to + ! zero in the matrix. + self%diag(:,:self%ncel-1) = - self%spr + self%diag(:,self%ncel) = - self%right_bound + self%diag(:,1) = self%diag(:,1) - self%left_bound + self%diag(:,2:) = self%diag(:,2:) - self%sub + + end subroutine make_tridiag_deriv_diag + + ! Sum two TriDiagOp objects into a new one; this is just the addition of + ! all the entries. + function add_tridiag_ops(op1, op2) result(new_op) + + type(TriDiagOp), intent(in) :: op1, op2 + type(TriDiagOp) :: new_op + + new_op = op1 + + call new_op%add(op2) + + end function add_tridiag_ops + + subroutine add_in_place_tridiag_ops(self, other) + + class(TriDiagOp), intent(inout) :: self + class(TriDiagOp), intent(in) :: other + + self%spr = self%spr + other%spr + self%sub = self%sub + other%sub + self%diag = self%diag + other%diag + + self%left_bound = self%left_bound + other%left_bound + self%right_bound = self%right_bound + other%right_bound + + end subroutine add_in_place_tridiag_ops + + ! Subtract two TriDiagOp objects. + function subtract_tridiag_ops(op1, op2) result(new_op) + + type(TriDiagOp), intent(in) :: op1, op2 + type(TriDiagOp) :: new_op + + new_op = op1 + + call new_op%subtract(op2) + + end function subtract_tridiag_ops + + ! Subtract two TriDiagOp objects. + subroutine subtract_in_place_tridiag_ops(self, other) + + class(TriDiagOp), intent(inout) :: self + class(TriDiagOp), intent(in) :: other + + self%spr = self%spr - other%spr + self%sub = self%sub - other%sub + self%diag = self%diag - other%diag + + self%left_bound = self%left_bound - other%left_bound + self%right_bound = self%right_bound - other%right_bound + + end subroutine subtract_in_place_tridiag_ops + + ! Equivalent to adding a multiple of the identity. + subroutine scalar_add_tridiag(self, constant) + + class(TriDiagOp), intent(inout) :: self + real(r8), intent(in) :: constant + + self%diag = self%diag + constant + + end subroutine scalar_add_tridiag + + ! Equivalent to adding the diagonal operator constructed from diag_array. + subroutine diagonal_add_tridiag(self, diag_array) + + class(TriDiagOp), intent(inout) :: self + real(r8), USE_CONTIGUOUS intent(in) :: diag_array(:,:) + + self%diag = self%diag + diag_array + + end subroutine diagonal_add_tridiag + + ! Multiply a scalar by an array. + subroutine scalar_lmult_tridiag(self, constant) + + class(TriDiagOp), intent(inout) :: self + real(r8), intent(in) :: constant + + self%spr = self%spr * constant + self%sub = self%sub * constant + self%diag = self%diag * constant + + self%left_bound = self%left_bound * constant + self%right_bound = self%right_bound * constant + + end subroutine scalar_lmult_tridiag + + ! Multiply in an array as if it contained the entries of a diagonal matrix + ! being multiplied from the left. + subroutine diagonal_lmult_tridiag(self, diag_array) + + class(TriDiagOp), intent(inout) :: self + real(r8), USE_CONTIGUOUS intent(in) :: diag_array(:,:) + + self%spr = self%spr * diag_array(:,:self%ncel-1) + self%sub = self%sub * diag_array(:,2:) + self%diag = self%diag * diag_array(:,:) + + self%left_bound = self%left_bound * diag_array(:,1) + self%right_bound = self%right_bound * diag_array(:,self%ncel) + + end subroutine diagonal_lmult_tridiag + + ! Decomposition constructor + ! + ! The equation to be solved later (with left_div) is: + ! - A(k)*q(k+1) + B(k)*q(k) - C(k)*q(k-1) = D(k) + ! + ! The solution (effectively via LU decomposition) has the form: + ! E(k) = C(k) / (B(k) - A(k)*E(k+1)) + ! F(k) = (D(k) + A(k)*F(k+1)) / (B(k) - A(k)*E(k+1)) + ! q(k) = E(k) * q(k-1) + F(k) + ! + ! Unlike Richtmyer and Morton, E and F are defined by iterating backward + ! down to level 1, and then q iterates forward. + ! + ! E can be calculated and stored now, without knowing D. + ! To calculate F later, we store A and the denominator. + function new_TriDiagDecomp(op, graft_decomp) result(decomp) + type(TriDiagOp), intent(in) :: op + type(TriDiagDecomp), intent(in), optional :: graft_decomp + + type(TriDiagDecomp) :: decomp + + integer :: k + + if (present(graft_decomp)) then + decomp%nsys = graft_decomp%nsys + decomp%ncel = graft_decomp%ncel + else + decomp%nsys = op%nsys + decomp%ncel = op%ncel + end if + + ! Simple allocation with no error checking. + allocate(decomp%ca(decomp%nsys,decomp%ncel)) + allocate(decomp%dnom(decomp%nsys,decomp%ncel)) + allocate(decomp%ze(decomp%nsys,decomp%ncel)) + + ! decomp%ca is simply the negative of the tridiagonal's superdiagonal. + decomp%ca(:,:op%ncel-1) = -op%spr + decomp%ca(:,op%ncel) = -op%right_bound + + if (present(graft_decomp)) then + ! Copy in graft_decomp beyond op%ncel. + decomp%ca(:,op%ncel+1:) = graft_decomp%ca(:,op%ncel+1:) + decomp%dnom(:,op%ncel+1:) = graft_decomp%dnom(:,op%ncel+1:) + decomp%ze(:,op%ncel+1:) = graft_decomp%ze(:,op%ncel+1:) + ! Fill in dnom edge value. + decomp%dnom(:,op%ncel) = 1._r8 / (op%diag(:,op%ncel) - & + decomp%ca(:,op%ncel)*decomp%ze(:,op%ncel+1)) + else + ! If no grafting, the edge value of dnom comes from the diagonal. + decomp%dnom(:,op%ncel) = 1._r8 / op%diag(:,op%ncel) + end if + + do k = op%ncel - 1, 1, -1 + decomp%ze(:,k+1) = - op%sub(:,k) * decomp%dnom(:,k+1) + decomp%dnom(:,k) = 1._r8 / & + (op%diag(:,k) - decomp%ca(:,k)*decomp%ze(:,k+1)) + end do + + ! Don't multiply edge level by denom, because we want to leave it up to + ! the BoundaryCond object to decide what this means in left_div. + decomp%ze(:,1) = -op%left_bound + + end function new_TriDiagDecomp + + ! Left-division (multiplication by inverse) using a decomposed operator. + ! + ! See the comment above for the constructor for a quick explanation of the + ! intermediate variables. The "q" argument is "D(k)" on input and "q(k)" on + ! output. + subroutine decomp_left_div(decomp, q, l_cond, r_cond) + + ! Decomposed matrix. + class(TriDiagDecomp), intent(in) :: decomp + ! Data to left-divide by the matrix. + real(r8), USE_CONTIGUOUS intent(inout) :: q(:,:) + ! Objects representing boundary conditions. + class(BoundaryCond), intent(in), optional :: l_cond, r_cond + + ! "F" from the equation above. + real(r8) :: zf(decomp%nsys,decomp%ncel) + + ! Level index. + integer :: k + + ! Include boundary conditions. + if (present(l_cond)) then + q(:,1) = q(:,1) + l_cond%apply_left(decomp%ze(:,1), q) + end if + + if (present(r_cond)) then + q(:,decomp%ncel) = q(:,decomp%ncel) + & + r_cond%apply_right(decomp%ca(:,decomp%ncel), q) + end if + + zf(:,decomp%ncel) = q(:,decomp%ncel) * decomp%dnom(:,decomp%ncel) + + do k = decomp%ncel - 1, 1, -1 + zf(:,k) = (q(:,k) + decomp%ca(:,k)*zf(:,k+1)) * decomp%dnom(:,k) + end do + + ! Perform back substitution + + q(:,1) = zf(:,1) + + do k = 2, decomp%ncel + q(:,k) = zf(:,k) + decomp%ze(:,k)*q(:,k-1) + end do + + end subroutine decomp_left_div + + ! Decomposition deallocation. + subroutine decomp_finalize(decomp) + class(TriDiagDecomp), intent(inout) :: decomp + + decomp%nsys = 0 + decomp%ncel = 0 + + if (allocated(decomp%ca)) deallocate(decomp%ca) + if (allocated(decomp%dnom)) deallocate(decomp%dnom) + if (allocated(decomp%ze)) deallocate(decomp%ze) + + end subroutine decomp_finalize + +end module linear_1d_operators + \ No newline at end of file diff --git a/to-be-ccpp-ized/utilities/shr_infnan_mod.F90.in b/to-be-ccpp-ized/utilities/shr_infnan_mod.F90.in new file mode 100644 index 00000000..992c46fc --- /dev/null +++ b/to-be-ccpp-ized/utilities/shr_infnan_mod.F90.in @@ -0,0 +1,406 @@ +! Flag representing compiler support of Fortran 2003's +! ieee_arithmetic intrinsic module. +#if defined CPRIBM || defined CPRPGI || defined CPRINTEL || defined CPRCRAY || defined CPRNAG +#define HAVE_IEEE_ARITHMETIC +#endif + +module shr_infnan_mod +!--------------------------------------------------------------------- +! Module to test for IEEE Inf and NaN values, which also provides a +! method of setting +/-Inf and signaling or quiet NaN. +! +! All functions are elemental, and thus work on arrays. +!--------------------------------------------------------------------- +! To test for these values, just call the corresponding function, e.g: +! +! var_is_nan = shr_infnan_isnan(x) +! +! You can also use it on arrays: +! +! array_contains_nan = any(shr_infnan_isnan(my_array)) +! +!--------------------------------------------------------------------- +! To generate these values, assign one of the provided derived-type +! variables to a real: +! +! use shr_infnan_mod, only: nan => shr_infnan_nan, & +! inf => shr_infnan_inf, & +! assignment(=) +! real(r4) :: my_nan +! real(r8) :: my_inf_array(2,2) +! my_nan = nan +! my_inf_array = inf +! +! Keep in mind that "shr_infnan_nan" and "shr_infnan_inf" cannot be +! passed to functions that expect real arguments. To pass a real +! NaN, you will have to use shr_infnan_nan to set a local real of +! the correct kind. +!--------------------------------------------------------------------- + +use shr_kind_mod, only: & + r4 => SHR_KIND_R4, & + r8 => SHR_KIND_R8 + +#ifdef HAVE_IEEE_ARITHMETIC + +! If we have IEEE_ARITHMETIC, the NaN test is provided for us. +use, intrinsic :: ieee_arithmetic, only: & + shr_infnan_isnan => ieee_is_nan + +#else + +! Integers of correct size for bit patterns below. +use shr_kind_mod, only: i4 => shr_kind_i4, i8 => shr_kind_i8 + +#endif + +implicit none +private +save + +! Test functions for NaN/Inf values. +public :: shr_infnan_isnan +public :: shr_infnan_isinf +public :: shr_infnan_isposinf +public :: shr_infnan_isneginf + +! Locally defined isnan. +#ifndef HAVE_IEEE_ARITHMETIC +interface shr_infnan_isnan + ! TYPE double,real + module procedure shr_infnan_isnan_{TYPE} +end interface +#endif + +interface shr_infnan_isinf + ! TYPE double,real + module procedure shr_infnan_isinf_{TYPE} +end interface + +interface shr_infnan_isposinf + ! TYPE double,real + module procedure shr_infnan_isposinf_{TYPE} +end interface + +interface shr_infnan_isneginf + ! TYPE double,real + module procedure shr_infnan_isneginf_{TYPE} +end interface + +! Derived types for generation of NaN/Inf +! Even though there's no reason to "use" the types directly, some compilers +! might have trouble with an object being used without its type. +public :: shr_infnan_nan_type +public :: shr_infnan_inf_type +public :: assignment(=) +public :: shr_infnan_to_r4 +public :: shr_infnan_to_r8 + +! Type representing Not A Number. +type :: shr_infnan_nan_type + logical :: quiet = .false. +end type shr_infnan_nan_type + +! Type representing +/-Infinity. +type :: shr_infnan_inf_type + logical :: positive = .true. +end type shr_infnan_inf_type + +! Allow assigning reals to NaN or Inf. +interface assignment(=) + ! TYPE double,real + ! DIMS 0,1,2,3,4,5,6,7 + module procedure set_nan_{DIMS}d_{TYPE} + ! TYPE double,real + ! DIMS 0,1,2,3,4,5,6,7 + module procedure set_inf_{DIMS}d_{TYPE} +end interface + +! Conversion functions. +interface shr_infnan_to_r8 + module procedure nan_r8 + module procedure inf_r8 +end interface + +interface shr_infnan_to_r4 + module procedure nan_r4 + module procedure inf_r4 +end interface + +! Initialize objects of NaN/Inf type for other modules to use. + +! Default NaN is signaling, but also provide snan and qnan to choose +! explicitly. +type(shr_infnan_nan_type), public, parameter :: shr_infnan_nan = & + shr_infnan_nan_type(.false.) +type(shr_infnan_nan_type), public, parameter :: shr_infnan_snan = & + shr_infnan_nan_type(.false.) +type(shr_infnan_nan_type), public, parameter :: shr_infnan_qnan = & + shr_infnan_nan_type(.true.) + +! Default Inf is positive, but provide posinf to go with neginf. +type(shr_infnan_inf_type), public, parameter :: shr_infnan_inf = & + shr_infnan_inf_type(.true.) +type(shr_infnan_inf_type), public, parameter :: shr_infnan_posinf = & + shr_infnan_inf_type(.true.) +type(shr_infnan_inf_type), public, parameter :: shr_infnan_neginf = & + shr_infnan_inf_type(.false.) + +! Bit patterns for implementation without ieee_arithmetic. +! Note that in order to satisfy gfortran's range check, we have to use +! ibset to set the sign bit from a BOZ pattern. +#ifndef HAVE_IEEE_ARITHMETIC +! Single precision. +integer(i4), parameter :: ssnan_pat = int(Z'7FA00000',i4) +integer(i4), parameter :: sqnan_pat = int(Z'7FC00000',i4) +integer(i4), parameter :: sposinf_pat = int(Z'7F800000',i4) +integer(i4), parameter :: sneginf_pat = ibset(sposinf_pat,bit_size(1_i4)-1) +! Double precision. +integer(i8), parameter :: dsnan_pat = int(Z'7FF4000000000000',i8) +integer(i8), parameter :: dqnan_pat = int(Z'7FF8000000000000',i8) +integer(i8), parameter :: dposinf_pat = int(Z'7FF0000000000000',i8) +integer(i8), parameter :: dneginf_pat = ibset(dposinf_pat,bit_size(1_i8)-1) +#endif + +contains + +!--------------------------------------------------------------------- +! TEST FUNCTIONS +!--------------------------------------------------------------------- +! The "isinf" function simply calls "isposinf" and "isneginf". +!--------------------------------------------------------------------- + +! TYPE double,real +elemental function shr_infnan_isinf_{TYPE}(x) result(isinf) + {VTYPE}, intent(in) :: x + logical :: isinf + + isinf = shr_infnan_isposinf(x) .or. shr_infnan_isneginf(x) + +end function shr_infnan_isinf_{TYPE} + +#ifdef HAVE_IEEE_ARITHMETIC + +!--------------------------------------------------------------------- +! The "isposinf" and "isneginf" functions get the IEEE class of a +! real, and test to see if the class is equal to ieee_positive_inf +! or ieee_negative_inf. +!--------------------------------------------------------------------- + +! TYPE double,real +elemental function shr_infnan_isposinf_{TYPE}(x) result(isposinf) + use, intrinsic :: ieee_arithmetic, only: & + ieee_class, & + ieee_positive_inf, & + operator(==) + {VTYPE}, intent(in) :: x + logical :: isposinf + + isposinf = (ieee_positive_inf == ieee_class(x)) + +end function shr_infnan_isposinf_{TYPE} + +! TYPE double,real +elemental function shr_infnan_isneginf_{TYPE}(x) result(isneginf) + use, intrinsic :: ieee_arithmetic, only: & + ieee_class, & + ieee_negative_inf, & + operator(==) + {VTYPE}, intent(in) :: x + logical :: isneginf + + isneginf = (ieee_negative_inf == ieee_class(x)) + +end function shr_infnan_isneginf_{TYPE} + +#else +! Don't have ieee_arithmetic. + +#ifdef CPRGNU +! NaN testing on gfortran. +! TYPE double,real +elemental function shr_infnan_isnan_{TYPE}(x) result(is_nan) + {VTYPE}, intent(in) :: x + logical :: is_nan + + is_nan = isnan(x) + +end function shr_infnan_isnan_{TYPE} +! End GNU section. +#endif + +!--------------------------------------------------------------------- +! The "isposinf" and "isneginf" functions just test against a known +! bit pattern if we don't have ieee_arithmetic. +!--------------------------------------------------------------------- + +! TYPE double,real +elemental function shr_infnan_isposinf_{TYPE}(x) result(isposinf) + {VTYPE}, intent(in) :: x + logical :: isposinf +#if ({ITYPE} == TYPEREAL) + integer(i4), parameter :: posinf_pat = sposinf_pat +#else + integer(i8), parameter :: posinf_pat = dposinf_pat +#endif + + isposinf = (x == transfer(posinf_pat,x)) + +end function shr_infnan_isposinf_{TYPE} + +! TYPE double,real +elemental function shr_infnan_isneginf_{TYPE}(x) result(isneginf) + {VTYPE}, intent(in) :: x + logical :: isneginf +#if ({ITYPE} == TYPEREAL) + integer(i4), parameter :: neginf_pat = sneginf_pat +#else + integer(i8), parameter :: neginf_pat = dneginf_pat +#endif + + isneginf = (x == transfer(neginf_pat,x)) + +end function shr_infnan_isneginf_{TYPE} + +! End ieee_arithmetic conditional. +#endif + +!--------------------------------------------------------------------- +! GENERATION FUNCTIONS +!--------------------------------------------------------------------- +! Two approaches for generation of NaN and Inf values: +! 1. With Fortran 2003, use the ieee_value intrinsic to get a value +! from the corresponding class. These are: +! - ieee_signaling_nan +! - ieee_quiet_nan +! - ieee_positive_inf +! - ieee_negative_inf +! 2. Without Fortran 2003, set the IEEE bit patterns directly. +! Use BOZ literals to get an integer with the correct bit +! pattern, then use "transfer" to transfer those bits into a +! real. +!--------------------------------------------------------------------- + +! TYPE double,real +! DIMS 0,1,2,3,4,5,6,7 +pure subroutine set_nan_{DIMS}d_{TYPE}(output, nan) +#ifdef HAVE_IEEE_ARITHMETIC + use, intrinsic :: ieee_arithmetic, only: & + ieee_signaling_nan, & + ieee_quiet_nan, & + ieee_value +#else +#if ({ITYPE} == TYPEREAL) + integer(i4), parameter :: snan_pat = ssnan_pat + integer(i4), parameter :: qnan_pat = sqnan_pat +#else + integer(i8), parameter :: snan_pat = dsnan_pat + integer(i8), parameter :: qnan_pat = dqnan_pat +#endif +#endif + {VTYPE}, intent(out) :: output{DIMSTR} + type(shr_infnan_nan_type), intent(in) :: nan + + ! Use scalar temporary for performance reasons, to reduce the cost of + ! the ieee_value call. + {VTYPE} :: tmp + +#ifdef HAVE_IEEE_ARITHMETIC + if (nan%quiet) then + tmp = ieee_value(tmp, ieee_quiet_nan) + else + tmp = ieee_value(tmp, ieee_signaling_nan) + end if +#else + if (nan%quiet) then + tmp = transfer(qnan_pat, tmp) + else + tmp = transfer(snan_pat, tmp) + end if +#endif + + output = tmp + +end subroutine set_nan_{DIMS}d_{TYPE} + +! TYPE double,real +! DIMS 0,1,2,3,4,5,6,7 +pure subroutine set_inf_{DIMS}d_{TYPE}(output, inf) +#ifdef HAVE_IEEE_ARITHMETIC + use, intrinsic :: ieee_arithmetic, only: & + ieee_positive_inf, & + ieee_negative_inf, & + ieee_value +#else +#if ({ITYPE} == TYPEREAL) + integer(i4), parameter :: posinf_pat = sposinf_pat + integer(i4), parameter :: neginf_pat = sneginf_pat +#else + integer(i8), parameter :: posinf_pat = dposinf_pat + integer(i8), parameter :: neginf_pat = dneginf_pat +#endif +#endif + {VTYPE}, intent(out) :: output{DIMSTR} + type(shr_infnan_inf_type), intent(in) :: inf + + ! Use scalar temporary for performance reasons, to reduce the cost of + ! the ieee_value call. + {VTYPE} :: tmp + +#ifdef HAVE_IEEE_ARITHMETIC + if (inf%positive) then + tmp = ieee_value(tmp,ieee_positive_inf) + else + tmp = ieee_value(tmp,ieee_negative_inf) + end if +#else + if (inf%positive) then + tmp = transfer(posinf_pat, tmp) + else + tmp = transfer(neginf_pat, tmp) + end if +#endif + + output = tmp + +end subroutine set_inf_{DIMS}d_{TYPE} + +!--------------------------------------------------------------------- +! CONVERSION INTERFACES. +!--------------------------------------------------------------------- +! Function methods to get reals from nan/inf types. +!--------------------------------------------------------------------- + +pure function nan_r8(nan) result(output) + class(shr_infnan_nan_type), intent(in) :: nan + real(r8) :: output + + output = nan + +end function nan_r8 + +pure function nan_r4(nan) result(output) + class(shr_infnan_nan_type), intent(in) :: nan + real(r4) :: output + + output = nan + +end function nan_r4 + +pure function inf_r8(inf) result(output) + class(shr_infnan_inf_type), intent(in) :: inf + real(r8) :: output + + output = inf + +end function inf_r8 + +pure function inf_r4(inf) result(output) + class(shr_infnan_inf_type), intent(in) :: inf + real(r4) :: output + + output = inf + +end function inf_r4 + +end module shr_infnan_mod diff --git a/to-be-ccpp-ized/utilities/shr_kind_mod.F90 b/to-be-ccpp-ized/utilities/shr_kind_mod.F90 new file mode 100644 index 00000000..9ab62a4e --- /dev/null +++ b/to-be-ccpp-ized/utilities/shr_kind_mod.F90 @@ -0,0 +1,21 @@ +MODULE shr_kind_mod + + !---------------------------------------------------------------------------- + ! precision/kind constants add data public + !---------------------------------------------------------------------------- + public + integer,parameter :: SHR_KIND_R8 = selected_real_kind(12) ! 8 byte real + integer,parameter :: SHR_KIND_R4 = selected_real_kind( 6) ! 4 byte real + integer,parameter :: SHR_KIND_RN = kind(1.0) ! native real + integer,parameter :: SHR_KIND_I8 = selected_int_kind (13) ! 8 byte integer + integer,parameter :: SHR_KIND_I4 = selected_int_kind ( 6) ! 4 byte integer + integer,parameter :: SHR_KIND_I2 = selected_int_kind ( 4) ! 2 byte integer + integer,parameter :: SHR_KIND_IN = kind(1) ! native integer + integer,parameter :: SHR_KIND_CS = 80 ! short char + integer,parameter :: SHR_KIND_CM = 160 ! mid-sized char + integer,parameter :: SHR_KIND_CL = 256 ! long char + integer,parameter :: SHR_KIND_CX = 512 ! extra-long char + integer,parameter :: SHR_KIND_CXX= 4096 ! extra-extra-long char + +END MODULE shr_kind_mod + \ No newline at end of file diff --git a/to-be-ccpp-ized/utilities/shr_log_mod.F90 b/to-be-ccpp-ized/utilities/shr_log_mod.F90 new file mode 100644 index 00000000..e028f497 --- /dev/null +++ b/to-be-ccpp-ized/utilities/shr_log_mod.F90 @@ -0,0 +1,121 @@ +!BOP =========================================================================== +! +! !MODULE: shr_log_mod -- variables and methods for logging +! +! !DESCRIPTION: +! Low-level shared variables for logging. +! +! Also, routines for generating log file messages. +! +! !INTERFACE: ------------------------------------------------------------------ + +module shr_log_mod + + ! !USES: + + use shr_kind_mod, only: shr_kind_in, shr_kind_cx + use shr_strconvert_mod, only: toString + + use, intrinsic :: iso_fortran_env, only: output_unit + + implicit none + private + + ! !PUBLIC TYPES: + + ! no public types + + ! !PUBLIC MEMBER FUNCTIONS: + + public :: shr_log_errMsg + public :: shr_log_OOBMsg + public :: shr_log_setLogUnit + public :: shr_log_getLogUnit + + ! !PUBLIC DATA MEMBERS: + + public :: shr_log_Level + public :: shr_log_Unit + + !EOP + + ! low-level shared variables for logging, these may not be parameters + integer(SHR_KIND_IN) :: shr_log_Level = 0 + integer(SHR_KIND_IN) :: shr_log_Unit = output_unit + + contains + + !=============================================================================== + !BOP =========================================================================== + ! + ! !IROUTINE: shr_log_errMsg -- Return an error message containing file & line info + ! + ! !DESCRIPTION: + ! Return an error message containing file & line info + ! \newline + ! errMsg = shr\_log\_errMsg(__FILE__, __LINE__) + ! + ! This is meant to be used when a routine expects a string argument for some message, + ! but you want to provide file and line information. + ! + ! However: Note that the performance of this function can be very bad. It is currently + ! maintained because it is used by old code, but you should probably avoid using this + ! in new code if possible. + ! + ! !REVISION HISTORY: + ! 2013-July-23 - Bill Sacks + ! + ! !INTERFACE: ------------------------------------------------------------------ + + pure function shr_log_errMsg(file, line) + + ! !INPUT/OUTPUT PARAMETERS: + + character(len=SHR_KIND_CX) :: shr_log_errMsg + character(len=*), intent(in) :: file + integer , intent(in) :: line + + !EOP + + shr_log_errMsg = 'ERROR in '//trim(file)//' at line '//toString(line) + + end function shr_log_errMsg + + ! Create a message for an out of bounds error. + pure function shr_log_OOBMsg(operation, bounds, idx) result(OOBMsg) + + ! A name for the operation being attempted when the bounds error + ! occurred. A string containing the subroutine name is ideal, but more + ! generic descriptions such as "read", "modify", or "insert" could be used. + character(len=*), intent(in) :: operation + + ! Upper and lower bounds allowed for the operation. + integer, intent(in) :: bounds(2) + + ! Index at which access was attempted. + integer, intent(in) :: idx + + ! Output message + character(len=:), allocatable :: OOBMsg + + allocate(OOBMsg, source=(operation//": "//toString(idx)//" not in range ["//& + toString(bounds(1))//", "//toString(bounds(2))//"].")) + + end function shr_log_OOBMsg + + subroutine shr_log_setLogUnit(unit) + integer, intent(in) :: unit + + shr_log_unit = unit + + end subroutine shr_log_setLogUnit + + subroutine shr_log_getLogUnit(unit) + integer, intent(out) :: unit + + unit = shr_log_unit + + end subroutine shr_log_getLogUnit + + end module shr_log_mod + \ No newline at end of file diff --git a/to-be-ccpp-ized/utilities/shr_strconvert_mod.F90 b/to-be-ccpp-ized/utilities/shr_strconvert_mod.F90 new file mode 100644 index 00000000..56ec172a --- /dev/null +++ b/to-be-ccpp-ized/utilities/shr_strconvert_mod.F90 @@ -0,0 +1,167 @@ +module shr_strconvert_mod + + ! This module defines toString, a generic function for creating character type + ! representations of data, as implemented for the most commonly used intrinsic + ! types: + ! + ! - 4 and 8 byte integer + ! - 4 and 8 byte real + ! - logical + ! + ! No toString implementation is provided for character input, but this may be + ! added if some use case arises. + ! + ! Currently, only scalar inputs are supported. The return type of this function + ! is character with deferred (allocatable) length. + ! + ! The functions for integers and reals allow an optional format_string argument, + ! which can be used to control the padding and precision of output as with any + ! write statement. However, the implementations internally must use a + ! preallocated buffer, so a format_string that significantly increases the size + ! of the output may cause a run-time error or undefined behavior in the program. + ! + ! Other modules may want to provide extensions of toString for their own derived + ! types. In this case there are two guidelines to observe: + ! + ! - It is preferable to have only one mandatory argument, which is the object to + ! produce a string from. There may be other formatting options, but the + ! implementation should do something sensible without these. + ! + ! - Since the main purpose of toString is to provide a human-readable + ! representation of a type, especially for documentation or debugging + ! purposes, refrain from printing large array components in their entirety + ! (instead consider printing only the shape, or statistics such as + ! min/mean/max for arrays of numbers). + + use shr_kind_mod, only: & + i4 => shr_kind_i4, & + i8 => shr_kind_i8, & + r4 => shr_kind_r4, & + r8 => shr_kind_r8, & + cs => shr_kind_cs + + use shr_infnan_mod, only: & + isnan => shr_infnan_isnan + + implicit none + private + + ! Human-readable representation of data. + public :: toString + + interface toString + module procedure i4ToString + module procedure i8ToString + module procedure r4ToString + module procedure r8ToString + module procedure logicalToString + end interface toString + + contains + + pure function i4ToString(input, format_string) result(string) + integer(i4), intent(in) :: input + character(len=*), intent(in), optional :: format_string + character(len=:), allocatable :: string + + character(len=cs) :: buffer + + if (present(format_string)) then + write(buffer, format_string) input + else + ! For most compilers, these two statements are equivalent to a format of + ! '(I0)', but that's not technically in the standard. + write(buffer, '(I11)') input + buffer = adjustl(buffer) + end if + + allocate(string, source=trim(buffer)) + + end function i4ToString + + pure function i8ToString(input, format_string) result(string) + integer(i8), intent(in) :: input + character(len=*), intent(in), optional :: format_string + character(len=:), allocatable :: string + + character(len=cs) :: buffer + + if (present(format_string)) then + write(buffer, format_string) input + else + ! For most compilers, these two statements are equivalent to a format of + ! '(I0)', but that's not technically in the standard. + write(buffer, '(I20)') input + buffer = adjustl(buffer) + end if + + allocate(string, source=trim(buffer)) + + end function i8ToString + + pure function r4ToString(input, format_string) result(string) + real(r4), intent(in) :: input + character(len=*), intent(in), optional :: format_string + character(len=:), allocatable :: string + + character(len=cs) :: buffer + + if (present(format_string)) then + write(buffer, format_string) input + else + write(buffer, '(ES15.8 E2)') input + buffer = adjustl(buffer) + ! Deal with the fact that the "+" sign is optional by simply adding it if + ! it is not present, so that the default format is standardized across + ! compilers. + ! Assumes that compilers do not treat the sign bit on NaN values specially. + if (.not. isnan(input) .and. all(buffer(1:1) /= ["-", "+"])) then + buffer = "+" // trim(buffer) + end if + end if + + allocate(string, source=trim(buffer)) + + end function r4ToString + + pure function r8ToString(input, format_string) result(string) + real(r8), intent(in) :: input + character(len=*), intent(in), optional :: format_string + character(len=:), allocatable :: string + + character(len=cs) :: buffer + + if (present(format_string)) then + write(buffer, format_string) input + else + write(buffer, '(ES24.16 E3)') input + buffer = adjustl(buffer) + ! Deal with the fact that the "+" sign is optional by simply adding it if + ! it is not present, so that the default format is standardized across + ! compilers. + ! Assumes that compilers do not treat the sign bit on NaN values specially. + if (.not. isnan(input) .and. all(buffer(1:1) /= ["-", "+"])) then + buffer = "+" // trim(buffer) + end if + end if + + allocate(string, source=trim(buffer)) + + end function r8ToString + + pure function logicalToString(input) result(string) + logical, intent(in) :: input + character(len=:), allocatable :: string + + ! We could use a write statement, but this is easier. + allocate(character(len=1) :: string) + if (input) then + string = "T" + else + string = "F" + end if + + end function logicalToString + +end module shr_strconvert_mod + \ No newline at end of file diff --git a/to-be-ccpp-ized/utilities/shr_sys_mod.F90 b/to-be-ccpp-ized/utilities/shr_sys_mod.F90 new file mode 100644 index 00000000..bf463d55 --- /dev/null +++ b/to-be-ccpp-ized/utilities/shr_sys_mod.F90 @@ -0,0 +1,332 @@ +!=============================================================================== +! SVN $Id: shr_sys_mod.F90 66411 2014-12-19 22:40:08Z santos@ucar.edu $ +! SVN $URL: https://svn-ccsm-models.cgd.ucar.edu/csm_share/trunk_tags/share3_150116/shr/shr_sys_mod.F90 $ +!=============================================================================== + +! Currently supported by all compilers +#define HAVE_GET_ENVIRONMENT +#define HAVE_SLEEP + +! Except this combination? +#if defined CPRPGI && defined CNL +#undef HAVE_GET_ENVIRONMENT +#endif + +#if defined CPRNAG +#define HAVE_EXECUTE +#endif + +MODULE shr_sys_mod + + use shr_kind_mod ! defines real & integer kinds + use shr_log_mod, only: s_loglev => shr_log_Level + use shr_log_mod, only: s_logunit => shr_log_Unit + use shr_abort_mod, only: shr_sys_abort => shr_abort_abort + use shr_abort_mod, only: shr_sys_backtrace => shr_abort_backtrace + +#ifdef CPRNAG + ! NAG does not provide these as intrinsics, but it does provide modules + ! that implement commonly used POSIX routines. + use f90_unix_dir, only: chdir + use f90_unix_proc, only: abort, sleep +#endif + + implicit none + +! PUBLIC: Public interfaces + + private + + public :: shr_sys_system ! make a system call + public :: shr_sys_chdir ! change current working dir + public :: shr_sys_getenv ! get an environment variable + public :: shr_sys_irtc ! returns real-time clock tick + public :: shr_sys_sleep ! have program sleep for a while + public :: shr_sys_flush ! flush an i/o buffer + + ! Imported from shr_abort_mod and republished with renames. Other code that wishes to + ! use these routines should use these shr_sys names rather than directly using the + ! routines from shr_abort_abort. (This is for consistency with older code, from when + ! these routines were defined in shr_sys_mod.) + public :: shr_sys_abort ! abort a program + public :: shr_sys_backtrace ! print a backtrace, if possible + +!=============================================================================== +CONTAINS +!=============================================================================== + +!=============================================================================== +!=============================================================================== + +SUBROUTINE shr_sys_system(str,rcode) + + IMPLICIT none + + !----- arguments --- + character(*) ,intent(in) :: str ! system/shell command string + integer(SHR_KIND_IN),intent(out) :: rcode ! function return error code + + !----- functions ----- +#if (defined LINUX && !defined CPRGNU) + integer(SHR_KIND_IN),external :: system ! function to envoke shell command +#endif + + !----- formats ----- + character(*),parameter :: subName = '(shr_sys_system) ' + character(*),parameter :: F00 = "('(shr_sys_system) ',4a)" + +!------------------------------------------------------------------------------- +! PURPOSE: an architecture independent system call +!------------------------------------------------------------------------------- + rcode = 0 +#ifdef HAVE_EXECUTE + call execute_command_line(str,exitstat=rcode) ! Intrinsic as of F2008 +#else +#if (defined AIX) + + call system(str,rcode) + +#elif (defined CPRGNU || defined LINUX) + + rcode = system(str) + +#else + + write(s_logunit,F00) 'ERROR: no implementation of system call for this architecture' + call shr_sys_abort(subName//'no implementation of system call for this architecture') +#endif +#endif + +END SUBROUTINE shr_sys_system + +!=============================================================================== +!=============================================================================== + +SUBROUTINE shr_sys_chdir(path, rcode) + + IMPLICIT none + + !----- arguments ----- + character(*) ,intent(in) :: path ! chdir to this dir + integer(SHR_KIND_IN),intent(out) :: rcode ! return code + + !----- local ----- + integer(SHR_KIND_IN) :: lenpath ! length of path +#if (defined AIX || (defined LINUX && !defined CPRGNU && !defined CPRNAG) || defined CPRINTEL) + integer(SHR_KIND_IN),external :: chdir ! AIX system call +#endif + + !----- formats ----- + character(*),parameter :: subName = '(shr_sys_chdir) ' + character(*),parameter :: F00 = "('(shr_sys_chdir) ',4a)" + +!------------------------------------------------------------------------------- +! PURPOSE: an architecture independent system call +!------------------------------------------------------------------------------- + rcode = 0 + lenpath=len_trim(path) + +#if (defined AIX) + + rcode = chdir(%ref(path(1:lenpath)//'\0')) + +#elif (defined Darwin || (defined LINUX && !defined CPRNAG)) + + rcode=chdir(path(1:lenpath)) + +#elif (defined CPRNAG) + + call chdir(path(1:lenpath), errno=rcode) + +#else + + write(s_logunit,F00) 'ERROR: no implementation of chdir for this architecture' + call shr_sys_abort(subname//'no implementation of chdir for this machine') + +#endif + +END SUBROUTINE shr_sys_chdir + +!=============================================================================== +!=============================================================================== + +SUBROUTINE shr_sys_getenv(name, val, rcode) + + IMPLICIT none + + !----- arguments ----- + character(*) ,intent(in) :: name ! env var name + character(*) ,intent(out) :: val ! env var value + integer(SHR_KIND_IN),intent(out) :: rcode ! return code + + !----- local ----- +#ifndef HAVE_GET_ENVIRONMENT + integer(SHR_KIND_IN) :: lenname ! length of env var name + integer(SHR_KIND_IN) :: lenval ! length of env var value + character(SHR_KIND_CL) :: tmpval ! temporary env var value +#endif + !----- formats ----- + character(*),parameter :: subName = '(shr_sys_getenv) ' + character(*),parameter :: F00 = "('(shr_sys_getenv) ',4a)" + +!------------------------------------------------------------------------------- +! PURPOSE: an architecture independent system call +!------------------------------------------------------------------------------- + +!$OMP master + + +#ifdef HAVE_GET_ENVIRONMENT + call get_environment_variable(name=name,value=val,status=rcode) ! Intrinsic in F2003 +#else + lenname=len_trim(name) +#if (defined AIX || defined LINUX) + + call getenv(trim(name),tmpval) + val=trim(tmpval) + rcode = 0 + if (len_trim(val) == 0 ) rcode = 1 + if (len_trim(val) > SHR_KIND_CL) rcode = 2 + +#else + + write(s_logunit,F00) 'ERROR: no implementation of getenv for this architecture' + call shr_sys_abort(subname//'no implementation of getenv for this machine') + +#endif +#endif +!$OMP end master + +END SUBROUTINE shr_sys_getenv + +!=============================================================================== +!=============================================================================== + +integer(SHR_KIND_I8) FUNCTION shr_sys_irtc( rate ) + + IMPLICIT none + + !----- arguments ----- + integer(SHR_KIND_I8), optional :: rate + + !----- local ----- + integer(SHR_KIND_IN) :: count + integer(SHR_KIND_IN) :: count_rate + integer(SHR_KIND_IN) :: count_max + + integer(SHR_KIND_IN),save :: last_count = -1 + integer(SHR_KIND_I8),save :: count_offset = 0 +!$OMP THREADPRIVATE (last_count, count_offset) + + !----- formats ----- + character(*),parameter :: subName = '(shr_sys_irtc) ' + character(*),parameter :: F00 = "('(shr_sys_irtc) ',4a)" + +!------------------------------------------------------------------------------- +! emulates Cray/SGI irtc function (returns clock tick since last reboot) +! +! This function is not intended to measure elapsed time between +! multi-threaded regions with different numbers of threads. However, +! use of the threadprivate declaration does guarantee accurate +! measurement per thread within a single multi-threaded region as +! long as the number of threads is not changed dynamically during +! execution within the multi-threaded region. +! +!------------------------------------------------------------------------------- + + call system_clock(count=count,count_rate=count_rate, count_max=count_max) + if ( present(rate) ) rate = count_rate + shr_sys_irtc = count + + !--- adjust for clock wrap-around --- + if ( last_count /= -1 ) then + if ( count < last_count ) count_offset = count_offset + count_max + end if + shr_sys_irtc = shr_sys_irtc + count_offset + last_count = count + +END FUNCTION shr_sys_irtc + +!=============================================================================== +!=============================================================================== + +SUBROUTINE shr_sys_sleep(sec) + + IMPLICIT none + + !----- arguments ----- + real (SHR_KIND_R8),intent(in) :: sec ! number of seconds to sleep + + !----- local ----- + integer(SHR_KIND_IN) :: isec ! integer number of seconds +#ifndef HAVE_SLEEP + integer(SHR_KIND_IN) :: rcode ! return code + character(90) :: str ! system call string +#endif + !----- formats ----- + character(*),parameter :: subName = '(shr_sys_sleep) ' + character(*),parameter :: F00 = "('(shr_sys_sleep) ',4a)" + character(*),parameter :: F10 = "('sleep ',i8 )" + +!------------------------------------------------------------------------------- +! PURPOSE: Sleep for approximately sec seconds +!------------------------------------------------------------------------------- + + isec = nint(sec) + + if (isec < 0) then + if (s_loglev > 0) write(s_logunit,F00) 'ERROR: seconds must be > 0, sec=',sec + else if (isec == 0) then + ! Don't consider this an error and don't call system sleep + else +#ifdef HAVE_SLEEP + call sleep(isec) +#else + write(str,FMT=F10) isec + call shr_sys_system( str, rcode ) +#endif + endif + +END SUBROUTINE shr_sys_sleep + +!=============================================================================== +!=============================================================================== + +SUBROUTINE shr_sys_flush(unit) + + IMPLICIT none + + !----- arguments ----- + integer(SHR_KIND_IN) :: unit ! flush output buffer for this unit + + !----- local ----- + !----- formats ----- + character(*),parameter :: subName = '(shr_sys_flush) ' + character(*),parameter :: F00 = "('(shr_sys_flush) ',4a)" + +!------------------------------------------------------------------------------- +! PURPOSE: an architecture independent system call +! +! This is probably no longer needed; the "flush" statement is supported by +! all compilers that CESM supports for years now. +! +!------------------------------------------------------------------------------- +!$OMP SINGLE + flush(unit) +!$OMP END SINGLE +! +! The following code was originally present, but there's an obvious issue. +! Since shr_sys_flush is usually used to flush output to a log, when it +! returns an error, does it do any good to print that error to the log? +! +! if (ierr > 0) then +! write(s_logunit,*) subname,' Flush reports error: ',ierr +! endif +! + +END SUBROUTINE shr_sys_flush + +!=============================================================================== +!=============================================================================== + +END MODULE shr_sys_mod diff --git a/to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion.F90 b/to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion.F90 new file mode 100644 index 00000000..4e9d87d5 --- /dev/null +++ b/to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion.F90 @@ -0,0 +1,136 @@ +module vertical_diffusion + + implicit none + private + save + + public :: fin_vol_solve + + contains + + ! Designed to solve the equation: + ! + ! w * dq/dt = d/dp (D q' - v q) + c q + ! + ! where q is a grid-cell average, and p is the vertical coordinate + ! (presumably pressure). + ! + ! In this function, coef_q_weight == w, coef_q_diff == D, + ! coef_q_adv == v, and coef_q == c. All these are optional; omitting a + ! coefficient is equivalent to setting the entire array to 0. + ! + ! coef_q_diff and coef_q_adv are defined at the level interfaces, while + ! coef_q and coef_q_weight are grid-cell averages. + + function fin_vol_solve(dt, p, toSolve, ncols, pver, coef_q, coef_q_diff, coef_q_adv, & + coef_q_weight, upper_bndry, lower_bndry, l_cond, r_cond) result(solution) + + use linear_1d_operators, only: & + zero_operator, & + diagonal_operator, & + diffusion_operator, & + advection_operator, & + BoundaryType, & + TriDiagDecomp, & + TriDiagOp, & + BoundaryCond, & + operator(+) + use shr_kind_mod, only: r8 => shr_kind_r8 + use coords_1d, only: Coords1D + + ! ---------------------- ! + ! Input-Output Arguments ! + ! ---------------------- ! + + ! Time step. + real(r8), intent(in) :: dt + ! Grid spacings. + type(Coords1D), intent(in) :: p + + ! Matrix to decomp from. + ! real(r8), intent(in) :: u(ncols,pver) + integer, intent(in) :: ncols + integer, intent(in) :: pver + real(r8), intent(in) :: toSolve(ncols,pver) + + ! Coefficients for diffusion and advection. + ! + ! The sizes must be consistent among all the coefficients that are + ! actually present, i.e. coef_q_diff and coef_q_adv should be one level + ! bigger than coef_q and coef_q_weight, and have the same column number. + real(r8), contiguous, intent(in), optional :: coef_q(:,:), & + coef_q_diff(:,:), coef_q_adv(:,:), coef_q_weight(:,:) + + ! Boundary conditions (optional, default to 0 flux through boundary). + class(BoundaryType), target, intent(in), optional :: & + upper_bndry, lower_bndry + + ! Objects representing boundary conditions. + class(BoundaryCond), intent(in), optional :: l_cond, r_cond + + real(r8) :: solution(ncols,pver) + + ! decomposition. + type(TriDiagDecomp) :: decomp + + ! --------------- ! + ! Local Variables ! + ! --------------- ! + + ! Operator objects. + type(TriDiagOp) :: add_term + type(TriDiagOp) :: net_operator + + ! ----------------------- ! + ! Main Computation Begins ! + ! ----------------------- ! + + ! A diffusion term is probably present, so start with that. Otherwise + ! start with an operator of all 0s. + + if (present(coef_q_diff)) then + net_operator = diffusion_operator(p, coef_q_diff, & + upper_bndry, lower_bndry) + else + net_operator = zero_operator(p%n, p%d) + end if + + ! Constant term (damping). + if (present(coef_q)) then + add_term = diagonal_operator(coef_q) + call net_operator%add(add_term) + end if + + ! Effective advection. + if (present(coef_q_adv)) then + add_term = advection_operator(p, coef_q_adv, & + upper_bndry, lower_bndry) + call net_operator%add(add_term) + end if + + ! We want I-dt*(w^-1)*A for a single time step, implicit method, where + ! A is the right-hand-side operator (i.e. what net_operator is now). + if (present(coef_q_weight)) then + call net_operator%lmult_as_diag(-dt/coef_q_weight) + else + call net_operator%lmult_as_diag(-dt) + end if + call net_operator%add_to_diag(1._r8) + + ! Decompose + decomp = TriDiagDecomp(net_operator) + solution = toSolve + + call net_operator%finalize() + call add_term%finalize() + + call decomp%left_div(solution(:ncols, :), l_cond=l_cond, r_cond=r_cond) + !tendency = tendency - toSolve + + ! Ensure local objects are deallocated. + call decomp%finalize() + + end function fin_vol_solve + +end module vertical_diffusion + \ No newline at end of file From f1e03e42cf85223aa921ab0ea8a71535c6444f96 Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Sun, 6 Oct 2024 20:41:32 -0600 Subject: [PATCH 04/14] Removing shr_ code. --- .../utilities/shr_infnan_mod.F90.in | 406 ------------------ to-be-ccpp-ized/utilities/shr_kind_mod.F90 | 21 - to-be-ccpp-ized/utilities/shr_log_mod.F90 | 121 ------ .../utilities/shr_strconvert_mod.F90 | 167 ------- to-be-ccpp-ized/utilities/shr_sys_mod.F90 | 332 -------------- 5 files changed, 1047 deletions(-) delete mode 100644 to-be-ccpp-ized/utilities/shr_infnan_mod.F90.in delete mode 100644 to-be-ccpp-ized/utilities/shr_kind_mod.F90 delete mode 100644 to-be-ccpp-ized/utilities/shr_log_mod.F90 delete mode 100644 to-be-ccpp-ized/utilities/shr_strconvert_mod.F90 delete mode 100644 to-be-ccpp-ized/utilities/shr_sys_mod.F90 diff --git a/to-be-ccpp-ized/utilities/shr_infnan_mod.F90.in b/to-be-ccpp-ized/utilities/shr_infnan_mod.F90.in deleted file mode 100644 index 992c46fc..00000000 --- a/to-be-ccpp-ized/utilities/shr_infnan_mod.F90.in +++ /dev/null @@ -1,406 +0,0 @@ -! Flag representing compiler support of Fortran 2003's -! ieee_arithmetic intrinsic module. -#if defined CPRIBM || defined CPRPGI || defined CPRINTEL || defined CPRCRAY || defined CPRNAG -#define HAVE_IEEE_ARITHMETIC -#endif - -module shr_infnan_mod -!--------------------------------------------------------------------- -! Module to test for IEEE Inf and NaN values, which also provides a -! method of setting +/-Inf and signaling or quiet NaN. -! -! All functions are elemental, and thus work on arrays. -!--------------------------------------------------------------------- -! To test for these values, just call the corresponding function, e.g: -! -! var_is_nan = shr_infnan_isnan(x) -! -! You can also use it on arrays: -! -! array_contains_nan = any(shr_infnan_isnan(my_array)) -! -!--------------------------------------------------------------------- -! To generate these values, assign one of the provided derived-type -! variables to a real: -! -! use shr_infnan_mod, only: nan => shr_infnan_nan, & -! inf => shr_infnan_inf, & -! assignment(=) -! real(r4) :: my_nan -! real(r8) :: my_inf_array(2,2) -! my_nan = nan -! my_inf_array = inf -! -! Keep in mind that "shr_infnan_nan" and "shr_infnan_inf" cannot be -! passed to functions that expect real arguments. To pass a real -! NaN, you will have to use shr_infnan_nan to set a local real of -! the correct kind. -!--------------------------------------------------------------------- - -use shr_kind_mod, only: & - r4 => SHR_KIND_R4, & - r8 => SHR_KIND_R8 - -#ifdef HAVE_IEEE_ARITHMETIC - -! If we have IEEE_ARITHMETIC, the NaN test is provided for us. -use, intrinsic :: ieee_arithmetic, only: & - shr_infnan_isnan => ieee_is_nan - -#else - -! Integers of correct size for bit patterns below. -use shr_kind_mod, only: i4 => shr_kind_i4, i8 => shr_kind_i8 - -#endif - -implicit none -private -save - -! Test functions for NaN/Inf values. -public :: shr_infnan_isnan -public :: shr_infnan_isinf -public :: shr_infnan_isposinf -public :: shr_infnan_isneginf - -! Locally defined isnan. -#ifndef HAVE_IEEE_ARITHMETIC -interface shr_infnan_isnan - ! TYPE double,real - module procedure shr_infnan_isnan_{TYPE} -end interface -#endif - -interface shr_infnan_isinf - ! TYPE double,real - module procedure shr_infnan_isinf_{TYPE} -end interface - -interface shr_infnan_isposinf - ! TYPE double,real - module procedure shr_infnan_isposinf_{TYPE} -end interface - -interface shr_infnan_isneginf - ! TYPE double,real - module procedure shr_infnan_isneginf_{TYPE} -end interface - -! Derived types for generation of NaN/Inf -! Even though there's no reason to "use" the types directly, some compilers -! might have trouble with an object being used without its type. -public :: shr_infnan_nan_type -public :: shr_infnan_inf_type -public :: assignment(=) -public :: shr_infnan_to_r4 -public :: shr_infnan_to_r8 - -! Type representing Not A Number. -type :: shr_infnan_nan_type - logical :: quiet = .false. -end type shr_infnan_nan_type - -! Type representing +/-Infinity. -type :: shr_infnan_inf_type - logical :: positive = .true. -end type shr_infnan_inf_type - -! Allow assigning reals to NaN or Inf. -interface assignment(=) - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_nan_{DIMS}d_{TYPE} - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_inf_{DIMS}d_{TYPE} -end interface - -! Conversion functions. -interface shr_infnan_to_r8 - module procedure nan_r8 - module procedure inf_r8 -end interface - -interface shr_infnan_to_r4 - module procedure nan_r4 - module procedure inf_r4 -end interface - -! Initialize objects of NaN/Inf type for other modules to use. - -! Default NaN is signaling, but also provide snan and qnan to choose -! explicitly. -type(shr_infnan_nan_type), public, parameter :: shr_infnan_nan = & - shr_infnan_nan_type(.false.) -type(shr_infnan_nan_type), public, parameter :: shr_infnan_snan = & - shr_infnan_nan_type(.false.) -type(shr_infnan_nan_type), public, parameter :: shr_infnan_qnan = & - shr_infnan_nan_type(.true.) - -! Default Inf is positive, but provide posinf to go with neginf. -type(shr_infnan_inf_type), public, parameter :: shr_infnan_inf = & - shr_infnan_inf_type(.true.) -type(shr_infnan_inf_type), public, parameter :: shr_infnan_posinf = & - shr_infnan_inf_type(.true.) -type(shr_infnan_inf_type), public, parameter :: shr_infnan_neginf = & - shr_infnan_inf_type(.false.) - -! Bit patterns for implementation without ieee_arithmetic. -! Note that in order to satisfy gfortran's range check, we have to use -! ibset to set the sign bit from a BOZ pattern. -#ifndef HAVE_IEEE_ARITHMETIC -! Single precision. -integer(i4), parameter :: ssnan_pat = int(Z'7FA00000',i4) -integer(i4), parameter :: sqnan_pat = int(Z'7FC00000',i4) -integer(i4), parameter :: sposinf_pat = int(Z'7F800000',i4) -integer(i4), parameter :: sneginf_pat = ibset(sposinf_pat,bit_size(1_i4)-1) -! Double precision. -integer(i8), parameter :: dsnan_pat = int(Z'7FF4000000000000',i8) -integer(i8), parameter :: dqnan_pat = int(Z'7FF8000000000000',i8) -integer(i8), parameter :: dposinf_pat = int(Z'7FF0000000000000',i8) -integer(i8), parameter :: dneginf_pat = ibset(dposinf_pat,bit_size(1_i8)-1) -#endif - -contains - -!--------------------------------------------------------------------- -! TEST FUNCTIONS -!--------------------------------------------------------------------- -! The "isinf" function simply calls "isposinf" and "isneginf". -!--------------------------------------------------------------------- - -! TYPE double,real -elemental function shr_infnan_isinf_{TYPE}(x) result(isinf) - {VTYPE}, intent(in) :: x - logical :: isinf - - isinf = shr_infnan_isposinf(x) .or. shr_infnan_isneginf(x) - -end function shr_infnan_isinf_{TYPE} - -#ifdef HAVE_IEEE_ARITHMETIC - -!--------------------------------------------------------------------- -! The "isposinf" and "isneginf" functions get the IEEE class of a -! real, and test to see if the class is equal to ieee_positive_inf -! or ieee_negative_inf. -!--------------------------------------------------------------------- - -! TYPE double,real -elemental function shr_infnan_isposinf_{TYPE}(x) result(isposinf) - use, intrinsic :: ieee_arithmetic, only: & - ieee_class, & - ieee_positive_inf, & - operator(==) - {VTYPE}, intent(in) :: x - logical :: isposinf - - isposinf = (ieee_positive_inf == ieee_class(x)) - -end function shr_infnan_isposinf_{TYPE} - -! TYPE double,real -elemental function shr_infnan_isneginf_{TYPE}(x) result(isneginf) - use, intrinsic :: ieee_arithmetic, only: & - ieee_class, & - ieee_negative_inf, & - operator(==) - {VTYPE}, intent(in) :: x - logical :: isneginf - - isneginf = (ieee_negative_inf == ieee_class(x)) - -end function shr_infnan_isneginf_{TYPE} - -#else -! Don't have ieee_arithmetic. - -#ifdef CPRGNU -! NaN testing on gfortran. -! TYPE double,real -elemental function shr_infnan_isnan_{TYPE}(x) result(is_nan) - {VTYPE}, intent(in) :: x - logical :: is_nan - - is_nan = isnan(x) - -end function shr_infnan_isnan_{TYPE} -! End GNU section. -#endif - -!--------------------------------------------------------------------- -! The "isposinf" and "isneginf" functions just test against a known -! bit pattern if we don't have ieee_arithmetic. -!--------------------------------------------------------------------- - -! TYPE double,real -elemental function shr_infnan_isposinf_{TYPE}(x) result(isposinf) - {VTYPE}, intent(in) :: x - logical :: isposinf -#if ({ITYPE} == TYPEREAL) - integer(i4), parameter :: posinf_pat = sposinf_pat -#else - integer(i8), parameter :: posinf_pat = dposinf_pat -#endif - - isposinf = (x == transfer(posinf_pat,x)) - -end function shr_infnan_isposinf_{TYPE} - -! TYPE double,real -elemental function shr_infnan_isneginf_{TYPE}(x) result(isneginf) - {VTYPE}, intent(in) :: x - logical :: isneginf -#if ({ITYPE} == TYPEREAL) - integer(i4), parameter :: neginf_pat = sneginf_pat -#else - integer(i8), parameter :: neginf_pat = dneginf_pat -#endif - - isneginf = (x == transfer(neginf_pat,x)) - -end function shr_infnan_isneginf_{TYPE} - -! End ieee_arithmetic conditional. -#endif - -!--------------------------------------------------------------------- -! GENERATION FUNCTIONS -!--------------------------------------------------------------------- -! Two approaches for generation of NaN and Inf values: -! 1. With Fortran 2003, use the ieee_value intrinsic to get a value -! from the corresponding class. These are: -! - ieee_signaling_nan -! - ieee_quiet_nan -! - ieee_positive_inf -! - ieee_negative_inf -! 2. Without Fortran 2003, set the IEEE bit patterns directly. -! Use BOZ literals to get an integer with the correct bit -! pattern, then use "transfer" to transfer those bits into a -! real. -!--------------------------------------------------------------------- - -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 -pure subroutine set_nan_{DIMS}d_{TYPE}(output, nan) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_signaling_nan, & - ieee_quiet_nan, & - ieee_value -#else -#if ({ITYPE} == TYPEREAL) - integer(i4), parameter :: snan_pat = ssnan_pat - integer(i4), parameter :: qnan_pat = sqnan_pat -#else - integer(i8), parameter :: snan_pat = dsnan_pat - integer(i8), parameter :: qnan_pat = dqnan_pat -#endif -#endif - {VTYPE}, intent(out) :: output{DIMSTR} - type(shr_infnan_nan_type), intent(in) :: nan - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - {VTYPE} :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (nan%quiet) then - tmp = ieee_value(tmp, ieee_quiet_nan) - else - tmp = ieee_value(tmp, ieee_signaling_nan) - end if -#else - if (nan%quiet) then - tmp = transfer(qnan_pat, tmp) - else - tmp = transfer(snan_pat, tmp) - end if -#endif - - output = tmp - -end subroutine set_nan_{DIMS}d_{TYPE} - -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 -pure subroutine set_inf_{DIMS}d_{TYPE}(output, inf) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_positive_inf, & - ieee_negative_inf, & - ieee_value -#else -#if ({ITYPE} == TYPEREAL) - integer(i4), parameter :: posinf_pat = sposinf_pat - integer(i4), parameter :: neginf_pat = sneginf_pat -#else - integer(i8), parameter :: posinf_pat = dposinf_pat - integer(i8), parameter :: neginf_pat = dneginf_pat -#endif -#endif - {VTYPE}, intent(out) :: output{DIMSTR} - type(shr_infnan_inf_type), intent(in) :: inf - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - {VTYPE} :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (inf%positive) then - tmp = ieee_value(tmp,ieee_positive_inf) - else - tmp = ieee_value(tmp,ieee_negative_inf) - end if -#else - if (inf%positive) then - tmp = transfer(posinf_pat, tmp) - else - tmp = transfer(neginf_pat, tmp) - end if -#endif - - output = tmp - -end subroutine set_inf_{DIMS}d_{TYPE} - -!--------------------------------------------------------------------- -! CONVERSION INTERFACES. -!--------------------------------------------------------------------- -! Function methods to get reals from nan/inf types. -!--------------------------------------------------------------------- - -pure function nan_r8(nan) result(output) - class(shr_infnan_nan_type), intent(in) :: nan - real(r8) :: output - - output = nan - -end function nan_r8 - -pure function nan_r4(nan) result(output) - class(shr_infnan_nan_type), intent(in) :: nan - real(r4) :: output - - output = nan - -end function nan_r4 - -pure function inf_r8(inf) result(output) - class(shr_infnan_inf_type), intent(in) :: inf - real(r8) :: output - - output = inf - -end function inf_r8 - -pure function inf_r4(inf) result(output) - class(shr_infnan_inf_type), intent(in) :: inf - real(r4) :: output - - output = inf - -end function inf_r4 - -end module shr_infnan_mod diff --git a/to-be-ccpp-ized/utilities/shr_kind_mod.F90 b/to-be-ccpp-ized/utilities/shr_kind_mod.F90 deleted file mode 100644 index 9ab62a4e..00000000 --- a/to-be-ccpp-ized/utilities/shr_kind_mod.F90 +++ /dev/null @@ -1,21 +0,0 @@ -MODULE shr_kind_mod - - !---------------------------------------------------------------------------- - ! precision/kind constants add data public - !---------------------------------------------------------------------------- - public - integer,parameter :: SHR_KIND_R8 = selected_real_kind(12) ! 8 byte real - integer,parameter :: SHR_KIND_R4 = selected_real_kind( 6) ! 4 byte real - integer,parameter :: SHR_KIND_RN = kind(1.0) ! native real - integer,parameter :: SHR_KIND_I8 = selected_int_kind (13) ! 8 byte integer - integer,parameter :: SHR_KIND_I4 = selected_int_kind ( 6) ! 4 byte integer - integer,parameter :: SHR_KIND_I2 = selected_int_kind ( 4) ! 2 byte integer - integer,parameter :: SHR_KIND_IN = kind(1) ! native integer - integer,parameter :: SHR_KIND_CS = 80 ! short char - integer,parameter :: SHR_KIND_CM = 160 ! mid-sized char - integer,parameter :: SHR_KIND_CL = 256 ! long char - integer,parameter :: SHR_KIND_CX = 512 ! extra-long char - integer,parameter :: SHR_KIND_CXX= 4096 ! extra-extra-long char - -END MODULE shr_kind_mod - \ No newline at end of file diff --git a/to-be-ccpp-ized/utilities/shr_log_mod.F90 b/to-be-ccpp-ized/utilities/shr_log_mod.F90 deleted file mode 100644 index e028f497..00000000 --- a/to-be-ccpp-ized/utilities/shr_log_mod.F90 +++ /dev/null @@ -1,121 +0,0 @@ -!BOP =========================================================================== -! -! !MODULE: shr_log_mod -- variables and methods for logging -! -! !DESCRIPTION: -! Low-level shared variables for logging. -! -! Also, routines for generating log file messages. -! -! !INTERFACE: ------------------------------------------------------------------ - -module shr_log_mod - - ! !USES: - - use shr_kind_mod, only: shr_kind_in, shr_kind_cx - use shr_strconvert_mod, only: toString - - use, intrinsic :: iso_fortran_env, only: output_unit - - implicit none - private - - ! !PUBLIC TYPES: - - ! no public types - - ! !PUBLIC MEMBER FUNCTIONS: - - public :: shr_log_errMsg - public :: shr_log_OOBMsg - public :: shr_log_setLogUnit - public :: shr_log_getLogUnit - - ! !PUBLIC DATA MEMBERS: - - public :: shr_log_Level - public :: shr_log_Unit - - !EOP - - ! low-level shared variables for logging, these may not be parameters - integer(SHR_KIND_IN) :: shr_log_Level = 0 - integer(SHR_KIND_IN) :: shr_log_Unit = output_unit - - contains - - !=============================================================================== - !BOP =========================================================================== - ! - ! !IROUTINE: shr_log_errMsg -- Return an error message containing file & line info - ! - ! !DESCRIPTION: - ! Return an error message containing file & line info - ! \newline - ! errMsg = shr\_log\_errMsg(__FILE__, __LINE__) - ! - ! This is meant to be used when a routine expects a string argument for some message, - ! but you want to provide file and line information. - ! - ! However: Note that the performance of this function can be very bad. It is currently - ! maintained because it is used by old code, but you should probably avoid using this - ! in new code if possible. - ! - ! !REVISION HISTORY: - ! 2013-July-23 - Bill Sacks - ! - ! !INTERFACE: ------------------------------------------------------------------ - - pure function shr_log_errMsg(file, line) - - ! !INPUT/OUTPUT PARAMETERS: - - character(len=SHR_KIND_CX) :: shr_log_errMsg - character(len=*), intent(in) :: file - integer , intent(in) :: line - - !EOP - - shr_log_errMsg = 'ERROR in '//trim(file)//' at line '//toString(line) - - end function shr_log_errMsg - - ! Create a message for an out of bounds error. - pure function shr_log_OOBMsg(operation, bounds, idx) result(OOBMsg) - - ! A name for the operation being attempted when the bounds error - ! occurred. A string containing the subroutine name is ideal, but more - ! generic descriptions such as "read", "modify", or "insert" could be used. - character(len=*), intent(in) :: operation - - ! Upper and lower bounds allowed for the operation. - integer, intent(in) :: bounds(2) - - ! Index at which access was attempted. - integer, intent(in) :: idx - - ! Output message - character(len=:), allocatable :: OOBMsg - - allocate(OOBMsg, source=(operation//": "//toString(idx)//" not in range ["//& - toString(bounds(1))//", "//toString(bounds(2))//"].")) - - end function shr_log_OOBMsg - - subroutine shr_log_setLogUnit(unit) - integer, intent(in) :: unit - - shr_log_unit = unit - - end subroutine shr_log_setLogUnit - - subroutine shr_log_getLogUnit(unit) - integer, intent(out) :: unit - - unit = shr_log_unit - - end subroutine shr_log_getLogUnit - - end module shr_log_mod - \ No newline at end of file diff --git a/to-be-ccpp-ized/utilities/shr_strconvert_mod.F90 b/to-be-ccpp-ized/utilities/shr_strconvert_mod.F90 deleted file mode 100644 index 56ec172a..00000000 --- a/to-be-ccpp-ized/utilities/shr_strconvert_mod.F90 +++ /dev/null @@ -1,167 +0,0 @@ -module shr_strconvert_mod - - ! This module defines toString, a generic function for creating character type - ! representations of data, as implemented for the most commonly used intrinsic - ! types: - ! - ! - 4 and 8 byte integer - ! - 4 and 8 byte real - ! - logical - ! - ! No toString implementation is provided for character input, but this may be - ! added if some use case arises. - ! - ! Currently, only scalar inputs are supported. The return type of this function - ! is character with deferred (allocatable) length. - ! - ! The functions for integers and reals allow an optional format_string argument, - ! which can be used to control the padding and precision of output as with any - ! write statement. However, the implementations internally must use a - ! preallocated buffer, so a format_string that significantly increases the size - ! of the output may cause a run-time error or undefined behavior in the program. - ! - ! Other modules may want to provide extensions of toString for their own derived - ! types. In this case there are two guidelines to observe: - ! - ! - It is preferable to have only one mandatory argument, which is the object to - ! produce a string from. There may be other formatting options, but the - ! implementation should do something sensible without these. - ! - ! - Since the main purpose of toString is to provide a human-readable - ! representation of a type, especially for documentation or debugging - ! purposes, refrain from printing large array components in their entirety - ! (instead consider printing only the shape, or statistics such as - ! min/mean/max for arrays of numbers). - - use shr_kind_mod, only: & - i4 => shr_kind_i4, & - i8 => shr_kind_i8, & - r4 => shr_kind_r4, & - r8 => shr_kind_r8, & - cs => shr_kind_cs - - use shr_infnan_mod, only: & - isnan => shr_infnan_isnan - - implicit none - private - - ! Human-readable representation of data. - public :: toString - - interface toString - module procedure i4ToString - module procedure i8ToString - module procedure r4ToString - module procedure r8ToString - module procedure logicalToString - end interface toString - - contains - - pure function i4ToString(input, format_string) result(string) - integer(i4), intent(in) :: input - character(len=*), intent(in), optional :: format_string - character(len=:), allocatable :: string - - character(len=cs) :: buffer - - if (present(format_string)) then - write(buffer, format_string) input - else - ! For most compilers, these two statements are equivalent to a format of - ! '(I0)', but that's not technically in the standard. - write(buffer, '(I11)') input - buffer = adjustl(buffer) - end if - - allocate(string, source=trim(buffer)) - - end function i4ToString - - pure function i8ToString(input, format_string) result(string) - integer(i8), intent(in) :: input - character(len=*), intent(in), optional :: format_string - character(len=:), allocatable :: string - - character(len=cs) :: buffer - - if (present(format_string)) then - write(buffer, format_string) input - else - ! For most compilers, these two statements are equivalent to a format of - ! '(I0)', but that's not technically in the standard. - write(buffer, '(I20)') input - buffer = adjustl(buffer) - end if - - allocate(string, source=trim(buffer)) - - end function i8ToString - - pure function r4ToString(input, format_string) result(string) - real(r4), intent(in) :: input - character(len=*), intent(in), optional :: format_string - character(len=:), allocatable :: string - - character(len=cs) :: buffer - - if (present(format_string)) then - write(buffer, format_string) input - else - write(buffer, '(ES15.8 E2)') input - buffer = adjustl(buffer) - ! Deal with the fact that the "+" sign is optional by simply adding it if - ! it is not present, so that the default format is standardized across - ! compilers. - ! Assumes that compilers do not treat the sign bit on NaN values specially. - if (.not. isnan(input) .and. all(buffer(1:1) /= ["-", "+"])) then - buffer = "+" // trim(buffer) - end if - end if - - allocate(string, source=trim(buffer)) - - end function r4ToString - - pure function r8ToString(input, format_string) result(string) - real(r8), intent(in) :: input - character(len=*), intent(in), optional :: format_string - character(len=:), allocatable :: string - - character(len=cs) :: buffer - - if (present(format_string)) then - write(buffer, format_string) input - else - write(buffer, '(ES24.16 E3)') input - buffer = adjustl(buffer) - ! Deal with the fact that the "+" sign is optional by simply adding it if - ! it is not present, so that the default format is standardized across - ! compilers. - ! Assumes that compilers do not treat the sign bit on NaN values specially. - if (.not. isnan(input) .and. all(buffer(1:1) /= ["-", "+"])) then - buffer = "+" // trim(buffer) - end if - end if - - allocate(string, source=trim(buffer)) - - end function r8ToString - - pure function logicalToString(input) result(string) - logical, intent(in) :: input - character(len=:), allocatable :: string - - ! We could use a write statement, but this is easier. - allocate(character(len=1) :: string) - if (input) then - string = "T" - else - string = "F" - end if - - end function logicalToString - -end module shr_strconvert_mod - \ No newline at end of file diff --git a/to-be-ccpp-ized/utilities/shr_sys_mod.F90 b/to-be-ccpp-ized/utilities/shr_sys_mod.F90 deleted file mode 100644 index bf463d55..00000000 --- a/to-be-ccpp-ized/utilities/shr_sys_mod.F90 +++ /dev/null @@ -1,332 +0,0 @@ -!=============================================================================== -! SVN $Id: shr_sys_mod.F90 66411 2014-12-19 22:40:08Z santos@ucar.edu $ -! SVN $URL: https://svn-ccsm-models.cgd.ucar.edu/csm_share/trunk_tags/share3_150116/shr/shr_sys_mod.F90 $ -!=============================================================================== - -! Currently supported by all compilers -#define HAVE_GET_ENVIRONMENT -#define HAVE_SLEEP - -! Except this combination? -#if defined CPRPGI && defined CNL -#undef HAVE_GET_ENVIRONMENT -#endif - -#if defined CPRNAG -#define HAVE_EXECUTE -#endif - -MODULE shr_sys_mod - - use shr_kind_mod ! defines real & integer kinds - use shr_log_mod, only: s_loglev => shr_log_Level - use shr_log_mod, only: s_logunit => shr_log_Unit - use shr_abort_mod, only: shr_sys_abort => shr_abort_abort - use shr_abort_mod, only: shr_sys_backtrace => shr_abort_backtrace - -#ifdef CPRNAG - ! NAG does not provide these as intrinsics, but it does provide modules - ! that implement commonly used POSIX routines. - use f90_unix_dir, only: chdir - use f90_unix_proc, only: abort, sleep -#endif - - implicit none - -! PUBLIC: Public interfaces - - private - - public :: shr_sys_system ! make a system call - public :: shr_sys_chdir ! change current working dir - public :: shr_sys_getenv ! get an environment variable - public :: shr_sys_irtc ! returns real-time clock tick - public :: shr_sys_sleep ! have program sleep for a while - public :: shr_sys_flush ! flush an i/o buffer - - ! Imported from shr_abort_mod and republished with renames. Other code that wishes to - ! use these routines should use these shr_sys names rather than directly using the - ! routines from shr_abort_abort. (This is for consistency with older code, from when - ! these routines were defined in shr_sys_mod.) - public :: shr_sys_abort ! abort a program - public :: shr_sys_backtrace ! print a backtrace, if possible - -!=============================================================================== -CONTAINS -!=============================================================================== - -!=============================================================================== -!=============================================================================== - -SUBROUTINE shr_sys_system(str,rcode) - - IMPLICIT none - - !----- arguments --- - character(*) ,intent(in) :: str ! system/shell command string - integer(SHR_KIND_IN),intent(out) :: rcode ! function return error code - - !----- functions ----- -#if (defined LINUX && !defined CPRGNU) - integer(SHR_KIND_IN),external :: system ! function to envoke shell command -#endif - - !----- formats ----- - character(*),parameter :: subName = '(shr_sys_system) ' - character(*),parameter :: F00 = "('(shr_sys_system) ',4a)" - -!------------------------------------------------------------------------------- -! PURPOSE: an architecture independent system call -!------------------------------------------------------------------------------- - rcode = 0 -#ifdef HAVE_EXECUTE - call execute_command_line(str,exitstat=rcode) ! Intrinsic as of F2008 -#else -#if (defined AIX) - - call system(str,rcode) - -#elif (defined CPRGNU || defined LINUX) - - rcode = system(str) - -#else - - write(s_logunit,F00) 'ERROR: no implementation of system call for this architecture' - call shr_sys_abort(subName//'no implementation of system call for this architecture') -#endif -#endif - -END SUBROUTINE shr_sys_system - -!=============================================================================== -!=============================================================================== - -SUBROUTINE shr_sys_chdir(path, rcode) - - IMPLICIT none - - !----- arguments ----- - character(*) ,intent(in) :: path ! chdir to this dir - integer(SHR_KIND_IN),intent(out) :: rcode ! return code - - !----- local ----- - integer(SHR_KIND_IN) :: lenpath ! length of path -#if (defined AIX || (defined LINUX && !defined CPRGNU && !defined CPRNAG) || defined CPRINTEL) - integer(SHR_KIND_IN),external :: chdir ! AIX system call -#endif - - !----- formats ----- - character(*),parameter :: subName = '(shr_sys_chdir) ' - character(*),parameter :: F00 = "('(shr_sys_chdir) ',4a)" - -!------------------------------------------------------------------------------- -! PURPOSE: an architecture independent system call -!------------------------------------------------------------------------------- - rcode = 0 - lenpath=len_trim(path) - -#if (defined AIX) - - rcode = chdir(%ref(path(1:lenpath)//'\0')) - -#elif (defined Darwin || (defined LINUX && !defined CPRNAG)) - - rcode=chdir(path(1:lenpath)) - -#elif (defined CPRNAG) - - call chdir(path(1:lenpath), errno=rcode) - -#else - - write(s_logunit,F00) 'ERROR: no implementation of chdir for this architecture' - call shr_sys_abort(subname//'no implementation of chdir for this machine') - -#endif - -END SUBROUTINE shr_sys_chdir - -!=============================================================================== -!=============================================================================== - -SUBROUTINE shr_sys_getenv(name, val, rcode) - - IMPLICIT none - - !----- arguments ----- - character(*) ,intent(in) :: name ! env var name - character(*) ,intent(out) :: val ! env var value - integer(SHR_KIND_IN),intent(out) :: rcode ! return code - - !----- local ----- -#ifndef HAVE_GET_ENVIRONMENT - integer(SHR_KIND_IN) :: lenname ! length of env var name - integer(SHR_KIND_IN) :: lenval ! length of env var value - character(SHR_KIND_CL) :: tmpval ! temporary env var value -#endif - !----- formats ----- - character(*),parameter :: subName = '(shr_sys_getenv) ' - character(*),parameter :: F00 = "('(shr_sys_getenv) ',4a)" - -!------------------------------------------------------------------------------- -! PURPOSE: an architecture independent system call -!------------------------------------------------------------------------------- - -!$OMP master - - -#ifdef HAVE_GET_ENVIRONMENT - call get_environment_variable(name=name,value=val,status=rcode) ! Intrinsic in F2003 -#else - lenname=len_trim(name) -#if (defined AIX || defined LINUX) - - call getenv(trim(name),tmpval) - val=trim(tmpval) - rcode = 0 - if (len_trim(val) == 0 ) rcode = 1 - if (len_trim(val) > SHR_KIND_CL) rcode = 2 - -#else - - write(s_logunit,F00) 'ERROR: no implementation of getenv for this architecture' - call shr_sys_abort(subname//'no implementation of getenv for this machine') - -#endif -#endif -!$OMP end master - -END SUBROUTINE shr_sys_getenv - -!=============================================================================== -!=============================================================================== - -integer(SHR_KIND_I8) FUNCTION shr_sys_irtc( rate ) - - IMPLICIT none - - !----- arguments ----- - integer(SHR_KIND_I8), optional :: rate - - !----- local ----- - integer(SHR_KIND_IN) :: count - integer(SHR_KIND_IN) :: count_rate - integer(SHR_KIND_IN) :: count_max - - integer(SHR_KIND_IN),save :: last_count = -1 - integer(SHR_KIND_I8),save :: count_offset = 0 -!$OMP THREADPRIVATE (last_count, count_offset) - - !----- formats ----- - character(*),parameter :: subName = '(shr_sys_irtc) ' - character(*),parameter :: F00 = "('(shr_sys_irtc) ',4a)" - -!------------------------------------------------------------------------------- -! emulates Cray/SGI irtc function (returns clock tick since last reboot) -! -! This function is not intended to measure elapsed time between -! multi-threaded regions with different numbers of threads. However, -! use of the threadprivate declaration does guarantee accurate -! measurement per thread within a single multi-threaded region as -! long as the number of threads is not changed dynamically during -! execution within the multi-threaded region. -! -!------------------------------------------------------------------------------- - - call system_clock(count=count,count_rate=count_rate, count_max=count_max) - if ( present(rate) ) rate = count_rate - shr_sys_irtc = count - - !--- adjust for clock wrap-around --- - if ( last_count /= -1 ) then - if ( count < last_count ) count_offset = count_offset + count_max - end if - shr_sys_irtc = shr_sys_irtc + count_offset - last_count = count - -END FUNCTION shr_sys_irtc - -!=============================================================================== -!=============================================================================== - -SUBROUTINE shr_sys_sleep(sec) - - IMPLICIT none - - !----- arguments ----- - real (SHR_KIND_R8),intent(in) :: sec ! number of seconds to sleep - - !----- local ----- - integer(SHR_KIND_IN) :: isec ! integer number of seconds -#ifndef HAVE_SLEEP - integer(SHR_KIND_IN) :: rcode ! return code - character(90) :: str ! system call string -#endif - !----- formats ----- - character(*),parameter :: subName = '(shr_sys_sleep) ' - character(*),parameter :: F00 = "('(shr_sys_sleep) ',4a)" - character(*),parameter :: F10 = "('sleep ',i8 )" - -!------------------------------------------------------------------------------- -! PURPOSE: Sleep for approximately sec seconds -!------------------------------------------------------------------------------- - - isec = nint(sec) - - if (isec < 0) then - if (s_loglev > 0) write(s_logunit,F00) 'ERROR: seconds must be > 0, sec=',sec - else if (isec == 0) then - ! Don't consider this an error and don't call system sleep - else -#ifdef HAVE_SLEEP - call sleep(isec) -#else - write(str,FMT=F10) isec - call shr_sys_system( str, rcode ) -#endif - endif - -END SUBROUTINE shr_sys_sleep - -!=============================================================================== -!=============================================================================== - -SUBROUTINE shr_sys_flush(unit) - - IMPLICIT none - - !----- arguments ----- - integer(SHR_KIND_IN) :: unit ! flush output buffer for this unit - - !----- local ----- - !----- formats ----- - character(*),parameter :: subName = '(shr_sys_flush) ' - character(*),parameter :: F00 = "('(shr_sys_flush) ',4a)" - -!------------------------------------------------------------------------------- -! PURPOSE: an architecture independent system call -! -! This is probably no longer needed; the "flush" statement is supported by -! all compilers that CESM supports for years now. -! -!------------------------------------------------------------------------------- -!$OMP SINGLE - flush(unit) -!$OMP END SINGLE -! -! The following code was originally present, but there's an obvious issue. -! Since shr_sys_flush is usually used to flush output to a log, when it -! returns an error, does it do any good to print that error to the log? -! -! if (ierr > 0) then -! write(s_logunit,*) subname,' Flush reports error: ',ierr -! endif -! - -END SUBROUTINE shr_sys_flush - -!=============================================================================== -!=============================================================================== - -END MODULE shr_sys_mod From f4231db52fcf266d82acad31014526ba9216a509 Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Sun, 6 Oct 2024 21:20:41 -0600 Subject: [PATCH 05/14] Updating module names to avoid conflicts. --- .../utilities/{coords_1d.F90 => atmos_phys_coords_1d.F90} | 4 ++-- ...d_operators.F90 => atmos_phys_linear_1d_operators.F90} | 2 +- ...rtical_diffusion.F90 => vertical_diffusion_solver.F90} | 8 ++++---- 3 files changed, 7 insertions(+), 7 deletions(-) rename to-be-ccpp-ized/utilities/{coords_1d.F90 => atmos_phys_coords_1d.F90} (98%) rename to-be-ccpp-ized/utilities/{linear_1d_operators.F90 => atmos_phys_linear_1d_operators.F90} (99%) rename to-be-ccpp-ized/utilities/vertical_diffusion/{vertical_diffusion.F90 => vertical_diffusion_solver.F90} (96%) diff --git a/to-be-ccpp-ized/utilities/coords_1d.F90 b/to-be-ccpp-ized/utilities/atmos_phys_coords_1d.F90 similarity index 98% rename from to-be-ccpp-ized/utilities/coords_1d.F90 rename to to-be-ccpp-ized/utilities/atmos_phys_coords_1d.F90 index 21475de0..30159518 100644 --- a/to-be-ccpp-ized/utilities/coords_1d.F90 +++ b/to-be-ccpp-ized/utilities/atmos_phys_coords_1d.F90 @@ -1,4 +1,4 @@ -module coords_1d +module atmos_phys_coords_1d ! This module defines the Coords1D type, which is intended to to cache ! commonly used information derived from a collection of sets of 1-D @@ -148,5 +148,5 @@ end subroutine guarded_deallocate end subroutine finalize -end module coords_1d +end module atmos_phys_coords_1d \ No newline at end of file diff --git a/to-be-ccpp-ized/utilities/linear_1d_operators.F90 b/to-be-ccpp-ized/utilities/atmos_phys_linear_1d_operators.F90 similarity index 99% rename from to-be-ccpp-ized/utilities/linear_1d_operators.F90 rename to to-be-ccpp-ized/utilities/atmos_phys_linear_1d_operators.F90 index 901200f4..abc1dad8 100644 --- a/to-be-ccpp-ized/utilities/linear_1d_operators.F90 +++ b/to-be-ccpp-ized/utilities/atmos_phys_linear_1d_operators.F90 @@ -41,7 +41,7 @@ module linear_1d_operators use shr_kind_mod, only: r8 => shr_kind_r8 use shr_log_mod, only: errMsg => shr_log_errMsg use shr_sys_mod, only: shr_sys_abort - use coords_1d, only: Coords1D + use atmos_phys_coords_1d, only: Coords1D implicit none private diff --git a/to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion.F90 b/to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 similarity index 96% rename from to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion.F90 rename to to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 index 4e9d87d5..145e5971 100644 --- a/to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion.F90 +++ b/to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 @@ -1,4 +1,4 @@ -module vertical_diffusion +module vertical_diffusion_solver implicit none private @@ -25,7 +25,7 @@ module vertical_diffusion function fin_vol_solve(dt, p, toSolve, ncols, pver, coef_q, coef_q_diff, coef_q_adv, & coef_q_weight, upper_bndry, lower_bndry, l_cond, r_cond) result(solution) - use linear_1d_operators, only: & + use atmos_phys_linear_1d_operators, only: & zero_operator, & diagonal_operator, & diffusion_operator, & @@ -36,7 +36,7 @@ function fin_vol_solve(dt, p, toSolve, ncols, pver, coef_q, coef_q_diff, coef_q_ BoundaryCond, & operator(+) use shr_kind_mod, only: r8 => shr_kind_r8 - use coords_1d, only: Coords1D + use atmos_phys_coords_1d, only: Coords1D ! ---------------------- ! ! Input-Output Arguments ! @@ -132,5 +132,5 @@ function fin_vol_solve(dt, p, toSolve, ncols, pver, coef_q, coef_q_diff, coef_q_ end function fin_vol_solve -end module vertical_diffusion +end module vertical_diffusion_solver \ No newline at end of file From 88f1dbbecafa9cb10be3c28a26c2f183487e134a Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Mon, 7 Oct 2024 09:00:52 -0600 Subject: [PATCH 06/14] Moving modules back to original names and removing host model copies instead. --- .../utilities/{atmos_phys_coords_1d.F90 => coords_1d.F90} | 4 ++-- ...s_phys_linear_1d_operators.F90 => linear_1d_operators.F90} | 0 .../vertical_diffusion/vertical_diffusion_solver.F90 | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) rename to-be-ccpp-ized/utilities/{atmos_phys_coords_1d.F90 => coords_1d.F90} (98%) rename to-be-ccpp-ized/utilities/{atmos_phys_linear_1d_operators.F90 => linear_1d_operators.F90} (100%) diff --git a/to-be-ccpp-ized/utilities/atmos_phys_coords_1d.F90 b/to-be-ccpp-ized/utilities/coords_1d.F90 similarity index 98% rename from to-be-ccpp-ized/utilities/atmos_phys_coords_1d.F90 rename to to-be-ccpp-ized/utilities/coords_1d.F90 index 30159518..21475de0 100644 --- a/to-be-ccpp-ized/utilities/atmos_phys_coords_1d.F90 +++ b/to-be-ccpp-ized/utilities/coords_1d.F90 @@ -1,4 +1,4 @@ -module atmos_phys_coords_1d +module coords_1d ! This module defines the Coords1D type, which is intended to to cache ! commonly used information derived from a collection of sets of 1-D @@ -148,5 +148,5 @@ end subroutine guarded_deallocate end subroutine finalize -end module atmos_phys_coords_1d +end module coords_1d \ No newline at end of file diff --git a/to-be-ccpp-ized/utilities/atmos_phys_linear_1d_operators.F90 b/to-be-ccpp-ized/utilities/linear_1d_operators.F90 similarity index 100% rename from to-be-ccpp-ized/utilities/atmos_phys_linear_1d_operators.F90 rename to to-be-ccpp-ized/utilities/linear_1d_operators.F90 diff --git a/to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 b/to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 index 145e5971..7b1023c2 100644 --- a/to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 +++ b/to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 @@ -36,7 +36,7 @@ function fin_vol_solve(dt, p, toSolve, ncols, pver, coef_q, coef_q_diff, coef_q_ BoundaryCond, & operator(+) use shr_kind_mod, only: r8 => shr_kind_r8 - use atmos_phys_coords_1d, only: Coords1D + use coords_1d, only: Coords1D ! ---------------------- ! ! Input-Output Arguments ! From 55db37b8abf02eff89f7b992fcad78653c156fbb Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Mon, 7 Oct 2024 10:49:57 -0600 Subject: [PATCH 07/14] Removing atmos_phys module prefix missed in previous commit. --- to-be-ccpp-ized/utilities/linear_1d_operators.F90 | 2 +- .../utilities/vertical_diffusion/vertical_diffusion_solver.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/to-be-ccpp-ized/utilities/linear_1d_operators.F90 b/to-be-ccpp-ized/utilities/linear_1d_operators.F90 index abc1dad8..901200f4 100644 --- a/to-be-ccpp-ized/utilities/linear_1d_operators.F90 +++ b/to-be-ccpp-ized/utilities/linear_1d_operators.F90 @@ -41,7 +41,7 @@ module linear_1d_operators use shr_kind_mod, only: r8 => shr_kind_r8 use shr_log_mod, only: errMsg => shr_log_errMsg use shr_sys_mod, only: shr_sys_abort - use atmos_phys_coords_1d, only: Coords1D + use coords_1d, only: Coords1D implicit none private diff --git a/to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 b/to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 index 7b1023c2..62b0495e 100644 --- a/to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 +++ b/to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 @@ -25,7 +25,7 @@ module vertical_diffusion_solver function fin_vol_solve(dt, p, toSolve, ncols, pver, coef_q, coef_q_diff, coef_q_adv, & coef_q_weight, upper_bndry, lower_bndry, l_cond, r_cond) result(solution) - use atmos_phys_linear_1d_operators, only: & + use linear_1d_operators, only: & zero_operator, & diagonal_operator, & diffusion_operator, & From 3d3fa2ae279c31fc53331b86006d9acc0875546b Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Tue, 15 Oct 2024 06:05:38 -0600 Subject: [PATCH 08/14] Removing un-needed save --- .../utilities/vertical_diffusion/vertical_diffusion_solver.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 b/to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 index 62b0495e..23641187 100644 --- a/to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 +++ b/to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 @@ -2,7 +2,6 @@ module vertical_diffusion_solver implicit none private - save public :: fin_vol_solve From ab9d1cc990ab1ab5e3f44148a934ca174cfbd42a Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Tue, 15 Oct 2024 06:20:47 -0600 Subject: [PATCH 09/14] Consolidating finalize calls for deallocating local objects. --- .../vertical_diffusion/vertical_diffusion_solver.F90 | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 b/to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 index 23641187..d4ce847b 100644 --- a/to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 +++ b/to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 @@ -120,14 +120,13 @@ function fin_vol_solve(dt, p, toSolve, ncols, pver, coef_q, coef_q_diff, coef_q_ decomp = TriDiagDecomp(net_operator) solution = toSolve - call net_operator%finalize() - call add_term%finalize() - call decomp%left_div(solution(:ncols, :), l_cond=l_cond, r_cond=r_cond) !tendency = tendency - toSolve ! Ensure local objects are deallocated. call decomp%finalize() + call net_operator%finalize() + call add_term%finalize() end function fin_vol_solve From 9ba5be1eb17b4f57a86dd5be0848854f7216e526 Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Mon, 21 Oct 2024 06:21:07 -0600 Subject: [PATCH 10/14] Moving code to be converted to new directory name. --- {to-be-ccpp-ized => to_be_ccppized}/utilities/coords_1d.F90 | 0 .../utilities/linear_1d_operators.F90 | 0 .../utilities/vertical_diffusion/vertical_diffusion_solver.F90 | 0 3 files changed, 0 insertions(+), 0 deletions(-) rename {to-be-ccpp-ized => to_be_ccppized}/utilities/coords_1d.F90 (100%) rename {to-be-ccpp-ized => to_be_ccppized}/utilities/linear_1d_operators.F90 (100%) rename {to-be-ccpp-ized => to_be_ccppized}/utilities/vertical_diffusion/vertical_diffusion_solver.F90 (100%) diff --git a/to-be-ccpp-ized/utilities/coords_1d.F90 b/to_be_ccppized/utilities/coords_1d.F90 similarity index 100% rename from to-be-ccpp-ized/utilities/coords_1d.F90 rename to to_be_ccppized/utilities/coords_1d.F90 diff --git a/to-be-ccpp-ized/utilities/linear_1d_operators.F90 b/to_be_ccppized/utilities/linear_1d_operators.F90 similarity index 100% rename from to-be-ccpp-ized/utilities/linear_1d_operators.F90 rename to to_be_ccppized/utilities/linear_1d_operators.F90 diff --git a/to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 b/to_be_ccppized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 similarity index 100% rename from to-be-ccpp-ized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 rename to to_be_ccppized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 From 8b1bd8676e3c7b8fb7547b4df1270320d74f6aab Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Mon, 21 Oct 2024 06:23:35 -0600 Subject: [PATCH 11/14] Moving vertical diffusion to similar folder as CAM. --- .../cam}/vertical_diffusion_solver.F90 | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename to_be_ccppized/{utilities/vertical_diffusion => physics/cam}/vertical_diffusion_solver.F90 (100%) diff --git a/to_be_ccppized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 b/to_be_ccppized/physics/cam/vertical_diffusion_solver.F90 similarity index 100% rename from to_be_ccppized/utilities/vertical_diffusion/vertical_diffusion_solver.F90 rename to to_be_ccppized/physics/cam/vertical_diffusion_solver.F90 From a1e05a5e8e30d36d2b36866a84568a791c37d679 Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Mon, 21 Oct 2024 16:59:18 -0600 Subject: [PATCH 12/14] Reorganizing files per PR comments. --- to_be_ccppized/{utilities => }/coords_1d.F90 | 0 to_be_ccppized/{utilities => }/linear_1d_operators.F90 | 0 to_be_ccppized/{physics/cam => }/vertical_diffusion_solver.F90 | 0 3 files changed, 0 insertions(+), 0 deletions(-) rename to_be_ccppized/{utilities => }/coords_1d.F90 (100%) rename to_be_ccppized/{utilities => }/linear_1d_operators.F90 (100%) rename to_be_ccppized/{physics/cam => }/vertical_diffusion_solver.F90 (100%) diff --git a/to_be_ccppized/utilities/coords_1d.F90 b/to_be_ccppized/coords_1d.F90 similarity index 100% rename from to_be_ccppized/utilities/coords_1d.F90 rename to to_be_ccppized/coords_1d.F90 diff --git a/to_be_ccppized/utilities/linear_1d_operators.F90 b/to_be_ccppized/linear_1d_operators.F90 similarity index 100% rename from to_be_ccppized/utilities/linear_1d_operators.F90 rename to to_be_ccppized/linear_1d_operators.F90 diff --git a/to_be_ccppized/physics/cam/vertical_diffusion_solver.F90 b/to_be_ccppized/vertical_diffusion_solver.F90 similarity index 100% rename from to_be_ccppized/physics/cam/vertical_diffusion_solver.F90 rename to to_be_ccppized/vertical_diffusion_solver.F90 From 396aaaca743ae6c6106f0a7e54133d34c819c804 Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Wed, 23 Oct 2024 07:03:04 -0600 Subject: [PATCH 13/14] Addressing PR comments. --- to_be_ccppized/vertical_diffusion_solver.F90 | 24 ++++++++------------ 1 file changed, 9 insertions(+), 15 deletions(-) diff --git a/to_be_ccppized/vertical_diffusion_solver.F90 b/to_be_ccppized/vertical_diffusion_solver.F90 index d4ce847b..1efdca06 100644 --- a/to_be_ccppized/vertical_diffusion_solver.F90 +++ b/to_be_ccppized/vertical_diffusion_solver.F90 @@ -46,8 +46,6 @@ function fin_vol_solve(dt, p, toSolve, ncols, pver, coef_q, coef_q_diff, coef_q_ ! Grid spacings. type(Coords1D), intent(in) :: p - ! Matrix to decomp from. - ! real(r8), intent(in) :: u(ncols,pver) integer, intent(in) :: ncols integer, intent(in) :: pver real(r8), intent(in) :: toSolve(ncols,pver) @@ -61,8 +59,7 @@ function fin_vol_solve(dt, p, toSolve, ncols, pver, coef_q, coef_q_diff, coef_q_ coef_q_diff(:,:), coef_q_adv(:,:), coef_q_weight(:,:) ! Boundary conditions (optional, default to 0 flux through boundary). - class(BoundaryType), target, intent(in), optional :: & - upper_bndry, lower_bndry + class(BoundaryType), intent(in), optional :: upper_bndry, lower_bndry ! Objects representing boundary conditions. class(BoundaryCond), intent(in), optional :: l_cond, r_cond @@ -88,31 +85,29 @@ function fin_vol_solve(dt, p, toSolve, ncols, pver, coef_q, coef_q_diff, coef_q_ ! start with an operator of all 0s. if (present(coef_q_diff)) then - net_operator = diffusion_operator(p, coef_q_diff, & - upper_bndry, lower_bndry) + net_operator = diffusion_operator(p, coef_q_diff, upper_bndry, lower_bndry) else - net_operator = zero_operator(p%n, p%d) + net_operator = zero_operator(p%n, p%d) end if ! Constant term (damping). if (present(coef_q)) then - add_term = diagonal_operator(coef_q) - call net_operator%add(add_term) + add_term = diagonal_operator(coef_q) + call net_operator%add(add_term) end if ! Effective advection. if (present(coef_q_adv)) then - add_term = advection_operator(p, coef_q_adv, & - upper_bndry, lower_bndry) - call net_operator%add(add_term) + add_term = advection_operator(p, coef_q_adv, upper_bndry, lower_bndry) + call net_operator%add(add_term) end if ! We want I-dt*(w^-1)*A for a single time step, implicit method, where ! A is the right-hand-side operator (i.e. what net_operator is now). if (present(coef_q_weight)) then - call net_operator%lmult_as_diag(-dt/coef_q_weight) + call net_operator%lmult_as_diag(-dt/coef_q_weight) else - call net_operator%lmult_as_diag(-dt) + call net_operator%lmult_as_diag(-dt) end if call net_operator%add_to_diag(1._r8) @@ -121,7 +116,6 @@ function fin_vol_solve(dt, p, toSolve, ncols, pver, coef_q, coef_q_diff, coef_q_ solution = toSolve call decomp%left_div(solution(:ncols, :), l_cond=l_cond, r_cond=r_cond) - !tendency = tendency - toSolve ! Ensure local objects are deallocated. call decomp%finalize() From 0c91251d3b7fccf41b1af57e4ef14d600f791767 Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Wed, 23 Oct 2024 07:05:06 -0600 Subject: [PATCH 14/14] Additional formatting changes. --- to_be_ccppized/vertical_diffusion_solver.F90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/to_be_ccppized/vertical_diffusion_solver.F90 b/to_be_ccppized/vertical_diffusion_solver.F90 index 1efdca06..c49d0c45 100644 --- a/to_be_ccppized/vertical_diffusion_solver.F90 +++ b/to_be_ccppized/vertical_diffusion_solver.F90 @@ -55,8 +55,10 @@ function fin_vol_solve(dt, p, toSolve, ncols, pver, coef_q, coef_q_diff, coef_q_ ! The sizes must be consistent among all the coefficients that are ! actually present, i.e. coef_q_diff and coef_q_adv should be one level ! bigger than coef_q and coef_q_weight, and have the same column number. - real(r8), contiguous, intent(in), optional :: coef_q(:,:), & - coef_q_diff(:,:), coef_q_adv(:,:), coef_q_weight(:,:) + real(r8), contiguous, intent(in), optional :: coef_q(:,:), & + coef_q_diff(:,:), & + coef_q_adv(:,:), & + coef_q_weight(:,:) ! Boundary conditions (optional, default to 0 flux through boundary). class(BoundaryType), intent(in), optional :: upper_bndry, lower_bndry