/[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.23 by adcroft, Wed Dec 9 16:11:52 1998 UTC revision 1.44 by jmc, Tue Oct 17 18:52:34 2006 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7  CStartOfInterface  CBOP
8    C     !ROUTINE: INI_CG2D
9    C     !INTERFACE:
10        SUBROUTINE INI_CG2D( myThid )        SUBROUTINE INI_CG2D( myThid )
 C     /==========================================================\  
 C     | SUBROUTINE INI_CG2D                                      |  
 C     | o Initialise 2d conjugate gradient solver operators.     |  
 C     |==========================================================|  
 C     | These arrays are purely a function of the basin geom.    |  
 C     | We set then here once and them use then repeatedly.      |  
 C     \==========================================================/  
       IMPLICIT NONE  
11    
12    C     !DESCRIPTION: \bv
13    C     *==========================================================*
14    C     | SUBROUTINE INI_CG2D                                      
15    C     | o Initialise 2d conjugate gradient solver operators.      
16    C     *==========================================================*
17    C     | These arrays are purely a function of the basin geom.    
18    C     | We set then here once and them use then repeatedly.      
19    C     *==========================================================*
20    C     \ev
21    
22    C     !USES:
23          IMPLICIT NONE
24  C     === Global variables ===  C     === Global variables ===
25  #include "SIZE.h"  #include "SIZE.h"
26  #include "EEPARAMS.h"  #include "EEPARAMS.h"
27  #include "PARAMS.h"  #include "PARAMS.h"
28  #include "GRID.h"  #include "GRID.h"
29  #include "DYNVARS.h"  #include "SURFACE.h"
30  #include "CG2D.h"  #include "CG2D.h"
31    #ifdef ALLOW_OBCS
32  #include "OBCS.h"  #include "OBCS.h"
33    #endif
34    
35    C     !INPUT/OUTPUT PARAMETERS:
36  C     === Routine arguments ===  C     === Routine arguments ===
37  C     myThid - Thread no. that called this routine.  C     myThid - Thread no. that called this routine.
38        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
39    
40    C     !LOCAL VARIABLES:
41  C     === Local variables ===  C     === Local variables ===
42  C     xG, yG - Global coordinate location.  C     xG, yG - Global coordinate location.
43  C     zG  C     zG
# Line 34  C     iG, jG - Global coordinate index Line 45  C     iG, jG - Global coordinate index
45  C     bi,bj  - Loop counters  C     bi,bj  - Loop counters
46  C     faceArea - Temporary used to hold cell face areas.  C     faceArea - Temporary used to hold cell face areas.
47  C     I,J,K  C     I,J,K
48  C     myNorm - Work variable used in clculating normalisation factor  C     myNorm - Work variable used in calculating normalisation factor
49    C     sumArea - Work variable used to compute the total Domain Area
50        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
51        INTEGER bi, bj        INTEGER bi, bj
52        INTEGER I,  J, K        INTEGER I,  J, K
53        _RL     faceArea        _RL     faceArea
54        _RS     myNorm        _RS     myNorm
55        _RL     aC, aCw, aCs        _RL     aC, aCw, aCs
56    CEOP
57    
58    C--   Initialize arrays in common blocs (CG2D.h) ; not really necessary
59    C     but safer when EXCH do not fill all the overlap regions.
60          DO bj=myByLo(myThid),myByHi(myThid)
61           DO bi=myBxLo(myThid),myBxHi(myThid)
62            DO J=1-OLy,sNy+OLy
63             DO I=1-OLx,sNx+OLx
64              aW2d(I,J,bi,bj) = 0. _d 0
65              aS2d(I,J,bi,bj) = 0. _d 0
66              pW(I,J,bi,bj) = 0. _d 0
67              pS(I,J,bi,bj) = 0. _d 0
68              pC(I,J,bi,bj) = 0. _d 0
69              cg2d_q(I,J,bi,bj) = 0. _d 0
70             ENDDO
71            ENDDO
72            DO J=1-1,sNy+1
73             DO I=1-1,sNx+1
74              cg2d_r(I,J,bi,bj) = 0. _d 0
75              cg2d_s(I,J,bi,bj) = 0. _d 0
76             ENDDO
77            ENDDO
78           ENDDO
79          ENDDO
80    
81  C--   Initialise laplace operator  C--   Initialise laplace operator
82  C     aW2d: integral in Z Ax/dX  C     aW2d: integral in Z Ax/dX
# Line 60  C     aS2d: integral in Z Ay/dY Line 96  C     aS2d: integral in Z Ay/dY
96             faceArea = _dyG(I,J,bi,bj)*drF(K)             faceArea = _dyG(I,J,bi,bj)*drF(K)
97       &               *_hFacW(I,J,K,bi,bj)       &               *_hFacW(I,J,K,bi,bj)
98             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +
99       &      gBaro*faceArea*recip_dxC(I,J,bi,bj)       &      implicSurfPress*implicDiv2DFlow
100         &           *faceArea*recip_dxC(I,J,bi,bj)
101             faceArea = _dxG(I,J,bi,bj)*drF(K)             faceArea = _dxG(I,J,bi,bj)*drF(K)
102       &               *_hFacS(I,J,K,bi,bj)       &               *_hFacS(I,J,K,bi,bj)
103             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +
104       &      gBaro*faceArea*recip_dyC(I,J,bi,bj)       &      implicSurfPress*implicDiv2DFlow
105         &           *faceArea*recip_dyC(I,J,bi,bj)
106            ENDDO            ENDDO
107           ENDDO           ENDDO
108          ENDDO          ENDDO
109          IF (openBoundaries) THEN  #ifdef ALLOW_OBCS
110            IF (useOBCS) THEN
111           DO I=1,sNx           DO I=1,sNx
112            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),bi,bj)=0.
113            IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj)+1,bi,bj)=0.            IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj)+1,bi,bj)=0.
# Line 82  C     aS2d: integral in Z Ay/dY Line 121  C     aS2d: integral in Z Ay/dY
121            IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj),J,bi,bj)=0.            IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj),J,bi,bj)=0.
122           ENDDO           ENDDO
123          ENDIF          ENDIF
124    #endif
125          DO J=1,sNy          DO J=1,sNy
126           DO I=1,sNx           DO I=1,sNx
127            myNorm = MAX(ABS(aW2d(I,J,bi,bj)),myNorm)            myNorm = MAX(ABS(aW2d(I,J,bi,bj)),myNorm)
# Line 90  C     aS2d: integral in Z Ay/dY Line 130  C     aS2d: integral in Z Ay/dY
130          ENDDO          ENDDO
131         ENDDO         ENDDO
132        ENDDO        ENDDO
133        cg2dNbuf(1,myThid) = myNorm        _GLOBAL_MAX_R4( myNorm, myThid )
134        _GLOBAL_MAX_R4( cg2dNbuf, myNorm, myThid )        IF ( myNorm .NE. 0. _d 0 ) THEN
135        IF ( cg2dNbuf(1,1) .NE. 0. _d 0 ) THEN         myNorm = 1. _d 0/myNorm
        myNorm = 1. _d 0/cg2dNbuf(1,1)  
