/[MITgcm]/MITgcm_contrib/dgoldberg/streamice/streamice_tridiag_solve.F
ViewVC logotype

Annotation of /MITgcm_contrib/dgoldberg/streamice/streamice_tridiag_solve.F

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


Revision 1.1 - (hide annotations) (download)
Tue May 28 22:32:39 2013 UTC (12 years, 1 month ago) by dgoldberg
Branch: MAIN
tridiagonal matrix solver for flowline case

1 dgoldberg 1.1 #include "STREAMICE_OPTIONS.h"
2    
3     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
4    
5     CBOP
6    
7     SUBROUTINE STREAMICE_TRIDIAG_SOLVE(
8     U cg_Uin, ! x-velocities
9     U cg_Vin, ! x-velocities
10     U cg_Bu, ! force in x dir
11     I A_uu, ! section of matrix that multiplies u and projects on u
12     I umask,
13     I myThid )
14     C /============================================================\
15     C | SUBROUTINE |
16     C | o |
17     C |============================================================|
18     C | |
19     C \============================================================/
20     IMPLICIT NONE
21    
22     #include "SIZE.h"
23     #include "EEPARAMS.h"
24     #include "PARAMS.h"
25     #include "STREAMICE.h"
26     #include "STREAMICE_CG.h"
27    
28     _RL cg_Uin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
29     _RL cg_Vin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
30     _RL cg_Bu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
31     _RL A_uu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1)
32     _RS umask (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
33     INTEGER myThid
34    
35     INTEGER iMin,iMax,i,j
36     _RL aMat(1:Nx)
37     _RL bMat(1:Nx)
38     _RL cMat(1:Nx)
39     _RL yMat(1:Nx)
40     _RL bet(1:Nx)
41     _RL tmpvar
42     INTEGER errCode
43    
44    
45     IF (nPx.gt.1 .or. nSx.gt.1) THEN
46     STOP 'must be serial for tridiag solve'
47     ENDIF
48    
49     errCode = 0
50    
51     imax = 0
52     iMin = 2
53     do i=imin,Nx
54     if (umask(i,1,1,1).eq.1.0) THEN
55    
56     aMat(i)=0.0
57     bmat(i)=0.0
58     cmat(i)=0.0
59     ymat(i)=0.0
60     do j=-1,1
61     aMat(i) = amat(i)+A_uu(i,1,1,1,-1,j)
62     bMat(i) = bmat(i)+A_uu(i,1,1,1,0,j)
63     cMat(i) = cmat(i)+A_uu(i,1,1,1,1,j)
64     enddo
65     yMat(i) = ymat(i)+cg_Bu(i,1,1,1)
66     else
67     iMax = i-1
68     exit
69     endif
70     enddo
71    
72     IF(imax.eq.0) THEN
73     imax=Nx
74     ENDIF
75    
76    
77     IF ( bMat(imin).NE.0. _d 0 ) THEN
78     bet(imin) = 1. _d 0 / bMat(imin)
79     ELSE
80     bet(imin) = 0. _d 0
81     errCode = 1
82     ENDIF
83    
84     DO i=imin+1,imax
85     tmpvar = bmat(i) - amat(i)*cmat(i-1)*bet(i-1)
86     IF ( tmpvar .NE. 0. _d 0 ) THEN
87     bet(i) = 1. _d 0 / tmpvar
88     ELSE
89     bet(i) = 0. _d 0
90     errCode = 1
91     ENDIF
92     ENDDO
93    
94    
95     ymat(imin) = ymat(imin)*bet(imin)
96    
97     DO i=imin+1,imax
98     ymat(i) = ( ymat(i)
99     & - amat(i)*ymat(i-1)
100     & )*bet(i)
101     ENDDO
102    
103    
104     DO i=imax-1,imin,-1
105     ymat(i) = ymat(i)
106     & - cmat(i)*bet(i)*ymat(i+1)
107     ENDDO
108    
109     DO j=1,sNy
110     DO i=imin,imax
111     cg_Uin (i,j,1,1) = ymat(i)
112     ENDDO
113     ENDDO
114    
115     DO j=1,sNy
116     DO i=1,sNx
117     cg_Vin (i,j,1,1) = 0. _d 0
118     ENDDO
119     ENDDO
120    
121     print *, "ERRORCODE", errcode
122    
123     RETURN
124     END

  ViewVC Help
Powered by ViewVC 1.1.22