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

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

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

revision 1.1 by dgoldberg, Tue May 28 22:32:39 2013 UTC revision 1.3 by dgoldberg, Wed Aug 27 19:29:14 2014 UTC
# Line 1  Line 1 
1    C $Header$
2    C $Name$
3    
4  #include "STREAMICE_OPTIONS.h"  #include "STREAMICE_OPTIONS.h"
5    
6  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7    
8  CBOP  CBOP
9    
10        SUBROUTINE STREAMICE_TRIDIAG_SOLVE(        SUBROUTINE STREAMICE_TRIDIAG_SOLVE(
11       U                               cg_Uin,     ! x-velocities       U                               cg_Uin,     ! x-velocities
12       U                               cg_Vin,     ! x-velocities       U                               cg_Vin,     ! x-velocities
13       U                               cg_Bu,      ! force in x dir       U                               cg_Bu,      ! force in x dir
# Line 12  CBOP Line 15  CBOP
15       I                               umask,       I                               umask,
16       I                               myThid )       I                               myThid )
17  C     /============================================================\  C     /============================================================\
18  C     | SUBROUTINE                                                 |    C     | SUBROUTINE                                                 |
19  C     | o                                                          |  C     | o                                                          |
20  C     |============================================================|  C     |============================================================|
21  C     |                                                            |  C     |                                                            |
# Line 32  C     \================================= Line 35  C     \=================================
35        _RS umask (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RS umask (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
36        INTEGER myThid        INTEGER myThid
37    
38        INTEGER iMin,iMax,i,j        INTEGER iMin,iMax,i,j,k
39        _RL aMat(1:Nx)        _RL aMat(1:Nx)
40        _RL bMat(1:Nx)        _RL bMat(1:Nx)
41        _RL cMat(1:Nx)        _RL cMat(1:Nx)
# Line 40  C     \================================= Line 43  C     \=================================
43        _RL bet(1:Nx)        _RL bet(1:Nx)
44        _RL tmpvar        _RL tmpvar
45        INTEGER errCode        INTEGER errCode
46    
47          
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    
61        IF (nPx.gt.1 .or. nSx.gt.1) THEN        IF (nPx.gt.1 .or. nSx.gt.1) THEN
62         STOP 'must be serial for tridiag solve'         STOP 'must be serial for tridiag solve'
63        ENDIF        ENDIF
# Line 58  C     \================================= Line 74  C     \=================================
74          cmat(i)=0.0          cmat(i)=0.0
75          ymat(i)=0.0          ymat(i)=0.0
76          do j=-1,1          do j=-1,1
77          aMat(i) = amat(i)+A_uu(i,1,1,1,-1,j)          do k=1,3
78          bMat(i) = bmat(i)+A_uu(i,1,1,1,0,j)          aMat(i) = amat(i)+A_uu(i,k,1,1,-1,j)
79          cMat(i) = cmat(i)+A_uu(i,1,1,1,1,j)          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          enddo
         yMat(i) = ymat(i)+cg_Bu(i,1,1,1)  
84         else         else
85          iMax = i-1          iMax = i-1
86          exit          exit
87         endif         endif
88        enddo        enddo
89          
90        IF(imax.eq.0) THEN        IF(imax.eq.0) THEN
91         imax=Nx         imax=Nx
92        ENDIF        ENDIF
   
93    
94        IF ( bMat(imin).NE.0. _d 0 ) THEN  
95          IF ( bMat(imin).NE.0. _d 0 ) THEN
96         bet(imin) = 1. _d 0 / bMat(imin)         bet(imin) = 1. _d 0 / bMat(imin)
97        ELSE        ELSE
98         bet(imin) = 0. _d 0         bet(imin) = 0. _d 0
99         errCode = 1         errCode = 1
100        ENDIF        ENDIF
101          
102        DO i=imin+1,imax        DO i=imin+1,imax
103         tmpvar = bmat(i) - amat(i)*cmat(i-1)*bet(i-1)         tmpvar = bmat(i) - amat(i)*cmat(i-1)*bet(i-1)
104         IF ( tmpvar .NE. 0. _d 0 ) THEN         IF ( tmpvar .NE. 0. _d 0 ) THEN
# Line 95  C     \================================= Line 113  C     \=================================
113        ymat(imin) = ymat(imin)*bet(imin)        ymat(imin) = ymat(imin)*bet(imin)
114    
115        DO i=imin+1,imax        DO i=imin+1,imax
116         ymat(i) = ( ymat(i)         ymat(i) = ( ymat(i)
117       &            - amat(i)*ymat(i-1)       &            - amat(i)*ymat(i-1)
118       &            )*bet(i)       &            )*bet(i)
119        ENDDO        ENDDO
120    
121    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22