/[MITgcm]/MITgcm/pkg/streamice/streamice_tridiag_solve.F
ViewVC logotype

Annotation of /MITgcm/pkg/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 Jun 21 20:49:51 2013 UTC (10 years, 11 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64o, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint64n, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65o, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, checkpoint64k, checkpoint65, checkpoint64j, checkpoint64m, checkpoint64l, HEAD
Changes since 1.1: +13 -10 lines
add CVS Header and Name

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

  ViewVC Help
Powered by ViewVC 1.1.22