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

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

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

revision 1.1 by cnh, Wed Apr 22 19:15:30 1998 UTC revision 1.11 by cnh, Thu May 28 03:34:52 1998 UTC
# Line 1  Line 1 
1  C $Id$  C $Header$
2    
3  #include "CPP_EEOPTIONS.h"  #include "CPP_EEOPTIONS.h"
4    
# Line 35  C     myNorm - Work variable used in clc Line 35  C     myNorm - Work variable used in clc
35        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
36        INTEGER bi, bj        INTEGER bi, bj
37        INTEGER I,  J, K        INTEGER I,  J, K
38        real    faceArea        _RL     faceArea
39        _RL     myNorm        _RL     myNorm
40          _RL     aC, aCw, aCs
41    
42  C--   Initialise laplace operator  C--   Initialise laplace operator
43  C     aW2d: integral in Z Ax/dX  C     aW2d: integral in Z Ax/dX
# Line 53  C     aS2d: integral in Z Ay/dY Line 54  C     aS2d: integral in Z Ay/dY
54          DO K=1,Nz          DO K=1,Nz
55           DO J=1,sNy           DO J=1,sNy
56            DO I=1,sNx            DO I=1,sNx
57             faceArea = dyG(I,J,bi,bj)*dzF(K)*HFacW(I,J,K,bi,bj)             faceArea = _dyG(I,J,bi,bj)*dzF(K)*_hFacW(I,J,K,bi,bj)
58             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +
59       &      gravity*faceArea*rDxC(I,J,bi,bj)       &      gBaro*faceArea*_rdxC(I,J,bi,bj)
60             faceArea = dxG(I,J,bi,bj)*dzF(K)*HFacS(I,J,K,bi,bj)             faceArea = _dxG(I,J,bi,bj)*dzF(K)*_hFacS(I,J,K,bi,bj)
61             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +
62       &      gravity*faceArea*rDyC(I,J,bi,bj)       &      gBaro*faceArea*_rdyC(I,J,bi,bj)
63            ENDDO            ENDDO
64           ENDDO           ENDDO
65          ENDDO          ENDDO
# Line 77  C     aS2d: integral in Z Ay/dY Line 78  C     aS2d: integral in Z Ay/dY
78        ELSE        ELSE
79         myNorm = 1. _d 0         myNorm = 1. _d 0
80        ENDIF        ENDIF
       _BEGIN_MASTER( myThid )  
