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

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

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

revision 1.1 by jmc, Wed Jan 7 21:19:37 2004 UTC revision 1.2 by jmc, Thu Jan 8 16:18:11 2004 UTC
# Line 30  C     == Global data == Line 30  C     == Global data ==
30  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
31  C     == Routine Arguments ==  C     == Routine Arguments ==
32  C     INPUT:  C     INPUT:
33  C     a5d :: 2nd  lower diagonal of the pentadiagonal matrix  C     iMin,iMax,jMin,jMax  :: computational domain
34  C     b5d :: 1rst lower diagonal of the pentadiagonal matrix  C     a5d    :: 2nd  lower diagonal of the pentadiagonal matrix
35  C     c5d :: main diagonal       of the pentadiagonal matrix  C     b5d    :: 1rst lower diagonal of the pentadiagonal matrix
36  C     d5d :: 1rst upper diagonal of the pentadiagonal matrix  C     c5d    :: main diagonal       of the pentadiagonal matrix
37  C     e5d :: 2nd  upper diagonal of the pentadiagonal matrix  C     d5d    :: 1rst upper diagonal of the pentadiagonal matrix
38  C     y5d :: Y vector (R.H.S.);  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:  C     OUTPUT:
43  C     y5d :: X = solution of A*X=Y  C     y5d    :: X = solution of A*X=Y
44  C     a5d,b5d,c5d,d5d,e5d :: modified to enable to find X' solution of  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  C                        A*X'=Y' without solving the full system again
46  C     errCode :: > 0 if singular matrix  C     errCode :: > 0 if singular matrix
# Line 58  CEOP Line 61  CEOP
61    
62        errCode = 0        errCode = 0
63    
64          DO k=1,Nr
65  C--   forward sweep (starting from top)  C--   forward sweep (starting from top)
66        k = 1         IF (k.EQ.1) THEN
67        DO j=jMin,jMax          DO j=jMin,jMax
68         DO i=iMin,iMax           DO i=iMin,iMax
69           IF ( c5d(i,j,k).NE.0. _d 0 ) THEN            IF ( c5d(i,j,k).NE.0. _d 0 ) THEN
70             c5d(i,j,k) = 1. _d 0 / c5d(i,j,k)             c5d(i,j,k) = 1. _d 0 / c5d(i,j,k)
71           ELSE            ELSE
72             c5d(i,j,k) = 0. _d 0             c5d(i,j,k) = 0. _d 0
73             errCode = 1             errCode = 1
74           ENDIF            ENDIF
75         ENDDO           ENDDO
76        ENDDO          ENDDO
77    
78        k = 2         ELSEIF (k.EQ.2) THEN
79        DO j=jMin,jMax          DO j=jMin,jMax
80         DO i=iMin,iMax           DO i=iMin,iMax
81  C--      [k] <- [k] - b_k/c_k-1 * [k-1]  C--       [k] <- [k] - b_k/c_k-1 * [k-1]
82           b5d(i,j,k) = b5d(i,j,k)*c5d(i,j,k-1)            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)            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)            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)            y5d(i,j,k,bi,bj) = y5d(i,j,k,bi,bj)
86       &                    - b5d(i,j,k)*y5d(i,j,k-1,bi,bj)       &                     - b5d(i,j,k)*y5d(i,j,k-1,bi,bj)
87           IF ( c5d(i,j,k).NE.0. _d 0 ) THEN            IF ( c5d(i,j,k).NE.0. _d 0 ) THEN
88             c5d(i,j,k) = 1. _d 0 / c5d(i,j,k)             c5d(i,j,k) = 1. _d 0 / c5d(i,j,k)
89           ELSE            ELSE
90             c5d(i,j,k) = 0. _d 0             c5d(i,j,k) = 0. _d 0
91             errCode = 1             errCode = 1
92           ENDIF            ENDIF
93         ENDDO           ENDDO
94        ENDDO          ENDDO
95    
96           ELSE
97  C--   Middle of forward sweep  C--   Middle of forward sweep
98        DO k=3,Nr          DO j=jMin,jMax
99         DO j=jMin,jMax           DO i=iMin,iMax
100          DO i=iMin,iMax  C--       [k] <- [k] - a_k/c_k-2 * [k-2]
101  C--      [k] <- [k] - a_k/c_k-2 * [k-2]            a5d(i,j,k) = a5d(i,j,k)*c5d(i,j,k-2)
102           a5d(i,j,k) = a5d(i,j,k)*c5d(i,j,k-2)            b5d(i,j,k) = b5d(i,j,k) - a5d(i,j,k)*d5d(i,j,k-2)
103           b5d(i,j,k) = b5d(i,j,k) - a5d(i,j,k)*d5d(i,j,k-2)            c5d(i,j,k) = c5d(i,j,k) - a5d(i,j,k)*e5d(i,j,k-2)
104           c5d(i,j,k) = c5d(i,j,k) - a5d(i,j,k)*e5d(i,j,k-2)  C--       [k] <- [k] - b_k/c_k-1 * [k-1]
105  C--      [k] <- [k] - b_k/c_k-1 * [k-1]            b5d(i,j,k) = b5d(i,j,k)*c5d(i,j,k-1)
106           b5d(i,j,k) = b5d(i,j,k)*c5d(i,j,k-1)            c5d(i,j,k) = c5d(i,j,k) - b5d(i,j,k)*d5d(i,j,k-1)
107           c5d(i,j,k) = c5d(i,j,k) - b5d(i,j,k)*d5d(i,j,k-1)            d5d(i,j,k) = d5d(i,j,k) - b5d(i,j,k)*e5d(i,j,k-1)
108           d5d(i,j,k) = d5d(i,j,k) - b5d(i,j,k)*e5d(i,j,k-1)            y5d(i,j,k,bi,bj) = y5d(i,j,k,bi,bj)
109           y5d(i,j,k,bi,bj) = y5d(i,j,k,bi,bj)       &                     - b5d(i,j,k)*y5d(i,j,k-1,bi,bj)
110       &                    - b5d(i,j,k)*y5d(i,j,k-1,bi,bj)       &                     - a5d(i,j,k)*y5d(i,j,k-2,bi,bj)
111       &                    - a5d(i,j,k)*y5d(i,j,k-2,bi,bj)            IF ( c5d(i,j,k).NE.0. _d 0 ) THEN
          IF ( c5d(i,j,k).NE.0. _d 0 ) THEN  
