input.nml
+&topography_nml
+ topog_file = "topography.data.nc",
+ water_file = "water.data.nc"
+/
+
+EOF
+
+# Run the test.
+
+test_expect_success "Test topography: r4_kind" '
+ mpirun -n 1 ./test_topography_r4
+'
+
+sync; rm -f *.nc
+
+test_expect_success "Test topography: r8_kind" '
+ mpirun -n 1 ./test_topography_r8
+'
+
+rm -f *.nc
+
+test_done
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
new file mode 100644
index 0000000000..95788eb795
--- /dev/null
+++ b/tridiagonal/include/tridiagonal.inc
@@ -0,0 +1,106 @@
+!***********************************************************************
+!* 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 .
+!***********************************************************************
+
+!> @addtogroup tridiagonal_mod
+!> @{
+
+!> @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.
+!!
+!! 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
From 03ee7de39ea5f19e29746aadd3e46f8eef669bb1 Mon Sep 17 00:00:00 2001
From: Ryan Mulhall <35538242+rem1776@users.noreply.github.com>
Date: Fri, 6 Oct 2023 11:16:21 -0400
Subject: [PATCH 05/10] fix: update macros and build options for mixed
precision (#1383)
---
.github/workflows/github_autotools_gnu.yml | 2 +-
configure.ac | 43 +-
diag_manager/diag_axis.F90 | 2 -
diag_manager/diag_data.F90 | 6 -
diag_manager/diag_manager.F90 | 2 -
diag_manager/diag_output.F90 | 2 -
diag_manager/diag_util.F90 | 2 -
include/fms_platform.h | 59 +-
mosaic/read_mosaic.c | 2 -
test_fms/mosaic2/Makefile.am | 8 +-
time_interp/include/time_interp_external.inc | 1424 ------------------
11 files changed, 25 insertions(+), 1527 deletions(-)
delete mode 100644 time_interp/include/time_interp_external.inc
diff --git a/.github/workflows/github_autotools_gnu.yml b/.github/workflows/github_autotools_gnu.yml
index b7008a0814..00c7f5d31c 100644
--- a/.github/workflows/github_autotools_gnu.yml
+++ b/.github/workflows/github_autotools_gnu.yml
@@ -9,7 +9,7 @@ jobs:
runs-on: ubuntu-latest
strategy:
matrix:
- conf-flag: [ --disable-openmp, --enable-mixed-mode, --disable-setting-flags, --with-mpi=no]
+ conf-flag: [ --disable-openmp, --disable-setting-flags, --with-mpi=no, --disable-r8-defaults]
input-flag: [--with-yaml, --enable-test-input=/home/unit_tests_input]
exclude:
- conf-flag: --with-mpi=no
diff --git a/configure.ac b/configure.ac
index ce24e3f2c2..65e3c0fe73 100644
--- a/configure.ac
+++ b/configure.ac
@@ -54,12 +54,6 @@ AS_IF([test x${CRAYPE_VERSION:+yes} = "xyes"],[
# Process user optons.
-AC_ARG_ENABLE([mixed-mode],
- [AS_HELP_STRING([--enable-mixed-mode],
- [Build using mixed mode. Enables both 64-bit and 32-bit reals in Fortran. This option will be ignored if --disable-fortran-flag-setting is also given.])])
-AS_IF([test ${enable_mixed_mode:-no} = no],
- [enable_mixed_mode=no],
- [enable_mixed_mode=yes])
AC_ARG_WITH([mpi],
[AS_HELP_STRING([--with-mpi],
[Build with MPI support. This option will be ignored if --disable-fortran-flag-setting is also given. (Default yes)])])
@@ -85,12 +79,6 @@ AS_IF([test ${enable_code_coverage:-no} = no],
[enable_code_coverage=no],
[enable_code_coverage=yes])
# individual mixed precision overload macros
-AC_ARG_ENABLE([overload-r4],
- [AS_HELP_STRING([--enable-overload-r4],
- [Enables the OVERLOAD_R4 macro to compile with 4 byte real routine overloads. (Default no)])])
-AS_IF([test ${enable_overload_r4:-no} = yes],
- [enable_overload_r4=yes],
- [enable_overload_r4=no])
AC_ARG_ENABLE([overload-c4],
[AS_HELP_STRING([--enable-overload-c4],
[Enables the OVERLOAD_C4 macro to compile with 4 byte complex routine overloads. (Default no)])])
@@ -117,6 +105,13 @@ AS_IF([test ${enable_deprecated_io:-no} = yes],
[enable_deprecated_io=yes],
[enable_deprecated_io=no])
+AC_ARG_ENABLE([r8-default],
+ [AS_HELP_STRING([--disable-r8-default],
+ [Disables the build from adding the 8 byte default real kind flag during compilation (default no)])])
+AS_IF([test ${enable_r8_default:-yes} = yes],
+ [enable_r8_default=yes],
+ [enable_r8_default=no])
+
# user enabled testing with input files
AC_MSG_CHECKING([whether to enable tests with input files])
AC_ARG_ENABLE([test-input],
@@ -282,7 +277,6 @@ AC_OPENMP()
AC_LANG_POP(Fortran)
# We passed all the tests. Set the required defines.
-AC_DEFINE([use_netCDF], [1], [This is required for the library to build])
if test $with_mpi = yes; then
AC_DEFINE([use_libMPI], [1], [This is required for the library to build])
fi
@@ -291,16 +285,18 @@ fi
if test $enable_deprecated_io = yes; then
#If the test pass, define use_deprecated_io macro and skip it's unit tests
AC_DEFINE([use_deprecated_io], [1], [This is required to use mpp_io and fms_io modules])
+ # this macro was only removed from non-deprecated code so still needs to be set
+ AC_DEFINE([use_netCDF], [1], [This has been removed elsewhere but is still required to build with deprecated io])
AM_CONDITIONAL([SKIP_DEPRECATED_IO_TESTS], true)
else
AM_CONDITIONAL([SKIP_DEPRECATED_IO_TESTS], false)
fi
-# disable mosaic unit tests if FMS is compiled in r4
-if test $enable_mixed_mode = yes ; then
- AM_CONDITIONAL([SKIP_MOSAIC_TESTS], true)
+# Builds with r8 default unless disable flag is given
+if test $enable_r8_default = yes; then
+ AM_CONDITIONAL([SKIP_MOSAIC_TESTS], false)
else
- AM_CONDITIONAL([SKIP_MOSAIC_TESTS], false)
+ AM_CONDITIONAL([SKIP_MOSAIC_TESTS], true)
fi
# Set any required compile flags. This will not be done if the user wants to
@@ -314,20 +310,13 @@ if test $enable_setting_flags = yes; then
# necessary fortran flags.
AC_FC_LINE_LENGTH([unlimited])
- # Will we build with default 64-bit reals in Fortran, or do mixed mode?
- if test $enable_mixed_mode = yes; then
- GX_FC_DEFAULT_REAL_KIND4_FLAG([dnl
- FCFLAGS="$FCFLAGS $FC_DEFAULT_REAL_KIND8_FLAG"])
- AC_DEFINE([OVERLOAD_R4], [1], [Set to overload the R4 Fortran routines])
- AC_DEFINE([OVERLOAD_R8], [1], [Set to overload the R8 Fortran routines])
- else
+ # Builds with r8 default unless disable flag is given
+ if test $enable_r8_default = yes; then
GX_FC_DEFAULT_REAL_KIND8_FLAG([dnl
FCFLAGS="$FCFLAGS $FC_DEFAULT_REAL_KIND8_FLAG"])
fi
+
# individual mixed precision overloads
- if test $enable_overload_r4 = yes; then
- AC_DEFINE([OVERLOAD_R4], [1], [Set to overload with the R4 Fortran routines])
- fi
if test $enable_overload_c4 = yes; then
AC_DEFINE([OVERLOAD_C4], [1], [Set to overload with the C4 Fortran routines])
fi
diff --git a/diag_manager/diag_axis.F90 b/diag_manager/diag_axis.F90
index 96221aee36..7eb5eaf686 100644
--- a/diag_manager/diag_axis.F90
+++ b/diag_manager/diag_axis.F90
@@ -40,9 +40,7 @@ MODULE diag_axis_mod
USE diag_data_mod, ONLY: diag_axis_type, max_subaxes, max_axes,&
& max_num_axis_sets, max_axis_attributes, debug_diag_manager,&
& first_send_data_call, diag_atttype
-#ifdef use_netCDF
USE netcdf, ONLY: NF90_INT, NF90_FLOAT, NF90_CHAR
-#endif
IMPLICIT NONE
diff --git a/diag_manager/diag_data.F90 b/diag_manager/diag_data.F90
index a1f5947098..1f443ce220 100644
--- a/diag_manager/diag_data.F90
+++ b/diag_manager/diag_data.F90
@@ -53,10 +53,8 @@ MODULE diag_data_mod
USE fms_mod, ONLY: WARNING, write_version_number
USE fms_diag_bbox_mod, ONLY: fmsDiagIbounds_type
-#ifdef use_netCDF
! NF90_FILL_REAL has value of 9.9692099683868690e+36.
USE netcdf, ONLY: NF_FILL_REAL => NF90_FILL_REAL
-#endif
use fms2_io_mod
IMPLICIT NONE
@@ -335,13 +333,9 @@ MODULE diag_data_mod
!
-#ifdef use_netCDF
REAL :: FILL_VALUE = NF_FILL_REAL !< Fill value used. Value will be NF90_FILL_REAL if using the
!! netCDF module, otherwise will be 9.9692099683868690e+36.
! from file /usr/local/include/netcdf.inc
-#else
- REAL :: FILL_VALUE = 9.9692099683868690e+36
-#endif
INTEGER :: pack_size = 1 !< 1 for double and 2 for float
diff --git a/diag_manager/diag_manager.F90 b/diag_manager/diag_manager.F90
index 92fdf0e122..55048c04f8 100644
--- a/diag_manager/diag_manager.F90
+++ b/diag_manager/diag_manager.F90
@@ -240,9 +240,7 @@ MODULE diag_manager_mod
USE fms_diag_fieldbuff_update_mod, ONLY: fieldbuff_update, fieldbuff_copy_missvals, &
& fieldbuff_copy_fieldvals
-#ifdef use_netCDF
USE netcdf, ONLY: NF90_INT, NF90_FLOAT, NF90_CHAR
-#endif
!----------
!ug support
diff --git a/diag_manager/diag_output.F90 b/diag_manager/diag_output.F90
index de8605d34b..095f8e659c 100644
--- a/diag_manager/diag_output.F90
+++ b/diag_manager/diag_output.F90
@@ -44,9 +44,7 @@ MODULE diag_output_mod
USE time_manager_mod, ONLY: get_calendar_type, valid_calendar_types
USE fms_mod, ONLY: error_mesg, mpp_pe, write_version_number, fms_error_handler, FATAL, note
-#ifdef use_netCDF
USE netcdf, ONLY: NF90_INT, NF90_FLOAT, NF90_CHAR
-#endif
use mpp_domains_mod, only: mpp_get_UG_io_domain
use mpp_domains_mod, only: mpp_get_UG_domain_npes
diff --git a/diag_manager/diag_util.F90 b/diag_manager/diag_util.F90
index 9956c2d9c4..23191e62e8 100644
--- a/diag_manager/diag_util.F90
+++ b/diag_manager/diag_util.F90
@@ -72,9 +72,7 @@ MODULE diag_util_mod
USE constants_mod, ONLY: SECONDS_PER_DAY, SECONDS_PER_HOUR, SECONDS_PER_MINUTE
USE fms2_io_mod
USE fms_diag_bbox_mod, ONLY: fmsDiagIbounds_type
-#ifdef use_netCDF
USE netcdf, ONLY: NF90_CHAR
-#endif
IMPLICIT NONE
PRIVATE
diff --git a/include/fms_platform.h b/include/fms_platform.h
index 3b86ac054f..fdd0bd4e41 100644
--- a/include/fms_platform.h
+++ b/include/fms_platform.h
@@ -46,32 +46,8 @@ use,intrinsic :: iso_c_binding, only: c_double,c_float,c_int64_t, &
!DEC$ MESSAGE:'Using 8-byte addressing'
#endif
-
-!Control "pure" functions.
-#ifdef NO_F95
-#define _PURE
-!DEC$ MESSAGE:'Not using pure routines.'
-#else
-#define _PURE pure
-!DEC$ MESSAGE:'Using pure routines.'
-#endif
-
-
-!Control array members of derived types.
-#ifdef NO_F2000
-#define _ALLOCATABLE pointer
-#define _NULL =>null()
-#define _ALLOCATED associated
-!DEC$ MESSAGE:'Using pointer derived type array members.'
-#else
-#define _ALLOCATABLE allocatable
-#define _NULL
-#define _ALLOCATED allocated
-!DEC$ MESSAGE:'Using allocatable derived type array members.'
-#endif
-
-
-!Control use of cray pointers.
+!Control use of cray pointers within mpp_peset
+!Other cray pointer usage in mpp routines is compiled regardless
#ifdef NO_CRAY_POINTERS
#undef use_CRI_pointers
!DEC$ MESSAGE:'Not using cray pointers.'
@@ -80,40 +56,13 @@ use,intrinsic :: iso_c_binding, only: c_double,c_float,c_int64_t, &
!DEC$ MESSAGE:'Using cray pointers.'
#endif
-
-!Control size of integers that will hold address values.
-!Appears for legacy reasons, but seems rather dangerous.
-#ifdef _32bits
-#define POINTER_KIND 4
-!DEC$ MESSAGE:'Using 4-byte addressing'
-#endif
-
-
-!If you do not want to use 64-bit integers.
+ !If you do not want to use 64-bit integers.
#ifdef no_8byte_integers
#define LONG_KIND INT_KIND
#endif
-
-!If you do not want to use 32-bit floats.
-#ifdef no_4byte_reals
-#define FLOAT_KIND DOUBLE_KIND
-#define NF_GET_VAR_REAL nf_get_var_double
-#define NF_GET_VARA_REAL nf_get_vara_double
-#define NF_GET_ATT_REAL nf_get_att_double
-#undef OVERLOAD_R4
-#undef OVERLOAD_C4
-#endif
-
-
!If you want to use quad-precision.
-! The NO_QUAD_PRECISION macro will be deprecated and removed at some future time.
-! Model code will rely solely upon the ENABLE_QUAD_PRECISION macro thereafer.
-#if defined(ENABLE_QUAD_PRECISION)
-#undef NO_QUAD_PRECISION
-#else
-#define NO_QUAD_PRECISION
-#undef QUAD_KIND
+#ifndef ENABLE_QUAD_PRECISION
#define QUAD_KIND DOUBLE_KIND
#endif
diff --git a/mosaic/read_mosaic.c b/mosaic/read_mosaic.c
index cbbe4f4b58..9fafad1f2b 100644
--- a/mosaic/read_mosaic.c
+++ b/mosaic/read_mosaic.c
@@ -23,9 +23,7 @@
#include "read_mosaic.h"
#include "constant.h"
#include "mosaic_util.h"
-#ifdef use_netCDF
#include
-#endif
/** \file
* \ingroup mosaic
diff --git a/test_fms/mosaic2/Makefile.am b/test_fms/mosaic2/Makefile.am
index ca0b4804d4..7f2d6143ab 100644
--- a/test_fms/mosaic2/Makefile.am
+++ b/test_fms/mosaic2/Makefile.am
@@ -47,13 +47,13 @@ test_grid2_r8_CPPFLAGS =-DTEST_MOS_KIND_=8 $(AM_CPPFLAGS)
# These files are also included in the distribution.
EXTRA_DIST = test_mosaic2.sh
-if SKIP_MOSAIC_TESTS
-TESTS_ENVIRONMENT = SKIP_TESTS="test_mosaic2.1 test_mosaic2.2 test_mosaic2.3 test_mosaic2.4"
-endif
-
# Run the test program.
TESTS = test_mosaic2.sh
+if SKIP_MOSAIC_TESTS
+ TESTS_ENVIRONMENT = SKIP_TESTS="test_mosaic2.1 test_mosaic2.2 test_mosaic2.3 test_mosaic2.4"
+endif
+
TEST_EXTENSIONS = .sh
SH_LOG_DRIVER = env AM_TAP_AWK='$(AWK)' $(SHELL) \
$(abs_top_srcdir)/test_fms/tap-driver.sh
diff --git a/time_interp/include/time_interp_external.inc b/time_interp/include/time_interp_external.inc
deleted file mode 100644
index 7c446f4c52..0000000000
--- a/time_interp/include/time_interp_external.inc
+++ /dev/null
@@ -1,1424 +0,0 @@
-!***********************************************************************
-!* 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 .
-!***********************************************************************
-!> @defgroup time_interp_external_mod time_interp_external_mod
-!> @ingroup time_interp
-!> @brief Perform I/O and time interpolation of external fields (contained in a file).
-!> @author M.J. Harrison
-!!
-!! Perform I/O and time interpolation for external fields.
-!! Uses udunits library to calculate calendar dates and
-!! convert units. Allows for reading data decomposed across
-!! model horizontal grid using optional domain2d argument
-!!
-!! data are defined over data domain for domain2d data
-!! (halo values are NOT updated by this module)
-
-!> @addtogroup time_interp_external_mod
-!> @{
-module time_interp_external_mod
-#ifdef use_deprecated_io
-#include
-!
-!M.J. Harrison
-!
-!Harper Simmons
-!
-!
-!
-
-!
-!
-!
-!
-!
-!
-! size of record dimension for internal buffer. This is useful for tuning i/o performance
-! particularly for large datasets (e.g. daily flux fields)
-!
-!
-
- use fms_mod, only : write_version_number
- use mpp_mod, only : mpp_error,FATAL,WARNING,mpp_pe, stdout, stdlog, NOTE
- use mpp_mod, only : input_nml_file
- use mpp_io_mod, only : mpp_open, mpp_get_atts, mpp_get_info, MPP_NETCDF, MPP_MULTI, MPP_SINGLE,&
- mpp_get_times, MPP_RDONLY, MPP_ASCII, default_axis,axistype,fieldtype,atttype, &
- mpp_get_axes, mpp_get_fields, mpp_read, default_field, mpp_close, &
- mpp_get_tavg_info, validtype, mpp_is_valid, mpp_get_file_name
- use time_manager_mod, only : time_type, get_date, set_date, operator ( >= ) , operator ( + ) , days_in_month, &
- operator( - ), operator ( / ) , days_in_year, increment_time, &
- set_time, get_time, operator( > ), get_calendar_type, NO_CALENDAR
- use get_cal_time_mod, only : get_cal_time
- use mpp_domains_mod, only : domain2d, mpp_get_compute_domain, mpp_get_data_domain, &
- mpp_get_global_domain, NULL_DOMAIN2D
- use time_interp_mod, only : time_interp, time_interp_init
- use axis_utils_mod, only : get_axis_cart, get_axis_modulo, get_axis_modulo_times
- use fms_mod, only : lowercase, open_namelist_file, check_nml_error, close_file
- use platform_mod, only: r8_kind
- use horiz_interp_mod, only : horiz_interp, horiz_interp_type
-
- implicit none
- private
-
-! Include variable "version" to be written to log file.
-#include
-
- integer, parameter, public :: NO_REGION=0, INSIDE_REGION=1, OUTSIDE_REGION=2
- integer, parameter, private :: modulo_year= 0001
- integer, parameter, private :: LINEAR_TIME_INTERP = 1 ! not used currently
- integer, parameter, public :: SUCCESS = 0, ERR_FIELD_NOT_FOUND = 1
- integer, private :: max_fields = 100, max_files= 40
- integer, private :: num_fields = 0, num_files=0
- ! denotes time intervals in file (interpreted from metadata)
- integer, private :: num_io_buffers = 2 ! set -1 to read all records from disk into memory
- logical, private :: module_initialized = .false.
- logical, private :: debug_this_module = .false.
-
- public init_external_field, time_interp_external, time_interp_external_init, &
- time_interp_external_exit, get_external_field_size, get_time_axis, get_external_field_missing
- public set_override_region, reset_src_data_region, get_external_field_axes
-
- private find_buf_index,&
- set_time_modulo
-
- !> @}
-
- !> @ingroup time_interp_external_mod
- type, private :: ext_fieldtype
- integer :: unit ! keep unit open when not reading all records
- character(len=128) :: name, units
- integer :: siz(4), ndim
- type(domain2d) :: domain
- type(axistype) :: axes(4)
- type(time_type), dimension(:), pointer :: time =>NULL() ! midpoint of time interval
- type(time_type), dimension(:), pointer :: start_time =>NULL(), end_time =>NULL()
- type(fieldtype) :: field ! mpp_io type
- type(time_type), dimension(:), pointer :: period =>NULL()
- logical :: modulo_time ! denote climatological time axis
- real, dimension(:,:,:,:), pointer :: data =>NULL() ! defined over data domain or global domain
- logical, dimension(:,:,:,:), pointer :: mask =>NULL() ! defined over data domain or global domain
- integer, dimension(:), pointer :: ibuf =>NULL() ! record numbers associated with buffers
- real, dimension(:,:,:,:), pointer :: src_data =>NULL() ! input data buffer
- type(validtype) :: valid ! data validator
- integer :: nbuf
- logical :: domain_present
- real(DOUBLE_KIND) :: slope, intercept
- integer :: isc,iec,jsc,jec
- type(time_type) :: modulo_time_beg, modulo_time_end
- logical :: have_modulo_times, correct_leap_year_inconsistency
- integer :: region_type
- integer :: is_region, ie_region, js_region, je_region
- integer :: is_src, ie_src, js_src, je_src
- integer :: tdim
- integer :: numwindows
- logical, dimension(:,:), pointer :: need_compute=>NULL()
- real :: missing ! missing value
- end type ext_fieldtype
-
- !> @ingroup time_interp_external_mod
- type, private :: filetype
- character(len=128) :: filename = ''
- integer :: unit = -1
- end type filetype
-
- !> Provide data from external file interpolated to current model time.
- !! Data may be local to current processor or global, depending on
- !! "init_external_field" flags. Uses @ref mpp_io_mod for I/O.
- !!
- !! @param index index of external field from previous call to init_external_field
- !! @param time target time for data
- !! @param [inout] data global or local data array
- !! @param interp time_interp_external defined interpolation method (optional). Currently
- !! this module only supports LINEAR_TIME_INTERP.
- !! @param verbose verbose flag for debugging (optional).
- !!
- !> @ingroup time_interp_external_mod
- interface time_interp_external
- module procedure time_interp_external_0d
- module procedure time_interp_external_2d
- module procedure time_interp_external_3d
- end interface
-
- !> @addtogroup time_interp_external_mod
- !> @{
-
- integer :: outunit
-
- type(ext_fieldtype), save, private, pointer :: field(:) => NULL()
- type(filetype), save, private, pointer :: opened_files(:) => NULL()
-!Balaji: really should use field%missing
- integer, private, parameter :: dk = DOUBLE_KIND
- real(DOUBLE_KIND), private, parameter :: time_interp_missing=-1e99_dk
- contains
-
-!
-!
-!
-! Initialize the time_interp_external module
-!
-!
- subroutine time_interp_external_init()
-
- integer :: ioun, io_status, logunit, ierr
-
- namelist /time_interp_external_nml/ num_io_buffers, debug_this_module, &
- max_fields, max_files
-
- ! open and read namelist
-
- if(module_initialized) return
-
- logunit = stdlog()
- outunit = stdout()
- call write_version_number("TIME_INTERP_EXTERNAL_MOD", version)
-
-#ifdef INTERNAL_FILE_NML
- read (input_nml_file, time_interp_external_nml, iostat=io_status)
- ierr = check_nml_error(io_status, 'time_interp_external_nml')
-#else
- ioun = open_namelist_file ()
- ierr=1; do while (ierr /= 0)
- read (ioun, nml=time_interp_external_nml, iostat=io_status, end=10)
- ierr = check_nml_error(io_status, 'time_interp_external_nml')
- enddo
-10 call close_file (ioun)
-#endif
-
- write(logunit,time_interp_external_nml)
- call realloc_fields(max_fields)
- call realloc_files(max_files)
-
- module_initialized = .true.
-
- call time_interp_init()
-
- return
-
- end subroutine time_interp_external_init
-! NAME="time_interp_external_init"
-
-
-!
-!
-!
-! initialize an external field. Buffer "num_io_buffers" (default=2) in memory to reduce memory allocations.
-! distributed reads are supported using the optional "domain" flag.
-! Units conversion via the optional "desired_units" flag using udunits_mod.
-!
-! Return integer id of field for future calls to time_interp_external.
-!
-!
-!
-!
-! filename
-!
-!
-! fieldname (in file)
-!
-!
-! mpp_io flag for format of file (optional). Currently only "MPP_NETCDF" supported
-!
-!
-! mpp_io flag for threading (optional). "MPP_SINGLE" means root pe reads global field and distributes to other PEs
-! "MPP_MULTI" means all PEs read data
-!
-!
-! domain flag (optional)
-!
-!
-! Target units for data (optional), e.g. convert from deg_K to deg_C.
-! Failure to convert using udunits will result in failure of this module.
-!
-!
-! verbose flag for debugging (optional).
-!
-!
-! MPP_IO axistype array for grid centers ordered X-Y-Z-T (optional).
-!
-!
-! array of axis lengths ordered X-Y-Z-T (optional).
-!
-
-
- !> Initialize an external field. Buffer "num_io_buffers" (default=2) in memory to reduce memory allocations.
- !! distributed reads are supported using the optional "domain" flag.
- !! Units conversion via the optional "desired_units" flag using udunits_mod.
- !!
- !> @return integer id of field for future calls to time_interp_external.
- !> @param file filename
- !> @param fieldname fieldname (in file)
- !> @param format mpp_io flag for format of file(optional). Currently only "MPP_NETCDF" supported
- !> @param threading mpp_io flag for threading (optional). "MPP_SINGLE" means root pe reads
- !! global field and distributes to other PEs. "MPP_MULTI" means all PEs read data
- !> @param domain domain flag (optional)
- !> @param desired_units Target units for data (optional), e.g. convert from deg_K to deg_C.
- !! Failure to convert using udunits will result in failure of this module.
- !> @param verbose verbose flag for debugging (optional).
- !> @param [out] axis_names List of axis names (optional).
- !> @param [inout] axis_sizes array of axis lengths ordered X-Y-Z-T (optional).
- function init_external_field(file,fieldname,format,threading,domain,desired_units,&
- verbose,axis_centers,axis_sizes,override,correct_leap_year_inconsistency,&
- permit_calendar_conversion,use_comp_domain,ierr, nwindows, ignore_axis_atts )
-
- character(len=*), intent(in) :: file,fieldname
- integer, intent(in), optional :: format, threading
- logical, intent(in), optional :: verbose
- character(len=*), intent(in), optional :: desired_units
- type(domain2d), intent(in), optional :: domain
- type(axistype), intent(inout), optional :: axis_centers(4)
- integer, intent(inout), optional :: axis_sizes(4)
- logical, intent(in), optional :: override, correct_leap_year_inconsistency,&
- permit_calendar_conversion,use_comp_domain
- integer, intent(out), optional :: ierr
- integer, intent(in), optional :: nwindows
- logical, optional :: ignore_axis_atts
- real :: missing
-
- integer :: init_external_field
-
- type(fieldtype), dimension(:), allocatable :: flds
- type(axistype), dimension(:), allocatable :: axes, fld_axes
- type(axistype) :: time_axis
- type(atttype), allocatable, dimension(:) :: global_atts
-
- real(DOUBLE_KIND) :: slope, intercept
- integer :: form, thread, fset, unit,ndim,nvar,natt,ntime,i,j
- integer :: iscomp,iecomp,jscomp,jecomp,isglobal,ieglobal,jsglobal,jeglobal
- integer :: isdata,iedata,jsdata,jedata, dxsize, dysize,dxsize_max,dysize_max
- logical :: verb, transpose_xy,use_comp_domain1
- real, dimension(:), allocatable :: tstamp, tstart, tend, tavg
- character(len=1) :: cart
- character(len=1), dimension(4) :: cart_dir
- character(len=128) :: units, fld_units
- character(len=128) :: name, msg, calendar_type, timebeg, timeend
- integer :: siz(4), siz_in(4), gxsize, gysize,gxsize_max, gysize_max
- type(time_type) :: tdiff
- integer :: yr, mon, day, hr, minu, sec
- integer :: len, nfile, nfields_orig, nbuf, nx,ny
- integer :: numwindows
- logical :: ignore_axatts
-
-
- if (.not. module_initialized) call mpp_error(FATAL,'Must call time_interp_external_init first')
- if(present(ierr)) ierr = SUCCESS
- ignore_axatts=.false.
- cart_dir(1)='X';cart_dir(2)='Y';cart_dir(3)='Z';cart_dir(4)='T'
- if(present(ignore_axis_atts)) ignore_axatts = ignore_axis_atts
- use_comp_domain1 = .false.
- if(PRESENT(use_comp_domain)) use_comp_domain1 = use_comp_domain
- form=MPP_NETCDF
- if (PRESENT(format)) form = format
- thread = MPP_MULTI
- if (PRESENT(threading)) thread = threading
- fset = MPP_SINGLE
- verb=.false.
- if (PRESENT(verbose)) verb=verbose
- if (debug_this_module) verb = .true.
- numwindows = 1
- if(present(nwindows)) numwindows = nwindows
-
- units = 'same'
- if (PRESENT(desired_units)) then
- units = desired_units
- call mpp_error(FATAL,'==> Unit conversion via time_interp_external &
- &has been temporarily deprecated. Previous versions of&
- &this module used udunits_mod to perform unit conversion.&
- & Udunits_mod is in the process of being replaced since &
- &there were portability issues associated with this code.&
- & Please remove the desired_units argument from calls to &
- &this routine.')
- endif
- nfile = 0
- do i=1,num_files
- if(trim(opened_files(i)%filename) == trim(file)) then
- nfile = i
- exit ! file is already opened
- endif
- enddo
- if(nfile == 0) then
- call mpp_open(unit,trim(file),MPP_RDONLY,form,threading=thread,&
- fileset=fset)
- num_files = num_files + 1
- if(num_files > max_files) then ! not enough space in the file table, reallocate it
- !--- z1l: For the case of multiple thread, realoc_files will cause memory leak.
- !--- If multiple threads are working on file A. One of the thread finished first and
- !--- begin to work on file B, the realloc_files will cause problem for
- !--- other threads are working on the file A.
- ! call realloc_files(2*size(opened_files))
- call mpp_error(FATAL, "time_interp_external: num_files is greater than max_files, "// &
- "increase time_interp_external_nml max_files")
- endif
- opened_files(num_files)%filename = trim(file)
- opened_files(num_files)%unit = unit
- else
- unit = opened_files(nfile)%unit
- endif
-
- call mpp_get_info(unit,ndim,nvar,natt,ntime)
-
- if (ntime < 1) then
- write(msg,'(a15,a,a58)') 'external field ',trim(fieldname),&
- ' does not have an associated record dimension (REQUIRED) '
- call mpp_error(FATAL,trim(msg))
- endif
- allocate(global_atts(natt))
- call mpp_get_atts(unit, global_atts)
- allocate(axes(ndim))
- call mpp_get_axes(unit, axes, time_axis)
- allocate(flds(nvar))
- call mpp_get_fields(unit,flds)
- allocate(tstamp(ntime),tstart(ntime),tend(ntime),tavg(ntime))
- call mpp_get_times(unit,tstamp)
- transpose_xy = .false.
- isdata=1; iedata=1; jsdata=1; jedata=1
- gxsize=1; gysize=1
- siz_in = 1
-
- if (PRESENT(domain)) then
- call mpp_get_compute_domain(domain,iscomp,iecomp,jscomp,jecomp)
- nx = iecomp-iscomp+1; ny = jecomp-jscomp+1
- call mpp_get_data_domain(domain,isdata,iedata,jsdata,jedata,dxsize,dxsize_max,dysize,dysize_max)
- call mpp_get_global_domain(domain,isglobal,ieglobal,jsglobal,jeglobal,gxsize,gxsize_max,gysize,gysize_max)
- elseif(use_comp_domain1) then
- call mpp_error(FATAL,"init_external_field:"//&
- " use_comp_domain=true but domain is not present")
- endif
-
- init_external_field = -1
- nfields_orig = num_fields
-
- do i=1,nvar
- call mpp_get_atts(flds(i),name=name,units=fld_units,ndim=ndim,siz=siz_in)
- call mpp_get_tavg_info(unit,flds(i),flds,tstamp,tstart,tend,tavg)
- call mpp_get_atts(flds(i),missing=missing)
- ! why does it convert case of the field name?
- if (trim(lowercase(name)) /= trim(lowercase(fieldname))) cycle
-
- if (verb) write(outunit,*) 'found field ',trim(fieldname), ' in file !!'
- num_fields = num_fields + 1
- if(num_fields > max_fields) then
- !--- z1l: For the case of multiple thread, realoc_fields will cause memory leak.
- !--- If multiple threads are working on field A. One of the thread finished first and
- !--- begin to work on field B, the realloc_files will cause problem for
- !--- other threads are working on the field A.
- !call realloc_fields(size(field)*2)
- call mpp_error(FATAL, "time_interp_external: num_fields is greater than max_fields, "// &
- "increase time_interp_external_nml max_fields")
- endif
-
- init_external_field = num_fields
- field(num_fields)%unit = unit
- field(num_fields)%name = trim(name)
- field(num_fields)%units = trim(fld_units)
- field(num_fields)%field = flds(i)
- field(num_fields)%isc = 1
- field(num_fields)%iec = 1
- field(num_fields)%jsc = 1
- field(num_fields)%jec = 1
- field(num_fields)%region_type = NO_REGION
- field(num_fields)%is_region = 0
- field(num_fields)%ie_region = -1
- field(num_fields)%js_region = 0
- field(num_fields)%je_region = -1
- if (PRESENT(domain)) then
- field(num_fields)%domain_present = .true.
- field(num_fields)%domain = domain
- field(num_fields)%isc=iscomp;field(num_fields)%iec = iecomp
- field(num_fields)%jsc=jscomp;field(num_fields)%jec = jecomp
- else
- field(num_fields)%domain_present = .false.
- endif
-
- call mpp_get_atts(flds(i),valid=field(num_fields)%valid )
- allocate(fld_axes(ndim))
- call mpp_get_atts(flds(i),axes=fld_axes)
- if (ndim > 4) call mpp_error(FATAL, &
- 'invalid array rank <=4d fields supported')
- field(num_fields)%siz = 1
- field(num_fields)%ndim = ndim
- field(num_fields)%tdim = 4
- field(num_fields)%missing = missing
- do j=1,field(num_fields)%ndim
- cart = 'N'
- call get_axis_cart(fld_axes(j), cart)
- call mpp_get_atts(fld_axes(j),len=len)
- if (cart == 'N' .and. .not. ignore_axatts) then
- write(msg,'(a,"/",a)') trim(file),trim(fieldname)
- call mpp_error(FATAL,'file/field '//trim(msg)// &
- ' couldnt recognize axis atts in time_interp_external')
- else if (cart == 'N' .and. ignore_axatts) then
- cart = cart_dir(j)
- endif
- select case (cart)
- case ('X')
- if (j.eq.2) transpose_xy = .true.
- if (.not.PRESENT(domain) .and. .not.PRESENT(override)) then
- isdata=1;iedata=len
- iscomp=1;iecomp=len
- gxsize = len
- dxsize = len
- field(num_fields)%isc=iscomp;field(num_fields)%iec=iecomp
- elseif (PRESENT(override)) then
- gxsize = len
- if (PRESENT(axis_sizes)) axis_sizes(1) = len
- endif
- field(num_fields)%axes(1) = fld_axes(j)
- if(use_comp_domain1) then
- field(num_fields)%siz(1) = nx
- else
- field(num_fields)%siz(1) = dxsize
- endif
- if (len /= gxsize) then
- write(msg,'(a,"/",a)') trim(file),trim(fieldname)
- call mpp_error(FATAL,'time_interp_ext, file/field '//trim(msg)//' x dim doesnt match model')
- endif
- case ('Y')
- field(num_fields)%axes(2) = fld_axes(j)
- if (.not.PRESENT(domain) .and. .not.PRESENT(override)) then
- jsdata=1;jedata=len
- jscomp=1;jecomp=len
- gysize = len
- dysize = len
- field(num_fields)%jsc=jscomp;field(num_fields)%jec=jecomp
- elseif (PRESENT(override)) then
- gysize = len
- if (PRESENT(axis_sizes)) axis_sizes(2) = len
- endif
- if(use_comp_domain1) then
- field(num_fields)%siz(2) = ny
- else
- field(num_fields)%siz(2) = dysize
- endif
- if (len /= gysize) then
- write(msg,'(a,"/",a)') trim(file),trim(fieldname)
- call mpp_error(FATAL,'time_interp_ext, file/field '//trim(msg)//' y dim doesnt match model')
- endif
- case ('Z')
- field(num_fields)%axes(3) = fld_axes(j)
- field(num_fields)%siz(3) = siz_in(3)
- case ('T')
- field(num_fields)%axes(4) = fld_axes(j)
- field(num_fields)%siz(4) = ntime
- field(num_fields)%tdim = j
- end select
- enddo
- siz = field(num_fields)%siz
-
- if (PRESENT(axis_centers)) then
- axis_centers = field(num_fields)%axes
- endif
-
- if (PRESENT(axis_sizes) .and. .not.PRESENT(override)) then
- axis_sizes = field(num_fields)%siz
- endif
-
- deallocate(fld_axes)
- if (verb) write(outunit,'(a,4i6)') 'field x,y,z,t local size= ',siz
- if (verb) write(outunit,*) 'field contains data in units = ',trim(field(num_fields)%units)
- if (transpose_xy) call mpp_error(FATAL,'axis ordering not supported')
- if (num_io_buffers .le. 1) call mpp_error(FATAL,'time_interp_ext:num_io_buffers should be at least 2')
- nbuf = min(num_io_buffers,siz(4))
-
- field(num_fields)%numwindows = numwindows
- allocate(field(num_fields)%need_compute(nbuf, numwindows))
- field(num_fields)%need_compute = .true.
-
- allocate(field(num_fields)%data(isdata:iedata,jsdata:jedata,siz(3),nbuf),&
- field(num_fields)%mask(isdata:iedata,jsdata:jedata,siz(3),nbuf) )
- field(num_fields)%mask = .false.
- field(num_fields)%data = 0.0
- slope=1.0;intercept=0.0
-! if (units /= 'same') call convert_units(trim(field(num_fields)%units),trim(units),slope,intercept)
-! if (verb.and.units /= 'same') then
-! write(outunit,*) 'attempting to convert data to units = ',trim(units)
-! write(outunit,'(a,f8.3,a,f8.3)') 'factor = ',slope,' offset= ',intercept
-! endif
- field(num_fields)%slope = slope
- field(num_fields)%intercept = intercept
- allocate(field(num_fields)%ibuf(nbuf))
- field(num_fields)%ibuf = -1
- field(num_fields)%nbuf = 0 ! initialize buffer number so that first reading fills data(:,:,:,1)
- if(PRESENT(override)) then
- field(num_fields)%is_src = 1
- field(num_fields)%ie_src = gxsize
- field(num_fields)%js_src = 1
- field(num_fields)%je_src = gysize
- allocate(field(num_fields)%src_data(gxsize,gysize,siz(3),nbuf))
- else
- field(num_fields)%is_src = isdata
- field(num_fields)%ie_src = iedata
- field(num_fields)%js_src = jsdata
- field(num_fields)%je_src = jedata
- allocate(field(num_fields)%src_data(isdata:iedata,jsdata:jedata,siz(3),nbuf))
- endif
-
- allocate(field(num_fields)%time(ntime))
- allocate(field(num_fields)%period(ntime))
- allocate(field(num_fields)%start_time(ntime))
- allocate(field(num_fields)%end_time(ntime))
-
- call mpp_get_atts(time_axis,units=units,calendar=calendar_type)
- do j=1,ntime
- field(num_fields)%time(j) = get_cal_time(tstamp(j),trim(units),trim(calendar_type), &
- & permit_calendar_conversion)
- field(num_fields)%start_time(j) = get_cal_time(tstart(j),trim(units),trim(calendar_type), &
- & permit_calendar_conversion)
- field(num_fields)%end_time(j) = get_cal_time( tend(j),trim(units),trim(calendar_type), &
- & permit_calendar_conversion)
- enddo
-
- if (field(num_fields)%modulo_time) then
- call set_time_modulo(field(num_fields)%Time)
- call set_time_modulo(field(num_fields)%start_time)
- call set_time_modulo(field(num_fields)%end_time)
- endif
-
- if(present(correct_leap_year_inconsistency)) then
- field(num_fields)%correct_leap_year_inconsistency = correct_leap_year_inconsistency
- else
- field(num_fields)%correct_leap_year_inconsistency = .false.
- endif
-
- if(get_axis_modulo_times(time_axis, timebeg, timeend)) then
- if(get_calendar_type() == NO_CALENDAR) then
- field(num_fields)%modulo_time_beg = set_time(timebeg)
- field(num_fields)%modulo_time_end = set_time(timeend)
- else
- field(num_fields)%modulo_time_beg = set_date(timebeg)
- field(num_fields)%modulo_time_end = set_date(timeend)
- endif
- field(num_fields)%have_modulo_times = .true.
- else
- field(num_fields)%have_modulo_times = .false.
- endif
- if(ntime == 1) then
- call mpp_error(NOTE, 'time_interp_external_mod: file '//trim(file)//' has only one time level')
- else
- do j= 1, ntime
- field(num_fields)%period(j) = field(num_fields)%end_time(j)-field(num_fields)%start_time(j)
- if (field(num_fields)%period(j) > set_time(0,0)) then
- call get_time(field(num_fields)%period(j), sec, day)
- sec = sec/2+mod(day,2)*43200
- day = day/2
- field(num_fields)%time(j) = field(num_fields)%start_time(j)+&
- set_time(sec,day)
- else
- if (j > 1 .and. j < ntime) then
- tdiff = field(num_fields)%time(j+1) - field(num_fields)%time(j-1)
- call get_time(tdiff, sec, day)
- sec = sec/2+mod(day,2)*43200
- day = day/2
- field(num_fields)%period(j) = set_time(sec,day)
- sec = sec/2+mod(day,2)*43200
- day = day/2
- field(num_fields)%start_time(j) = field(num_fields)%time(j) - set_time(sec,day)
- field(num_fields)%end_time(j) = field(num_fields)%time(j) + set_time(sec,day)
- elseif ( j == 1) then
- tdiff = field(num_fields)%time(2) - field(num_fields)%time(1)
- call get_time(tdiff, sec, day)
- field(num_fields)%period(j) = set_time(sec,day)
- sec = sec/2+mod(day,2)*43200
- day = day/2
- field(num_fields)%start_time(j) = field(num_fields)%time(j) - set_time(sec,day)
- field(num_fields)%end_time(j) = field(num_fields)%time(j) + set_time(sec,day)
- else
- tdiff = field(num_fields)%time(ntime) - field(num_fields)%time(ntime-1)
- call get_time(tdiff, sec, day)
- field(num_fields)%period(j) = set_time(sec,day)
- sec = sec/2+mod(day,2)*43200
- day = day/2
- field(num_fields)%start_time(j) = field(num_fields)%time(j) - set_time(sec,day)
- field(num_fields)%end_time(j) = field(num_fields)%time(j) + set_time(sec,day)
- endif
- endif
- enddo
- endif
-
- do j=1,ntime-1
- if (field(num_fields)%time(j) >= field(num_fields)%time(j+1)) then
- write(msg,'(A,i20)') "times not monotonically increasing. Filename: " &
- //TRIM(file)//" field: "//TRIM(fieldname)//" timeslice: ", j
- call mpp_error(FATAL, TRIM(msg))
- endif
- enddo
-
- field(num_fields)%modulo_time = get_axis_modulo(time_axis)
-
- if (verb) then
- if (field(num_fields)%modulo_time) write(outunit,*) 'data are being treated as modulo in time'
- do j= 1, ntime
- write(outunit,*) 'time index, ', j
- call get_date(field(num_fields)%start_time(j),yr,mon,day,hr,minu,sec)
- write(outunit,'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') &
- 'start time: yyyy/mm/dd hh:mm:ss= ',yr,'/',mon,'/',day,hr,':',minu,':',sec
- call get_date(field(num_fields)%time(j),yr,mon,day,hr,minu,sec)
- write(outunit,'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') &
- 'mid time: yyyy/mm/dd hh:mm:ss= ',yr,'/',mon,'/',day,hr,':',minu,':',sec
- call get_date(field(num_fields)%end_time(j),yr,mon,day,hr,minu,sec)
- write(outunit,'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') &
- 'end time: yyyy/mm/dd hh:mm:ss= ',yr,'/',mon,'/',day,hr,':',minu,':',sec
- enddo
- end if
-
- enddo
-
- if (num_fields == nfields_orig) then
- if (present(ierr)) then
- ierr = ERR_FIELD_NOT_FOUND
- else
- call mpp_error(FATAL,'external field "'//trim(fieldname)//'" not found in file "'//trim(file)//'"')
- endif
- endif
-
- deallocate(global_atts)
- deallocate(axes)
- deallocate(flds)
- deallocate(tstamp, tstart, tend, tavg)
-
- return
-
- end function init_external_field
-
-! NAME="init_external_field"
-
-
- !> @brief 2D time interpolation for @ref time_interp_external
- subroutine time_interp_external_2d(index, time, data_in, interp, verbose,horz_interp, mask_out, &
- is_in, ie_in, js_in, je_in, window_id)
-
- integer, intent(in) :: index
- type(time_type), intent(in) :: time
- real, dimension(:,:), intent(inout) :: data_in
- integer, intent(in), optional :: interp
- logical, intent(in), optional :: verbose
- type(horiz_interp_type),intent(in), optional :: horz_interp
- logical, dimension(:,:), intent(out), optional :: mask_out ! set to true where output data is valid
- integer, intent(in), optional :: is_in, ie_in, js_in, je_in
- integer, intent(in), optional :: window_id
-
- real , dimension(size(data_in,1), size(data_in,2), 1) :: data_out
- logical, dimension(size(data_in,1), size(data_in,2), 1) :: mask3d
-
- data_out(:,:,1) = data_in(:,:) ! fill initial values for the portions of array that are not touched by 3d routine
- call time_interp_external_3d(index, time, data_out, interp, verbose, horz_interp, mask3d, &
- is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id)
- data_in(:,:) = data_out(:,:,1)
- if (PRESENT(mask_out)) mask_out(:,:) = mask3d(:,:,1)
-
- return
- end subroutine time_interp_external_2d
-
-!
-!
-!
-! Provide data from external file interpolated to current model time.
-! Data may be local to current processor or global, depending on
-! "init_external_field" flags.
-!
-!
-!
-! index of external field from previous call to init_external_field
-!
-!
-! target time for data
-!
-!
-! global or local data array
-!
-!
-! time_interp_external defined interpolation method (optional). Currently this module only supports
-! LINEAR_TIME_INTERP.
-!
-!
-! verbose flag for debugging (optional).
-!
-
- !> @brief 3D time interpolation for @ref time_interp_external
- subroutine time_interp_external_3d(index, time, data, interp,verbose,horz_interp, mask_out, is_in, ie_in, &
- & js_in, je_in, window_id)
-
- integer, intent(in) :: index
- type(time_type), intent(in) :: time
- real, dimension(:,:,:), intent(inout) :: data
- integer, intent(in), optional :: interp
- logical, intent(in), optional :: verbose
- type(horiz_interp_type), intent(in), optional :: horz_interp
- logical, dimension(:,:,:), intent(out), optional :: mask_out ! set to true where output data is valid
- integer, intent(in), optional :: is_in, ie_in, js_in, je_in
- integer, intent(in), optional :: window_id
-
- integer :: nx, ny, nz, interp_method, t1, t2
- integer :: i1, i2, isc, iec, jsc, jec, mod_time
- integer :: yy, mm, dd, hh, min, ss
- character(len=256) :: err_msg, filename
-
- integer :: isw, iew, jsw, jew, nxw, nyw
- ! these are boundaries of the updated portion of the "data" argument
- ! they are calculated using sizes of the "data" and isc,iec,jsc,jsc
- ! fileds from respective input field, to center the updated portion
- ! in the output array
-
- real :: w1,w2
- logical :: verb
- character(len=16) :: message1, message2
-
- nx = size(data,1)
- ny = size(data,2)
- nz = size(data,3)
-
- interp_method = LINEAR_TIME_INTERP
- if (PRESENT(interp)) interp_method = interp
- verb=.false.
- if (PRESENT(verbose)) verb=verbose
- if (debug_this_module) verb = .true.
-
- if (index < 1.or.index > num_fields) &
- call mpp_error(FATAL, &
- & 'invalid index in call to time_interp_ext -- field was not initialized or failed to initialize')
-
- isc=field(index)%isc;iec=field(index)%iec
- jsc=field(index)%jsc;jec=field(index)%jec
-
- if( field(index)%numwindows == 1 ) then
- nxw = iec-isc+1
- nyw = jec-jsc+1
- else
- if( .not. present(is_in) .or. .not. present(ie_in) .or. .not. present(js_in) .or. .not. present(je_in) ) then
- call mpp_error(FATAL, 'time_interp_external: is_in, ie_in, js_in and je_in must be present ' // &
- 'when numwindows > 1, field='//trim(field(index)%name))
- endif
- nxw = ie_in - is_in + 1
- nyw = je_in - js_in + 1
- isc = isc + is_in - 1
- iec = isc + ie_in - is_in
- jsc = jsc + js_in - 1
- jec = jsc + je_in - js_in
- endif
-
- isw = (nx-nxw)/2+1; iew = isw+nxw-1
- jsw = (ny-nyw)/2+1; jew = jsw+nyw-1
-
- if (nx < nxw .or. ny < nyw .or. nz < field(index)%siz(3)) then
- write(message1,'(i6,2i5)') nx,ny,nz
- call mpp_error(FATAL,'field '//trim(field(index)%name)//' Array size mismatch in time_interp_external.'// &
- ' Array "data" is too small. shape(data)='//message1)
- endif
- if(PRESENT(mask_out)) then
- if (size(mask_out,1) /= nx .or. size(mask_out,2) /= ny .or. size(mask_out,3) /= nz) then
- write(message1,'(i6,2i5)') nx,ny,nz
- write(message2,'(i6,2i5)') size(mask_out,1),size(mask_out,2),size(mask_out,3)
- call mpp_error(FATAL,'field '//trim(field(index)%name)//' array size mismatch in time_interp_external.'// &
- ' Shape of array "mask_out" does not match that of array "data".'// &
- ' shape(data)='//message1//' shape(mask_out)='//message2)
- endif
- endif
-
- if (field(index)%siz(4) == 1) then
- ! only one record in the file => time-independent field
- call load_record(field(index),1,horz_interp, is_in, ie_in ,js_in, je_in,window_id)
- i1 = find_buf_index(1,field(index)%ibuf)
- if( field(index)%region_type == NO_REGION ) then
- where(field(index)%mask(isc:iec,jsc:jec,:,i1))
- data(isw:iew,jsw:jew,:) = field(index)%data(isc:iec,jsc:jec,:,i1)
- elsewhere
-! data(isw:iew,jsw:jew,:) = time_interp_missing !field(index)%missing? Balaji
- data(isw:iew,jsw:jew,:) = field(index)%missing
- end where
- else
- where(field(index)%mask(isc:iec,jsc:jec,:,i1))
- data(isw:iew,jsw:jew,:) = field(index)%data(isc:iec,jsc:jec,:,i1)
- end where
- endif
- if(PRESENT(mask_out)) &
- mask_out(isw:iew,jsw:jew,:) = field(index)%mask(isc:iec,jsc:jec,:,i1)
- else
- if(field(index)%have_modulo_times) then
- call time_interp(time,field(index)%modulo_time_beg, field(index)%modulo_time_end, field(index)%time(:), &
- w2, t1, t2, field(index)%correct_leap_year_inconsistency, err_msg=err_msg)
- if(err_msg .NE. '') then
- filename = mpp_get_file_name(field(index)%unit)
- call mpp_error(FATAL,"time_interp_external 1: "//trim(err_msg)//&
- ",file="//trim(filename)//",field="//trim(field(index)%name) )
- endif
- else
- if(field(index)%modulo_time) then
- mod_time=1
- else
- mod_time=0
- endif
- call time_interp(time,field(index)%time(:),w2,t1,t2,modtime=mod_time, err_msg=err_msg)
- if(err_msg .NE. '') then
- filename = mpp_get_file_name(field(index)%unit)
- call mpp_error(FATAL,"time_interp_external 2: "//trim(err_msg)//&
- ",file="//trim(filename)//",field="//trim(field(index)%name) )
- endif
- endif
- w1 = 1.0-w2
- if (verb) then
- call get_date(time,yy,mm,dd,hh,min,ss)
- write(outunit,'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') &
- 'target time yyyy/mm/dd hh:mm:ss= ',yy,'/',mm,'/',dd,hh,':',min,':',ss
- write(outunit,*) 't1, t2, w1, w2= ', t1, t2, w1, w2
- endif
-
- call load_record(field(index),t1,horz_interp, is_in, ie_in ,js_in, je_in, window_id)
- call load_record(field(index),t2,horz_interp, is_in, ie_in ,js_in, je_in, window_id)
- i1 = find_buf_index(t1,field(index)%ibuf)
- i2 = find_buf_index(t2,field(index)%ibuf)
- if(i1<0.or.i2<0) &
- call mpp_error(FATAL,'time_interp_external : records were not loaded correctly in memory')
-
- if (verb) then
- write(outunit,*) 'ibuf= ',field(index)%ibuf
- write(outunit,*) 'i1,i2= ',i1, i2
- endif
-
- if( field(index)%region_type == NO_REGION ) then
- where(field(index)%mask(isc:iec,jsc:jec,:,i1).and.field(index)%mask(isc:iec,jsc:jec,:,i2))
- data(isw:iew,jsw:jew,:) = field(index)%data(isc:iec,jsc:jec,:,i1)*w1 + &
- field(index)%data(isc:iec,jsc:jec,:,i2)*w2
- elsewhere
-! data(isw:iew,jsw:jew,:) = time_interp_missing !field(index)%missing? Balaji
- data(isw:iew,jsw:jew,:) = field(index)%missing
- end where
- else
- where(field(index)%mask(isc:iec,jsc:jec,:,i1).and.field(index)%mask(isc:iec,jsc:jec,:,i2))
- data(isw:iew,jsw:jew,:) = field(index)%data(isc:iec,jsc:jec,:,i1)*w1 + &
- field(index)%data(isc:iec,jsc:jec,:,i2)*w2
- end where
- endif
- if(PRESENT(mask_out)) &
- mask_out(isw:iew,jsw:jew,:) = &
- field(index)%mask(isc:iec,jsc:jec,:,i1).and.&
- field(index)%mask(isc:iec,jsc:jec,:,i2)
- endif
-
- end subroutine time_interp_external_3d
-! NAME="time_interp_external"
-
- !> @brief Scalar time interpolation for @ref time_interp_external
- subroutine time_interp_external_0d(index, time, data, verbose)
-
- integer, intent(in) :: index
- type(time_type), intent(in) :: time
- real, intent(inout) :: data
- logical, intent(in), optional :: verbose
-
- integer :: t1, t2
- integer :: i1, i2, mod_time
- integer :: yy, mm, dd, hh, min, ss
- character(len=256) :: err_msg, filename
-
- real :: w1,w2
- logical :: verb
-
- verb=.false.
- if (PRESENT(verbose)) verb=verbose
- if (debug_this_module) verb = .true.
-
- if (index < 1.or.index > num_fields) &
- call mpp_error(FATAL, &
- & 'invalid index in call to time_interp_ext -- field was not initialized or failed to initialize')
-
- if (field(index)%siz(4) == 1) then
- ! only one record in the file => time-independent field
- call load_record_0d(field(index),1)
- i1 = find_buf_index(1,field(index)%ibuf)
- data = field(index)%data(1,1,1,i1)
- else
- if(field(index)%have_modulo_times) then
- call time_interp(time,field(index)%modulo_time_beg, field(index)%modulo_time_end, field(index)%time(:), &
- w2, t1, t2, field(index)%correct_leap_year_inconsistency, err_msg=err_msg)
- if(err_msg .NE. '') then
- filename = mpp_get_file_name(field(index)%unit)
- call mpp_error(FATAL,"time_interp_external 3:"//trim(err_msg)//&
- ",file="//trim(filename)//",field="//trim(field(index)%name) )
- endif
- else
- if(field(index)%modulo_time) then
- mod_time=1
- else
- mod_time=0
- endif
- call time_interp(time,field(index)%time(:),w2,t1,t2,modtime=mod_time, err_msg=err_msg)
- if(err_msg .NE. '') then
- filename = mpp_get_file_name(field(index)%unit)
- call mpp_error(FATAL,"time_interp_external 4:"//trim(err_msg)// &
- ",file="//trim(filename)//",field="//trim(field(index)%name) )
- endif
- endif
- w1 = 1.0-w2
- if (verb) then
- call get_date(time,yy,mm,dd,hh,min,ss)
- write(outunit,'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') &
- 'target time yyyy/mm/dd hh:mm:ss= ',yy,'/',mm,'/',dd,hh,':',min,':',ss
- write(outunit,*) 't1, t2, w1, w2= ', t1, t2, w1, w2
- endif
- call load_record_0d(field(index),t1)
- call load_record_0d(field(index),t2)
- i1 = find_buf_index(t1,field(index)%ibuf)
- i2 = find_buf_index(t2,field(index)%ibuf)
-
- if(i1<0.or.i2<0) &
- call mpp_error(FATAL,'time_interp_external : records were not loaded correctly in memory')
- data = field(index)%data(1,1,1,i1)*w1 + field(index)%data(1,1,1,i2)*w2
- if (verb) then
- write(outunit,*) 'ibuf= ',field(index)%ibuf
- write(outunit,*) 'i1,i2= ',i1, i2
- endif
- endif
-
- end subroutine time_interp_external_0d
-
- subroutine set_time_modulo(Time)
-
- type(time_type), intent(inout), dimension(:) :: Time
-
- integer :: ntime, n
- integer :: yr, mon, dy, hr, minu, sec
-
- ntime = size(Time(:))
-
- do n = 1, ntime
- call get_date(Time(n), yr, mon, dy, hr, minu, sec)
- yr = modulo_year
- Time(n) = set_date(yr, mon, dy, hr, minu, sec)
- enddo
-
-
- end subroutine set_time_modulo
-
-! ============================================================================
-! load specified record from file
-subroutine load_record(field, rec, interp, is_in, ie_in, js_in, je_in, window_id_in)
- type(ext_fieldtype), intent(inout) :: field
- integer , intent(in) :: rec ! record number
- type(horiz_interp_type), intent(in), optional :: interp
- integer, intent(in), optional :: is_in, ie_in, js_in, je_in
- integer, intent(in), optional :: window_id_in
-
- ! ---- local vars
- integer :: ib ! index in the array of input buffers
- integer :: isw,iew,jsw,jew ! boundaries of the domain on each window
- integer :: is_region, ie_region, js_region, je_region, i, j
- integer :: start(4), nread(4)
- logical :: need_compute
- real :: mask_in(size(field%src_data,1),size(field%src_data,2),size(field%src_data,3))
- real, allocatable :: mask_out(:,:,:)
- integer :: window_id
-
- window_id = 1
- if( PRESENT(window_id_in) ) window_id = window_id_in
- need_compute = .true.
-
-!$OMP CRITICAL
- ib = find_buf_index(rec,field%ibuf)
-
- if(ib>0) then
- !--- do nothing
- need_compute = .false.
- else
- ! calculate current buffer number in round-robin fasion
- field%nbuf = field%nbuf + 1
- if(field%nbuf > size(field%data,4).or.field%nbuf <= 0) field%nbuf = 1
- ib = field%nbuf
- field%ibuf(ib) = rec
- field%need_compute(ib,:) = .true.
-
- if (field%domain_present .and. .not.PRESENT(interp)) then
- if (debug_this_module) write(outunit,*) 'reading record with domain for field ',trim(field%name)
- call mpp_read(field%unit,field%field,field%domain,field%src_data(:,:,:,ib),rec)
- else
- if (debug_this_module) write(outunit,*) 'reading record without domain for field ',trim(field%name)
- start = 1; nread = 1
- start(1) = field%is_src; nread(1) = field%ie_src - field%is_src + 1
- start(2) = field%js_src; nread(2) = field%je_src - field%js_src + 1
- start(3) = 1; nread(3) = size(field%src_data,3)
- start(field%tdim) = rec; nread(field%tdim) = 1
- call mpp_read(field%unit,field%field,field%src_data(:,:,:,ib),start,nread)
- endif
- endif
-!$OMP END CRITICAL
- isw=field%isc;iew=field%iec
- jsw=field%jsc;jew=field%jec
-
- if( field%numwindows > 1) then
- if( .NOT. PRESENT(is_in) .OR. .NOT. PRESENT(ie_in) .OR. .NOT. PRESENT(js_in) .OR. .NOT. PRESENT(je_in) ) then
- call mpp_error(FATAL, &
- & 'time_interp_external(load_record): is_in, ie_in, js_in, je_in must be present when numwindows>1')
- endif
- isw = isw + is_in - 1
- iew = isw + ie_in - is_in
- jsw = jsw + js_in - 1
- jew = jsw + je_in - js_in
- endif
-
- ! interpolate to target grid
-
- need_compute = field%need_compute(ib, window_id)
- if(need_compute) then
- if(PRESENT(interp)) then
- is_region = field%is_region; ie_region = field%ie_region
- js_region = field%js_region; je_region = field%je_region
- mask_in = 0.0
- where (mpp_is_valid(field%src_data(:,:,:,ib), field%valid)) mask_in = 1.0
- if ( field%region_type .NE. NO_REGION ) then
- if( ANY(mask_in == 0.0) ) then
- call mpp_error(FATAL, "time_interp_external: mask_in should be all 1 when region_type is not NO_REGION")
- endif
- if( field%region_type == OUTSIDE_REGION) then
- do j = js_region, je_region
- do i = is_region, ie_region
- mask_in(i,j,:) = 0.0
- enddo
- enddo
- else ! field%region_choice == INSIDE_REGION
- do j = 1, size(mask_in,2)
- do i = 1, size(mask_in,1)
- if( jje_region .OR. iie_region ) mask_in(i,j,:) = 0.0
- enddo
- enddo
- endif
- endif
- allocate(mask_out(isw:iew,jsw:jew, size(field%src_data,3)))
- call horiz_interp(interp,field%src_data(:,:,:,ib),field%data(isw:iew,jsw:jew,:,ib), &
- mask_in=mask_in, &
- mask_out=mask_out)
-
- field%mask(isw:iew,jsw:jew,:,ib) = mask_out(isw:iew,jsw:jew,:) > 0
- deallocate(mask_out)
- else
- if ( field%region_type .NE. NO_REGION ) then
- call mpp_error(FATAL, "time_interp_external: region_type should be NO_REGION when interp is not present")
- endif
- field%data(isw:iew,jsw:jew,:,ib) = field%src_data(isw:iew,jsw:jew,:,ib)
- field%mask(isw:iew,jsw:jew,:,ib) = mpp_is_valid(field%data(isw:iew,jsw:jew,:,ib),field%valid)
- endif
- ! convert units
- where(field%mask(isw:iew,jsw:jew,:,ib)) field%data(isw:iew,jsw:jew,:,ib) = &
- field%data(isw:iew,jsw:jew,:,ib)*field%slope + field%intercept
- field%need_compute(ib, window_id) = .false.
- endif
-
-end subroutine load_record
-
-
-subroutine load_record_0d(field, rec)
- type(ext_fieldtype), intent(inout) :: field
- integer , intent(in) :: rec ! record number
- ! ---- local vars
- integer :: ib ! index in the array of input buffers
- integer :: start(4), nread(4)
-
- ib = find_buf_index(rec,field%ibuf)
-
- if(ib>0) then
- return
- else
- ! calculate current buffer number in round-robin fasion
- field%nbuf = field%nbuf + 1
- if(field%nbuf > size(field%data,4).or.field%nbuf <= 0) field%nbuf = 1
- ib = field%nbuf
- field%ibuf(ib) = rec
-
- if (debug_this_module) write(outunit,*) 'reading record without domain for field ',trim(field%name)
- start = 1; nread = 1
- start(3) = 1; nread(3) = size(field%src_data,3)
- start(field%tdim) = rec; nread(field%tdim) = 1
- call mpp_read(field%unit,field%field,field%src_data(:,:,:,ib),start,nread)
- if ( field%region_type .NE. NO_REGION ) then
- call mpp_error(FATAL, "time_interp_external: region_type should be NO_REGION when field is scalar")
- endif
- field%data(1,1,:,ib) = field%src_data(1,1,:,ib)
- field%mask(1,1,:,ib) = mpp_is_valid(field%data(1,1,:,ib),field%valid)
- ! convert units
- where(field%mask(1,1,:,ib)) field%data(1,1,:,ib) = &
- field%data(1,1,:,ib)*field%slope + field%intercept
- endif
-
-end subroutine load_record_0d
-
-! ============================================================================
-subroutine reset_src_data_region(index, is, ie, js, je)
- integer, intent(in) :: index
- integer, intent(in) :: is, ie, js, je
- integer :: nk, nbuf
-
- if( is == field(index)%is_src .AND. ie == field(index)%ie_src .AND. &
- js == field(index)%js_src .AND. ie == field(index)%je_src ) return
-
- if( .NOT. ASSOCIATED(field(index)%src_data) ) call mpp_error(FATAL, &
- "time_interp_external: field(index)%src_data is not associated")
- nk = size(field(index)%src_data,3)
- nbuf = size(field(index)%src_data,4)
- deallocate(field(index)%src_data)
- allocate(field(index)%src_data(is:ie,js:je,nk,nbuf))
- field(index)%is_src = is
- field(index)%ie_src = ie
- field(index)%js_src = js
- field(index)%je_src = je
-
-
-end subroutine reset_src_data_region
-
-! ============================================================================
-subroutine set_override_region(index, region_type, is_region, ie_region, js_region, je_region)
- integer, intent(in) :: index, region_type
- integer, intent(in) :: is_region, ie_region, js_region, je_region
-
- field(index)%region_type = region_type
- field(index)%is_region = is_region
- field(index)%ie_region = ie_region
- field(index)%js_region = js_region
- field(index)%je_region = je_region
-
- return
-
-end subroutine set_override_region
-
-! ============================================================================
-! reallocates array of fields, increasing its size
-subroutine realloc_files(n)
- integer, intent(in) :: n ! new size
-
- type(filetype), pointer :: ptr(:)
- integer :: i
-
- if (associated(opened_files)) then
- if (n <= size(opened_files)) return ! do nothing, if requested size no more than current
- endif
-
- allocate(ptr(n))
- do i = 1, size(ptr)
- ptr(i)%filename = ''
- ptr(i)%unit = -1
- enddo
-
- if (associated(opened_files))then
- ptr(1:size(opened_files)) = opened_files(:)
- deallocate(opened_files)
- endif
- opened_files => ptr
-
-end subroutine realloc_files
-
-! ============================================================================
-! reallocates array of fields,increasing its size
-subroutine realloc_fields(n)
- integer, intent(in) :: n ! new size
-
- type(ext_fieldtype), pointer :: ptr(:)
- integer :: i, ier
-
- if (associated(field)) then
- if (n <= size(field)) return ! do nothing if requested size no more then current
- endif
-
- allocate(ptr(n))
- do i=1,size(ptr)
- ptr(i)%unit=-1
- ptr(i)%name=''
- ptr(i)%units=''
- ptr(i)%siz=-1
- ptr(i)%ndim=-1
- ptr(i)%domain = NULL_DOMAIN2D
- ptr(i)%axes(:) = default_axis
- if (ASSOCIATED(ptr(i)%time)) DEALLOCATE(ptr(i)%time, stat=ier)
- if (ASSOCIATED(ptr(i)%start_time)) DEALLOCATE(ptr(i)%start_time, stat=ier)
- if (ASSOCIATED(ptr(i)%end_time)) DEALLOCATE(ptr(i)%end_time, stat=ier)
- ptr(i)%field = default_field
- if (ASSOCIATED(ptr(i)%period)) DEALLOCATE(ptr(i)%period, stat=ier)
- ptr(i)%modulo_time=.false.
- if (ASSOCIATED(ptr(i)%data)) DEALLOCATE(ptr(i)%data, stat=ier)
- if (ASSOCIATED(ptr(i)%ibuf)) DEALLOCATE(ptr(i)%ibuf, stat=ier)
- if (ASSOCIATED(ptr(i)%src_data)) DEALLOCATE(ptr(i)%src_data, stat=ier)
- ptr(i)%nbuf=-1
- ptr(i)%domain_present=.false.
- ptr(i)%slope=1.0
- ptr(i)%intercept=0.0
- ptr(i)%isc=-1;ptr(i)%iec=-1
- ptr(i)%jsc=-1;ptr(i)%jec=-1
- enddo
- if (associated(field)) then
- ptr(1:size(field)) = field(:)
- deallocate(field)
- endif
- field=>ptr
-
-end subroutine realloc_fields
-
-
- function find_buf_index(indx,buf)
- integer :: indx
- integer, dimension(:) :: buf
- integer :: find_buf_index
-
- integer :: nbuf, i
-
- nbuf = size(buf(:))
-
- find_buf_index = -1
-
- do i=1,nbuf
- if (buf(i) == indx) then
- find_buf_index = i
- exit
- endif
- enddo
-
- end function find_buf_index
-
-!
-!
-!
-! return size of field after call to init_external_field.
-! Ordering is X/Y/Z/T.
-! This call only makes sense for non-distributed reads.
-!
-!
-!
-! returned from previous call to init_external_field.
-!
-
- function get_external_field_size(index)
-
- integer :: index
- integer :: get_external_field_size(4)
-
- if (index .lt. 1 .or. index .gt. num_fields) &
- call mpp_error(FATAL,'invalid index in call to get_external_field_size')
-
-
- get_external_field_size(1) = field(index)%siz(1)
- get_external_field_size(2) = field(index)%siz(2)
- get_external_field_size(3) = field(index)%siz(3)
- get_external_field_size(4) = field(index)%siz(4)
-
- end function get_external_field_size
-! NAME="get_external_field_size"
-
-
-!
-!
-!
-! return missing value
-!
-!
-!
-! returned from previous call to init_external_field.
-!
-
- function get_external_field_missing(index)
-
- integer :: index
- real :: get_external_field_missing
-
- if (index .lt. 1 .or. index .gt. num_fields) &
- call mpp_error(FATAL,'invalid index in call to get_external_field_size')
-
-
-! call mpp_get_atts(field(index)%field,missing=missing)
- get_external_field_missing = field(index)%missing
-
- end function get_external_field_missing
-! NAME="get_external_field_missing"
-
-!
-!
-!
-! return field axes after call to init_external_field.
-! Ordering is X/Y/Z/T.
-!
-!
-!
-! returned from previous call to init_external_field.
-!
-
-
- function get_external_field_axes(index)
-
- integer :: index
- type(axistype), dimension(4) :: get_external_field_axes
-
- if (index .lt. 1 .or. index .gt. num_fields) &
- call mpp_error(FATAL,'invalid index in call to get_external_field_size')
-
-
- get_external_field_axes(1) = field(index)%axes(1)
- get_external_field_axes(2) = field(index)%axes(2)
- get_external_field_axes(3) = field(index)%axes(3)
- get_external_field_axes(4) = field(index)%axes(4)
-
- end function get_external_field_axes
-! NAME="get_external_field_axes"
-
-! ===========================================================================
-subroutine get_time_axis(index, time)
- integer , intent(in) :: index ! field id
- type(time_type), intent(out) :: time(:) ! array of time values to be filled
-
- integer :: n ! size of the data to be assigned
-
- if (index < 1.or.index > num_fields) &
- call mpp_error(FATAL,'invalid index in call to get_time_axis')
-
- n = min(size(time),size(field(index)%time))
-
- time(1:n) = field(index)%time(1:n)
-end subroutine
-
-!
-!
-!
-! exit time_interp_external_mod. Close all open files and
-! release storage
-!
-
- subroutine time_interp_external_exit()
-
- integer :: i,j
-!
-! release storage arrays
-!
- do i=1,num_fields
- deallocate(field(i)%time,field(i)%start_time,field(i)%end_time,&
- field(i)%period,field(i)%data,field(i)%mask,field(i)%ibuf)
- if (ASSOCIATED(field(i)%src_data)) deallocate(field(i)%src_data)
- do j=1,4
- field(i)%axes(j) = default_axis
- enddo
- field(i)%domain = NULL_DOMAIN2D
- field(i)%field = default_field
- field(i)%nbuf = 0
- field(i)%slope = 0.
- field(i)%intercept = 0.
- enddo
-
- deallocate(field)
- deallocate(opened_files)
-
- num_fields = 0
-
- module_initialized = .false.
-
- end subroutine time_interp_external_exit
-! NAME="time_interp_external_exit"
-#endif
-end module time_interp_external_mod
-!> @}
-! close documentation grouping
From 06b94a7f574e7794684b8584391744ded68e2989 Mon Sep 17 00:00:00 2001
From: Ryan Mulhall <35538242+rem1776@users.noreply.github.com>
Date: Tue, 10 Oct 2023 11:47:07 -0400
Subject: [PATCH 06/10] fix: add allocatable and pure macros back into
platform.h (#1386)
---
include/fms_platform.h | 22 ++++++++++++++++++++++
1 file changed, 22 insertions(+)
diff --git a/include/fms_platform.h b/include/fms_platform.h
index fdd0bd4e41..30feb9f73b 100644
--- a/include/fms_platform.h
+++ b/include/fms_platform.h
@@ -46,6 +46,28 @@ use,intrinsic :: iso_c_binding, only: c_double,c_float,c_int64_t, &
!DEC$ MESSAGE:'Using 8-byte addressing'
#endif
+!Control "pure" functions.
+#ifdef NO_F95
+#define _PURE
+!DEC$ MESSAGE:'Not using pure routines.'
+#else
+#define _PURE pure
+!DEC$ MESSAGE:'Using pure routines.'
+#endif
+
+!Control array members of derived types.
+#ifdef NO_F2000
+#define _ALLOCATABLE pointer
+#define _NULL =>null()
+#define _ALLOCATED associated
+!DEC$ MESSAGE:'Using pointer derived type array members.'
+#else
+#define _ALLOCATABLE allocatable
+#define _NULL
+#define _ALLOCATED allocated
+!DEC$ MESSAGE:'Using allocatable derived type array members.'
+#endif
+
!Control use of cray pointers within mpp_peset
!Other cray pointer usage in mpp routines is compiled regardless
#ifdef NO_CRAY_POINTERS
From 527a42f3ff36d75aac68b65757b35f73a74146bb Mon Sep 17 00:00:00 2001
From: Ryan Mulhall <35538242+rem1776@users.noreply.github.com>
Date: Fri, 27 Oct 2023 10:44:25 -0400
Subject: [PATCH 07/10] chore: updates the various versions and changelog
(#1395)
---
CHANGELOG.md | 36 ++++++++++++++++++++++++++++++++++++
CMakeLists.txt | 2 +-
configure.ac | 2 +-
libFMS/Makefile.am | 2 +-
4 files changed, 39 insertions(+), 3 deletions(-)
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2c616e647d..782f940443 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,42 @@ and this project uses `yyyy.rr[.pp]`, where `yyyy` is the year a patch is releas
`rr` is a sequential release number (starting from `01`), and an optional two-digit
sequential patch number (starting from `01`).
+## [2023.03] - 2023-10-27
+### Known Issues
+- GCC 9 and below as well as GCC 11.1.0 are unsupported due to compilation issues. See prior releases for more details.
+- `NO_QUAD_PRECISION` macro is no longer set by FMS, the `ENABLE_QUAD_PRECISION` macro has replaced prior usage of `NO_QUAD_PRECISION`. `-DENABLE_QUAD_PRECISION` should be set if quad precision is to be used, otherwise FMS will not use quad precision reals where applicable.
+
+### Added
+- UNIT_TESTS: New unit tests have been created or and existing ones expanded on for any modules utilizing mixed precision support.
+
+### Changed
+- MIXED PRECISION: Most subroutines and functions in FMS have been updated to simultaneously accept both 4 byte and 8 byte reals as arguments. This deprecates the `--enable-mixed-mode` option, which enabled similar functionality but was limited to certain directories and was not enabled by default. To facilitate easier testing of these code changes, the CMake precision options for default real size were left in (along with an equivalent `--disable-r8-default` flag for autotools). The resulting libraries will support mixed-precision real kinds regardless of default real size. It should also be noted that many routines that accept real arguments have been moved to include files along with headers in order to be compiled with both kinds. Most module level variables were explicitly declared as r8_kind for these updates.
+- Some type/module changes were made to facilitate mixed precision support. They are **intended** to have minimal impact to other codebases:
+ - COUPLER_TYPES: In coupler_types.F90, `coupler_nd_field_type` and `coupler_nd_values_type` have been renamed to indicate real kind value: `coupler_nd_real4/8_field_type` and `coupler_nd_real4/8_values_type`. The `bc` field within `coupler_nd_bc_type` was modified to use r8_kind within the value and field types, and an additional field added `bc_r4` to use r4_kind values.
+ - TRIDIAGONAL: Module state between r4 and r8 calls are distinct (ie. subsequent calls will only be affected by calls of the same precision). This behaviour can be changed via the `save_both_kinds` optional argument to `tri_invert`.
+- CODE_STYLE: has been updated to reflect the formatting used for the mixed precision support updates.
+
+### Fixed
+- DIAG_MANAGER: Tile number (ie. tileX) will now be added to filenames for sub-regional diagnostics.
+- MPP: Bug affecting non-intel compilers coming from uninitialized pointer in the `nest_domain_type`
+- MPP: Bug fix for unallocated field causing seg faults in `mpp_check_field`
+- FMS2_IO: Fixed segfault occuring from use of cray pointer remapping along with mpp_scatter/gather
+- TEST_FMS: Added various fixes for different compilers within test programs for fms2_io, mpp, diag_manager, parser, and sat_vapor_pres.
+- INTERPOLATOR: Deallocates fields in the type that were previously left out in `interpolator_end`
+
+### Removed
+- CPP MACROS:
+ - `no_4byte_reals` was removed and will not set any additional macros if used. `no_8byte_integers` is still functional.
+ - `NO_QUAD_PRECISION` was removed. It was conditionally set if ENABLE_QUAD_PRECISION was undefined. ENABLE_QUAD_PRECISION should be used in model components instead (logic is flipped)
+ - `use_netCDF` was set by autotools previously but wasn't consistently used in the code. FMS should always be compiled with netcdf installed so this was removed with the exception of its use in deprecated IO modules.
+- DRIFTERS: The drifters subdirectory has been deprecated. It will only be compiled if using the `-Duse_drifters` CPP flag.
+
+### Tag Commit Hashes
+- 2023.03-beta1 06b94a7f574e7794684b8584391744ded68e2989
+- 2023.03-alpha3 b25a7c52a27dfd52edc10bc0ebe12776af0f03df
+- 2023.03-alpha2 9983ce308e62e9f7215b04c227cebd30fd75e784
+- 2023.03-alpha1 a46bd94fd8dd1f6f021501e29179003ff28180ec
+
## [2023.02] - 2023-07-27
### Known Issues
- GCC 11.1.0 is unsupported due to compilation issues with select type. The issue is resolved in later GCC releases.
diff --git a/CMakeLists.txt b/CMakeLists.txt
index a70abe14da..89c7eb329d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -26,7 +26,7 @@ set(CMAKE_Fortran_FLAGS_DEBUG)
# Define the CMake project
project(FMS
- VERSION 2023.02.0
+ VERSION 2023.03.0
DESCRIPTION "GFDL FMS Library"
HOMEPAGE_URL "https://www.gfdl.noaa.gov/fms"
LANGUAGES C Fortran)
diff --git a/configure.ac b/configure.ac
index 65e3c0fe73..5a9aff17c7 100644
--- a/configure.ac
+++ b/configure.ac
@@ -25,7 +25,7 @@ AC_PREREQ([2.69])
# Initialize with name, version, and support email address.
AC_INIT([GFDL FMS Library],
- [2023.02.00-dev],
+ [2023.03.00],
[gfdl.climate.model.info@noaa.gov],
[FMS],
[https://www.github.com/NOAA-GFDL/FMS])
diff --git a/libFMS/Makefile.am b/libFMS/Makefile.am
index db57f86562..9605216504 100644
--- a/libFMS/Makefile.am
+++ b/libFMS/Makefile.am
@@ -28,7 +28,7 @@ lib_LTLIBRARIES = libFMS.la
# These linker flags specify libtool version info.
# See http://www.gnu.org/software/libtool/manual/libtool.html#Libtool-versioning
# for information regarding incrementing `-version-info`.
-libFMS_la_LDFLAGS = -version-info 16:0:0
+libFMS_la_LDFLAGS = -version-info 17:0:0
# Add the convenience libraries to the FMS library.
libFMS_la_LIBADD = $(top_builddir)/platform/libplatform.la
From 74338de3ee45eccdc8f2b045d096ed08a7f0012e Mon Sep 17 00:00:00 2001
From: "github-actions[bot]"
<41898282+github-actions[bot]@users.noreply.github.com>
Date: Fri, 27 Oct 2023 16:08:19 -0400
Subject: [PATCH 08/10] chore: add dev to version number (#1400)
---
configure.ac | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/configure.ac b/configure.ac
index 5a9aff17c7..3b242c0016 100644
--- a/configure.ac
+++ b/configure.ac
@@ -25,7 +25,7 @@ AC_PREREQ([2.69])
# Initialize with name, version, and support email address.
AC_INIT([GFDL FMS Library],
- [2023.03.00],
+ [2023.03.00-dev],
[gfdl.climate.model.info@noaa.gov],
[FMS],
[https://www.github.com/NOAA-GFDL/FMS])
From c70718d10bc91a27f588a63b14b7bceb292a9ed0 Mon Sep 17 00:00:00 2001
From: Ryan Mulhall <35538242+rem1776@users.noreply.github.com>
Date: Wed, 1 Nov 2023 15:05:19 -0400
Subject: [PATCH 09/10] fix: allow for arguments within mpi launcher override
(#1401)
---
test_fms/test-lib.sh.in | 4 +++-
1 file changed, 3 insertions(+), 1 deletion(-)
diff --git a/test_fms/test-lib.sh.in b/test_fms/test-lib.sh.in
index b983b48d84..9be57a630a 100644
--- a/test_fms/test-lib.sh.in
+++ b/test_fms/test-lib.sh.in
@@ -96,8 +96,10 @@ mpirun () {
# Set the name of the mpi launcher for use in test scripts.
local mpi_launcher='@MPI_LAUNCHER@'
local oversubscribe='@OVERSUBSCRIBE@'
+ # need to strip off any args that may be included with MPI_LAUNCHER arg for check below to work
+ local mpi_cmd="`echo $mpi_launcher | awk '{print $1;}'`"
# Check if running with MPI: if so, the mpi_launcher will point to a command
- command -v "$mpi_launcher" 2>&1 > /dev/null
+ command -v "$mpi_cmd" 2>&1 > /dev/null
if test $? -eq 0
then
# use `command` to keep from reusing this function
From b0dbc305b585a244f19880640e0d9ed508455ed8 Mon Sep 17 00:00:00 2001
From: Bennett Chang <98476720+bcc2761@users.noreply.github.com>
Date: Thu, 9 Nov 2023 10:09:19 -0500
Subject: [PATCH 10/10] feat: data_override namelist flag to enable YAML data
tables (#1382)
---
data_override/include/data_override.inc | 39 ++++++++++-----
test_fms/data_override/Makefile.am | 2 +
test_fms/data_override/test_data_override2.sh | 49 +++++++++++++++++++
.../data_override/test_data_override_init.F90 | 29 +++++++++++
4 files changed, 107 insertions(+), 12 deletions(-)
create mode 100644 test_fms/data_override/test_data_override_init.F90
diff --git a/data_override/include/data_override.inc b/data_override/include/data_override.inc
index 663d2b0fcf..6d76c0d537 100644
--- a/data_override/include/data_override.inc
+++ b/data_override/include/data_override.inc
@@ -103,11 +103,7 @@ real(FMS_DATA_OVERRIDE_KIND_) :: min_glo_lon_lnd, max_glo_lon_lnd
real(FMS_DATA_OVERRIDE_KIND_) :: min_glo_lon_ice, max_glo_lon_ice
integer :: num_fields = 0 !< number of fields in override_array already processed
-#ifdef use_yaml
type(data_type), dimension(:), allocatable :: data_table !< user-provided data table
-#else
-type(data_type), dimension(max_table) :: data_table !< user-provided data table
-#endif
type(data_type) :: default_table
type(override_type), dimension(max_array) :: override_array !< to store processed fields
@@ -118,8 +114,9 @@ logical :: reproduce_null_char_bug = .false.
!! to reproduce the mpp_io bug where lat/lon_bnd were
!! not read correctly if null characters are present in
!! the netcdf file
+logical :: use_data_table_yaml = .false.
-namelist /data_override_nml/ debug_data_override, grid_center_bug, reproduce_null_char_bug
+namelist /data_override_nml/ debug_data_override, grid_center_bug, reproduce_null_char_bug, use_data_table_yaml
public :: DATA_OVERRIDE_INIT_IMPL_, DATA_OVERRIDE_UNSET_ATM_, DATA_OVERRIDE_UNSET_OCN_, &
& DATA_OVERRIDE_UNSET_LND_, DATA_OVERRIDE_UNSET_ICE_, DATA_OVERRIDE_0D_, &
@@ -166,6 +163,12 @@ if (grid_center_bug) then
"that is no longer supported. Please remove this namelist variable.")
endif
+if (use_data_table_yaml) then
+ call mpp_error(NOTE, "You are using YAML.")
+else
+ call mpp_error(NOTE, "You are using the legacy table.")
+end if
+
atm_on = PRESENT(Atm_domain_in)
ocn_on = PRESENT(Ocean_domain_in)
lnd_on = PRESENT(Land_domain_in)
@@ -197,12 +200,25 @@ endif
default_table%interpol_method = 'bilinear'
#ifdef use_yaml
- call read_table_yaml(data_table)
+ if (use_data_table_yaml) then
+ call read_table_yaml(data_table)
+ else
+ allocate(data_table(max_table))
+ do i = 1, max_table
+ data_table(i) = default_table
+ enddo
+ call read_table(data_table)
+ end if
#else
- do i = 1,max_table
- data_table(i) = default_table
- enddo
- call read_table(data_table)
+ if (use_data_table_yaml) then
+ call mpp_error(FATAL, "compilation error, need to compile with `-Duse_yaml`")
+ else
+ allocate(data_table(max_table))
+ do i = 1, max_table
+ data_table(i) = default_table
+ enddo
+ call read_table(data_table)
+ end if
#endif
! Initialize override array
@@ -330,7 +346,6 @@ function count_ne_1(in_1, in_2, in_3)
count_ne_1 = .not.(in_1.NEQV.in_2.NEQV.in_3) .OR. (in_1.AND.in_2.AND.in_3)
end function count_ne_1
-#ifndef use_yaml
subroutine read_table(data_table)
type(data_type), dimension(max_table), intent(inout) :: data_table
@@ -475,7 +490,7 @@ subroutine read_table(data_table)
if(io_status/=0) call mpp_error(FATAL, 'data_override_mod: Error in closing file data_table')
end subroutine read_table
-#else
+#ifdef use_yaml
subroutine read_table_yaml(data_table)
type(data_type), dimension(:), allocatable, intent(out) :: data_table
diff --git a/test_fms/data_override/Makefile.am b/test_fms/data_override/Makefile.am
index f5c956c446..366c773800 100644
--- a/test_fms/data_override/Makefile.am
+++ b/test_fms/data_override/Makefile.am
@@ -30,6 +30,7 @@ LDADD = $(top_builddir)/libFMS/libFMS.la
# Build this test program.
check_PROGRAMS = \
+ test_data_override_init \
test_get_grid_v1_r4 \
test_get_grid_v1_r8 \
test_data_override_r4 \
@@ -38,6 +39,7 @@ check_PROGRAMS = \
test_data_override_ongrid_r8
# This is the source code for the test.
+test_data_override_init_SOURCES = test_data_override_init.F90
test_data_override_r4_SOURCES = test_data_override.F90
test_data_override_r8_SOURCES = test_data_override.F90
diff --git a/test_fms/data_override/test_data_override2.sh b/test_fms/data_override/test_data_override2.sh
index 35546b41d3..064b0511d4 100755
--- a/test_fms/data_override/test_data_override2.sh
+++ b/test_fms/data_override/test_data_override2.sh
@@ -117,4 +117,53 @@ fi
done
+# data_override with the default table (not setting namelist)
+cat <<_EOF > data_table
+"ICE", "sst_obs", "SST", "INPUT/sst_ice_clim.nc", .false., 300.0
+_EOF
+
+test_expect_success "data_override_init with the default table" '
+ mpirun -n 1 ./test_data_override_init
+'
+# data_override with yaml table (setting namelist to .True.)
+cat <<_EOF > input.nml
+&data_override_nml
+use_data_table_yaml=.true.
+/
+_EOF
+
+cat <<_EOF > data_table.yaml
+data_table:
+ - gridname : OCN
+ fieldname_code : runoff
+ fieldname_file : runoff
+ file_name : INPUT/runoff.daitren.clim.1440x1080.v20180328.nc
+ interpol_method : none
+ factor : 1.0
+_EOF
+
+if [ ! -z $parser_skip ]; then
+ test_expect_failure "data_override_init with the yaml table" '
+ mpirun -n 1 ./test_data_override_init
+ '
+else
+ test_expect_success "data_override_init with the yaml table" '
+ mpirun -n 1 ./test_data_override_init
+ '
+fi
+#data_override with default table (setting namelist to .True.)
+cat <<_EOF > data_table
+"ICE", "sst_obs", "SST", "INPUT/sst_ice_clim.nc", .true., 300.0
+_EOF
+
+cat <<_EOF > input.nml
+&data_override_nml
+use_data_table_yaml=.false.
+/
+_EOF
+
+test_expect_success "data_override_init with the default table" '
+ mpirun -n 1 ./test_data_override_init
+'
+
test_done
diff --git a/test_fms/data_override/test_data_override_init.F90 b/test_fms/data_override/test_data_override_init.F90
new file mode 100644
index 0000000000..dceec5aca3
--- /dev/null
+++ b/test_fms/data_override/test_data_override_init.F90
@@ -0,0 +1,29 @@
+!***********************************************************************
+!* 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 .
+!***********************************************************************
+
+program test_data_override_init
+
+ use fms_mod, only: fms_init, fms_end
+ use data_override_mod
+
+ call fms_init()
+ call data_override_init
+ call fms_end()
+
+end program test_data_override_init