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 |