pkg/streamice/streamice_tridiag_solve.F

5ca83cd8f7 Dani*0001 #include "STREAMICE_OPTIONS.h"
                0002
                0003 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0004
                0005 CBOP
                0006
887f427c62 Jean*0007       SUBROUTINE STREAMICE_TRIDIAG_SOLVE(
5ca83cd8f7 Dani*0008      U                               cg_Uin,     ! x-velocities
                0009      U                               cg_Vin,     ! x-velocities
                0010      U                               cg_Bu,      ! force in x dir
                0011      I                               A_uu,       ! section of matrix that multiplies u and projects on u
                0012      I                               umask,
                0013      I                               myThid )
                0014 C     /============================================================\
887f427c62 Jean*0015 C     | SUBROUTINE                                                 |
5ca83cd8f7 Dani*0016 C     | o                                                          |
                0017 C     |============================================================|
                0018 C     |                                                            |
                0019 C     \============================================================/
                0020       IMPLICIT NONE
                0021
                0022 #include "SIZE.h"
                0023 #include "EEPARAMS.h"
                0024 #include "PARAMS.h"
                0025 #include "STREAMICE.h"
                0026 #include "STREAMICE_CG.h"
                0027
                0028       _RL cg_Uin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0029       _RL cg_Vin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0030       _RL cg_Bu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0031       _RL A_uu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1)
                0032       _RS umask (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0033       INTEGER myThid
                0034
                0035       INTEGER iMin,iMax,i,j,k
                0036       _RL aMat(1:Nx)
                0037       _RL bMat(1:Nx)
                0038       _RL cMat(1:Nx)
                0039       _RL yMat(1:Nx)
                0040       _RL bet(1:Nx)
                0041       _RL tmpvar
                0042       INTEGER errCode
887f427c62 Jean*0043
5ca83cd8f7 Dani*0044
                0045 !      CALL WRITE_FLD_XY_RL ("taud_tri","",cg_Bu,0,mythid)
                0046 !      CALL WRITE_FLD_XY_RL ("A_m1m1","",A_uu(:,:,:,:,-1,-1),0,mythid)
                0047 !      CALL WRITE_FLD_XY_RL ("A_m1_0","",A_uu(:,:,:,:,-1,0),0,mythid)
                0048 !      CALL WRITE_FLD_XY_RL ("A_m1p1","",A_uu(:,:,:,:,-1,1),0,mythid)
                0049 !      CALL WRITE_FLD_XY_RL ("A_0_m1","",A_uu(:,:,:,:,0,-1),0,mythid)
                0050 !      CALL WRITE_FLD_XY_RL ("A_0_0","",A_uu(:,:,:,:,0,0),0,mythid)
                0051 !      CALL WRITE_FLD_XY_RL ("A_0_p1","",A_uu(:,:,:,:,0,1),0,mythid)
                0052 !      CALL WRITE_FLD_XY_RL ("A_p1m1","",A_uu(:,:,:,:,1,-1),0,mythid)
                0053 !      CALL WRITE_FLD_XY_RL ("A_p1_0","",A_uu(:,:,:,:,1,0),0,mythid)
                0054 !      CALL WRITE_FLD_XY_RL ("A_p1p1","",A_uu(:,:,:,:,1,1),0,mythid)
                0055
                0056
887f427c62 Jean*0057
5ca83cd8f7 Dani*0058       IF (nPx.gt.1 .or. nSx.gt.1) THEN
                0059        STOP 'must be serial for tridiag solve'
                0060       ENDIF
                0061
                0062       errCode = 0
                0063
                0064       imax = 0
                0065       iMin = 2
                0066       do i=imin,Nx
                0067        if (umask(i,1,1,1).eq.1.0) THEN
                0068
                0069         aMat(i)=0.0
                0070         bmat(i)=0.0
                0071         cmat(i)=0.0
                0072         ymat(i)=0.0
                0073         do j=-1,1
                0074         do k=1,3
                0075         aMat(i) = amat(i)+A_uu(i,k,1,1,-1,j)
                0076         bMat(i) = bmat(i)+A_uu(i,k,1,1,0,j)
                0077         cMat(i) = cmat(i)+A_uu(i,k,1,1,1,j)
                0078         enddo
                0079         yMat(i) = ymat(i)+cg_Bu(i,j+2,1,1)
                0080         enddo
                0081        else
                0082         iMax = i-1
                0083         exit
                0084        endif
                0085       enddo
887f427c62 Jean*0086
5ca83cd8f7 Dani*0087       IF(imax.eq.0) THEN
                0088        imax=Nx
                0089       ENDIF
                0090
887f427c62 Jean*0091
                0092       IF ( bMat(imin).NE.0. _d 0 ) THEN
5ca83cd8f7 Dani*0093        bet(imin) = 1. _d 0 / bMat(imin)
                0094       ELSE
                0095        bet(imin) = 0. _d 0
                0096        errCode = 1
                0097       ENDIF
887f427c62 Jean*0098
5ca83cd8f7 Dani*0099       DO i=imin+1,imax
                0100        tmpvar = bmat(i) - amat(i)*cmat(i-1)*bet(i-1)
                0101        IF ( tmpvar .NE. 0. _d 0 ) THEN
                0102         bet(i) = 1. _d 0 / tmpvar
                0103        ELSE
                0104         bet(i) = 0. _d 0
                0105         errCode = 1
                0106        ENDIF
                0107       ENDDO
                0108
                0109
                0110       ymat(imin) = ymat(imin)*bet(imin)
                0111
                0112       DO i=imin+1,imax
887f427c62 Jean*0113        ymat(i) = ( ymat(i)
5ca83cd8f7 Dani*0114      &            - amat(i)*ymat(i-1)
887f427c62 Jean*0115      &            )*bet(i)
5ca83cd8f7 Dani*0116       ENDDO
                0117
                0118
                0119       DO i=imax-1,imin,-1
                0120           ymat(i) = ymat(i)
                0121      &     - cmat(i)*bet(i)*ymat(i+1)
                0122       ENDDO
                0123
                0124       DO j=1,sNy
                0125        DO i=imin,imax
                0126         cg_Uin (i,j,1,1) = ymat(i)
                0127        ENDDO
                0128       ENDDO
                0129
                0130       DO j=1,sNy
                0131        DO i=1,sNx
                0132         cg_Vin (i,j,1,1) = 0. _d 0
                0133        ENDDO
                0134       ENDDO
                0135
                0136       print *, "ERRORCODE", errcode
                0137
                0138       RETURN
                0139       END