diff --git a/CMakeLists.txt b/CMakeLists.txt
index 5216e603fd..9084d59bb1 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -321,6 +321,7 @@ foreach(kind ${kinds})
field_manager/include
time_interp/include
tracer_manager/include
+ tridiagonal/include
interpolator/include
coupler/include
data_override/include)
diff --git a/Makefile.am b/Makefile.am
index ffb12344ea..dd1d27696d 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -35,8 +35,8 @@ endif
# Make targets will be run in each subdirectory. Order is significant.
SUBDIRS = \
platform \
- tridiagonal \
mpp \
+ tridiagonal \
constants \
constants4 \
memutils \
diff --git a/configure.ac b/configure.ac
index caf53f7de5..50a3b0f647 100644
--- a/configure.ac
+++ b/configure.ac
@@ -503,13 +503,13 @@ AC_CONFIG_FILES([
test_fms/coupler/Makefile
test_fms/parser/Makefile
test_fms/string_utils/Makefile
+ test_fms/tridiagonal/Makefile
test_fms/sat_vapor_pres/Makefile
test_fms/diag_integral/Makefile
test_fms/tracer_manager/Makefile
test_fms/random_numbers/Makefile
test_fms/topography/Makefile
test_fms/column_diagnostics/Makefile
-
FMS.pc
])
diff --git a/test_fms/Makefile.am b/test_fms/Makefile.am
index ea877c412e..3dd225360b 100644
--- a/test_fms/Makefile.am
+++ b/test_fms/Makefile.am
@@ -28,7 +28,7 @@ ACLOCAL_AMFLAGS = -I m4
SUBDIRS = astronomy coupler diag_manager data_override exchange monin_obukhov drifters \
mosaic interpolator fms mpp mpp_io time_interp time_manager horiz_interp topography \
field_manager axis_utils affinity fms2_io parser string_utils sat_vapor_pres tracer_manager \
-random_numbers diag_integral column_diagnostics
+random_numbers diag_integral column_diagnostics tridiagonal
# testing utility scripts to distribute
diff --git a/test_fms/tridiagonal/Makefile.am b/test_fms/tridiagonal/Makefile.am
new file mode 100644
index 0000000000..211869202b
--- /dev/null
+++ b/test_fms/tridiagonal/Makefile.am
@@ -0,0 +1,52 @@
+#***********************************************************************
+#* GNU Lesser General Public License
+#*
+#* This file is part of the GFDL Flexible Modeling System (FMS).
+#*
+#* FMS is free software: you can redistribute it and/or modify it under
+#* the terms of the GNU Lesser General Public License as published by
+#* the Free Software Foundation, either version 3 of the License, or (at
+#* your option) any later version.
+#*
+#* FMS is distributed in the hope that it will be useful, but WITHOUT
+#* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+#* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+#* for more details.
+#*
+#* You should have received a copy of the GNU Lesser General Public
+#* License along with FMS. If not, see .
+#***********************************************************************
+
+# This is an automake file for the test_fms/tridiagonal directory of the FMS
+# package.
+
+# Ryan Mulhall
+
+# Find the .mod directory
+AM_CPPFLAGS = -I$(top_srcdir)/include -I$(MODDIR)
+
+# Link to the FMS library.
+LDADD = $(top_builddir)/libFMS/libFMS.la
+
+# Build this test program.
+check_PROGRAMS = test_tridiagonal_r4 test_tridiagonal_r8
+
+# compiles test file with both kind sizes via macro
+test_tridiagonal_r4_SOURCES=test_tridiagonal.F90
+test_tridiagonal_r8_SOURCES=test_tridiagonal.F90
+
+test_tridiagonal_r4_CPPFLAGS=-DTRID_REAL_TYPE=tridiag_r4 -DTEST_TRIDIAG_REAL=r4_kind -I$(MODDIR)
+test_tridiagonal_r8_CPPFLAGS=-DTRID_REAL_TYPE=tridiag_r8 -DTEST_TRIDIAG_REAL=r8_kind -I$(MODDIR)
+
+# Run the test program.
+TESTS = test_tridiagonal.sh
+
+TEST_EXTENSIONS = .sh
+SH_LOG_DRIVER = env AM_TAP_AWK='$(AWK)' $(SHELL) \
+ $(abs_top_srcdir)/test_fms/tap-driver.sh
+
+# These files will be included in the distribution.
+EXTRA_DIST = test_tridiagonal.sh
+
+# Clean up
+CLEANFILES = *.nml *.out* *.dpi *.spi *.dyn *.spl *.o test_tridiagonal4 test_tridiagonal8 test_tridiagonal
diff --git a/test_fms/tridiagonal/test_tridiagonal.F90 b/test_fms/tridiagonal/test_tridiagonal.F90
new file mode 100644
index 0000000000..18200a8c77
--- /dev/null
+++ b/test_fms/tridiagonal/test_tridiagonal.F90
@@ -0,0 +1,173 @@
+!***********************************************************************
+!* GNU Lesser General Public License
+!*
+!* This file is part of the GFDL Flexible Modeling System (FMS).
+!*
+!* FMS is free software: you can redistribute it and/or modify it under
+!* the terms of the GNU Lesser General Public License as published by
+!* the Free Software Foundation, either version 3 of the License, or (at
+!* your option) any later version.
+!*
+!* FMS is distributed in the hope that it will be useful, but WITHOUT
+!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+!* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+!* for more details.
+!*
+!* You should have received a copy of the GNU Lesser General Public
+!* License along with FMS. If not, see .
+!***********************************************************************
+#ifndef TEST_TRIDIAG_KIND
+#define TEST_TRIDIAG_KIND 8
+#endif
+
+!> Tests the tridiagonal module routines (tri_invert and close_tridiagonal)
+!! Tests reals with the kind value set above,
+program test_tridiagonal
+
+ use tridiagonal_mod
+ use platform_mod
+ use mpp_mod
+ use fms_mod
+
+ implicit none
+
+ integer, parameter :: IN_LEN = 8 !< length of input arrays
+ integer, parameter :: kindl = TEST_TRIDIAG_KIND !< kind value for all reals in this test
+ !! set by TEST_TRIDIAG_KIND cpp macro
+ real(TEST_TRIDIAG_KIND), allocatable :: d(:,:,:), x(:,:,:), ref_array(:,:,:)
+ real(TEST_TRIDIAG_KIND), allocatable :: a(:,:,:), b(:,:,:), c(:,:,:)
+ real(r4_kind), allocatable :: d_r4(:,:,:), x_r4(:,:,:)
+ real(r4_kind), allocatable :: a_r4(:,:,:), b_r4(:,:,:), c_r4(:,:,:)
+ real(r8_kind), allocatable :: d_r8(:,:,:), x_r8(:,:,:)
+ real(r8_kind), allocatable :: a_r8(:,:,:), b_r8(:,:,:), c_r8(:,:,:)
+ integer :: i, end, ierr, io
+ real(TEST_TRIDIAG_KIND) :: k
+ ! nml
+ logical :: do_error_check = .false.
+ namelist / test_tridiagonal_nml/ do_error_check
+
+ call mpp_init
+
+ read (input_nml_file, test_tridiagonal_nml, iostat=io)
+ ierr = check_nml_error (io, 'test_tridiagonal_nml')
+
+ ! allocate input and output arrays
+ allocate(d(1,1,IN_LEN))
+ allocate(a(1,1,IN_LEN))
+ allocate(b(1,1,IN_LEN))
+ allocate(c(1,1,IN_LEN))
+ allocate(x(1,1,IN_LEN))
+
+ !! simple test with only 1 coeff
+ a = 0.0_kindl
+ b = 1.0_kindl
+ c = 0.0_kindl
+ d = 5.0_kindl
+ call tri_invert(x, d, a, b, c)
+ if(any(x .ne. 5.0_kindl)) call mpp_error(FATAL, "test_tridiagonal: invalid results for 1 coefficient check")
+ !! check with stored data arrays
+ d = -5.0_kindl
+ call tri_invert(x, d)
+ if(any(x .ne. -5.0_kindl)) call mpp_error(FATAL, "test_tridiagonal: invalid results for 1 coefficient check")
+
+ ! test with a,b,c
+ ! 0.5x(n-2) + x(n-1) + 0.5x(n) = 1
+ !
+ ! x(n) = k * [4, 1, 3, 2, 2, 3, 1, 4]
+ ! k * [8 , 1, 7, 2, 6, .. ] = k *(-n/2 + ((n%2)*arr_length/2))
+ a = 0.5_kindl
+ b = 1.0_kindl
+ c = 0.5_kindl
+ d = 1.0_kindl
+ call tri_invert(x, d, a, b, c)
+ ! set up reference answers
+ k = 1.0_kindl/(IN_LEN+1.0_kindl) * 2.0_kindl
+ allocate(ref_array(1,1,IN_LEN))
+ do i=1, IN_LEN/2
+ end=IN_LEN-i+1
+ if(mod(i, 2) .eq. 1) then
+ ref_array(1,1,i) = real(-(i/2) + (mod(i,2)* IN_LEN/2), kindl)
+ ref_array(1,1,end) = real(-(i/2) + (mod(i,2)* IN_LEN/2), kindl)
+ else
+ ref_array(1,1,i) = real(i/2, kindl)
+ ref_array(1,1,end) = real(i/2, kindl)
+ endif
+ enddo
+ ref_array = ref_array * k
+ ! check
+ do i=1, IN_LEN
+ if(ABS(x(1,1,i) - ref_array(1,1,i)) .gt. 0.1e-12_kindl) then
+ print *, i, x(1,1,i), ref_array(1,1,i)
+ call mpp_error(FATAL, "test_tridiagonal: failed reference check for tri_invert")
+ endif
+ enddo
+ !! check with stored data arrays
+ d = -1.0_kindl
+ ref_array = ref_array * -1.0_kindl
+ call tri_invert(x, d)
+ do i=1, IN_LEN
+ if(ABS(x(1,1,i) - ref_array(1,1,i)) .gt. 0.1e-12_kindl) then
+ print *, i, x(1,1,i), ref_array(1,1,i)
+ call mpp_error(FATAL, "test_tridiagonal: failed reference check for tri_invert with saved values")
+ endif
+ enddo
+ call close_tridiagonal()
+
+ !! tests for module state across kinds
+ !! default keeps stored values separate depending on kind
+ !! store_both_kinds argument can be specified to store both r4 and r8 kinds
+ if(kindl .eq. r8_kind) then
+ allocate(a_r4(1,1,IN_LEN), b_r4(1,1,IN_LEN), c_r4(1,1,IN_LEN))
+ allocate(d_r4(1,1,IN_LEN), x_r4(1,1,IN_LEN))
+ a_r4 = 0.0_r4_kind; b_r4 = 1.0_r4_kind; c_r4 = 0.0_r4_kind
+ d_r4 = 5.0_r4_kind; x_r4 = 0.0_r4_kind
+ a = 0.0_kindl; b = 2.0_kindl; c = 0.0_kindl
+ d = 5.0_kindl
+ ! default, module variables distinct per kind
+ call tri_invert(x_r4, d_r4, a_r4, b_r4, c_r4)
+ ! conditionally errors here for calling with unallocated a/b/c for kind
+ if( do_error_check ) call tri_invert(x, d)
+ call tri_invert(x, d, a, b, c)
+ ! check both values are correct from prior state
+ call tri_invert(x_r4, d_r4)
+ call tri_invert(x, d)
+ if(any(x_r4 .ne. 5.0_r4_kind)) call mpp_error(FATAL, "test_tridiagonal: invalid r4 kind result")
+ if(any(x .ne. 2.5_r8_kind)) call mpp_error(FATAL, "test_tridiagonal: invalid r8 kind result")
+ call close_tridiagonal()
+ ! run with storing for both kinds
+ call tri_invert(x_r4, d_r4, a_r4, b_r4, c_r4, store_both_kinds=.true.)
+ call tri_invert(x_r4, d_r4)
+ call tri_invert(x, d)
+ if(any(x_r4 .ne. 5.0_r4_kind)) call mpp_error(FATAL, "test_tridiagonal: invalid r4 kind result")
+ if(any(x .ne. 5.0_r8_kind)) call mpp_error(FATAL, "test_tridiagonal: invalid r8 kind result")
+ else
+ allocate(a_r8(1,1,IN_LEN), b_r8(1,1,IN_LEN), c_r8(1,1,IN_LEN))
+ allocate(d_r8(1,1,IN_LEN), x_r8(1,1,IN_LEN))
+ a_r8 = 0.0_r8_kind; b_r8 = 1.0_r8_kind; c_r8 = 0.0_r8_kind
+ d_r8 = 5.0_r8_kind; x_r8 = 0.0_r8_kind
+ a = 0.0_kindl; b = 2.0_kindl; c = 0.0_kindl
+ d = 5.0_kindl
+ ! default, module variables distinct per kind
+ call tri_invert(x_r8, d_r8, a_r8, b_r8, c_r8)
+ ! conditionally errors here for calling with unallocated a/b/c for kind
+ if( do_error_check ) call tri_invert(x, d)
+ call tri_invert(x, d, a, b, c)
+ ! check both values are correct from prior state
+ call tri_invert(x_r8, d_r8)
+ call tri_invert(x, d)
+ if(any(x_r8 .ne. 5.0_r8_kind)) call mpp_error(FATAL, "test_tridiagonal: invalid r8 kind result")
+ if(any(x .ne. 2.5_r8_kind)) call mpp_error(FATAL, "test_tridiagonal: invalid r8 kind result")
+ call close_tridiagonal()
+ ! run with storing for both kinds
+ call tri_invert(x_r8, d_r8, a_r8, b_r8, c_r8, store_both_kinds=.true.)
+ call tri_invert(x_r8, d_r8)
+ call tri_invert(x, d)
+ if(any(x_r8 .ne. 5.0_r8_kind)) call mpp_error(FATAL, "test_tridiagonal: invalid r8 kind result")
+ if(any(x .ne. 5.0_r8_kind)) call mpp_error(FATAL, "test_tridiagonal: invalid r8 kind result")
+ endif
+
+ call close_tridiagonal()
+
+ call mpp_exit
+
+end program
\ No newline at end of file
diff --git a/test_fms/tridiagonal/test_tridiagonal.sh b/test_fms/tridiagonal/test_tridiagonal.sh
new file mode 100755
index 0000000000..4be1fa80b1
--- /dev/null
+++ b/test_fms/tridiagonal/test_tridiagonal.sh
@@ -0,0 +1,51 @@
+#!/bin/sh
+
+#***********************************************************************
+#* GNU Lesser General Public License
+#*
+#* This file is part of the GFDL Flexible Modeling System (FMS).
+#*
+#* FMS is free software: you can redistribute it and/or modify it under
+#* the terms of the GNU Lesser General Public License as published by
+#* the Free Software Foundation, either version 3 of the License, or (at
+#* your option) any later version.
+#*
+#* FMS is distributed in the hope that it will be useful, but WITHOUT
+#* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+#* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+#* for more details.
+#*
+#* You should have received a copy of the GNU Lesser General Public
+#* License along with FMS. If not, see .
+#***********************************************************************
+
+# This is part of the GFDL FMS package. This is a shell script to
+# execute tests in the test_fms/time_manager directory.
+
+# Ryan Mulhall 9/2023
+
+# Set common test settings.
+. ../test-lib.sh
+
+rm -f input.nml && touch input.nml
+
+test_expect_success "test tridiagonal functionality 32 bit reals" '
+ mpirun -n 1 ./test_tridiagonal_r4
+'
+test_expect_success "test tridiagonal functionality 64 bit reals" '
+ mpirun -n 1 ./test_tridiagonal_r8
+'
+# tries to call without a,b,c args provided or previously set
+cat <<_EOF > input.nml
+&test_tridiagonal_nml
+do_error_check = .true.
+/
+_EOF
+test_expect_failure "error out if passed in incorrect real size (r4_kind)" '
+ mpirun -n 1 ./test_tridiagonal_r4
+'
+test_expect_failure "error out if passed in incorrect real size (r8_kind)" '
+ mpirun -n 1 ./test_tridiagonal_r8
+'
+
+test_done
diff --git a/tridiagonal/Makefile.am b/tridiagonal/Makefile.am
index 177ca904a5..d8b90c409b 100644
--- a/tridiagonal/Makefile.am
+++ b/tridiagonal/Makefile.am
@@ -23,14 +23,17 @@
# Ed Hartnett 2/22/19
# Include .h and .mod files.
-AM_CPPFLAGS = -I$(top_srcdir)/include
+AM_CPPFLAGS = -I$(top_srcdir)/include -I$(top_srcdir)/tridiagonal/include
AM_FCFLAGS = $(FC_MODINC). $(FC_MODOUT)$(MODDIR)
# Build this uninstalled convenience library.
noinst_LTLIBRARIES = libtridiagonal.la
# The convenience library depends on its source.
-libtridiagonal_la_SOURCES = tridiagonal.F90
+libtridiagonal_la_SOURCES = tridiagonal.F90 \
+ include/tridiagonal.inc \
+ include/tridiagonal_r4.fh \
+ include/tridiagonal_r8.fh
# Mod file depends on its o file, is built and then installed.
tridiagonal.lo: tridiagonal_mod.$(FC_MODEXT)
diff --git a/tridiagonal/include/tridiagonal.inc b/tridiagonal/include/tridiagonal.inc
index c22f99c4ee..95788eb795 100644
--- a/tridiagonal/include/tridiagonal.inc
+++ b/tridiagonal/include/tridiagonal.inc
@@ -16,58 +16,9 @@
!* You should have received a copy of the GNU Lesser General Public
!* License along with FMS. If not, see .
!***********************************************************************
-!> @defgroup tridiagonal_mod tridiagonal_mod
-!> @ingroup tridiagonal
-!> @brief Solves a tridiagonal system of equations.
-!!
-!> The following schematic represents the system of equations solved,
-!! where X is the solution.
-!!
-!! | B(1) A(1) 0 0 ....... 0 | |X(1)| |D(1)|
-!! | C(2) B(2) A(2) 0 ....... 0 | |X(2)| |D(2)|
-!! | 0 C(3) B(3) A(3) 0 ....... 0 | | .. | | .. |
-!! | .......................................... | | .. | = | .. |
-!! | .......................................... | | .. | | .. |
-!! | C(N-2) B(N-2) A(N-2) 0 | | .. | | .. |
-!! | 0 C(N-1) B(N-1) A(N-1)| | .. | | .. |
-!! | 0 0 C(N) B(N) | |X(N)| |D(N)|
-!!
-!!
-!! To solve this system
-!!
-!! call tri_invert(X,D,A,B,C)
-!!
-!! real, intent(out), dimension(:,:,:) :: X
-!! real, intent(in), dimension(:,:,:) :: D
-!! real, optional, dimension(:,:,:) :: A,B,C
-!!
-!! For simplicity (?), A and C are assumed to be dimensioned the same size
-!! as B, D, and X, although any input values for A(N) and C(1) are ignored.
-!! (some checks are needed here)
-!!
-!! If A is not present, it is assumed that the matrix (A,B.C) has not been changed
-!! since the last call to tri_invert.
-!!
-!! To release memory,
-!!
-!! call close_tridiagonal
-!!
-!! Arguments A, B, and C are optional, and are saved as module variables
-!! if one recalls tri_invert without changing (A,B,C)
!> @addtogroup tridiagonal_mod
!> @{
-module tridiagonal_mod
-
-!--------------------------------------------------------------------------
-real, private, allocatable, dimension(:,:,:) :: e,g,cc
-real, private, allocatable, dimension(:,:) :: bb
-logical, private :: init_tridiagonal = .false.
-!--------------------------------------------------------------------------
-
-contains
-
-!--------------------------------------------------------------------------
!> @brief Sets up and solves the tridiagonal system of equations
!!
@@ -76,104 +27,80 @@ contains
!! There are no checks to make sure the sizes agree.
!!
!! The value of A(N) is modified on output, and B and C are unchanged.
-subroutine tri_invert(x,d,a,b,c)
-
-implicit none
-
-real, intent(out), dimension(:,:,:) :: x !< Solution to the tridiagonal system of equations
-real, intent(in), dimension(:,:,:) :: d !< The right-hand side term, see the schematic above.
-real, optional, dimension(:,:,:) :: a,b,c !< Left hand side terms(see schematic above).
- !! If not provided, values from last call are used
-
-real, dimension(size(x,1),size(x,2),size(x,3)) :: f
-integer :: k
-
-if(present(a)) then
-
- !< Check if module variables are allocated
- !$OMP SINGLE
- init_tridiagonal = .true.
- if(allocated(e)) deallocate(e)
- if(allocated(g)) deallocate(g)
- if(allocated(bb)) deallocate(bb)
- if(allocated(cc)) deallocate(cc)
- allocate(e (size(x,1),size(x,2),size(x,3)))
- allocate(g (size(x,1),size(x,2),size(x,3)))
- allocate(bb(size(x,1),size(x,2)))
- allocate(cc(size(x,1),size(x,2),size(x,3)))
- !$OMP END SINGLE !< There is an implicit barrier.
-
- e(:,:,1) = - a(:,:,1)/b(:,:,1)
- a(:,:,size(x,3)) = 0.0
-
- do k= 2,size(x,3)
- g(:,:,k) = 1.0/(b(:,:,k)+c(:,:,k)*e(:,:,k-1))
- e(:,:,k) = - a(:,:,k)*g(:,:,k)
- end do
- cc = c
- bb = 1.0/b(:,:,1)
-
-end if
-
-! if(.not.init_tridiagonal) error
-
-f(:,:,1) = d(:,:,1)*bb
-do k= 2, size(x,3)
- f(:,:,k) = (d(:,:,k) - cc(:,:,k)*f(:,:,k-1))*g(:,:,k)
-end do
-
-x(:,:,size(x,3)) = f(:,:,size(x,3))
-do k = size(x,3)-1,1,-1
- x(:,:,k) = e(:,:,k)*x(:,:,k+1)+f(:,:,k)
-end do
-
-return
-end subroutine tri_invert
-
-!-----------------------------------------------------------------
-
-!> @brief Releases memory used by the solver
-subroutine close_tridiagonal
-
- implicit none
-
- !< Check if module variables are allocated
- !$OMP SINGLE
- if(allocated(e)) deallocate(e)
- if(allocated(g)) deallocate(g)
- if(allocated(bb)) deallocate(bb)
- if(allocated(cc)) deallocate(cc)
- !$OMP END SINGLE !< There is an implicit barrier.
-
-return
-end subroutine close_tridiagonal
-
-!----------------------------------------------------------------
-
-end module tridiagonal_mod
-
-!
-
-!
-! Optional arguments A,B,C have no intent declaration,
-! so the default intent is inout. The value of A(N) is modified
-! on output, and B and C are unchanged.
-!
-!
-! The following private allocatable arrays save the relevant information
-! if one recalls tri_invert without changing (A,B,C):
-!
-! allocate ( e (size(x,1), size(x,2), size(x,3)) )
-! allocate ( g (size(x,1), size(x,2), size(x,3)) )
-! allocate ( cc (size(x,1), size(x,2), size(x,3)) )
-! allocate ( bb (size(x,1), size(x,2)) )
-!
-! This storage is deallocated when close_tridiagonal is called.
-!
-!
-! Maybe a cleaner version?
-!
-
-!
-!> @}
-! close documentation grouping
+!!
+!! For mixed precision, this routine uses the kind size macro(FMS_TRID_KIND_) to determine
+!! which module variables are used/stored. This means a,b, and c values will only be stored for calls
+!! of the same real kind value unless store_both_kinds is present and .true..
+subroutine TRI_INVERT_(x,d,a,b,c, store_both_kinds)
+
+ real(FMS_TRID_KIND_), intent(out), dimension(:,:,:) :: x !< Solution to the tridiagonal system of equations
+ real(FMS_TRID_KIND_), intent(in), dimension(:,:,:) :: d !< The right-hand side term, see the schematic above.
+ real(FMS_TRID_KIND_), optional, dimension(:,:,:) :: a,b,c !< Left hand side terms(see schematic on module page).
+ !! If not provided, values from last call are used
+ logical, optional :: store_both_kinds !< Will save module state
+ !! variables for both kind types in order to be used in
+ !! subsequent calls with either kind.
+
+ real(FMS_TRID_KIND_), dimension(size(x,1),size(x,2),size(x,3)) :: f
+ integer, parameter :: kindl = FMS_TRID_KIND_
+
+ integer :: k
+
+ if(present(a)) then
+ !$OMP SINGLE
+ INIT_VAR = .true.
+ if(allocated(TRID_REAL_TYPE%e)) deallocate(TRID_REAL_TYPE%e)
+ if(allocated(TRID_REAL_TYPE%g)) deallocate(TRID_REAL_TYPE%g)
+ if(allocated(TRID_REAL_TYPE%bb)) deallocate(TRID_REAL_TYPE%bb)
+ if(allocated(TRID_REAL_TYPE%cc)) deallocate(TRID_REAL_TYPE%cc)
+ allocate(TRID_REAL_TYPE%e (size(x,1),size(x,2),size(x,3)))
+ allocate(TRID_REAL_TYPE%g (size(x,1),size(x,2),size(x,3)))
+ allocate(TRID_REAL_TYPE%bb(size(x,1),size(x,2)))
+ allocate(TRID_REAL_TYPE%cc(size(x,1),size(x,2),size(x,3)))
+ !$OMP END SINGLE
+
+ TRID_REAL_TYPE%e(:,:,1) = - a(:,:,1) / b(:,:,1)
+ a(:,:,size(x,3)) = 0.0_kindl
+
+ do k= 2,size(x,3)
+ TRID_REAL_TYPE%g(:,:,k) = 1.0_kindl/(b(:,:,k)+c(:,:,k)*TRID_REAL_TYPE%e(:,:,k-1))
+ TRID_REAL_TYPE%e(:,:,k) = - a(:,:,k)* TRID_REAL_TYPE%g(:,:,k)
+ end do
+ TRID_REAL_TYPE%cc = c
+ TRID_REAL_TYPE%bb = 1.0_kindl/b(:,:,1)
+
+ end if
+
+ if(.not.INIT_VAR) call mpp_error(FATAL, 'tri_invert: a,b,and c args not provided or previously calculated.')
+
+ f(:,:,1) = d(:,:,1)*TRID_REAL_TYPE%bb
+ do k= 2, size(x,3)
+ f(:,:,k) = (d(:,:,k) - TRID_REAL_TYPE%cc(:,:,k)*f(:,:,k-1))*TRID_REAL_TYPE%g(:,:,k)
+ end do
+
+ x(:,:,size(x,3)) = f(:,:,size(x,3))
+ do k = size(x,3)-1,1,-1
+ x(:,:,k) = TRID_REAL_TYPE%e(:,:,k)*x(:,:,k+1)+f(:,:,k)
+ end do
+
+ ! stores both kind values for subsequent calculations if running with option
+ if( present(store_both_kinds)) then
+ if( store_both_kinds ) then
+ if( FMS_TRID_KIND_ .eq. r8_kind) then
+ tridiag_r4%e = real(TRID_REAL_TYPE%e, r4_kind)
+ tridiag_r4%g = real(TRID_REAL_TYPE%g, r4_kind)
+ tridiag_r4%cc = real(TRID_REAL_TYPE%cc, r4_kind)
+ tridiag_r4%bb = real(TRID_REAL_TYPE%bb, r4_kind)
+ init_tridiagonal_r4 = .true.
+ else
+ tridiag_r8%e = real(TRID_REAL_TYPE%e, r8_kind)
+ tridiag_r8%g = real(TRID_REAL_TYPE%g, r8_kind)
+ tridiag_r8%cc = real(TRID_REAL_TYPE%cc, r8_kind)
+ tridiag_r8%bb = real(TRID_REAL_TYPE%bb, r8_kind)
+ init_tridiagonal_r8 = .true.
+ endif
+ endif
+ endif
+
+ return
+end subroutine TRI_INVERT_
\ No newline at end of file
diff --git a/tridiagonal/include/tridiagonal_r4.fh b/tridiagonal/include/tridiagonal_r4.fh
new file mode 100644
index 0000000000..09e0ad57ac
--- /dev/null
+++ b/tridiagonal/include/tridiagonal_r4.fh
@@ -0,0 +1,32 @@
+!***********************************************************************
+!* GNU Lesser General Public License
+!*
+!* This file is part of the GFDL Flexible Modeling System (FMS).
+!*
+!* FMS is free software: you can redistribute it and/or modify it under
+!* the terms of the GNU Lesser General Public License as published by
+!* the Free Software Foundation, either version 3 of the License, or (at
+!* your option) any later version.
+!*
+!* FMS is distributed in the hope that it will be useful, but WITHOUT
+!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+!* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+!* for more details.
+!*
+!* You should have received a copy of the GNU Lesser General Public
+!* License along with FMS. If not, see .
+!***********************************************************************
+
+#undef FMS_TRID_KIND_
+#define FMS_TRID_KIND_ r4_kind
+
+#undef TRID_REAL_TYPE
+#define TRID_REAL_TYPE tridiag_r4
+
+#undef TRI_INVERT_
+#define TRI_INVERT_ tri_invert_r4
+
+#undef INIT_VAR
+#define INIT_VAR init_tridiagonal_r4
+
+#include "tridiagonal.inc"
\ No newline at end of file
diff --git a/tridiagonal/include/tridiagonal_r8.fh b/tridiagonal/include/tridiagonal_r8.fh
new file mode 100644
index 0000000000..b007941723
--- /dev/null
+++ b/tridiagonal/include/tridiagonal_r8.fh
@@ -0,0 +1,32 @@
+!***********************************************************************
+!* GNU Lesser General Public License
+!*
+!* This file is part of the GFDL Flexible Modeling System (FMS).
+!*
+!* FMS is free software: you can redistribute it and/or modify it under
+!* the terms of the GNU Lesser General Public License as published by
+!* the Free Software Foundation, either version 3 of the License, or (at
+!* your option) any later version.
+!*
+!* FMS is distributed in the hope that it will be useful, but WITHOUT
+!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+!* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+!* for more details.
+!*
+!* You should have received a copy of the GNU Lesser General Public
+!* License along with FMS. If not, see .
+!***********************************************************************
+
+#undef FMS_TRID_KIND_
+#define FMS_TRID_KIND_ r8_kind
+
+#undef TRID_REAL_TYPE
+#define TRID_REAL_TYPE tridiag_r8
+
+#undef TRI_INVERT_
+#define TRI_INVERT_ tri_invert_r8
+
+#undef INIT_VAR
+#define INIT_VAR init_tridiagonal_r8
+
+#include "tridiagonal.inc"
\ No newline at end of file
diff --git a/tridiagonal/tridiagonal.F90 b/tridiagonal/tridiagonal.F90
index c22f99c4ee..e34feb4d92 100644
--- a/tridiagonal/tridiagonal.F90
+++ b/tridiagonal/tridiagonal.F90
@@ -52,128 +52,85 @@
!!
!! call close_tridiagonal
!!
+!!
+!!
!! Arguments A, B, and C are optional, and are saved as module variables
!! if one recalls tri_invert without changing (A,B,C)
+!!
+!! @note
+!! Optional arguments A,B,C have no intent declaration,
+!! so the default intent is inout. The value of A(N) is modified
+!! on output, and B and C are unchanged.
+!!
+!! The following private allocatable arrays save the relevant information
+!! if one recalls tri_invert without changing (A,B,C):
+!!
+!! allocate ( e (size(x,1), size(x,2), size(x,3)) )
+!! allocate ( g (size(x,1), size(x,2), size(x,3)) )
+!! allocate ( cc (size(x,1), size(x,2), size(x,3)) )
+!! allocate ( bb (size(x,1), size(x,2)) )
+!!
+!! This storage is deallocated when close_tridiagonal is called.
!> @addtogroup tridiagonal_mod
!> @{
module tridiagonal_mod
-!--------------------------------------------------------------------------
-real, private, allocatable, dimension(:,:,:) :: e,g,cc
-real, private, allocatable, dimension(:,:) :: bb
-logical, private :: init_tridiagonal = .false.
-!--------------------------------------------------------------------------
-
-contains
-
-!--------------------------------------------------------------------------
-
-!> @brief Sets up and solves the tridiagonal system of equations
-!!
-!> For simplicity, A and C are assumed to be dimensioned the same size
-!! as B, D, and X, although any input values for A(N) and C(1) are ignored.
-!! There are no checks to make sure the sizes agree.
-!!
-!! The value of A(N) is modified on output, and B and C are unchanged.
-subroutine tri_invert(x,d,a,b,c)
-
-implicit none
-
-real, intent(out), dimension(:,:,:) :: x !< Solution to the tridiagonal system of equations
-real, intent(in), dimension(:,:,:) :: d !< The right-hand side term, see the schematic above.
-real, optional, dimension(:,:,:) :: a,b,c !< Left hand side terms(see schematic above).
- !! If not provided, values from last call are used
-
-real, dimension(size(x,1),size(x,2),size(x,3)) :: f
-integer :: k
-
-if(present(a)) then
-
- !< Check if module variables are allocated
- !$OMP SINGLE
- init_tridiagonal = .true.
- if(allocated(e)) deallocate(e)
- if(allocated(g)) deallocate(g)
- if(allocated(bb)) deallocate(bb)
- if(allocated(cc)) deallocate(cc)
- allocate(e (size(x,1),size(x,2),size(x,3)))
- allocate(g (size(x,1),size(x,2),size(x,3)))
- allocate(bb(size(x,1),size(x,2)))
- allocate(cc(size(x,1),size(x,2),size(x,3)))
- !$OMP END SINGLE !< There is an implicit barrier.
-
- e(:,:,1) = - a(:,:,1)/b(:,:,1)
- a(:,:,size(x,3)) = 0.0
-
- do k= 2,size(x,3)
- g(:,:,k) = 1.0/(b(:,:,k)+c(:,:,k)*e(:,:,k-1))
- e(:,:,k) = - a(:,:,k)*g(:,:,k)
- end do
- cc = c
- bb = 1.0/b(:,:,1)
-
-end if
-
-! if(.not.init_tridiagonal) error
-
-f(:,:,1) = d(:,:,1)*bb
-do k= 2, size(x,3)
- f(:,:,k) = (d(:,:,k) - cc(:,:,k)*f(:,:,k-1))*g(:,:,k)
-end do
-
-x(:,:,size(x,3)) = f(:,:,size(x,3))
-do k = size(x,3)-1,1,-1
- x(:,:,k) = e(:,:,k)*x(:,:,k+1)+f(:,:,k)
-end do
-
-return
-end subroutine tri_invert
-
-!-----------------------------------------------------------------
-
-!> @brief Releases memory used by the solver
-subroutine close_tridiagonal
-
- implicit none
-
- !< Check if module variables are allocated
- !$OMP SINGLE
- if(allocated(e)) deallocate(e)
- if(allocated(g)) deallocate(g)
- if(allocated(bb)) deallocate(bb)
- if(allocated(cc)) deallocate(cc)
- !$OMP END SINGLE !< There is an implicit barrier.
-
-return
-end subroutine close_tridiagonal
-
-!----------------------------------------------------------------
+ use platform_mod, only: r4_kind, r8_kind
+ use mpp_mod, only: mpp_error, FATAL
+ implicit none
+
+ type :: tridiag_reals_r4
+ real(r4_kind), private, allocatable, dimension(:,:,:) :: e, g, cc
+ real(r4_kind), private, allocatable, dimension(:,:) :: bb
+ end type
+
+ type :: tridiag_reals_r8
+ real(r8_kind), private, allocatable, dimension(:,:,:) :: e, g, cc
+ real(r8_kind), private, allocatable, dimension(:,:) :: bb
+ end type
+
+ type(tridiag_reals_r4) :: tridiag_r4 !< holds reals stored from r4_kind calls to tri_invert
+ type(tridiag_reals_r8) :: tridiag_r8 !< holds reals stored from r8_kind calls to tri_invert
+
+ logical, private :: init_tridiagonal_r4 = .false. !< true when fields in tridiag_r4 are allocated
+ logical, private :: init_tridiagonal_r8 = .false. !< true when fields in tridiag_r8 are allocated
+
+ !> Interface to solve tridiagonal systems of equations for either kind value.
+ !! Module level variables will be deallocated and allocated for every
+ !! Since this relies on the state of module variables (unless A,B,C are specified)
+ !! the values stored are distinct for each kind call unless the added optional argument store_both_kinds
+ !! is true.
+ interface tri_invert
+ module procedure tri_invert_r4
+ module procedure tri_invert_r8
+ end interface
+
+ public :: tri_invert
+
+ contains
+
+ !> @brief Releases memory used by the solver
+ subroutine close_tridiagonal
+ if(.not. init_tridiagonal_r4 .and. .not. init_tridiagonal_r8) return
+ !$OMP SINGLE
+ if(allocated(tridiag_r4%e)) deallocate(tridiag_r4%e)
+ if(allocated(tridiag_r4%g)) deallocate(tridiag_r4%g)
+ if(allocated(tridiag_r4%cc)) deallocate(tridiag_r4%cc)
+ if(allocated(tridiag_r4%bb)) deallocate(tridiag_r4%bb)
+ if(allocated(tridiag_r8%e)) deallocate(tridiag_r8%e)
+ if(allocated(tridiag_r8%g)) deallocate(tridiag_r8%g)
+ if(allocated(tridiag_r8%cc)) deallocate(tridiag_r8%cc)
+ if(allocated(tridiag_r8%bb)) deallocate(tridiag_r8%bb)
+ init_tridiagonal_r4 = .false.; init_tridiagonal_r8 = .false.
+ !$OMP END SINGLE
+ return
+ end subroutine close_tridiagonal
+
+#include "tridiagonal_r4.fh"
+#include "tridiagonal_r8.fh"
end module tridiagonal_mod
-!
-
-!
-! Optional arguments A,B,C have no intent declaration,
-! so the default intent is inout. The value of A(N) is modified
-! on output, and B and C are unchanged.
-!
-!
-! The following private allocatable arrays save the relevant information
-! if one recalls tri_invert without changing (A,B,C):
-!
-! allocate ( e (size(x,1), size(x,2), size(x,3)) )
-! allocate ( g (size(x,1), size(x,2), size(x,3)) )
-! allocate ( cc (size(x,1), size(x,2), size(x,3)) )
-! allocate ( bb (size(x,1), size(x,2)) )
-!
-! This storage is deallocated when close_tridiagonal is called.
-!
-!
-! Maybe a cleaner version?
-!
-
-!
!> @}
-! close documentation grouping
+! close documentation grouping
\ No newline at end of file