81         cg2dNorm = myNorm         cg2dNorm = myNorm
82          _BEGIN_MASTER( myThid )
83  CcnhDebugStarts  CcnhDebugStarts
84         WRITE(msgBuf,*) ' CG2D normalisation factor = ', cg2dNorm         WRITE(msgBuf,'(A,F)') '// CG2D normalisation factor = ', cg2dNorm
85           CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
86           WRITE(msgBuf,*) '                               '
87         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
88  CcnhDebugEnds  CcnhDebugEnds
89        _END_MASTER( myThid )        _END_MASTER( myThid )
# Line 108  C     CALL PLOT_FIELD_XYR8( aS2d, 'AS2D Line 111  C     CALL PLOT_FIELD_XYR8( aS2d, 'AS2D
111  CcnhDebugEnds  CcnhDebugEnds
112    
113  C--   Initialise preconditioner  C--   Initialise preconditioner
114    C     Note. 20th May 1998
115    C           I made a weird discovery! In the model paper we argue
116    C           for the form of the preconditioner used here ( see
117    C           A Finite-volume, Incompressible Navier-Stokes Model
118    C           ...., Marshall et. al ). The algebra gives a simple
119    C           0.5 factor for the averaging of ac and aCw to get a
120    C           symmettric pre-conditioner. By using a factor of 0.51
121    C           i.e. scaling the off-diagonal terms in the
122    C           preconditioner down slightly I managed to get the
123    C           number of iterations for convergence in a test case to
124    C           drop form 192 -> 134! Need to investigate this further!
125    C           For now I have introduced a parameter cg2dpcOffDFac which
126    C           defaults to 0.51 but can be set at runtime.
127        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
128         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
129          DO J=1,sNy          DO J=1,sNy
130           DO I=1,sNx           DO I=1,sNx
131            pC(I,J,bi,bj) = 1. _d 0            pC(I,J,bi,bj) = 1. _d 0
132            IF (            aC = -(
133       &     aW2d(I,J,bi,bj) + aW2d(I+1,J,bi,bj)       &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)
134       &    +aS2d(I,J,bi,bj) + aS2D(I,J+1,bi,bj)       &    +aS2d(I,J,bi,bj) + aS2D(I  ,J+1,bi,bj)
135       &     .EQ. 0.       &    +freeSurfFac*myNorm*
136       &       )  pC(I,J,bi,bj) = 0. _d 0       &     zA(I,J,bi,bj)/deltaTMom/deltaTMom
137            pW(I,J,bi,bj) = 0.       &    )
138            pS(I,J,bi,bj) = 0.            aCs = -(
139         &     aW2d(I,J-1,bi,bj) + aW2d(I+1,J-1,bi,bj)
140         &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)
141         &    +freeSurfFac*myNorm*
142         &     zA(I,J-1,bi,bj)/deltaTMom/deltaTMom
143         &    )
144              aCw = -(
145         &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)
146         &    +aS2d(I-1,J,bi,bj) + aS2d(I-1,J+1,bi,bj)
147         &    +freeSurfFac*myNorm*
148         &     zA(I-1,J,bi,bj)/deltaTMom/deltaTMom
149         &    )
150              IF ( aC .EQ. 0. ) THEN
151                pC(I,J,bi,bj) = 0. _d 0
152              ELSE
153               pC(I,J,bi,bj) =  1. _d 0 / aC
154              ENDIF
155              IF ( aC + aCw .EQ. 0. ) THEN
156               pW(I,J,bi,bj) = 0.
157              ELSE
158               pW(I,J,bi,bj) =
159         &     -aW2d(I  ,J  ,bi,bj)/((cg2dpcOffDFac *(aCw+aC))**2 )
160              ENDIF
161              IF ( aC + aCs .EQ. 0. ) THEN
162               pS(I,J,bi,bj) = 0.
163              ELSE
164               pS(I,J,bi,bj) =
165         &     -aS2d(I  ,J  ,bi,bj)/((cg2dpcOffDFac *(aCs+aC))**2 )
166              ENDIF
167    C         pC(I,J,bi,bj) = 1.
168    C         pW(I,J,bi,bj) = 0.
169    C         pS(I,J,bi,bj) = 0.
170           ENDDO           ENDDO
171          ENDDO          ENDDO
172         ENDDO         ENDDO
# Line 129  C--   Update overlap regions Line 176  C--   Update overlap regions
176        _EXCH_XY_R4(pW, myThid)        _EXCH_XY_R4(pW, myThid)
177        _EXCH_XY_R4(pS, myThid)        _EXCH_XY_R4(pS, myThid)
178    
179  C--   Set default values for initial guess  C--   Set default values for initial guess and RHS
180        DO bj=myByLo(myThid),myByHi(myThid)        IF ( startTime .EQ. 0 ) THEN
181         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
182          DO J=1,sNy          DO bi=myBxLo(myThid),myBxHi(myThid)
183           DO I=1,sNx           DO J=1,sNy
184            cg2d_x(I,J,bi,bj) = 0. _d 0            DO I=1,sNx
185               cg2d_x(I,J,bi,bj) = 0. _d 0
186               cg2d_b(I,J,bi,bj) = 0. _d 0
187              ENDDO
188           ENDDO           ENDDO
189          ENDDO          ENDDO
190         ENDDO         ENDDO
191        ENDDO  C--    Update overlap regions
192  C--   Update overlap regions         _EXCH_XY_R8(cg2d_x, myThid)
193        _EXCH_XY_R8(cg2d_x, myThid)         _EXCH_XY_R8(cg2d_b, myThid)
194          ENDIF
195    
196        RETURN        RETURN
197        END        END
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.22