112             c5d(i,j,k) = 1. _d 0 / c5d(i,j,k)             c5d(i,j,k) = 1. _d 0 / c5d(i,j,k)
113           ELSE            ELSE
114             c5d(i,j,k) = 0. _d 0             c5d(i,j,k) = 0. _d 0
115             errCode = 1             errCode = 1
116           ENDIF            ENDIF
117             ENDDO
118          ENDDO          ENDDO
119         ENDDO  C-      end if k= .. ; end of k loop
120           ENDIF
121        ENDDO        ENDDO
122    
123  C--   Backward sweep (starting from bottom)  C--   Backward sweep (starting from bottom)
124        k = Nr        DO k=Nr,1,-1
125        DO j=jMin,jMax         IF (k.EQ.Nr) THEN
126         DO i=iMin,iMax          DO j=jMin,jMax
127           y5d(i,j,k,bi,bj) =   y5d(i,j,k,bi,bj)*c5d(i,j,k)           DO i=iMin,iMax
128         ENDDO            y5d(i,j,k,bi,bj) =   y5d(i,j,k,bi,bj)*c5d(i,j,k)
129        ENDDO           ENDDO
130        k = Nr-1          ENDDO
131        DO j=jMin,jMax         ELSEIF (k.EQ.Nr-1) THEN
132         DO i=iMin,iMax          DO j=jMin,jMax
133           y5d(i,j,k,bi,bj) = ( y5d(i,j,k,bi,bj)           DO i=iMin,iMax
134       &                      - d5d(i,j,k)*y5d(i,j,k+1,bi,bj)            y5d(i,j,k,bi,bj) = ( y5d(i,j,k,bi,bj)
135       &                      )*c5d(i,j,k)       &                       - d5d(i,j,k)*y5d(i,j,k+1,bi,bj)
136         ENDDO       &                       )*c5d(i,j,k)
137        ENDDO           ENDDO
138        DO k=Nr-2,1,-1          ENDDO
139         DO j=jMin,jMax         ELSE
140          DO i=iMin,iMax          DO j=jMin,jMax
141           y5d(i,j,k,bi,bj) = ( y5d(i,j,k,bi,bj)           DO i=iMin,iMax
142       &                      - d5d(i,j,k)*y5d(i,j,k+1,bi,bj)            y5d(i,j,k,bi,bj) = ( y5d(i,j,k,bi,bj)
143       &                      - e5d(i,j,k)*y5d(i,j,k+2,bi,bj)       &                       - d5d(i,j,k)*y5d(i,j,k+1,bi,bj)
144       &                      )*c5d(i,j,k)       &                       - e5d(i,j,k)*y5d(i,j,k+2,bi,bj)
145         &                       )*c5d(i,j,k)
146             ENDDO
147          ENDDO          ENDDO
148         ENDDO  C-      end if k= .. ; end of k loop
149           ENDIF
150        ENDDO        ENDDO
151    
152        RETURN        RETURN
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22