/[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.2 - (hide annotations) (download)
Fri May 31 17:33:03 2013 UTC (12 years, 1 month ago) by dgoldberg
Branch: MAIN
Changes since 1.1: +20 -5 lines
fixing bugs

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 dgoldberg 1.2 INTEGER iMin,iMax,i,j,k
36 dgoldberg 1.1 _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 dgoldberg 1.2
45     ! CALL WRITE_FLD_XY_RL ("taud_tri","",cg_Bu,0,mythid)
46     ! CALL WRITE_FLD_XY_RL ("A_m1m1","",A_uu(:,:,:,:,-1,-1),0,mythid)
47     ! CALL WRITE_FLD_XY_RL ("A_m1_0","",A_uu(:,:,:,:,-1,0),0,mythid)
48     ! CALL WRITE_FLD_XY_RL ("A_m1p1","",A_uu(:,:,:,:,-1,1),0,mythid)
49     ! CALL WRITE_FLD_XY_RL ("A_0_m1","",A_uu(:,:,:,:,0,-1),0,mythid)
50     ! CALL WRITE_FLD_XY_RL ("A_0_0","",A_uu(:,:,:,:,0,0),0,mythid)
51     ! CALL WRITE_FLD_XY_RL ("A_0_p1","",A_uu(:,:,:,:,0,1),0,mythid)
52     ! CALL WRITE_FLD_XY_RL ("A_p1m1","",A_uu(:,:,:,:,1,-1),0,mythid)
53     ! CALL WRITE_FLD_XY_RL ("A_p1_0","",A_uu(:,:,:,:,1,0),0,mythid)
54     ! CALL WRITE_FLD_XY_RL ("A_p1p1","",A_uu(:,:,:,:,1,1),0,mythid)
55    
56    
57 dgoldberg 1.1
58     IF (nPx.gt.1 .or. nSx.gt.1) THEN
59     STOP 'must be serial for tridiag solve'
60     ENDIF
61    
62     errCode = 0
63    
64     imax = 0
65     iMin = 2
66     do i=imin,Nx
67     if (umask(i,1,1,1).eq.1.0) THEN
68    
69     aMat(i)=0.0
70     bmat(i)=0.0
71     cmat(i)=0.0
72     ymat(i)=0.0
73     do j=-1,1
74 dgoldberg 1.2 do k=1,3
75     aMat(i) = amat(i)+A_uu(i,k,1,1,-1,j)
76     bMat(i) = bmat(i)+A_uu(i,k,1,1,0,j)
77     cMat(i) = cmat(i)+A_uu(i,k,1,1,1,j)
78     enddo
79     yMat(i) = ymat(i)+cg_Bu(i,j+2,1,1)
80 dgoldberg 1.1 enddo
81     else
82     iMax = i-1
83     exit
84     endif
85     enddo
86    
87     IF(imax.eq.0) THEN
88     imax=Nx
89     ENDIF
90    
91    
92     IF ( bMat(imin).NE.0. _d 0 ) THEN
93     bet(imin) = 1. _d 0 / bMat(imin)
94     ELSE
95     bet(imin) = 0. _d 0
96     errCode = 1
97     ENDIF
98    
99     DO i=imin+1,imax
100     tmpvar = bmat(i) - amat(i)*cmat(i-1)*bet(i-1)
101     IF ( tmpvar .NE. 0. _d 0 ) THEN
102     bet(i) = 1. _d 0 / tmpvar
103     ELSE
104     bet(i) = 0. _d 0
105     errCode = 1
106     ENDIF
107     ENDDO
108    
109    
110     ymat(imin) = ymat(imin)*bet(imin)
111    
112     DO i=imin+1,imax
113     ymat(i) = ( ymat(i)
114     & - amat(i)*ymat(i-1)
115     & )*bet(i)
116     ENDDO
117    
118    
119     DO i=imax-1,imin,-1
120     ymat(i) = ymat(i)
121     & - cmat(i)*bet(i)*ymat(i+1)
122     ENDDO
123    
124     DO j=1,sNy
125     DO i=imin,imax
126     cg_Uin (i,j,1,1) = ymat(i)
127     ENDDO
128     ENDDO
129    
130     DO j=1,sNy
131     DO i=1,sNx
132     cg_Vin (i,j,1,1) = 0. _d 0
133     ENDDO
134     ENDDO
135    
136     print *, "ERRORCODE", errcode
137    
138     RETURN
139     END

  ViewVC Help
Powered by ViewVC 1.1.22