/[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.7 by jmc, Tue Dec 5 05:25:08 2006 UTC revision 1.8 by jmc, Wed May 18 23:37:38 2011 UTC
# Line 29  C     === Global variables === Line 29  C     === Global variables ===
29  #include "GRID.h"  #include "GRID.h"
30  #include "SURFACE.h"  #include "SURFACE.h"
31  #include "CG2D.h"  #include "CG2D.h"
 #ifdef ALLOW_OBCS  
 #include "OBCS.h"  
 #endif  
32    
33  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
34  C     === Routine arguments ===  C     === Routine arguments ===
35  C     myTime - Current time in simulation  C     myTime :: Current time in simulation
36  C     myIter - Current iteration number in simulation  C     myIter :: Current iteration number in simulation
37  C     myThid - Thread number for this instance of the routine.  C     myThid :: Thread number for this instance of the routine.
38        _RL myTime        _RL myTime
39        INTEGER myIter        INTEGER myIter
40        INTEGER myThid        INTEGER myThid
# Line 48  C-- Note : compared to "INI_CG2D", no ne Line 45  C-- Note : compared to "INI_CG2D", no ne
45  C   the solver norn=malisation factor of the solver tolerance  C   the solver norn=malisation factor of the solver tolerance
46  C     === Local variables ===  C     === Local variables ===
47  C     bi,bj  :: tile indices  C     bi,bj  :: tile indices
48  C     I,J,K  :: Loop counters  C     i,j,k  :: Loop counters
49  C     faceArea :: Temporary used to hold cell face areas.  C     faceArea :: Temporary used to hold cell face areas.
50        INTEGER bi, bj        INTEGER bi, bj
51        INTEGER I, J, K, ks        INTEGER i, j, k, ks
52        _RL     faceArea        _RL     faceArea
53        _RL     pW_tmp, pS_tmp        _RL     pW_tmp, pS_tmp
54        LOGICAL updatePreCond        LOGICAL updatePreCond
# Line 70  C     aW2d: integral in Z Ax/dX Line 67  C     aW2d: integral in Z Ax/dX
67  C     aS2d: integral in Z Ay/dY  C     aS2d: integral in Z Ay/dY
68        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
69         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
70          DO J=1,sNy+1          DO j=1,sNy+1
71           DO I=1,sNx+1           DO i=1,sNx+1
72            aW2d(I,J,bi,bj) = 0. _d 0            aW2d(i,j,bi,bj) = 0. _d 0
73            aS2d(I,J,bi,bj) = 0. _d 0            aS2d(i,j,bi,bj) = 0. _d 0
74           ENDDO           ENDDO
75          ENDDO          ENDDO
76          DO K=1,Nr          DO k=1,Nr
77           DO J=1,sNy+1           DO j=1,sNy+1
78            DO I=1,sNx+1            DO i=1,sNx+1
79  C  deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect  C  deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect
80             faceArea = _dyG(I,J,bi,bj)*drF(K)             faceArea = _dyG(i,j,bi,bj)*drF(k)
81       &               *_hFacW(I,J,K,bi,bj)       &               *_hFacW(i,j,k,bi,bj)
82             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj)             aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
83       &              + faceArea*recip_dxC(I,J,bi,bj)       &              + faceArea*recip_dxC(i,j,bi,bj)
84             faceArea = _dxG(I,J,bi,bj)*drF(K)             faceArea = _dxG(i,j,bi,bj)*drF(k)
85       &               *_hFacS(I,J,K,bi,bj)       &               *_hFacS(i,j,k,bi,bj)
86             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj)             aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
87       &              + faceArea*recip_dyC(I,J,bi,bj)       &              + faceArea*recip_dyC(i,j,bi,bj)
88            ENDDO            ENDDO
89           ENDDO           ENDDO
90          ENDDO          ENDDO
91            DO j=1,sNy+1
92             DO i=1,sNx+1
93              aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)*cg2dNorm
94         &                     *implicSurfPress*implicDiv2DFlow
95  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
96          IF (useOBCS) THEN       &                  *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj)
 C  Note: would need loop from 1 to sNx+1 (and below, from 1 to sNy+1)  
 C    to get the same solver-matrix as from INI_CG2D,  
 C    since aS2d & aW2d are not exchanged here (but in INI_CG2D, they are).  
          DO I=1,sNx  
           IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj),bi,bj)=0.  
           IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj)+1,bi,bj)=0.  
           IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj)+1,bi,bj)=0.  
           IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj),bi,bj)=0.  
          ENDDO  
          DO J=1,sNy  
           IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj),J,bi,bj)=0.  
           IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj)+1,J,bi,bj)=0.  
           IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj)+1,J,bi,bj)=0.  
           IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj),J,bi,bj)=0.  
          ENDDO  
         ENDIF  
97  #endif  #endif
98          DO J=1,sNy+1            aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)*cg2dNorm
          DO I=1,sNx+1  
           aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj)*cg2dNorm  
      &                     *implicSurfPress*implicDiv2DFlow  
           aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj)*cg2dNorm  
