tridiagonal/tridiagonal.F90

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