/[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.51 by heimbach, Fri Jun 15 20:13:16 2012 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"
 #include "OBCS.h"  
31    
32    C     !INPUT/OUTPUT PARAMETERS:
33  C     === Routine arguments ===  C     === Routine arguments ===
34  C     myThid - Thread no. that called this routine.  C     myThid - Thread no. that called this routine.
35        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
36    
37    C     !LOCAL VARIABLES:
38  C     === Local variables ===  C     === Local variables ===
39  C     xG, yG - Global coordinate location.  C     bi,bj  :: tile indices
40  C     zG  C     i,j,k  :: Loop counters
41  C     iG, jG - Global coordinate index  C     faceArea :: Temporary used to hold cell face areas.
42  C     bi,bj  - Loop counters  C     myNorm   :: Work variable used in calculating normalisation factor
 C     faceArea - Temporary used to hold cell face areas.  
 C     I,J,K  
 C     myNorm - Work variable used in clculating normalisation factor  
43        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
44        INTEGER bi, bj        INTEGER bi, bj
45        INTEGER I,  J, K        INTEGER i, j, k, ks
46        _RL     faceArea        _RL     faceArea
47        _RS     myNorm        _RS     myNorm
48        _RL     aC, aCw, aCs        _RS     aC, aCw, aCs
49    CEOP
50    
51    C--   Initialize arrays in common blocs (CG2D.h) ; not really necessary
52    C     but safer when EXCH do not fill all the overlap regions.
53          DO bj=myByLo(myThid),myByHi(myThid)
54           DO bi=myBxLo(myThid),myBxHi(myThid)
55            DO j=1-OLy,sNy+OLy
56             DO i=1-OLx,sNx+OLx
57              aW2d(i,j,bi,bj) = 0. _d 0
58              aS2d(i,j,bi,bj) = 0. _d 0
59              aC2d(i,j,bi,bj) = 0. _d 0
60              pW(i,j,bi,bj) = 0. _d 0
61              pS(i,j,bi,bj) = 0. _d 0
62              pC(i,j,bi,bj) = 0. _d 0
63             ENDDO
64            ENDDO
65            DO j=1-1,sNy+1
66             DO i=1-1,sNx+1
67              cg2d_q(i,j,bi,bj) = 0. _d 0
68              cg2d_r(i,j,bi,bj) = 0. _d 0
69              cg2d_s(i,j,bi,bj) = 0. _d 0
70    #ifdef ALLOW_CG2D_NSA
71              cg2d_z(i,j,bi,bj) = 0. _d 0
72    #endif /* ALLOW_CG2D_NSA */
73    #ifdef ALLOW_SRCG
74              cg2d_y(i,j,bi,bj) = 0. _d 0
75              cg2d_v(i,j,bi,bj) = 0. _d 0
76    #endif /*  ALLOW_SRCG */
77             ENDDO
78            ENDDO
79           ENDDO
80          ENDDO
81    
82    C--   Init. scalars
83          cg2dNorm = 0. _d 0
84          cg2dNormaliseRHS = .FALSE.
85          cg2dtolerance = 0. _d 0
86    
87  C--   Initialise laplace operator  C--   Initialise laplace operator
88  C     aW2d: integral in Z Ax/dX  C     aW2d: integral in Z Ax/dX
# Line 48  C     aS2d: integral in Z Ay/dY Line 90  C     aS2d: integral in Z Ay/dY
90        myNorm = 0. _d 0        myNorm = 0. _d 0
91        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
92         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
93          DO J=1,sNy          DO j=1,sNy
94           DO I=1,sNx           DO i=1,sNx
95            aW2d(I,J,bi,bj) = 0. _d 0            aW2d(i,j,bi,bj) = 0. _d 0
96            aS2d(I,J,bi,bj) = 0. _d 0            aS2d(i,j,bi,bj) = 0. _d 0
97           ENDDO           ENDDO
98          ENDDO          ENDDO
99          DO K=1,Nr          DO k=1,Nr
100           DO J=1,sNy           DO j=1,sNy
101            DO I=1,sNx            DO i=1,sNx
102             faceArea = _dyG(I,J,bi,bj)*drF(K)  C  deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect
103       &               *_hFacW(I,J,K,bi,bj)             faceArea = _dyG(i,j,bi,bj)*drF(k)
104             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +       &               *_hFacW(i,j,k,bi,bj)
105       &      gBaro*faceArea*recip_dxC(I,J,bi,bj)             aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
106             faceArea = _dxG(I,J,bi,bj)*drF(K)       &              + implicSurfPress*implicDiv2DFlow
107       &               *_hFacS(I,J,K,bi,bj)       &               *faceArea*recip_dxC(i,j,bi,bj)
108             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +             faceArea = _dxG(i,j,bi,bj)*drF(k)
109       &      gBaro*faceArea*recip_dyC(I,J,bi,bj)       &               *_hFacS(i,j,k,bi,bj)
110               aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
111         &              + implicSurfPress*implicDiv2DFlow
112         &               *faceArea*recip_dyC(i,j,bi,bj)
113            ENDDO            ENDDO
114           ENDDO           ENDDO
115          ENDDO          ENDDO
116          IF (openBoundaries) THEN          DO j=1,sNy
117           DO I=1,sNx           DO i=1,sNx
118            IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj),bi,bj)=0.  #ifdef ALLOW_OBCS
119            IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj)+1,bi,bj)=0.             aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
120            IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj)+1,bi,bj)=0.       &                   *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj)
121            IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj),bi,bj)=0.             aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
122           ENDDO       &                   *maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj)
123           DO J=1,sNy  #endif /* ALLOW_OBCS */
124            IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj),J,bi,bj)=0.             myNorm = MAX(ABS(aW2d(i,j,bi,bj)),myNorm)
125            IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj)+1,J,bi,bj)=0.             myNorm = MAX(ABS(aS2d(i,j,bi,bj)),myNorm)
           IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj)+1,J,bi,bj)=0.  
           IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj),J,bi,bj)=0.  
          ENDDO  
         ENDIF  
         DO J=1,sNy  
          DO I=1,sNx  
           myNorm = MAX(ABS(aW2d(I,J,bi,bj)),myNorm)  
           myNorm = MAX(ABS(aS2d(I,J,bi,bj)),myNorm)  
