/[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.30 by jmc, Tue Feb 20 15:06:21 2001 UTC
# Line 1  Line 1 
1  C $Id$  C $Header$
2    C $Name$
3    
4  #include "CPP_EEOPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6  CStartOfInterface  CStartOfInterface
7        SUBROUTINE INI_CG2D( myThid )        SUBROUTINE INI_CG2D( myThid )
# Line 11  C     |================================= Line 12  C     |=================================
12  C     | These arrays are purely a function of the basin geom.    |  C     | These arrays are purely a function of the basin geom.    |
13  C     | We set then here once and them use then repeatedly.      |  C     | We set then here once and them use then repeatedly.      |
14  C     \==========================================================/  C     \==========================================================/
15          IMPLICIT NONE
16    
17  C     === Global variables ===  C     === Global variables ===
18  #include "SIZE.h"  #include "SIZE.h"
19  #include "EEPARAMS.h"  #include "EEPARAMS.h"
20  #include "PARAMS.h"  #include "PARAMS.h"
21  #include "GRID.h"  #include "GRID.h"
22    #include "DYNVARS.h"
23  #include "CG2D.h"  #include "CG2D.h"
24    #ifdef ALLOW_OBCS
25    #include "OBCS.h"
26    #endif
27    
28  C     === Routine arguments ===  C     === Routine arguments ===
29  C     myThid - Thread no. that called this routine.  C     myThid - Thread no. that called this routine.
# Line 35  C     myNorm - Work variable used in clc Line 41  C     myNorm - Work variable used in clc
41        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
42        INTEGER bi, bj        INTEGER bi, bj
43        INTEGER I,  J, K        INTEGER I,  J, K
44        real    faceArea        _RL     faceArea
45        _RL     myNorm        _RS     myNorm
46          _RL     aC, aCw, aCs
47    
48  C--   Initialise laplace operator  C--   Initialise laplace operator
49  C     aW2d: integral in Z Ax/dX  C     aW2d: integral in Z Ax/dX
# Line 50  C     aS2d: integral in Z Ay/dY Line 57  C     aS2d: integral in Z Ay/dY
57            aS2d(I,J,bi,bj) = 0. _d 0            aS2d(I,J,bi,bj) = 0. _d 0
58           ENDDO           ENDDO
59          ENDDO          ENDDO
60          DO K=1,Nz          DO K=1,Nr
61           DO J=1,sNy           DO J=1,sNy
62            DO I=1,sNx            DO I=1,sNx
63             faceArea = dyG(I,J,bi,bj)*dzF(K)*HFacW(I,J,K,bi,bj)             faceArea = _dyG(I,J,bi,bj)*drF(K)
64         &               *_hFacW(I,J,K,bi,bj)
65             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +
66       &      gravity*faceArea*rDxC(I,J,bi,bj)       &      implicSurfPress*implicDiv2DFlow*
67             faceArea = dxG(I,J,bi,bj)*dzF(K)*HFacS(I,J,K,bi,bj)       &      gBaro*faceArea*recip_dxC(I,J,bi,bj)
68               faceArea = _dxG(I,J,bi,bj)*drF(K)
69         &               *_hFacS(I,J,K,bi,bj)
70             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +
71       &      gravity*faceArea*rDyC(I,J,bi,bj)       &      implicSurfPress*implicDiv2DFlow*
72         &      gBaro*faceArea*recip_dyC(I,J,bi,bj)
73            ENDDO            ENDDO
74           ENDDO           ENDDO
75          ENDDO          ENDDO
76    #ifdef ALLOW_OBCS
77            IF (useOBCS) THEN
78             DO I=1,sNx
79              IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj),bi,bj)=0.
80              IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj)+1,bi,bj)=0.
81              IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj)+1,bi,bj)=0.
82              IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj),bi,bj)=0.
83             ENDDO
84             DO J=1,sNy
85              IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj),J,bi,bj)=0.
86              IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj)+1,J,bi,bj)=0.
87              IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj)+1,J,bi,bj)=0.
88              IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj),J,bi,bj)=0.
89             ENDDO
90            ENDIF
91    #endif
92          DO J=1,sNy          DO J=1,sNy
93           DO I=1,sNx           DO I=1,sNx
94            myNorm = MAX(ABS(aW2d(I,J,bi,bj)),myNorm)            myNorm = MAX(ABS(aW2d(I,J,bi,bj)),myNorm)
# Line 70  C     aS2d: integral in Z Ay/dY Line 97  C     aS2d: integral in Z Ay/dY
97          ENDDO          ENDDO
98         ENDDO         ENDDO
99        ENDDO        ENDDO
100        cg2dNbuf(1,myThid) = myNorm        _GLOBAL_MAX_R4( myNorm, myThid )
101        _GLOBAL_MAX_R8( cg2dNbuf, myNorm, myThid )        IF ( myNorm .NE. 0. _d 0 ) THEN
102        IF ( cg2dNbuf(1,1) .NE. 0. _d 0 ) THEN         myNorm = 1. _d 0/myNorm
        myNorm = 1. _d 0/cg2dNbuf(1,1)  
