/[MITgcm]/MITgcm/model/src/solve_pentadiagonal.F
ViewVC logotype

Contents of /MITgcm/model/src/solve_pentadiagonal.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.2 - (show annotations) (download)
Thu Jan 8 16:18:11 2004 UTC (20 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint58l_post, checkpoint57t_post, checkpoint57o_post, checkpoint52l_pre, hrcube4, checkpoint58e_post, mitgcm_mapl_00, checkpoint52n_post, checkpoint52j_post, checkpoint53d_post, checkpoint58u_post, checkpoint58w_post, checkpoint54a_pre, checkpoint57m_post, checkpoint55c_post, checkpoint54e_post, checkpoint57s_post, checkpoint54a_post, checkpoint53c_post, checkpoint57k_post, checkpoint55d_pre, checkpoint57d_post, checkpoint57g_post, checkpoint60, checkpoint61, checkpoint62, checkpoint57b_post, checkpoint57c_pre, checkpoint58r_post, checkpoint55j_post, checkpoint56b_post, checkpoint57i_post, checkpoint57y_post, hrcube_1, checkpoint57e_post, checkpoint52l_post, checkpoint55h_post, checkpoint58n_post, checkpoint53b_post, checkpoint58x_post, checkpoint52k_post, checkpoint57g_pre, checkpoint54b_post, checkpoint53b_pre, checkpoint55b_post, checkpoint58t_post, checkpoint58h_post, checkpoint54d_post, checkpoint56c_post, checkpoint52m_post, checkpoint57y_pre, checkpoint55, checkpoint53a_post, checkpoint57f_pre, checkpoint57a_post, checkpoint54, checkpoint58q_post, checkpoint54f_post, checkpoint57v_post, checkpoint59q, checkpoint59p, checkpoint55g_post, checkpoint59r, checkpoint58j_post, checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59a, checkpoint55f_post, checkpoint59c, checkpoint59b, checkpoint59m, checkpoint59l, checkpoint59o, checkpoint59n, checkpoint59i, checkpoint59h, checkpoint59k, checkpoint59j, checkpoint57r_post, checkpoint59, checkpoint58, checkpoint57a_pre, checkpoint55i_post, checkpoint57, checkpoint56, checkpoint53, eckpoint57e_pre, checkpoint57h_done, checkpoint58f_post, checkpoint53g_post, checkpoint52f_post, checkpoint57x_post, checkpoint57n_post, checkpoint58d_post, checkpoint58c_post, checkpoint57w_post, checkpoint57p_post, checkpint57u_post, checkpoint57f_post, checkpoint58a_post, checkpoint58i_post, checkpoint57q_post, checkpoint58g_post, hrcube5, checkpoint58o_post, checkpoint57z_post, checkpoint62c, checkpoint62b, checkpoint62a, checkpoint57c_post, checkpoint58y_post, checkpoint55e_post, checkpoint58k_post, checkpoint52i_post, checkpoint52j_pre, checkpoint58v_post, checkpoint53f_post, checkpoint55a_post, checkpoint53d_pre, checkpoint54c_post, checkpoint58s_post, checkpoint61f, checkpoint61g, checkpoint61d, checkpoint61e, checkpoint61b, checkpoint61c, checkpoint58p_post, checkpoint61a, checkpoint61n, checkpoint61o, checkpoint61l, checkpoint61m, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i, checkpoint61v, checkpoint61w, checkpoint61t, checkpoint61u, checkpoint61r, checkpoint61s, checkpoint61p, checkpoint61q, checkpoint57j_post, checkpoint61z, checkpoint61x, checkpoint61y, checkpoint58b_post, checkpoint57h_pre, checkpoint58m_post, checkpoint57l_post, checkpoint52i_pre, checkpoint52h_pre, checkpoint52f_pre, checkpoint57h_post, hrcube_2, hrcube_3, checkpoint56a_post, checkpoint55d_post
Changes since 1.1: +79 -70 lines 
make it safer in case Nr=1 or Nr=2
1 C $Header: /u/gcmpack/MITgcm/model/src/solve_pentadiagonal.F,v 1.1 2004/01/07 21:19:37 jmc Exp $
2 C $Name: $
3
4 #include "CPP_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: SOLVE_PENTADIAGONAL
8 C !INTERFACE:
9 SUBROUTINE SOLVE_PENTADIAGONAL(
10 I iMin,iMax, jMin,jMax,
11 U a5d, b5d, c5d, d5d, e5d,
12 U y5d,
13 O errCode,
14 I bi, bj, myThid )
15 C !DESCRIPTION: \bv
16 C *==========================================================*
17 C | S/R SOLVE_PENTADIAGONAL
18 C | o Solve a penta-diagonal system A*X=Y (dimension Nr)
19 C *==========================================================*
20 C | o Used to solve implicitly vertical advection & diffusion
21 C *==========================================================*
22 C \ev
23
24 C !USES:
25 IMPLICIT NONE
26 C == Global data ==
27 #include "SIZE.h"
28 #include "EEPARAMS.h"
29
30 C !INPUT/OUTPUT PARAMETERS:
31 C == Routine Arguments ==
32 C INPUT:
33 C iMin,iMax,jMin,jMax :: computational domain
34 C a5d :: 2nd lower diagonal of the pentadiagonal matrix
35 C b5d :: 1rst lower diagonal of the pentadiagonal matrix
36 C c5d :: main diagonal of the pentadiagonal matrix
37 C d5d :: 1rst upper diagonal of the pentadiagonal matrix
38 C e5d :: 2nd upper diagonal of the pentadiagonal matrix
39 C y5d :: Y vector (R.H.S.);
40 C bi,bj :: tile indices
41 C myThid :: thread number
42 C OUTPUT:
43 C y5d :: X = solution of A*X=Y
44 C a5d,b5d,c5d,d5d,e5d :: modified to enable to find X' solution of
45 C A*X'=Y' without solving the full system again
46 C errCode :: > 0 if singular matrix
47 INTEGER iMin,iMax,jMin,jMax
48 _RL a5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
49 _RL b5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
50 _RL c5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
51 _RL d5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
52 _RL e5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
53 _RL y5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
54 INTEGER errCode
55 INTEGER bi, bj, myThid
56
57 C !LOCAL VARIABLES:
58 C == Local variables ==
59 INTEGER i,j,k
60 CEOP
61
62 errCode = 0
63
64 DO k=1,Nr
65 C-- forward sweep (starting from top)
66 IF (k.EQ.1) THEN
67 DO j=jMin,jMax
68 DO i=iMin,iMax
69 IF ( c5d(i,j,k).NE.0. _d 0 ) THEN
70 c5d(i,j,k) = 1. _d 0 / c5d(i,j,k)
71 ELSE
72 c5d(i,j,k) = 0. _d 0
73 errCode = 1
74 ENDIF
75 ENDDO
76 ENDDO
77
78 ELSEIF (k.EQ.2) THEN
79 DO j=jMin,jMax
80 DO i=iMin,iMax
81 C-- [k] <- [k] - b_k/c_k-1 * [k-1]
82 b5d(i,j,k) = b5d(i,j,k)*c5d(i,j,k-1)
83 c5d(i,j,k) = c5d(i,j,k) - b5d(i,j,k)*d5d(i,j,k-1)
84 d5d(i,j,k) = d5d(i,j,k) - b5d(i,j,k)*e5d(i,j,k-1)
85 y5d(i,j,k,bi,bj) = y5d(i,j,k,bi,bj)
86 & - b5d(i,j,k)*y5d(i,j,k-1,bi,bj)
87 IF ( c5d(i,j,k).NE.0. _d 0 ) THEN
88 c5d(i,j,k) = 1. _d 0 / c5d(i,j,k)
89 ELSE
90 c5d(i,j,k) = 0. _d 0
91 errCode = 1
92 ENDIF
93 ENDDO
94 ENDDO
95
96 ELSE
97 C-- Middle of forward sweep
98 DO j=jMin,jMax
99 DO i=iMin,iMax
100 C-- [k] <- [k] - a_k/c_k-2 * [k-2]
101 a5d(i,j,k) = a5d(i,j,k)*c5d(i,j,k-2)
102 b5d(i,j,k) = b5d(i,j,k) - a5d(i,j,k)*d5d(i,j,k-2)
103 c5d(i,j,k) = c5d(i,j,k) - a5d(i,j,k)*e5d(i,j,k-2)
104 C-- [k] <- [k] - b_k/c_k-1 * [k-1]
105 b5d(i,j,k) = b5d(i,j,k)*c5d(i,j,k-1)
106 c5d(i,j,k) = c5d(i,j,k) - b5d(i,j,k)*d5d(i,j,k-1)
107 d5d(i,j,k) = d5d(i,j,k) - b5d(i,j,k)*e5d(i,j,k-1)
108 y5d(i,j,k,bi,bj) = y5d(i,j,k,bi,bj)
109 & - b5d(i,j,k)*y5d(i,j,k-1,bi,bj)
110 & - a5d(i,j,k)*y5d(i,j,k-2,bi,bj)
111 IF ( c5d(i,j,k).NE.0. _d 0 ) THEN
112 c5d(i,j,k) = 1. _d 0 / c5d(i,j,k)
113 ELSE
114 c5d(i,j,k) = 0. _d 0
115 errCode = 1
116 ENDIF
117 ENDDO
118 ENDDO
119 C- end if k= .. ; end of k loop
120 ENDIF
121 ENDDO
122
123 C-- Backward sweep (starting from bottom)
124 DO k=Nr,1,-1
125 IF (k.EQ.Nr) THEN
126 DO j=jMin,jMax
127 DO i=iMin,iMax
128 y5d(i,j,k,bi,bj) = y5d(i,j,k,bi,bj)*c5d(i,j,k)
129 ENDDO
130 ENDDO
131 ELSEIF (k.EQ.Nr-1) THEN
132 DO j=jMin,jMax
133 DO i=iMin,iMax
134 y5d(i,j,k,bi,bj) = ( y5d(i,j,k,bi,bj)
135 & - d5d(i,j,k)*y5d(i,j,k+1,bi,bj)
136 & )*c5d(i,j,k)
137 ENDDO
138 ENDDO
139 ELSE
140 DO j=jMin,jMax
141 DO i=iMin,iMax
142 y5d(i,j,k,bi,bj) = ( y5d(i,j,k,bi,bj)
143 & - d5d(i,j,k)*y5d(i,j,k+1,bi,bj)
144 & - e5d(i,j,k)*y5d(i,j,k+2,bi,bj)
145 & )*c5d(i,j,k)
146 ENDDO
147 ENDDO
148 C- end if k= .. ; end of k loop
149 ENDIF
150 ENDDO
151
152 RETURN
153 END
  ViewVC Help
Powered by ViewVC 1.1.22