126           ENDDO           ENDDO
127          ENDDO          ENDDO
128         ENDDO         ENDDO
129        ENDDO        ENDDO
130        cg2dNbuf(1,myThid) = myNorm        _GLOBAL_MAX_RS( myNorm, myThid )
131        _GLOBAL_MAX_R4( cg2dNbuf, myNorm, myThid )        IF ( myNorm .NE. 0. _d 0 ) THEN
132        IF ( cg2dNbuf(1,1) .NE. 0. _d 0 ) THEN         myNorm = 1. _d 0/myNorm
        myNorm = 1. _d 0/cg2dNbuf(1,1)  
133        ELSE        ELSE
134         myNorm = 1. _d 0         myNorm = 1. _d 0
135        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 )  
136        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
137         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
138          DO J=1,sNy          DO j=1,sNy
139           DO I=1,sNx           DO i=1,sNx
140            aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj)*myNorm            aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)*myNorm
141            aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj)*myNorm            aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)*myNorm
142           ENDDO           ENDDO
143          ENDDO          ENDDO
144         ENDDO         ENDDO
145        ENDDO        ENDDO
146          
147  C--   Update overlap regions  C--   Update overlap regions
148  CcnhDebugStarts  CcnhDebugStarts
149  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )
150  C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )
151  CcnhDebugEnds  CcnhDebugEnds
152        _EXCH_XY_R4(aW2d, myThid)        CALL EXCH_UV_XY_RS( aW2d, aS2d, .FALSE., myThid )
       _EXCH_XY_R4(aS2d, myThid)  