99       &                     *implicSurfPress*implicDiv2DFlow       &                     *implicSurfPress*implicDiv2DFlow
100    #ifdef ALLOW_OBCS
101         &                  *maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj)
102    #endif
103           ENDDO           ENDDO
104          ENDDO          ENDDO
105  C--   compute matrix main diagonal :  C--   compute matrix main diagonal :
106          IF ( deepAtmosphere ) THEN          IF ( deepAtmosphere ) THEN
107           DO J=1,sNy           DO j=1,sNy
108            DO I=1,sNx            DO i=1,sNx
109             ks = ksurfC(I,J,bi,bj)             ks = ksurfC(i,j,bi,bj)
110             aC2d(I,J,bi,bj) = -(             aC2d(i,j,bi,bj) = -(
111       &       aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)       &       aW2d(i,j,bi,bj) + aW2d(i+1, j ,bi,bj)
112       &      +aS2d(I,J,bi,bj) + aS2d(I  ,J+1,bi,bj)       &      +aS2d(i,j,bi,bj) + aS2d( i ,j+1,bi,bj)
113       &      +freeSurfFac*cg2dNorm*recip_Bo(I,J,bi,bj)*deepFac2F(ks)       &      +freeSurfFac*cg2dNorm*recip_Bo(i,j,bi,bj)*deepFac2F(ks)
114       &                  *rA(I,J,bi,bj)/deltaTMom/deltaTfreesurf       &                  *rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
115       &                        )       &                        )
116            ENDDO            ENDDO
117           ENDDO           ENDDO
118          ELSE          ELSE
119           DO J=1,sNy           DO j=1,sNy
120            DO I=1,sNx            DO i=1,sNx
121             aC2d(I,J,bi,bj) = -(             aC2d(i,j,bi,bj) = -(
122       &       aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)       &       aW2d(i,j,bi,bj) + aW2d(i+1, j ,bi,bj)
123       &      +aS2d(I,J,bi,bj) + aS2d(I  ,J+1,bi,bj)       &      +aS2d(i,j,bi,bj) + aS2d( i ,j+1,bi,bj)
124       &      +freeSurfFac*cg2dNorm*recip_Bo(I,J,bi,bj)       &      +freeSurfFac*cg2dNorm*recip_Bo(i,j,bi,bj)
125       &                  *rA(I,J,bi,bj)/deltaTMom/deltaTfreesurf       &                  *rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
126       &                        )       &                        )
127            ENDDO            ENDDO
128           ENDDO           ENDDO
# Line 154  C--   Update overlap regions Line 138  C--   Update overlap regions
138  C--   Initialise preconditioner  C--   Initialise preconditioner
139        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
140         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
141          DO J=1,sNy+1          DO j=1,sNy+1
142           DO I=1,sNx+1           DO i=1,sNx+1
143            IF ( aC2d(I,J,bi,bj) .EQ. 0. ) THEN            IF ( aC2d(i,j,bi,bj) .EQ. 0. ) THEN
144              pC(I,J,bi,bj) = 1. _d 0              pC(i,j,bi,bj) = 1. _d 0
145            ELSE            ELSE
146             pC(I,J,bi,bj) =  1. _d 0 / aC2d(I,J,bi,bj)             pC(i,j,bi,bj) =  1. _d 0 / aC2d(i,j,bi,bj)
147            ENDIF            ENDIF
148            pW_tmp = aC2d(I,J,bi,bj)+aC2d(I-1,J,bi,bj)            pW_tmp = aC2d(i,j,bi,bj)+aC2d(i-1,j,bi,bj)
149            IF ( pW_tmp .EQ. 0. ) THEN            IF ( pW_tmp .EQ. 0. ) THEN
150             pW(I,J,bi,bj) = 0.             pW(i,j,bi,bj) = 0.
151            ELSE            ELSE
152             pW(I,J,bi,bj) =             pW(i,j,bi,bj) =
153       &     -aW2d(I,J,bi,bj)/((cg2dpcOffDFac *pW_tmp)**2 )       &     -aW2d(i,j,bi,bj)/((cg2dpcOffDFac *pW_tmp)**2 )
154            ENDIF            ENDIF
155            pS_tmp = aC2d(I,J,bi,bj)+aC2d(I,J-1,bi,bj)            pS_tmp = aC2d(i,j,bi,bj)+aC2d(i,j-1,bi,bj)
156            IF ( pS_tmp .EQ. 0. ) THEN            IF ( pS_tmp .EQ. 0. ) THEN
157             pS(I,J,bi,bj) = 0.             pS(i,j,bi,bj) = 0.
158            ELSE            ELSE
159             pS(I,J,bi,bj) =             pS(i,j,bi,bj) =
160       &     -aS2d(I,J,bi,bj)/((cg2dpcOffDFac *pS_tmp)**2 )       &     -aS2d(i,j,bi,bj)/((cg2dpcOffDFac *pS_tmp)**2 )
161            ENDIF            ENDIF
162           ENDDO           ENDDO
163          ENDDO          ENDDO
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22