35 real(FMS_TRID_KIND_),
intent(out),
dimension(:,:,:) :: x
36 real(FMS_TRID_KIND_),
intent(in),
dimension(:,:,:) :: d
37 real(FMS_TRID_KIND_),
optional,
dimension(:,:,:) :: a,b,c
39 logical,
optional :: store_both_kinds
43 real(FMS_TRID_KIND_),
dimension(size(x,1),size(x,2),size(x,3)) :: f
44 integer,
parameter :: kindl = fms_trid_kind_
51 if(
allocated(trid_real_type%e))
deallocate(trid_real_type%e)
52 if(
allocated(trid_real_type%g))
deallocate(trid_real_type%g)
53 if(
allocated(trid_real_type%bb))
deallocate(trid_real_type%bb)
54 if(
allocated(trid_real_type%cc))
deallocate(trid_real_type%cc)
55 allocate(trid_real_type%e (
size(x,1),
size(x,2),
size(x,3)))
56 allocate(trid_real_type%g (
size(x,1),
size(x,2),
size(x,3)))
57 allocate(trid_real_type%bb(
size(x,1),
size(x,2)))
58 allocate(trid_real_type%cc(
size(x,1),
size(x,2),
size(x,3)))
61 trid_real_type%e(:,:,1) = - a(:,:,1) / b(:,:,1)
62 a(:,:,
size(x,3)) = 0.0_kindl
65 trid_real_type%g(:,:,k) = 1.0_kindl/(b(:,:,k)+c(:,:,k)*trid_real_type%e(:,:,k-1))
66 trid_real_type%e(:,:,k) = - a(:,:,k)* trid_real_type%g(:,:,k)
69 trid_real_type%bb = 1.0_kindl/b(:,:,1)
73 if(.not.init_var)
call mpp_error(fatal,
'tri_invert: a,b,and c args not provided or previously calculated.')
75 f(:,:,1) = d(:,:,1)*trid_real_type%bb
77 f(:,:,k) = (d(:,:,k) - trid_real_type%cc(:,:,k)*f(:,:,k-1))*trid_real_type%g(:,:,k)
80 x(:,:,
size(x,3)) = f(:,:,
size(x,3))
81 do k =
size(x,3)-1,1,-1
82 x(:,:,k) = trid_real_type%e(:,:,k)*x(:,:,k+1)+f(:,:,k)
86 if(
present(store_both_kinds))
then
87 if( store_both_kinds )
then
88 if( fms_trid_kind_ .eq. r8_kind)
then
89 tridiag_r4%e = real(trid_real_type%e, r4_kind)
90 tridiag_r4%g = real(trid_real_type%g, r4_kind)
91 tridiag_r4%cc = real(trid_real_type%cc, r4_kind)
92 tridiag_r4%bb = real(trid_real_type%bb, r4_kind)
93 init_tridiagonal_r4 = .true.
95 tridiag_r8%e = real(trid_real_type%e, r8_kind)
96 tridiag_r8%g = real(trid_real_type%g, r8_kind)
97 tridiag_r8%cc = real(trid_real_type%cc, r8_kind)
98 tridiag_r8%bb = real(trid_real_type%bb, r8_kind)
99 init_tridiagonal_r8 = .true.
105 end subroutine TRI_INVERT_
subroutine tri_invert_(x, d, a, b, c, store_both_kinds)
Sets up and solves the tridiagonal system of equations.