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