/[MITgcm]/MITgcm/model/src/update_cg2d.F
ViewVC logotype

Diff of /MITgcm/model/src/update_cg2d.F

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

revision 1.2 by cnh, Wed Sep 26 18:09:16 2001 UTC revision 1.3 by jmc, Sun Feb 10 20:04:51 2002 UTC
# Line 51  C     faceArea - Temporary used to hold Line 51  C     faceArea - Temporary used to hold
51        INTEGER bi, bj        INTEGER bi, bj
52        INTEGER I, J, K        INTEGER I, J, K
53        _RL     faceArea        _RL     faceArea
54        _RL     aC, aCw, aCs        _RL     pW_tmp, pS_tmp
55  CEOP  CEOP
56    
57  C--   Initialise laplace operator  C--   Initialise laplace operator
# Line 59  C     aW2d: integral in Z Ax/dX Line 59  C     aW2d: integral in Z Ax/dX
59  C     aS2d: integral in Z Ay/dY  C     aS2d: integral in Z Ay/dY
60        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
61         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
62          DO J=1,sNy          DO J=1,sNy+1
63           DO I=1,sNx           DO I=1,sNx+1
64            aW2d(I,J,bi,bj) = 0. _d 0            aW2d(I,J,bi,bj) = 0. _d 0
65            aS2d(I,J,bi,bj) = 0. _d 0            aS2d(I,J,bi,bj) = 0. _d 0
66           ENDDO           ENDDO
67          ENDDO          ENDDO
68          DO K=1,Nr          DO K=1,Nr
69           DO J=1,sNy           DO J=1,sNy+1
70            DO I=1,sNx            DO I=1,sNx+1
71             faceArea = _dyG(I,J,bi,bj)*drF(K)             faceArea = _dyG(I,J,bi,bj)*drF(K)
72       &               *_hFacW(I,J,K,bi,bj)       &               *_hFacW(I,J,K,bi,bj)
73             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +
# Line 95  C     aS2d: integral in Z Ay/dY Line 95  C     aS2d: integral in Z Ay/dY
95           ENDDO           ENDDO
96          ENDIF          ENDIF
97  #endif  #endif
98          DO J=1,sNy          DO J=1,sNy+1
99           DO I=1,sNx           DO I=1,sNx+1
100            aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj)*cg2dNorm            aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj)*cg2dNorm
101       &                     *implicSurfPress*implicDiv2DFlow       &                     *implicSurfPress*implicDiv2DFlow
102            aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj)*cg2dNorm            aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj)*cg2dNorm
103       &                     *implicSurfPress*implicDiv2DFlow       &                     *implicSurfPress*implicDiv2DFlow
104           ENDDO           ENDDO
105          ENDDO          ENDDO
106    C--   Start to compute preconditioner matrix (use cg2d_q as temporary array)
107            DO J=1,sNy
108             DO I=1,sNx
109              cg2d_q(I,J,bi,bj) = -(
110         &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)
111         &    +aS2d(I,J,bi,bj) + aS2d(I  ,J+1,bi,bj)
112         &    +freeSurfFac*cg2dNorm*recip_Bo(I,J,bi,bj)*
113         &     rA(I,J,bi,bj)/deltaTMom/deltaTMom
114         &    )
115             ENDDO
116            ENDDO
117         ENDDO         ENDDO
118        ENDDO        ENDDO
119                
120  C--   Update overlap regions  C--   Update overlap regions
121  c     _EXCH_XY_R4(aW2d, myThid)        _EXCH_XY_R8(cg2d_q, myThid)
 c     _EXCH_XY_R4(aS2d, myThid)  
       CALL EXCH_UV_XY_RS(aW2d,aS2d,.FALSE.,myThid)  
122    
123  C--   Initialise preconditioner  C--   Initialise preconditioner
124        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
125         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
126          DO J=1,sNy          DO J=1,sNy+1
127           DO I=1,sNx           DO I=1,sNx+1
128            pC(I,J,bi,bj) = 1. _d 0            IF ( cg2d_q(I,J,bi,bj) .EQ. 0. ) THEN
           aC = -(  
      &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)  
      &    +aS2d(I,J,bi,bj) + aS2d(I  ,J+1,bi,bj)  
      &    +freeSurfFac*cg2dNorm*recip_Bo(I,J,bi,bj)*  
      &     rA(I,J,bi,bj)/deltaTMom/deltaTMom  
      &    )  
           aCs = -(  
      &     aW2d(I,J-1,bi,bj) + aW2d(I+1,J-1,bi,bj)  
      &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)  
      &    +freeSurfFac*cg2dNorm*recip_Bo(I,J-1,bi,bj)*  
      &     rA(I,J-1,bi,bj)/deltaTMom/deltaTMom  
      &    )  
           aCw = -(  
      &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)  
      &    +aS2d(I-1,J,bi,bj) + aS2d(I-1,J+1,bi,bj)  
      &    +freeSurfFac*cg2dNorm*recip_Bo(I-1,J,bi,bj)*  
      &     rA(I-1,J,bi,bj)/deltaTMom/deltaTMom  
      &    )  
           IF ( aC .EQ. 0. ) THEN  
129              pC(I,J,bi,bj) = 1. _d 0              pC(I,J,bi,bj) = 1. _d 0
130            ELSE            ELSE
131             pC(I,J,bi,bj) =  1. _d 0 / aC             pC(I,J,bi,bj) =  1. _d 0 / cg2d_q(I,J,bi,bj)
132            ENDIF            ENDIF
133            IF ( aC + aCw .EQ. 0. ) THEN            pW_tmp = cg2d_q(I,J,bi,bj)+cg2d_q(I-1,J,bi,bj)
134              IF ( pW_tmp .EQ. 0. ) THEN
135             pW(I,J,bi,bj) = 0.             pW(I,J,bi,bj) = 0.
136            ELSE            ELSE
137             pW(I,J,bi,bj) =             pW(I,J,bi,bj) =
138       &     -aW2d(I  ,J  ,bi,bj)/((cg2dpcOffDFac *(aCw+aC))**2 )       &     -aW2d(I,J,bi,bj)/((cg2dpcOffDFac *pW_tmp)**2 )
139            ENDIF            ENDIF
140            IF ( aC + aCs .EQ. 0. ) THEN            pS_tmp = cg2d_q(I,J,bi,bj)+cg2d_q(I,J-1,bi,bj)
141              IF ( pS_tmp .EQ. 0. ) THEN
142             pS(I,J,bi,bj) = 0.             pS(I,J,bi,bj) = 0.
143            ELSE            ELSE
144             pS(I,J,bi,bj) =             pS(I,J,bi,bj) =
145       &     -aS2d(I  ,J  ,bi,bj)/((cg2dpcOffDFac *(aCs+aC))**2 )       &     -aS2d(I,J,bi,bj)/((cg2dpcOffDFac *pS_tmp)**2 )
146            ENDIF            ENDIF
147           ENDDO           ENDDO
148          ENDDO          ENDDO
149         ENDDO         ENDDO
150        ENDDO        ENDDO
 C--   Update overlap regions  
       _EXCH_XY_R4(pC, myThid)  
 c     _EXCH_XY_R4(pW, myThid)  
 c     _EXCH_XY_R4(pS, myThid)  
       CALL EXCH_UV_XY_RS(pW,pS,.FALSE.,myThid)  
151    
152  #endif /* NONLIN_FRSURF */  #endif /* NONLIN_FRSURF */
153    
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22