FMS  2025.04
Flexible Modeling System
tridiagonal.F90
1 !***********************************************************************
2 !* Apache License 2.0
3 !*
4 !* This file is part of the GFDL Flexible Modeling System (FMS).
5 !*
6 !* Licensed under the Apache License, Version 2.0 (the "License");
7 !* you may not use this file except in compliance with the License.
8 !* You may obtain a copy of the License at
9 !*
10 !* http://www.apache.org/licenses/LICENSE-2.0
11 !*
12 !* FMS is distributed in the hope that it will be useful, but WITHOUT
13 !* WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied;
14 !* without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
15 !* PARTICULAR PURPOSE. See the License for the specific language
16 !* governing permissions and limitations under the License.
17 !***********************************************************************
18 !> @defgroup tridiagonal_mod tridiagonal_mod
19 !> @ingroup tridiagonal
20 !> @brief Solves a tridiagonal system of equations.
21 !!
22 !> The following schematic represents the system of equations solved,
23 !! where X is the solution.
24 !! <PRE>
25 !! | B(1) A(1) 0 0 ....... 0 | |X(1)| |D(1)|
26 !! | C(2) B(2) A(2) 0 ....... 0 | |X(2)| |D(2)|
27 !! | 0 C(3) B(3) A(3) 0 ....... 0 | | .. | | .. |
28 !! | .......................................... | | .. | = | .. |
29 !! | .......................................... | | .. | | .. |
30 !! | C(N-2) B(N-2) A(N-2) 0 | | .. | | .. |
31 !! | 0 C(N-1) B(N-1) A(N-1)| | .. | | .. |
32 !! | 0 0 C(N) B(N) | |X(N)| |D(N)|
33 !!
34 !! </PRE>
35 !! To solve this system
36 !! <PRE>
37 !! call tri_invert(X,D,A,B,C)
38 !!
39 !! real, intent(out), dimension(:,:,:) :: X
40 !! real, intent(in), dimension(:,:,:) :: D
41 !! real, optional, dimension(:,:,:) :: A,B,C
42 !! </PRE>
43 !! For simplicity (?), A and C are assumed to be dimensioned the same size
44 !! as B, D, and X, although any input values for A(N) and C(1) are ignored.
45 !! (some checks are needed here)
46 !!
47 !! If A is not present, it is assumed that the matrix (A,B.C) has not been changed
48 !! since the last call to tri_invert.
49 !!
50 !! To release memory,
51 !! <PRE>
52 !! call close_tridiagonal
53 !! </PRE>
54 !!
55 !!
56 !! Arguments A, B, and C are optional, and are saved as module variables
57 !! if one recalls tri_invert without changing (A,B,C)
58 !!
59 !! @note
60 !! Optional arguments A,B,C have no intent declaration,
61 !! so the default intent is inout. The value of A(N) is modified
62 !! on output, and B and C are unchanged.
63 !!
64 !! The following private allocatable arrays save the relevant information
65 !! if one recalls tri_invert without changing (A,B,C):
66 !! <PRE>
67 !! allocate ( e (size(x,1), size(x,2), size(x,3)) )
68 !! allocate ( g (size(x,1), size(x,2), size(x,3)) )
69 !! allocate ( cc (size(x,1), size(x,2), size(x,3)) )
70 !! allocate ( bb (size(x,1), size(x,2)) )
71 !! </PRE>
72 !! This storage is deallocated when close_tridiagonal is called.
73 
74 !> @addtogroup tridiagonal_mod
75 !> @{
76 module tridiagonal_mod
77 
78  use platform_mod, only: r4_kind, r8_kind
79  use mpp_mod, only: mpp_error, fatal
80  implicit none
81 
83  real(r4_kind), private, allocatable, dimension(:,:,:) :: e, g, cc
84  real(r4_kind), private, allocatable, dimension(:,:) :: bb
85  end type
86 
88  real(r8_kind), private, allocatable, dimension(:,:,:) :: e, g, cc
89  real(r8_kind), private, allocatable, dimension(:,:) :: bb
90  end type
91 
92  type(tridiag_reals_r4) :: tridiag_r4 !< holds reals stored from r4_kind calls to tri_invert
93  type(tridiag_reals_r8) :: tridiag_r8 !< holds reals stored from r8_kind calls to tri_invert
94 
95  logical, private :: init_tridiagonal_r4 = .false. !< true when fields in tridiag_r4 are allocated
96  logical, private :: init_tridiagonal_r8 = .false. !< true when fields in tridiag_r8 are allocated
97 
98  !> Interface to solve tridiagonal systems of equations for either kind value.
99  !! Module level variables will be deallocated and allocated for every
100  !! Since this relies on the state of module variables (unless A,B,C are specified)
101  !! the values stored are distinct for each kind call unless the added optional argument store_both_kinds
102  !! is true.
103  interface tri_invert
104  module procedure tri_invert_r4
105  module procedure tri_invert_r8
106  end interface
107 
108  public :: tri_invert
109 
110  contains
111 
112  !> @brief Releases memory used by the solver
113  subroutine close_tridiagonal
114  if(.not. init_tridiagonal_r4 .and. .not. init_tridiagonal_r8) return
115  !$OMP SINGLE
116  if(allocated(tridiag_r4%e)) deallocate(tridiag_r4%e)
117  if(allocated(tridiag_r4%g)) deallocate(tridiag_r4%g)
118  if(allocated(tridiag_r4%cc)) deallocate(tridiag_r4%cc)
119  if(allocated(tridiag_r4%bb)) deallocate(tridiag_r4%bb)
120  if(allocated(tridiag_r8%e)) deallocate(tridiag_r8%e)
121  if(allocated(tridiag_r8%g)) deallocate(tridiag_r8%g)
122  if(allocated(tridiag_r8%cc)) deallocate(tridiag_r8%cc)
123  if(allocated(tridiag_r8%bb)) deallocate(tridiag_r8%bb)
124  init_tridiagonal_r4 = .false.; init_tridiagonal_r8 = .false.
125  !$OMP END SINGLE
126  return
127  end subroutine close_tridiagonal
128 
129 #include "tridiagonal_r4.fh"
130 #include "tridiagonal_r8.fh"
131 
132 end module tridiagonal_mod
133 
134 !> @}
135 ! close documentation grouping
Error handler.
Definition: mpp.F90:381
logical, private init_tridiagonal_r4
true when fields in tridiag_r4 are allocated
Definition: tridiagonal.F90:95
type(tridiag_reals_r4) tridiag_r4
holds reals stored from r4_kind calls to tri_invert
Definition: tridiagonal.F90:92
subroutine close_tridiagonal
Releases memory used by the solver.
logical, private init_tridiagonal_r8
true when fields in tridiag_r8 are allocated
Definition: tridiagonal.F90:96
type(tridiag_reals_r8) tridiag_r8
holds reals stored from r8_kind calls to tri_invert
Definition: tridiagonal.F90:93
Interface to solve tridiagonal systems of equations for either kind value. Module level variables wil...