/[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.36 by jmc, Fri Jul 6 21:39:37 2001 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_EEOPTIONS.h"  #include "CPP_OPTIONS.h"
5    
 CStartOfInterface  
6        SUBROUTINE INI_CG2D( myThid )        SUBROUTINE INI_CG2D( myThid )
7  C     /==========================================================\  C     /==========================================================\
8  C     | SUBROUTINE INI_CG2D                                      |  C     | SUBROUTINE INI_CG2D                                      |
# Line 11  C     |================================= Line 11  C     |=================================
11  C     | These arrays are purely a function of the basin geom.    |  C     | These arrays are purely a function of the basin geom.    |
12  C     | We set then here once and them use then repeatedly.      |  C     | We set then here once and them use then repeatedly.      |
13  C     \==========================================================/  C     \==========================================================/
14          IMPLICIT NONE
15    
16  C     === Global variables ===  C     === Global variables ===
17  #include "SIZE.h"  #include "SIZE.h"
18  #include "EEPARAMS.h"  #include "EEPARAMS.h"
19  #include "PARAMS.h"  #include "PARAMS.h"
20  #include "GRID.h"  #include "GRID.h"
21    c#include "DYNVARS.h"
22    #include "SURFACE.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.
30        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
31    
32  C     === Local variables ===  C     === Local variables ===
33  C     xG, yG - Global coordinate location.  C     xG, yG - Global coordinate location.
# Line 31  C     iG, jG - Global coordinate index Line 36  C     iG, jG - Global coordinate index
36  C     bi,bj  - Loop counters  C     bi,bj  - Loop counters
37  C     faceArea - Temporary used to hold cell face areas.  C     faceArea - Temporary used to hold cell face areas.
38  C     I,J,K  C     I,J,K
39  C     myNorm - Work variable used in clculating normalisation factor  C     myNorm - Work variable used in calculating normalisation factor
40    C     sumArea - Work variable used to compute the total Domain Area
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     sumArea
47          _RL     aC, aCw, aCs
48    
49  C--   Initialise laplace operator  C--   Initialise laplace operator
50  C     aW2d: integral in Z Ax/dX  C     aW2d: integral in Z Ax/dX
# Line 50  C     aS2d: integral in Z Ay/dY Line 58  C     aS2d: integral in Z Ay/dY
58            aS2d(I,J,bi,bj) = 0. _d 0            aS2d(I,J,bi,bj) = 0. _d 0
59           ENDDO           ENDDO
60          ENDDO          ENDDO
61          DO K=1,Nz          DO K=1,Nr
62           DO J=1,sNy           DO J=1,sNy
63            DO I=1,sNx            DO I=1,sNx
64             faceArea = dyG(I,J,bi,bj)*dzF(K)*HFacW(I,J,K,bi,bj)             faceArea = _dyG(I,J,bi,bj)*drF(K)
65         &               *_hFacW(I,J,K,bi,bj)
66             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +
67       &      gravity*faceArea*rDxC(I,J,bi,bj)       &      implicSurfPress*implicDiv2DFlow
68             faceArea = dxG(I,J,bi,bj)*dzF(K)*HFacS(I,J,K,bi,bj)       &           *faceArea*recip_dxC(I,J,bi,bj)
69               faceArea = _dxG(I,J,bi,bj)*drF(K)
70         &               *_hFacS(I,J,K,bi,bj)
71             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +
72       &      gravity*faceArea*rDyC(I,J,bi,bj)       &      implicSurfPress*implicDiv2DFlow
73         &           *faceArea*recip_dyC(I,J,bi,bj)
74            ENDDO            ENDDO
75           ENDDO           ENDDO
76          ENDDO          ENDDO
77    #ifdef ALLOW_OBCS
78            IF (useOBCS) THEN
79             DO I=1,sNx
80              IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj),bi,bj)=0.
81              IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj)+1,bi,bj)=0.
82              IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj)+1,bi,bj)=0.
83              IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj),bi,bj)=0.
84             ENDDO
85             DO J=1,sNy
86              IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj),J,bi,bj)=0.
87              IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj)+1,J,bi,bj)=0.
88              IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj)+1,J,bi,bj)=0.
89              IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj),J,bi,bj)=0.
90             ENDDO
91            ENDIF
92    #endif
93          DO J=1,sNy          DO J=1,sNy
94           DO I=1,sNx           DO I=1,sNx
95            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 98  C     aS2d: integral in Z Ay/dY
98          ENDDO          ENDDO
99         ENDDO         ENDDO
100        ENDDO        ENDDO
101        cg2dNbuf(1,myThid) = myNorm        _GLOBAL_MAX_R4( myNorm, myThid )
102        _GLOBAL_MAX_R8( cg2dNbuf, myNorm, myThid )        IF ( myNorm .NE. 0. _d 0 ) THEN
103        IF ( cg2dNbuf(1,1) .NE. 0. _d 0 ) THEN         myNorm = 1. _d 0/myNorm
        myNorm = 1. _d 0/cg2dNbuf(1,1)  