136        ELSE        ELSE
137         myNorm = 1. _d 0         myNorm = 1. _d 0
138        ENDIF        ENDIF
        cg2dNorm = myNorm  
       _BEGIN_MASTER( myThid )  
 CcnhDebugStarts  
        WRITE(msgBuf,'(A,E40.25)') '// CG2D normalisation factor = ',  
      &               cg2dNorm  
        CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)  
        WRITE(msgBuf,*) '                               '  
        CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)  
 CcnhDebugEnds  
       _END_MASTER( myThid )  
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          DO J=1,sNy
# Line 123  CcnhDebugStarts Line 152  CcnhDebugStarts
152  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )
153  C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )
154  CcnhDebugEnds  CcnhDebugEnds
155        _EXCH_XY_R4(aW2d, myThid)  c     _EXCH_XY_R4(aW2d, myThid)
156        _EXCH_XY_R4(aS2d, myThid)  c     _EXCH_XY_R4(aS2d, myThid)
157          CALL EXCH_UV_XY_RS(aW2d,aS2d,.FALSE.,myThid)
158  CcnhDebugStarts  CcnhDebugStarts
159        CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )
160        CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )
161  CcnhDebugEnds  CcnhDebugEnds
162    
163          _BEGIN_MASTER(myThid)
164    C-- set global parameter in common block:
165          cg2dNorm = myNorm
166    C-- Define the solver tolerance in the appropriate Unit :
167          cg2dNormaliseRHS = cg2dTargetResWunit.LE.0.
168          IF (cg2dNormaliseRHS) THEN
169    C-  when using a normalisation of RHS, tolerance has no unit => no conversion
170            cg2dTolerance = cg2dTargetResidual
171          ELSE
172    C-  convert Target-Residual (in W unit) to cg2d-solver residual unit [m^2/s^2]
173            cg2dTolerance = cg2dNorm * cg2dTargetResWunit
174         &                           * globalArea / deltaTmom
175          ENDIF
176          _END_MASTER(myThid)
177    
178    CcnhDebugStarts
179          _BEGIN_MASTER( myThid )
180           WRITE(msgBuf,'(2A,1PE23.16)') 'INI_CG2D: ',
181         &      'CG2D normalisation factor = ', cg2dNorm
182           CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
183           IF (.NOT.cg2dNormaliseRHS) THEN
184            WRITE(msgBuf,'(2A,1PE22.15,A,1PE16.10,A)') 'INI_CG2D: ',
185         &      'cg2dTolerance =', cg2dTolerance, ' (Area=',globalArea,')'
186            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
187           ENDIF
188           WRITE(msgBuf,*) ' '
189           CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
190          _END_MASTER( myThid )
191    CcnhDebugEnds
192        
193  C--   Initialise preconditioner  C--   Initialise preconditioner
194  C     Note. 20th May 1998  C     Note. 20th May 1998
195  C           I made a weird discovery! In the model paper we argue  C           I made a weird discovery! In the model paper we argue
# Line 151  C           defaults to 0.51 but can be Line 211  C           defaults to 0.51 but can be
211            pC(I,J,bi,bj) = 1. _d 0            pC(I,J,bi,bj) = 1. _d 0
212            aC = -(            aC = -(
213       &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)       &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)
214       &    +aS2d(I,J,bi,bj) + aS2D(I  ,J+1,bi,bj)       &    +aS2d(I,J,bi,bj) + aS2d(I  ,J+1,bi,bj)
215       &    +freeSurfFac*myNorm* horiVertRatio*       &    +freeSurfFac*myNorm*recip_Bo(I,J,bi,bj)*
216       &     rA(I,J,bi,bj)/deltaTMom/deltaTMom       &     rA(I,J,bi,bj)/deltaTMom/deltaTfreesurf
217       &    )       &    )
218            aCs = -(            aCs = -(
219       &     aW2d(I,J-1,bi,bj) + aW2d(I+1,J-1,bi,bj)       &     aW2d(I,J-1,bi,bj) + aW2d(I+1,J-1,bi,bj)
220       &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)       &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)
221       &    +freeSurfFac*myNorm* horiVertRatio*       &    +freeSurfFac*myNorm*recip_Bo(I,J-1,bi,bj)*
222       &     rA(I,J-1,bi,bj)/deltaTMom/deltaTMom       &     rA(I,J-1,bi,bj)/deltaTMom/deltaTfreesurf
223       &    )       &    )
224            aCw = -(            aCw = -(
225       &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)       &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)
226       &    +aS2d(I-1,J,bi,bj) + aS2d(I-1,J+1,bi,bj)       &    +aS2d(I-1,J,bi,bj) + aS2d(I-1,J+1,bi,bj)
227       &    +freeSurfFac*myNorm* horiVertRatio*       &    +freeSurfFac*myNorm*recip_Bo(I-1,J,bi,bj)*
228       &     rA(I-1,J,bi,bj)/deltaTMom/deltaTMom       &     rA(I-1,J,bi,bj)/deltaTMom/deltaTfreesurf
229       &    )       &    )
230            IF ( aC .EQ. 0. ) THEN            IF ( aC .EQ. 0. ) THEN
231              pC(I,J,bi,bj) = 0. _d 0              pC(I,J,bi,bj) = 1. _d 0
232            ELSE            ELSE
233             pC(I,J,bi,bj) =  1. _d 0 / aC             pC(I,J,bi,bj) =  1. _d 0 / aC
234            ENDIF            ENDIF
# Line 193  C         pS(I,J,bi,bj) = 0. Line 253  C         pS(I,J,bi,bj) = 0.
253        ENDDO        ENDDO
254  C--   Update overlap regions  C--   Update overlap regions
255        _EXCH_XY_R4(pC, myThid)        _EXCH_XY_R4(pC, myThid)
256        _EXCH_XY_R4(pW, myThid)  c     _EXCH_XY_R4(pW, myThid)
257        _EXCH_XY_R4(pS, myThid)  c     _EXCH_XY_R4(pS, myThid)
258          CALL EXCH_UV_XY_RS(pW,pS,.FALSE.,myThid)
259  CcnhDebugStarts  CcnhDebugStarts
260        CALL PLOT_FIELD_XYRS( pC, 'pC   INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( pC, 'pC   INI_CG2D.2' , 1, myThid )
261        CALL PLOT_FIELD_XYRS( pW, 'pW   INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( pW, 'pW   INI_CG2D.2' , 1, myThid )
262        CALL PLOT_FIELD_XYRS( pS, 'pS   INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( pS, 'pS   INI_CG2D.2' , 1, myThid )
263  CcnhDebugEnds  CcnhDebugEnds
264    
 C--   Set default values for initial guess and RHS  
       IF ( startTime .EQ. 0 ) THEN  
        DO bj=myByLo(myThid),myByHi(myThid)  
         DO bi=myBxLo(myThid),myBxHi(myThid)  
          DO J=1,sNy  
           DO I=1,sNx  
            cg2d_x(I,J,bi,bj) = 0. _d 0  
            cg2d_b(I,J,bi,bj) = 0. _d 0  
 #ifdef INCLUDE_CD_CODE  
            cg2d_xNM1(I,J,bi,bj) = 0. _d 0  
 #endif  
           ENDDO  
          ENDDO  
         ENDDO  
        ENDDO  
 C--    Update overlap regions  
        _EXCH_XY_R8(cg2d_x, myThid)  
        _EXCH_XY_R8(cg2d_b, myThid)  
 #ifdef INCLUDE_CD_CODE  
        _EXCH_XY_R8(cg2d_xNM1, myThid)  
 #endif  
       ENDIF  
   
265        RETURN        RETURN
266        END        END
Legend:
Removed from v.1.23  
changed lines
  Added in v.1.44
  ViewVC Help
Powered by ViewVC 1.1.22