Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-02 18:37:06 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
1eaf9121ed Jean*0001 #include "CPP_OPTIONS.h"
                0002 
                0003 CBOP
                0004 C     !ROUTINE: SOLVE_UV_TRIDIAGO
                0005 C     !INTERFACE:
                0006       SUBROUTINE SOLVE_UV_TRIDIAGO(
                0007      I                     kSize, ols, solve4u, solve4v,
                0008      I                     aU, bU, cU, rhsU,
                0009      I                     aV, bV, cV, rhsV,
                0010      O                     uFld, vFld,
                0011      O                     errCode,
                0012      I                     subIter, myIter, myThid )
                0013 C     !DESCRIPTION: \bv
                0014 C     *==========================================================*
                0015 C     | S/R SOLVE_UV_TRIDIAGO
                0016 C     | o Solve a pair of tri-diagonal system along X and Y lines
                0017 C     |   (in X-dir for uFld and in Y-dir for vFld)
                0018 C     *==========================================================*
                0019 C     | o Used, e.g., in linear part of seaice LSR solver
                0020 C     *==========================================================*
                0021 C     \ev
                0022 
                0023 C     !USES:
                0024       IMPLICIT NONE
                0025 C     == Global data ==
                0026 #include "SIZE.h"
                0027 #include "EEPARAMS.h"
                0028 
                0029 C     !INPUT/OUTPUT PARAMETERS:
                0030 C     == Routine Arguments ==
                0031 C     kSize    :: size in 3rd dimension
                0032 C     ols      :: size of overlap (of input arg array)
                0033 C     solve4u  :: logical flag, do solve for u-component if true
                0034 C     solve4v  :: logical flag, do solve for v-component if true
                0035 C     aU,bU,cU :: u-matrix (lower diagonal, main diagonal & upper diagonal)
                0036 C     rhsU     :: RHS vector (u-component)
                0037 C     aV,bV,cV :: v-matrix (lower diagonal, main diagonal & upper diagonal)
                0038 C     rhsV     :: RHS vector (v-component)
                0039 C     uFld     :: X = solution of: A_u * X = rhsU
                0040 C     vFld     :: X = solution of: A_v * X = rhsV
                0041 C     errCode  :: > 0 if singular matrix
                0042 C     subIter  :: current sub-iteration number
                0043 C     myIter   :: current iteration number
                0044 C     myThid   :: my Thread Id number
                0045       INTEGER kSize, ols
                0046       LOGICAL solve4u, solve4v
                0047       _RL  aU (1-ols:sNx+ols,1-ols:sNy+ols,kSize,nSx,nSy)
                0048       _RL  bU (1-ols:sNx+ols,1-ols:sNy+ols,kSize,nSx,nSy)
                0049       _RL  cU (1-ols:sNx+ols,1-ols:sNy+ols,kSize,nSx,nSy)
                0050       _RL rhsU(1-ols:sNx+ols,1-ols:sNy+ols,kSize,nSx,nSy)
                0051       _RL  aV (1-ols:sNx+ols,1-ols:sNy+ols,kSize,nSx,nSy)
                0052       _RL  bV (1-ols:sNx+ols,1-ols:sNy+ols,kSize,nSx,nSy)
                0053       _RL  cV (1-ols:sNx+ols,1-ols:sNy+ols,kSize,nSx,nSy)
                0054       _RL rhsV(1-ols:sNx+ols,1-ols:sNy+ols,kSize,nSx,nSy)
                0055       _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
                0056       _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
                0057       INTEGER errCode
                0058       INTEGER subIter, myIter, myThid
                0059 
                0060 C     !SHARED LOCAL VARIABLES:
                0061 C     aTu, cTu, yTu :: tile edges coeff and RHS for u-component
                0062 C     aTv, cTv, yTv :: tile edges coeff and RHS for v-component
                0063       COMMON /SOLVE_UV_3DIAG_LOCAL/
                0064      &  aTu, cTu, yTu, aTv, cTv, yTv
                0065       _RL aTu(2,1:sNy,nSx,nSy)
                0066       _RL cTu(2,1:sNy,nSx,nSy)
                0067       _RL yTu(2,1:sNy,nSx,nSy)
                0068       _RL aTv(2,1:sNx,nSx,nSy)
                0069       _RL cTv(2,1:sNx,nSx,nSy)
                0070       _RL yTv(2,1:sNx,nSx,nSy)
                0071 
                0072 C     !LOCAL VARIABLES:
                0073 C     == Local variables ==
                0074       INTEGER bi, bj, bm, bp
                0075       INTEGER i,j,k
                0076       INTEGER ii, im, ip
                0077       INTEGER jj, jm, jp
                0078       _RL tmpVar
                0079       _RL uTmp1, uTmp2, vTmp1, vTmp2
                0080       _RL alpU(1:sNx,1:sNy,nSx,nSy)
                0081       _RL gamU(1:sNx,1:sNy,nSx,nSy)
                0082       _RL yy_U(1:sNx,1:sNy,nSx,nSy)
                0083       _RL alpV(1:sNx,1:sNy,nSx,nSy)
                0084       _RL gamV(1:sNx,1:sNy,nSx,nSy)
                0085       _RL yy_V(1:sNx,1:sNy,nSx,nSy)
                0086 CEOP
                0087 
                0088       errCode = 0
                0089       IF ( .NOT.solve4u .AND. .NOT.solve4v ) RETURN
                0090 
                0091 C--   outside loop on level number k
                0092       DO k = 1,kSize
                0093 
                0094        IF ( solve4u ) THEN
                0095         DO bj=myByLo(myThid),myByHi(myThid)
                0096          DO bi=myBxLo(myThid),myBxHi(myThid)
                0097 
                0098 C--   work on local copy:
                0099           DO j= 1,sNy
                0100            DO i= 1,sNx
                0101             alpU(i,j,bi,bj) = aU(i,j,k,bi,bj)
                0102             gamU(i,j,bi,bj) = cU(i,j,k,bi,bj)
                0103             yy_U(i,j,bi,bj) = rhsU(i,j,k,bi,bj)
                0104            ENDDO
                0105           ENDDO
                0106 
                0107 C--   Beginning of forward sweep (i=1)
                0108           i = 1
                0109           DO j= 1,sNy
                0110 C-    normalise row [1] ( 1 on main diagonal)
                0111             tmpVar = bU(i,j,k,bi,bj)
                0112             IF ( tmpVar.NE.0. _d 0 ) THEN
                0113              tmpVar = 1. _d 0 / tmpVar
                0114             ELSE
                0115              tmpVar = 0. _d 0
                0116              errCode = 1
                0117             ENDIF
                0118             gamU(i,j,bi,bj) = gamU(i,j,bi,bj)*tmpVar
                0119             alpU(i,j,bi,bj) = alpU(i,j,bi,bj)*tmpVar
                0120             yy_U(i,j,bi,bj) = yy_U(i,j,bi,bj)*tmpVar
                0121           ENDDO
                0122 
                0123 C--   Middle of forward sweep (i=2:sNx)
                0124           DO j= 1,sNy
                0125            DO i= 2,sNx
                0126             im = i-1
                0127 C-    update row [i] <-- [i] - alp_i * [i-1] and normalise (main diagonal = 1)
                0128             tmpVar = bU(i,j,k,bi,bj) - alpU(i,j,bi,bj)*gamU(im,j,bi,bj)
                0129             IF ( tmpVar.NE.0. _d 0 ) THEN
                0130              tmpVar = 1. _d 0 / tmpVar
                0131             ELSE
                0132              tmpVar = 0. _d 0
                0133              errCode = 1
                0134             ENDIF
                0135             yy_U(i,j,bi,bj) = ( yy_U(i,j,bi,bj)
                0136      &                        - alpU(i,j,bi,bj)*yy_U(im,j,bi,bj)
                0137      &                        )*tmpVar
                0138             gamU(i,j,bi,bj) =   gamU(i,j,bi,bj)*tmpVar
                0139             alpU(i,j,bi,bj) = - alpU(i,j,bi,bj)*alpU(im,j,bi,bj)*tmpVar
                0140            ENDDO
                0141           ENDDO
                0142 
                0143 C--   Backward sweep (i=sNx-1:-1:1)
                0144           DO j= 1,sNy
                0145            DO ii= 1,sNx-1
                0146             i = sNx - ii
                0147             ip = i+1
                0148 C-    update row [i] <-- [i] - gam_i * [i+1]
                0149             yy_U(i,j,bi,bj) =  yy_U(i,j,bi,bj)
                0150      &                       - gamU(i,j,bi,bj)*yy_U(ip,j,bi,bj)
                0151             alpU(i,j,bi,bj) =  alpU(i,j,bi,bj)
                0152      &                       - gamU(i,j,bi,bj)*alpU(ip,j,bi,bj)
                0153             gamU(i,j,bi,bj) = -gamU(i,j,bi,bj)*gamU(ip,j,bi,bj)
                0154            ENDDO
                0155           ENDDO
                0156 
                0157 C--    At this stage, the 3-diagonal system is reduced to Identity with two
                0158 C      more columns (alp & gam) corresponding to unknow X(i=0) and X(i=sNx+1):
                0159 C                                       X_0
                0160 C         alp  1 0    ...    0 0 gam    X_1        Y_1
                0161 C         alp  0 1    ...    0 0 gam    X_2        Y_2
                0162 C
                0163 C          .   . .    ...    . .  .      .          .
                0164 C       (  .   . .    ...    . .  .  )(  .   ) = (  .   )
                0165 C          .   . .    ...    . .  .      .          .
                0166 C
                0167 C         alp  0 0    ...    1 0 gam    X_n-1      Y_n-1
                0168 C         alp  0 0    ...    0 1 gam    X_n        Y_n
                0169 C                                       X_n+1
                0170 C-----
                0171 
                0172 C--   Store tile edges coeff: (1) <--> i=1 ; (2) <--> i=sNx
                0173           DO j= 1,sNy
                0174             aTu(1,j,bi,bj) = alpU( 1, j,bi,bj)
                0175             cTu(1,j,bi,bj) = gamU( 1, j,bi,bj)
                0176             yTu(1,j,bi,bj) = yy_U( 1, j,bi,bj)
                0177             aTu(2,j,bi,bj) = alpU(sNx,j,bi,bj)
                0178             cTu(2,j,bi,bj) = gamU(sNx,j,bi,bj)
                0179             yTu(2,j,bi,bj) = yy_U(sNx,j,bi,bj)
                0180           ENDDO
                0181 
                0182 C     end bi,bj-loops
                0183          ENDDO
                0184         ENDDO
                0185 
                0186 C--   Solve for tile edges values
                0187         IF ( nPx*nPy.GT.1 .OR. useCubedSphereExchange ) THEN
                0188           STOP 'ABNORMAL END: S/R SOLVE_UV_TRIDIAGO: missing code'
                0189         ENDIF
                0190         _BARRIER
                0191         _BEGIN_MASTER(myThid)
                0192         DO bj=1,nSy
                0193          DO j=1,sNy
                0194 
                0195           DO bi=2,nSx
                0196            bm = bi-1
                0197 C-    update row [1,bi] <- [1,bi] - a(1,bi)*[2,bi-1] (& normalise diag)
                0198             tmpVar = oneRL - aTu(1,j,bi,bj)*cTu(2,j,bm,bj)
                0199             IF ( tmpVar.NE.0. _d 0 ) THEN
                0200              tmpVar = 1. _d 0 / tmpVar
                0201             ELSE
                0202              tmpVar = 0. _d 0
                0203              errCode = 1
                0204             ENDIF
                0205             yTu(1,j,bi,bj) = ( yTu(1,j,bi,bj)
                0206      &                       - aTu(1,j,bi,bj)*yTu(2,j,bm,bj)
                0207      &                       )*tmpVar
                0208             cTu(1,j,bi,bj) =   cTu(1,j,bi,bj)*tmpVar
                0209             aTu(1,j,bi,bj) = - aTu(1,j,bi,bj)*aTu(2,j,bm,bj)*tmpVar
                0210 
                0211 C-    update row [2,bi] <- [2,bi] - a(2,bi)*[2,bi-1] + a(2,bi)*c(2,bi-1)*[1,bi]
                0212             tmpVar = aTu(2,j,bi,bj)*cTu(2,j,bm,bj)
                0213             yTu(2,j,bi,bj) =  yTu(2,j,bi,bj)
                0214      &                      - aTu(2,j,bi,bj)*yTu(2,j,bm,bj)
                0215      &                      + tmpVar*yTu(1,j,bi,bj)
                0216             cTu(2,j,bi,bj) =  cTu(2,j,bi,bj)
                0217      &                      + tmpVar*cTu(1,j,bi,bj)
                0218             aTu(2,j,bi,bj) = -aTu(2,j,bi,bj)*aTu(2,j,bm,bj)
                0219      &                      + tmpVar*aTu(1,j,bi,bj)
                0220           ENDDO
                0221 
                0222           DO bi=nSx-1,1,-1
                0223            bp = bi+1
                0224            DO i=1,2
                0225 C-    update row [1,bi] <- [1,bi] - c(1,bi)*[1,bi+1]
                0226 C-    update row [2,bi] <- [2,bi] - c(2,bi)*[1,bi+1]
                0227             aTu(i,j,bi,bj) =  aTu(i,j,bi,bj)
                0228      &                      - cTu(i,j,bi,bj)*aTu(1,j,bp,bj)
                0229             yTu(i,j,bi,bj) =  yTu(i,j,bi,bj)
                0230      &                      - cTu(i,j,bi,bj)*yTu(1,j,bp,bj)
                0231             cTu(i,j,bi,bj) = -cTu(i,j,bi,bj)*cTu(1,j,bp,bj)
                0232            ENDDO
                0233           ENDDO
                0234 
                0235 C--  periodic in X:  X_0 <=> X_Nx and X_(N+1) <=> X_1 ;
                0236 C    find the value at the 2 opposite location (i=1 and i=Nx)
                0237           bm = 1
                0238           bp = nSx
                0239           cTu(1,j,bm,bj) = oneRL + cTu(1,j,bm,bj)
                0240           aTu(2,j,bp,bj) = oneRL + aTu(2,j,bp,bj)
                0241           tmpVar = cTu(1,j,bm,bj) * aTu(2,j,bp,bj)
                0242      &           - aTu(1,j,bm,bj) * cTu(2,j,bp,bj)
                0243           IF ( tmpVar.NE.0. _d 0 ) THEN
                0244              tmpVar = 1. _d 0 / tmpVar
                0245           ELSE
                0246              tmpVar = 0. _d 0
                0247              errCode = 1
                0248           ENDIF
                0249           uTmp1 = ( aTu(2,j,bp,bj) * yTu(1,j,bm,bj)
                0250      &            - aTu(1,j,bm,bj) * yTu(2,j,bp,bj)
                0251      &            )*tmpVar
                0252           uTmp2 = ( cTu(1,j,bm,bj) * yTu(2,j,bp,bj)
                0253      &            - cTu(2,j,bp,bj) * yTu(1,j,bm,bj)
                0254      &            )*tmpVar
                0255 
                0256 C-    finalise tile-edges solution (put into RHS "yTu"):
                0257           DO bi=1,nSx
                0258            DO i=1,2
                0259             IF ( bi+i .EQ.2 ) THEN
                0260              yTu(i,j,bi,bj) = uTmp1
                0261             ELSEIF ( bi+i .EQ. nSx+2 ) THEN
                0262              yTu(i,j,bi,bj) = uTmp2
                0263             ELSE
                0264              yTu(i,j,bi,bj) = yTu(i,j,bi,bj)
                0265      &                      - aTu(i,j,bi,bj) * uTmp2
                0266      &                      - cTu(i,j,bi,bj) * uTmp1
                0267             ENDIF
                0268            ENDDO
                0269           ENDDO
                0270 
                0271          ENDDO
                0272         ENDDO
                0273         _END_MASTER(myThid)
                0274         _BARRIER
                0275 
                0276         DO bj=myByLo(myThid),myByHi(myThid)
                0277          DO bi=myBxLo(myThid),myBxHi(myThid)
                0278           bm = 1 + MOD(bi-2+nSx,nSx)
                0279           bp = 1 + MOD(bi-0+nSx,nSx)
                0280           DO j= 1,sNy
                0281            DO i= 1,sNx
                0282             uFld(i,j,k,bi,bj) = yy_U(i,j,bi,bj)
                0283      &                      - alpU(i,j,bi,bj) * yTu(2,j,bm,bj)
                0284      &                      - gamU(i,j,bi,bj) * yTu(1,j,bp,bj)
                0285            ENDDO
                0286           ENDDO
                0287          ENDDO
                0288         ENDDO
                0289 
                0290 C     end solve for uFld
                0291        ENDIF
                0292 
                0293        IF ( solve4v ) THEN
                0294         DO bj=myByLo(myThid),myByHi(myThid)
                0295          DO bi=myBxLo(myThid),myBxHi(myThid)
                0296 
                0297 C--   work on local copy:
                0298           DO j= 1,sNy
                0299            DO i= 1,sNx
                0300             alpV(i,j,bi,bj) = aV(i,j,k,bi,bj)
                0301             gamV(i,j,bi,bj) = cV(i,j,k,bi,bj)
                0302             yy_V(i,j,bi,bj) = rhsV(i,j,k,bi,bj)
                0303            ENDDO
                0304           ENDDO
                0305 
                0306 C--   Beginning of forward sweep (j=1)
                0307           j = 1
                0308           DO i= 1,sNx
                0309 C-    normalise row [1] ( 1 on main diagonal)
                0310             tmpVar = bV(i,j,k,bi,bj)
                0311             IF ( tmpVar.NE.0. _d 0 ) THEN
                0312              tmpVar = 1. _d 0 / tmpVar
                0313             ELSE
                0314              tmpVar = 0. _d 0
                0315              errCode = 1
                0316             ENDIF
                0317             gamV(i,j,bi,bj) = gamV(i,j,bi,bj)*tmpVar
                0318             alpV(i,j,bi,bj) = alpV(i,j,bi,bj)*tmpVar
                0319             yy_V(i,j,bi,bj) = yy_V(i,j,bi,bj)*tmpVar
                0320           ENDDO
                0321 
                0322 C--   Middle of forward sweep (j=2:sNy)
                0323           DO i= 1,sNx
                0324            DO j= 2,sNy
                0325             jm = j-1
                0326 C-    update row [j] <-- [j] - alp_j * [j-1] and normalise (main diagonal = 1)
                0327             tmpVar = bV(i,j,k,bi,bj) - alpV(i,j,bi,bj)*gamV(i,jm,bi,bj)
                0328             IF ( tmpVar.NE.0. _d 0 ) THEN
                0329              tmpVar = 1. _d 0 / tmpVar
                0330             ELSE
                0331              tmpVar = 0. _d 0
                0332              errCode = 1
                0333             ENDIF
                0334             yy_V(i,j,bi,bj) = ( yy_V(i,j,bi,bj)
                0335      &                        - alpV(i,j,bi,bj)*yy_V(i,jm,bi,bj)
                0336      &                        )*tmpVar
                0337             gamV(i,j,bi,bj) =   gamV(i,j,bi,bj)*tmpVar
                0338             alpV(i,j,bi,bj) = - alpV(i,j,bi,bj)*alpV(i,jm,bi,bj)*tmpVar
                0339            ENDDO
                0340           ENDDO
                0341 
                0342 C--   Backward sweep (j=sNy-1:-1:1)
                0343           DO i= 1,sNx
                0344            DO jj= 1,sNy-1
                0345             j = sNy - jj
                0346             jp = j+1
                0347 C-    update row [j] <-- [j] - gam_j * [j+1]
                0348             yy_V(i,j,bi,bj) =  yy_V(i,j,bi,bj)
                0349      &                       - gamV(i,j,bi,bj)*yy_V(i,jp,bi,bj)
                0350             alpV(i,j,bi,bj) =  alpV(i,j,bi,bj)
                0351      &                       - gamV(i,j,bi,bj)*alpV(i,jp,bi,bj)
                0352             gamV(i,j,bi,bj) = -gamV(i,j,bi,bj)*gamV(i,jp,bi,bj)
                0353            ENDDO
                0354           ENDDO
                0355 
                0356 C--    At this stage, the 3-diagonal system is reduced to Identity with two
                0357 C      more columns (alp & gam) corresponding to unknow X(j=0) and X(j=sNy+1)
                0358 
                0359 C--   Store tile edges coeff: (1) <--> j=1 ; (2) <--> j=sNy
                0360           DO i= 1,sNx
                0361             aTv(1,i,bi,bj) = alpV(i, 1, bi,bj)
                0362             cTv(1,i,bi,bj) = gamV(i, 1, bi,bj)
                0363             yTv(1,i,bi,bj) = yy_V(i, 1, bi,bj)
                0364             aTv(2,i,bi,bj) = alpV(i,sNy,bi,bj)
                0365             cTv(2,i,bi,bj) = gamV(i,sNy,bi,bj)
                0366             yTv(2,i,bi,bj) = yy_V(i,sNy,bi,bj)
                0367           ENDDO
                0368 
                0369 C     end bi,bj-loops
                0370          ENDDO
                0371         ENDDO
                0372 
                0373 C--   Solve for tile edges values
                0374         IF ( nPx*nPy.GT.1 .OR. useCubedSphereExchange ) THEN
                0375          STOP 'ABNORMAL END: S/R SOLVE_UV_TRIDIAGO: missing code'
                0376         ENDIF
                0377         _BARRIER
                0378         _BEGIN_MASTER(myThid)
                0379         DO bi=1,nSx
                0380          DO i=1,sNx
                0381 
                0382           DO bj=2,nSy
                0383            bm = bj-1
                0384 C-    update row [1,bj] <- [1,bj] - a(1,bj)*[2,bj-1] (& normalise diag)
                0385             tmpVar = oneRL - aTv(1,i,bi,bj)*cTv(2,i,bi,bm)
                0386             IF ( tmpVar.NE.0. _d 0 ) THEN
                0387              tmpVar = 1. _d 0 / tmpVar
                0388             ELSE
                0389              tmpVar = 0. _d 0
                0390              errCode = 1
                0391             ENDIF
                0392             yTv(1,i,bi,bj) = ( yTv(1,i,bi,bj)
                0393      &                       - aTv(1,i,bi,bj)*yTv(2,i,bi,bm)
                0394      &                       )*tmpVar
                0395             cTv(1,i,bi,bj) =   cTv(1,i,bi,bj)*tmpVar
                0396             aTv(1,i,bi,bj) = - aTv(1,i,bi,bj)*aTv(2,i,bi,bm)*tmpVar
                0397 
                0398 C-    update row [2,bj] <- [2,bj] - a(2,bj)*[2,bj-1] + a(2,bj)*c(2,bj-1)*[1,bj]
                0399             tmpVar = aTv(2,i,bi,bj)*cTv(2,i,bi,bm)
                0400             yTv(2,i,bi,bj) =  yTv(2,i,bi,bj)
                0401      &                      - aTv(2,i,bi,bj)*yTv(2,i,bi,bm)
                0402      &                      + tmpVar*yTv(1,i,bi,bj)
                0403             cTv(2,i,bi,bj) =  cTv(2,i,bi,bj)
                0404      &                      + tmpVar*cTv(1,i,bi,bj)
                0405             aTv(2,i,bi,bj) = -aTv(2,i,bi,bj)*aTv(2,i,bi,bm)
                0406      &                      + tmpVar*aTv(1,i,bi,bj)
                0407           ENDDO
                0408 
                0409           DO bj=nSy-1,1,-1
                0410            bp = bj+1
                0411            DO j=1,2
                0412 C-    update row [1,bj] <- [1,bj] - c(1,bj)*[1,bj+1]
                0413 C-    update row [2,bj] <- [2,bj] - c(2,bj)*[1,bj+1]
                0414             aTv(j,i,bi,bj) =  aTv(j,i,bi,bj)
                0415      &                      - cTv(j,i,bi,bj)*aTv(1,i,bi,bp)
                0416             yTv(j,i,bi,bj) =  yTv(j,i,bi,bj)
                0417      &                      - cTv(j,i,bi,bj)*yTv(1,i,bi,bp)
                0418             cTv(j,i,bi,bj) = -cTv(j,i,bi,bj)*cTv(1,i,bi,bp)
                0419            ENDDO
                0420           ENDDO
                0421 
                0422 C--  periodic in Y:  X_0 <=> X_Ny and X_(N+1) <=> X_1 ;
                0423 C    find the value at the 2 opposite location (j=1 and j=Ny)
                0424           bm = 1
                0425           bp = nSy
                0426           cTv(1,i,bi,bm) = oneRL + cTv(1,i,bi,bm)
                0427           aTv(2,i,bi,bp) = oneRL + aTv(2,i,bi,bp)
                0428           tmpVar = cTv(1,i,bi,bm) * aTv(2,i,bi,bp)
                0429      &           - aTv(1,i,bi,bm) * cTv(2,i,bi,bp)
                0430           IF ( tmpVar.NE.0. _d 0 ) THEN
                0431              tmpVar = 1. _d 0 / tmpVar
                0432           ELSE
                0433              tmpVar = 0. _d 0
                0434              errCode = 1
                0435           ENDIF
                0436           vTmp1 = ( aTv(2,i,bi,bp) * yTv(1,i,bi,bm)
                0437      &            - aTv(1,i,bi,bm) * yTv(2,i,bi,bp)
                0438      &            )*tmpVar
                0439           vTmp2 = ( cTv(1,i,bi,bm) * yTv(2,i,bi,bp)
                0440      &            - cTv(2,i,bi,bp) * yTv(1,i,bi,bm)
                0441      &            )*tmpVar
                0442 
                0443 C-    finalise tile-edges solution (put into RHS "yTv"):
                0444           DO bj=1,nSy
                0445            DO j=1,2
                0446             IF ( bj+j .EQ.2 ) THEN
                0447              yTv(j,i,bi,bj) = vTmp1
                0448             ELSEIF ( bj+j .EQ. nSy+2 ) THEN
                0449              yTv(j,i,bi,bj) = vTmp2
                0450             ELSE
                0451              yTv(j,i,bi,bj) = yTv(j,i,bi,bj)
                0452      &                      - aTv(j,i,bi,bj) * vTmp2
                0453      &                      - cTv(j,i,bi,bj) * vTmp1
                0454             ENDIF
                0455            ENDDO
                0456           ENDDO
                0457 
                0458          ENDDO
                0459         ENDDO
                0460         _END_MASTER(myThid)
                0461         _BARRIER
                0462 
                0463         DO bj=myByLo(myThid),myByHi(myThid)
                0464          DO bi=myBxLo(myThid),myBxHi(myThid)
                0465           bm = 1 + MOD(bj-2+nSy,nSy)
                0466           bp = 1 + MOD(bj-0+nSy,nSy)
                0467           DO j= 1,sNy
                0468            DO i= 1,sNx
                0469             vFld(i,j,k,bi,bj) = yy_V(i,j,bi,bj)
                0470      &                      - alpV(i,j,bi,bj) * yTv(2,i,bi,bm)
                0471      &                      - gamV(i,j,bi,bj) * yTv(1,i,bi,bp)
                0472            ENDDO
                0473           ENDDO
                0474          ENDDO
                0475         ENDDO
                0476 
                0477 C     end solve for vFld
                0478        ENDIF
                0479 
                0480 C     end k-loop
                0481       ENDDO
                0482 
                0483       RETURN
                0484       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!