104        ELSE        ELSE
105         myNorm = 1. _d 0         myNorm = 1. _d 0
106        ENDIF        ENDIF
       _BEGIN_MASTER( myThid )  
107         cg2dNorm = myNorm         cg2dNorm = myNorm
108          _BEGIN_MASTER( myThid )
109  CcnhDebugStarts  CcnhDebugStarts
110         WRITE(msgBuf,'(A,F)') '// CG2D normalisation factor = ', cg2dNorm         WRITE(msgBuf,'(A,E40.25)') '// CG2D normalisation factor = ',
111         &               cg2dNorm
112         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
113         WRITE(msgBuf,*) '                               '         WRITE(msgBuf,*) '                               '
114         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
# Line 99  CcnhDebugEnds Line 127  CcnhDebugEnds
127                
128  C--   Update overlap regions  C--   Update overlap regions
129  CcnhDebugStarts  CcnhDebugStarts
130  C     CALL PLOT_FIELD_XYR8( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )
131  C     CALL PLOT_FIELD_XYR8( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )
132  CcnhDebugEnds  CcnhDebugEnds
133        _EXCH_XY_R4(aW2d, myThid)  c     _EXCH_XY_R4(aW2d, myThid)
134        _EXCH_XY_R4(aS2d, myThid)  c     _EXCH_XY_R4(aS2d, myThid)
135          CALL EXCH_UV_XY_RS(aW2d,aS2d,.FALSE.,myThid)
136  CcnhDebugStarts  CcnhDebugStarts
137  C     CALL PLOT_FIELD_XYR8( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )
138  C     CALL PLOT_FIELD_XYR8( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )
139  CcnhDebugEnds  CcnhDebugEnds
140    
141    C-- Define the solver tolerance in the appropriate Unit :
142          cg2dNormaliseRHS = cg2dTargetResWunit.LE.0
143          IF (cg2dNormaliseRHS) THEN
144    C-  when using a normalisation of RHS, tolerance has no unit => no conversion
145            cg2dTolerance = cg2dTargetResidual
146          ELSE
147    C-  compute the total Area of the domain :
148            sumArea = 0.
149            DO bj=myByLo(myThid),myByHi(myThid)
150             DO bi=myBxLo(myThid),myBxHi(myThid)
151              DO j=1,sNy
152               DO i=1,sNx
153                 IF (Ro_surf(i,j,bi,bj).GT.R_low(i,j,bi,bj)) THEN
154                   sumArea = sumArea + rA(i,j,bi,bj)
155                 ENDIF
156               ENDDO
157              ENDDO
158             ENDDO
159            ENDDO
160    c       WRITE(*,*) ' mythid, sumArea = ', mythid, sumArea
161            _GLOBAL_SUM_R8( sumArea, myThid )
162    C-  convert Target-Residual (in W unit) to cg2d-solver residual unit [m^2/s^2]
163            cg2dTolerance = cg2dNorm * cg2dTargetResWunit
164         &                           * sumArea / deltaTMom
165            WRITE(*,'(2A,1P2E22.14)') ' ini_cg2d: ',
166         &          'sumArea,cg2dTolerance =', sumArea,cg2dTolerance
167          ENDIF
168        
169  C--   Initialise preconditioner  C--   Initialise preconditioner
170    C     Note. 20th May 1998
171    C           I made a weird discovery! In the model paper we argue
172    C           for the form of the preconditioner used here ( see
173    C           A Finite-volume, Incompressible Navier-Stokes Model
174    C           ...., Marshall et. al ). The algebra gives a simple
175    C           0.5 factor for the averaging of ac and aCw to get a
176    C           symmettric pre-conditioner. By using a factor of 0.51
177    C           i.e. scaling the off-diagonal terms in the
178    C           preconditioner down slightly I managed to get the
179    C           number of iterations for convergence in a test case to
180    C           drop form 192 -> 134! Need to investigate this further!
181    C           For now I have introduced a parameter cg2dpcOffDFac which
182    C           defaults to 0.51 but can be set at runtime.
183        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
184         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
185          DO J=1,sNy          DO J=1,sNy
186           DO I=1,sNx           DO I=1,sNx
187            pC(I,J,bi,bj) = 1. _d 0            pC(I,J,bi,bj) = 1. _d 0
188            IF (            aC = -(
189       &     aW2d(I,J,bi,bj) + aW2d(I+1,J,bi,bj)       &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)
190       &    +aS2d(I,J,bi,bj) + aS2D(I,J+1,bi,bj)       &    +aS2d(I,J,bi,bj) + aS2d(I  ,J+1,bi,bj)
191       &     .EQ. 0.       &    +freeSurfFac*myNorm*recip_Bo(I,J,bi,bj)*
192       &       )  pC(I,J,bi,bj) = 0. _d 0       &     rA(I,J,bi,bj)/deltaTMom/deltaTMom
193            pW(I,J,bi,bj) = 0.       &    )
194            pS(I,J,bi,bj) = 0.            aCs = -(
195         &     aW2d(I,J-1,bi,bj) + aW2d(I+1,J-1,bi,bj)
196         &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)
197         &    +freeSurfFac*myNorm*recip_Bo(I,J-1,bi,bj)*
198         &     rA(I,J-1,bi,bj)/deltaTMom/deltaTMom
199         &    )
200              aCw = -(
201         &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)
202         &    +aS2d(I-1,J,bi,bj) + aS2d(I-1,J+1,bi,bj)
203         &    +freeSurfFac*myNorm*recip_Bo(I-1,J,bi,bj)*
204         &     rA(I-1,J,bi,bj)/deltaTMom/deltaTMom
205         &    )
206              IF ( aC .EQ. 0. ) THEN
207                pC(I,J,bi,bj) = 1. _d 0
208              ELSE
209               pC(I,J,bi,bj) =  1. _d 0 / aC
210              ENDIF
211              IF ( aC + aCw .EQ. 0. ) THEN
212               pW(I,J,bi,bj) = 0.
213              ELSE
214               pW(I,J,bi,bj) =
215         &     -aW2d(I  ,J  ,bi,bj)/((cg2dpcOffDFac *(aCw+aC))**2 )
216              ENDIF
217              IF ( aC + aCs .EQ. 0. ) THEN
218               pS(I,J,bi,bj) = 0.
219              ELSE
220               pS(I,J,bi,bj) =
221         &     -aS2d(I  ,J  ,bi,bj)/((cg2dpcOffDFac *(aCs+aC))**2 )
222              ENDIF
223    C         pC(I,J,bi,bj) = 1.
224    C         pW(I,J,bi,bj) = 0.
225    C         pS(I,J,bi,bj) = 0.
226           ENDDO           ENDDO
227          ENDDO          ENDDO
228         ENDDO         ENDDO
229        ENDDO        ENDDO
230  C--   Update overlap regions  C--   Update overlap regions
231        _EXCH_XY_R4(pC, myThid)        _EXCH_XY_R4(pC, myThid)
232        _EXCH_XY_R4(pW, myThid)  c     _EXCH_XY_R4(pW, myThid)
233        _EXCH_XY_R4(pS, myThid)  c     _EXCH_XY_R4(pS, myThid)
234          CALL EXCH_UV_XY_RS(pW,pS,.FALSE.,myThid)
235  C--   Set default values for initial guess  CcnhDebugStarts
236        DO bj=myByLo(myThid),myByHi(myThid)  C     CALL PLOT_FIELD_XYRS( pC, 'pC   INI_CG2D.2' , 1, myThid )
237         DO bi=myBxLo(myThid),myBxHi(myThid)  C     CALL PLOT_FIELD_XYRS( pW, 'pW   INI_CG2D.2' , 1, myThid )
238          DO J=1,sNy  C     CALL PLOT_FIELD_XYRS( pS, 'pS   INI_CG2D.2' , 1, myThid )
239           DO I=1,sNx  CcnhDebugEnds
           cg2d_x(I,J,bi,bj) = 0. _d 0  
          ENDDO  
         ENDDO  
        ENDDO  
       ENDDO  
 C--   Update overlap regions  
       _EXCH_XY_R8(cg2d_x, myThid)  
240    
241        RETURN        RETURN
242        END        END
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.36
  ViewVC Help
Powered by ViewVC 1.1.22