103        ELSE        ELSE
104         myNorm = 1. _d 0         myNorm = 1. _d 0
105        ENDIF        ENDIF
       _BEGIN_MASTER( myThid )  
106         cg2dNorm = myNorm         cg2dNorm = myNorm
107          _BEGIN_MASTER( myThid )
108  CcnhDebugStarts  CcnhDebugStarts
109         WRITE(msgBuf,*) ' CG2D normalisation factor = ', cg2dNorm         WRITE(msgBuf,'(A,E40.25)') '// CG2D normalisation factor = ',
110         &               cg2dNorm
111           CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
112           WRITE(msgBuf,*) '                               '
113         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
114  CcnhDebugEnds  CcnhDebugEnds
115        _END_MASTER( myThid )        _END_MASTER( myThid )
# Line 97  CcnhDebugEnds Line 126  CcnhDebugEnds
126                
127  C--   Update overlap regions  C--   Update overlap regions
128  CcnhDebugStarts  CcnhDebugStarts
129  C     CALL PLOT_FIELD_XYR8( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )
130  C     CALL PLOT_FIELD_XYR8( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )
131  CcnhDebugEnds  CcnhDebugEnds
132        _EXCH_XY_R4(aW2d, myThid)        _EXCH_XY_R4(aW2d, myThid)
133        _EXCH_XY_R4(aS2d, myThid)        _EXCH_XY_R4(aS2d, myThid)
134  CcnhDebugStarts  CcnhDebugStarts
135  C     CALL PLOT_FIELD_XYR8( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )
136  C     CALL PLOT_FIELD_XYR8( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )
137  CcnhDebugEnds  CcnhDebugEnds
138    
139  C--   Initialise preconditioner  C--   Initialise preconditioner
140    C     Note. 20th May 1998
141    C           I made a weird discovery! In the model paper we argue
142    C           for the form of the preconditioner used here ( see
143    C           A Finite-volume, Incompressible Navier-Stokes Model
144    C           ...., Marshall et. al ). The algebra gives a simple
145    C           0.5 factor for the averaging of ac and aCw to get a
146    C           symmettric pre-conditioner. By using a factor of 0.51
147    C           i.e. scaling the off-diagonal terms in the
148    C           preconditioner down slightly I managed to get the
149    C           number of iterations for convergence in a test case to
150    C           drop form 192 -> 134! Need to investigate this further!
151    C           For now I have introduced a parameter cg2dpcOffDFac which
152    C           defaults to 0.51 but can be set at runtime.
153        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
154         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
155          DO J=1,sNy          DO J=1,sNy
156           DO I=1,sNx           DO I=1,sNx
157            pC(I,J,bi,bj) = 1. _d 0            pC(I,J,bi,bj) = 1. _d 0
158            IF (            aC = -(
159       &     aW2d(I,J,bi,bj) + aW2d(I+1,J,bi,bj)       &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)
160       &    +aS2d(I,J,bi,bj) + aS2D(I,J+1,bi,bj)       &    +aS2d(I,J,bi,bj) + aS2d(I  ,J+1,bi,bj)
161       &     .EQ. 0.       &    +freeSurfFac*myNorm* horiVertRatio*
162       &       )  pC(I,J,bi,bj) = 0. _d 0       &     rA(I,J,bi,bj)/deltaTMom/deltaTMom
163            pW(I,J,bi,bj) = 0.       &    )
164            pS(I,J,bi,bj) = 0.            aCs = -(
165         &     aW2d(I,J-1,bi,bj) + aW2d(I+1,J-1,bi,bj)
166         &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)
167         &    +freeSurfFac*myNorm* horiVertRatio*
168         &     rA(I,J-1,bi,bj)/deltaTMom/deltaTMom
169         &    )
170              aCw = -(
171         &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)
172         &    +aS2d(I-1,J,bi,bj) + aS2d(I-1,J+1,bi,bj)
173         &    +freeSurfFac*myNorm* horiVertRatio*
174         &     rA(I-1,J,bi,bj)/deltaTMom/deltaTMom
175         &    )
176              IF ( aC .EQ. 0. ) THEN
177                pC(I,J,bi,bj) = 1. _d 0
178              ELSE
179               pC(I,J,bi,bj) =  1. _d 0 / aC
180              ENDIF
181              IF ( aC + aCw .EQ. 0. ) THEN
182               pW(I,J,bi,bj) = 0.
183              ELSE
184               pW(I,J,bi,bj) =
185         &     -aW2d(I  ,J  ,bi,bj)/((cg2dpcOffDFac *(aCw+aC))**2 )
186              ENDIF
187              IF ( aC + aCs .EQ. 0. ) THEN
188               pS(I,J,bi,bj) = 0.
189              ELSE
190               pS(I,J,bi,bj) =
191         &     -aS2d(I  ,J  ,bi,bj)/((cg2dpcOffDFac *(aCs+aC))**2 )
192              ENDIF
193    C         pC(I,J,bi,bj) = 1.
194    C         pW(I,J,bi,bj) = 0.
195    C         pS(I,J,bi,bj) = 0.
196           ENDDO           ENDDO
197          ENDDO          ENDDO
198         ENDDO         ENDDO
# Line 128  C--   Update overlap regions Line 201  C--   Update overlap regions
201        _EXCH_XY_R4(pC, myThid)        _EXCH_XY_R4(pC, myThid)
202        _EXCH_XY_R4(pW, myThid)        _EXCH_XY_R4(pW, myThid)
203        _EXCH_XY_R4(pS, myThid)        _EXCH_XY_R4(pS, myThid)
204    CcnhDebugStarts
205    C     CALL PLOT_FIELD_XYRS( pC, 'pC   INI_CG2D.2' , 1, myThid )
206    C     CALL PLOT_FIELD_XYRS( pW, 'pW   INI_CG2D.2' , 1, myThid )
207    C     CALL PLOT_FIELD_XYRS( pS, 'pS   INI_CG2D.2' , 1, myThid )
208    CcnhDebugEnds
209    
210  C--   Set default values for initial guess  C--   Set default values for initial guess and RHS
211        DO bj=myByLo(myThid),myByHi(myThid)        IF ( startTime .EQ. 0 ) THEN
212         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
213          DO J=1,sNy          DO bi=myBxLo(myThid),myBxHi(myThid)
214           DO I=1,sNx           DO J=1,sNy
215            cg2d_x(I,J,bi,bj) = 0. _d 0            DO I=1,sNx
216               cg2d_x(I,J,bi,bj) = 0. _d 0
217               cg2d_b(I,J,bi,bj) = 0. _d 0
218    #ifdef INCLUDE_CD_CODE
219               cg2d_xNM1(I,J,bi,bj) = 0. _d 0
220    #endif
221              ENDDO
222           ENDDO           ENDDO
223          ENDDO          ENDDO
224         ENDDO         ENDDO
225        ENDDO  C--    Update overlap regions
226  C--   Update overlap regions         _EXCH_XY_R8(cg2d_x, myThid)
227        _EXCH_XY_R8(cg2d_x, myThid)         _EXCH_XY_R8(cg2d_b, myThid)
228    #ifdef INCLUDE_CD_CODE
229           _EXCH_XY_R8(cg2d_xNM1, myThid)
230    #endif
231          ENDIF
232    
233        RETURN        RETURN
234        END        END
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.30
  ViewVC Help
Powered by ViewVC 1.1.22