/[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.3 by cnh, Sun Apr 26 23:41:54 1998 UTC revision 1.19 by cnh, Sun Sep 6 17:35:20 1998 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    
3  #include "CPP_EEOPTIONS.h"  #include "CPP_OPTIONS.h"
4    
5  CStartOfInterface  CStartOfInterface
6        SUBROUTINE INI_CG2D( myThid )        SUBROUTINE INI_CG2D( myThid )
# Line 17  C     === Global variables === Line 17  C     === Global variables ===
17  #include "EEPARAMS.h"  #include "EEPARAMS.h"
18  #include "PARAMS.h"  #include "PARAMS.h"
19  #include "GRID.h"  #include "GRID.h"
20    #include "DYNVARS.h"
21  #include "CG2D.h"  #include "CG2D.h"
22    
23  C     === Routine arguments ===  C     === Routine arguments ===
# Line 35  C     myNorm - Work variable used in clc Line 36  C     myNorm - Work variable used in clc
36        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
37        INTEGER bi, bj        INTEGER bi, bj
38        INTEGER I,  J, K        INTEGER I,  J, K
39        real    faceArea        _RL     faceArea
40        _RL     myNorm        _RS     myNorm
41          _RL     aC, aCw, aCs
42    
43  C--   Initialise laplace operator  C--   Initialise laplace operator
44  C     aW2d: integral in Z Ax/dX  C     aW2d: integral in Z Ax/dX
# Line 50  C     aS2d: integral in Z Ay/dY Line 52  C     aS2d: integral in Z Ay/dY
52            aS2d(I,J,bi,bj) = 0. _d 0            aS2d(I,J,bi,bj) = 0. _d 0
53           ENDDO           ENDDO
54          ENDDO          ENDDO
55          DO K=1,Nz          DO K=1,Nr
56           DO J=1,sNy           DO J=1,sNy
57            DO I=1,sNx            DO I=1,sNx
58             faceArea = dyG(I,J,bi,bj)*dzF(K)*HFacW(I,J,K,bi,bj)             faceArea = _dyG(I,J,bi,bj)*drF(K)*_hFacW(I,J,K,bi,bj)
59             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +
60       &      gravity*faceArea*rDxC(I,J,bi,bj)       &      gBaro*faceArea*recip_dxC(I,J,bi,bj)
61             faceArea = dxG(I,J,bi,bj)*dzF(K)*HFacS(I,J,K,bi,bj)             faceArea = _dxG(I,J,bi,bj)*drF(K)*_hFacS(I,J,K,bi,bj)
62             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +
63       &      gravity*faceArea*rDyC(I,J,bi,bj)       &      gBaro*faceArea*recip_dyC(I,J,bi,bj)
64            ENDDO            ENDDO
65           ENDDO           ENDDO
66          ENDDO          ENDDO
# Line 71  C     aS2d: integral in Z Ay/dY Line 73  C     aS2d: integral in Z Ay/dY
73         ENDDO         ENDDO
74        ENDDO        ENDDO
75        cg2dNbuf(1,myThid) = myNorm        cg2dNbuf(1,myThid) = myNorm
76        _GLOBAL_MAX_R8( cg2dNbuf, myNorm, myThid )        _GLOBAL_MAX_R4( cg2dNbuf, myNorm, myThid )
77        IF ( cg2dNbuf(1,1) .NE. 0. _d 0 ) THEN        IF ( cg2dNbuf(1,1) .NE. 0. _d 0 ) THEN
78         myNorm = 1. _d 0/cg2dNbuf(1,1)         myNorm = 1. _d 0/cg2dNbuf(1,1)
79        ELSE        ELSE
80         myNorm = 1. _d 0         myNorm = 1. _d 0
81        ENDIF        ENDIF
       _BEGIN_MASTER( myThid )  
82         cg2dNorm = myNorm         cg2dNorm = myNorm
83          _BEGIN_MASTER( myThid )
84  CcnhDebugStarts  CcnhDebugStarts
85         WRITE(msgBuf,'(A,F)') '// CG2D normalisation factor = ', cg2dNorm         WRITE(msgBuf,'(A,E40.25)') '// CG2D normalisation factor = ', cg2dNorm
86         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
87         WRITE(msgBuf,*) '                               '         WRITE(msgBuf,*) '                               '
88         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
# Line 99  CcnhDebugEnds Line 101  CcnhDebugEnds
101                
102  C--   Update overlap regions  C--   Update overlap regions
103  CcnhDebugStarts  CcnhDebugStarts
104  C     CALL PLOT_FIELD_XYR8( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )
105  C     CALL PLOT_FIELD_XYR8( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )
106  CcnhDebugEnds  CcnhDebugEnds
107        _EXCH_XY_R4(aW2d, myThid)        _EXCH_XY_R4(aW2d, myThid)
108        _EXCH_XY_R4(aS2d, myThid)        _EXCH_XY_R4(aS2d, myThid)
109  CcnhDebugStarts  CcnhDebugStarts
110  C     CALL PLOT_FIELD_XYR8( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )        CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )
111  C     CALL PLOT_FIELD_XYR8( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )        CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )
112  CcnhDebugEnds  CcnhDebugEnds
113    
114  C--   Initialise preconditioner  C--   Initialise preconditioner
115    C     Note. 20th May 1998
116    C           I made a weird discovery! In the model paper we argue
117    C           for the form of the preconditioner used here ( see
118    C           A Finite-volume, Incompressible Navier-Stokes Model
119    C           ...., Marshall et. al ). The algebra gives a simple
120    C           0.5 factor for the averaging of ac and aCw to get a
121    C           symmettric pre-conditioner. By using a factor of 0.51
122    C           i.e. scaling the off-diagonal terms in the
123    C           preconditioner down slightly I managed to get the
124    C           number of iterations for convergence in a test case to
125    C           drop form 192 -> 134! Need to investigate this further!
126    C           For now I have introduced a parameter cg2dpcOffDFac which
127    C           defaults to 0.51 but can be set at runtime.
128        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
129         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
130          DO J=1,sNy          DO J=1,sNy
131           DO I=1,sNx           DO I=1,sNx
132            pC(I,J,bi,bj) = 1. _d 0            pC(I,J,bi,bj) = 1. _d 0
133            IF (            aC = -(
134       &     aW2d(I,J,bi,bj) + aW2d(I+1,J,bi,bj)       &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)
135       &    +aS2d(I,J,bi,bj) + aS2D(I,J+1,bi,bj)       &    +aS2d(I,J,bi,bj) + aS2D(I  ,J+1,bi,bj)
136       &     .EQ. 0.       &    +freeSurfFac*myNorm* horiVertRatio*
137       &       )  pC(I,J,bi,bj) = 0. _d 0       &     rA(I,J,bi,bj)/deltaTMom/deltaTMom
138            pW(I,J,bi,bj) = 0.       &    )
139            pS(I,J,bi,bj) = 0.            aCs = -(
140         &     aW2d(I,J-1,bi,bj) + aW2d(I+1,J-1,bi,bj)
141         &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)
142         &    +freeSurfFac*myNorm* horiVertRatio*
143         &     rA(I,J-1,bi,bj)/deltaTMom/deltaTMom
144         &    )
145              aCw = -(
146         &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)
147         &    +aS2d(I-1,J,bi,bj) + aS2d(I-1,J+1,bi,bj)
148         &    +freeSurfFac*myNorm* horiVertRatio*
149         &     rA(I-1,J,bi,bj)/deltaTMom/deltaTMom
150         &    )
151              IF ( aC .EQ. 0. ) THEN
152                pC(I,J,bi,bj) = 0. _d 0
153              ELSE
154               pC(I,J,bi,bj) =  1. _d 0 / aC
155              ENDIF
156              IF ( aC + aCw .EQ. 0. ) THEN
157               pW(I,J,bi,bj) = 0.
158              ELSE
159               pW(I,J,bi,bj) =
160         &     -aW2d(I  ,J  ,bi,bj)/((cg2dpcOffDFac *(aCw+aC))**2 )
161              ENDIF
162              IF ( aC + aCs .EQ. 0. ) THEN
163               pS(I,J,bi,bj) = 0.
164              ELSE
165               pS(I,J,bi,bj) =
166         &     -aS2d(I  ,J  ,bi,bj)/((cg2dpcOffDFac *(aCs+aC))**2 )
167              ENDIF
168    C         pC(I,J,bi,bj) = 1.
169    C         pW(I,J,bi,bj) = 0.
170    C         pS(I,J,bi,bj) = 0.
171           ENDDO           ENDDO
172          ENDDO          ENDDO
173         ENDDO         ENDDO
# Line 130  C--   Update overlap regions Line 176  C--   Update overlap regions
176        _EXCH_XY_R4(pC, myThid)        _EXCH_XY_R4(pC, myThid)
177        _EXCH_XY_R4(pW, myThid)        _EXCH_XY_R4(pW, myThid)
178        _EXCH_XY_R4(pS, myThid)        _EXCH_XY_R4(pS, myThid)
179    CcnhDebugStarts
180          CALL PLOT_FIELD_XYRS( pC, 'pC   INI_CG2D.2' , 1, myThid )
181          CALL PLOT_FIELD_XYRS( pW, 'pW   INI_CG2D.2' , 1, myThid )
182          CALL PLOT_FIELD_XYRS( pS, 'pS   INI_CG2D.2' , 1, myThid )
183    CcnhDebugEnds
184    
185  C--   Set default values for initial guess  C--   Set default values for initial guess and RHS
186        DO bj=myByLo(myThid),myByHi(myThid)        IF ( startTime .EQ. 0 ) THEN
187         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
188          DO J=1,sNy          DO bi=myBxLo(myThid),myBxHi(myThid)
189           DO I=1,sNx           DO J=1,sNy
190            cg2d_x(I,J,bi,bj) = 0. _d 0            DO I=1,sNx
191               cg2d_x(I,J,bi,bj) = 0. _d 0
192               cg2d_b(I,J,bi,bj) = 0. _d 0
193    #ifdef ALLOW_CD
194               cg2d_xNM1(I,J,bi,bj) = 0. _d 0
195    #endif
196              ENDDO
197           ENDDO           ENDDO
198          ENDDO          ENDDO
199         ENDDO         ENDDO
200        ENDDO  C--    Update overlap regions
201  C--   Update overlap regions         _EXCH_XY_R8(cg2d_x, myThid)
202        _EXCH_XY_R8(cg2d_x, myThid)         _EXCH_XY_R8(cg2d_b, myThid)
203    #ifdef ALLOW_CD
204           _EXCH_XY_R8(cg2d_xNM1, myThid)
205    #endif
206          ENDIF
207    
208        RETURN        RETURN
209        END        END
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.19
  ViewVC Help
Powered by ViewVC 1.1.22