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

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

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

Revision 1.2 - (hide 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 jmc 1.2 C $Header: /u/gcmpack/MITgcm/model/src/solve_pentadiagonal.F,v 1.1 2004/01/07 21:19:37 jmc Exp $
2 jmc 1.1 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 jmc 1.2 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 jmc 1.1 C OUTPUT:
43 jmc 1.2 C y5d :: X = solution of A*X=Y
44 jmc 1.1 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 jmc 1.2 DO k=1,Nr
65 jmc 1.1 C-- forward sweep (starting from top)
66 jmc 1.2 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 jmc 1.1 c5d(i,j,k) = 1. _d 0 / c5d(i,j,k)
71 jmc 1.2 ELSE
72 jmc 1.1 c5d(i,j,k) = 0. _d 0
73     errCode = 1
74 jmc 1.2 ENDIF
75     ENDDO
76     ENDDO
77 jmc 1.1
78 jmc 1.2 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 jmc 1.1 c5d(i,j,k) = 1. _d 0 / c5d(i,j,k)
89 jmc 1.2 ELSE
90 jmc 1.1 c5d(i,j,k) = 0. _d 0
91     errCode = 1
92 jmc 1.2 ENDIF
93     ENDDO
94     ENDDO
95 jmc 1.1
96 jmc 1.2 ELSE
97 jmc 1.1 C-- Middle of forward sweep
98 jmc 1.2 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 jmc 1.1 c5d(i,j,k) = 1. _d 0 / c5d(i,j,k)
113 jmc 1.2 ELSE
114 jmc 1.1 c5d(i,j,k) = 0. _d 0
115     errCode = 1
116 jmc 1.2 ENDIF
117     ENDDO
118 jmc 1.1 ENDDO
119 jmc 1.2 C- end if k= .. ; end of k loop
120     ENDIF
121 jmc 1.1 ENDDO
122    
123     C-- Backward sweep (starting from bottom)
124 jmc 1.2 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 jmc 1.1 ENDDO
148 jmc 1.2 C- end if k= .. ; end of k loop
149     ENDIF
150 jmc 1.1 ENDDO
151    
152     RETURN
153     END
  ViewVC Help
Powered by ViewVC 1.1.22