153  CcnhDebugStarts  CcnhDebugStarts
154        CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )
155        CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )
156    CcnhDebugEnds
157    
158          _BEGIN_MASTER(myThid)
159    C-- set global parameter in common block:
160          cg2dNorm = myNorm
161    C-- Define the solver tolerance in the appropriate Unit :
162          cg2dNormaliseRHS = cg2dTargetResWunit.LE.0.
163          IF (cg2dNormaliseRHS) THEN
164    C-  when using a normalisation of RHS, tolerance has no unit => no conversion
165            cg2dTolerance = cg2dTargetResidual
166          ELSE
167    C-  convert Target-Residual (in W unit) to cg2d-solver residual unit [m^2/s^2]
168            cg2dTolerance = cg2dNorm * cg2dTargetResWunit
169         &                           * globalArea / deltaTmom
170          ENDIF
171          _END_MASTER(myThid)
172    
173    CcnhDebugStarts
174          _BEGIN_MASTER( myThid )
175           WRITE(msgBuf,'(2A,1PE23.16)') 'INI_CG2D: ',
176         &      'CG2D normalisation factor = ', cg2dNorm
177           CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
178           IF (.NOT.cg2dNormaliseRHS) THEN
179            WRITE(msgBuf,'(2A,1PE22.15,A,1PE16.10,A)') 'INI_CG2D: ',
180         &      'cg2dTolerance =', cg2dTolerance, ' (Area=',globalArea,')'
181            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
182           ENDIF
183           WRITE(msgBuf,*) ' '
184           CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
185          _END_MASTER( myThid )
186  CcnhDebugEnds  CcnhDebugEnds
187    
188  C--   Initialise preconditioner  C--   Initialise preconditioner
# Line 146  C           For now I have introduced a Line 201  C           For now I have introduced a
201  C           defaults to 0.51 but can be set at runtime.  C           defaults to 0.51 but can be set at runtime.
202        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
203         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
204          DO J=1,sNy  C-    calculate and store solver main diagonal :
205           DO I=1,sNx          DO j=0,sNy+1
206            pC(I,J,bi,bj) = 1. _d 0           DO i=0,sNx+1
207            aC = -(             ks = ksurfC(i,j,bi,bj)
208       &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)             aC2d(i,j,bi,bj) = -(
209       &    +aS2d(I,J,bi,bj) + aS2D(I  ,J+1,bi,bj)       &       aW2d(i,j,bi,bj) + aW2d(i+1,j, bi,bj)
210       &    +freeSurfFac*myNorm* horiVertRatio*       &      +aS2d(i,j,bi,bj) + aS2d( i,j+1,bi,bj)
211       &     rA(I,J,bi,bj)/deltaTMom/deltaTMom       &      +freeSurfFac*myNorm*recip_Bo(i,j,bi,bj)*deepFac2F(ks)
212       &    )       &                  *rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
213            aCs = -(       &                        )
214       &     aW2d(I,J-1,bi,bj) + aW2d(I+1,J-1,bi,bj)           ENDDO
215       &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)          ENDDO
216       &    +freeSurfFac*myNorm* horiVertRatio*          DO j=1,sNy
217       &     rA(I,J-1,bi,bj)/deltaTMom/deltaTMom           DO i=1,sNx
218       &    )             aC  = aC2d( i, j, bi,bj)
219            aCw = -(             aCs = aC2d( i,j-1,bi,bj)
220       &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)             aCw = aC2d(i-1,j, bi,bj)
221       &    +aS2d(I-1,J,bi,bj) + aS2d(I-1,J+1,bi,bj)             IF ( aC .EQ. 0. ) THEN
222       &    +freeSurfFac*myNorm* horiVertRatio*              pC(i,j,bi,bj) = 1. _d 0
223       &     rA(I-1,J,bi,bj)/deltaTMom/deltaTMom             ELSE
224       &    )              pC(i,j,bi,bj) =  1. _d 0 / aC
225            IF ( aC .EQ. 0. ) THEN             ENDIF
226              pC(I,J,bi,bj) = 0. _d 0             IF ( aC + aCw .EQ. 0. ) THEN
227            ELSE              pW(i,j,bi,bj) = 0.
228             pC(I,J,bi,bj) =  1. _d 0 / aC             ELSE
229            ENDIF              pW(i,j,bi,bj) =
230            IF ( aC + aCw .EQ. 0. ) THEN       &         -aW2d(i,j,bi,bj)/((cg2dpcOffDFac *(aCw+aC))**2 )
231             pW(I,J,bi,bj) = 0.             ENDIF
232            ELSE             IF ( aC + aCs .EQ. 0. ) THEN
233             pW(I,J,bi,bj) =              pS(i,j,bi,bj) = 0.
234       &     -aW2d(I  ,J  ,bi,bj)/((cg2dpcOffDFac *(aCw+aC))**2 )             ELSE
235            ENDIF              pS(i,j,bi,bj) =
236            IF ( aC + aCs .EQ. 0. ) THEN       &         -aS2d(i,j,bi,bj)/((cg2dpcOffDFac *(aCs+aC))**2 )
237             pS(I,J,bi,bj) = 0.             ENDIF
238            ELSE  C          pC(i,j,bi,bj) = 1.
239             pS(I,J,bi,bj) =  C          pW(i,j,bi,bj) = 0.
240       &     -aS2d(I  ,J  ,bi,bj)/((cg2dpcOffDFac *(aCs+aC))**2 )  C          pS(i,j,bi,bj) = 0.
           ENDIF  
 C         pC(I,J,bi,bj) = 1.  
 C         pW(I,J,bi,bj) = 0.  
 C         pS(I,J,bi,bj) = 0.  
241           ENDDO           ENDDO
242          ENDDO          ENDDO
243         ENDDO         ENDDO
244        ENDDO        ENDDO
245  C--   Update overlap regions  C--   Update overlap regions
246        _EXCH_XY_R4(pC, myThid)        CALL EXCH_XY_RS( pC, myThid )
247        _EXCH_XY_R4(pW, myThid)        CALL EXCH_UV_XY_RS( pW, pS, .FALSE., myThid )
       _EXCH_XY_R4(pS, myThid)  
248  CcnhDebugStarts  CcnhDebugStarts
249        CALL PLOT_FIELD_XYRS( pC, 'pC   INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( pC, 'pC   INI_CG2D.2' , 1, myThid )
250        CALL PLOT_FIELD_XYRS( pW, 'pW   INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( pW, 'pW   INI_CG2D.2' , 1, myThid )
251        CALL PLOT_FIELD_XYRS( pS, 'pS   INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( pS, 'pS   INI_CG2D.2' , 1, myThid )
252  CcnhDebugEnds  CcnhDebugEnds
253    
 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  
   
254        RETURN        RETURN
255        END        END
Legend:
Removed from v.1.23  
changed lines
  Added in v.1.51
  ViewVC Help
Powered by ViewVC 1.1.22