/[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.37.6.2 by heimbach, Tue Jun 24 23:05:29 2003 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$  C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7  CBOP  CBOP
# Line 10  C     !INTERFACE: Line 11  C     !INTERFACE:
11    
12  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
13  C     *==========================================================*  C     *==========================================================*
14  C     | SUBROUTINE INI_CG2D                                        C     | SUBROUTINE INI_CG2D
15  C     | o Initialise 2d conjugate gradient solver operators.        C     | o Initialise 2d conjugate gradient solver operators.
16  C     *==========================================================*  C     *==========================================================*
17  C     | These arrays are purely a function of the basin geom.      C     | These arrays are purely a function of the basin geom.
18  C     | We set then here once and them use then repeatedly.        C     | We set then here once and them use then repeatedly.
19  C     *==========================================================*  C     *==========================================================*
20  C     \ev  C     \ev
21    
# Line 25  C     === Global variables === Line 26  C     === Global variables ===
26  #include "EEPARAMS.h"  #include "EEPARAMS.h"
27  #include "PARAMS.h"  #include "PARAMS.h"
28  #include "GRID.h"  #include "GRID.h"
 c#include "DYNVARS.h"  
29  #include "SURFACE.h"  #include "SURFACE.h"
30  #include "CG2D.h"  #include "CG2D.h"
 #ifdef ALLOW_OBCS  
 #include "OBCS.h"  
 #endif  
31    
32  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
33  C     === Routine arguments ===  C     === Routine arguments ===
# Line 39  C     myThid - Thread no. that called th Line 36  C     myThid - Thread no. that called th
36    
37  C     !LOCAL VARIABLES:  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 calculating normalisation factor  
 C     sumArea - Work variable used to compute the total Domain Area  
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     sumArea        _RS     aC, aCw, aCs
       _RL     aC, aCw, aCs  
49  CEOP  CEOP
50    
51  C--   Initialize arrays in common blocs (CG2D.h) ; not really necessary  C--   Initialize arrays in common blocs (CG2D.h) ; not really necessary
52  C     but safer when EXCH do not fill all the overlap regions.  C     but safer when EXCH do not fill all the overlap regions.
53        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
54         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
55          DO J=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
56           DO I=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
57            aW2d(I,J,bi,bj) = 0. _d 0            aW2d(i,j,bi,bj) = 0. _d 0
58            aS2d(I,J,bi,bj) = 0. _d 0            aS2d(i,j,bi,bj) = 0. _d 0
59            pW(I,J,bi,bj) = 0. _d 0            aC2d(i,j,bi,bj) = 0. _d 0
60            pS(I,J,bi,bj) = 0. _d 0            pW(i,j,bi,bj) = 0. _d 0
61            pC(I,J,bi,bj) = 0. _d 0            pS(i,j,bi,bj) = 0. _d 0
62            cg2d_q(I,J,bi,bj) = 0. _d 0            pC(i,j,bi,bj) = 0. _d 0
63           ENDDO           ENDDO
64          ENDDO          ENDDO
65          DO J=1-1,sNy+1          DO j=1-1,sNy+1
66           DO I=1-1,sNx+1           DO i=1-1,sNx+1
67            cg2d_r(I,J,bi,bj) = 0. _d 0            cg2d_q(i,j,bi,bj) = 0. _d 0
68            cg2d_s(I,J,bi,bj) = 0. _d 0            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           ENDDO
78          ENDDO          ENDDO
79         ENDDO         ENDDO
80        ENDDO        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
89  C     aS2d: integral in Z Ay/dY  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       &      implicSurfPress*implicDiv2DFlow             aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
106       &           *faceArea*recip_dxC(I,J,bi,bj)       &              + implicSurfPress*implicDiv2DFlow
107             faceArea = _dxG(I,J,bi,bj)*drF(K)       &               *faceArea*recip_dxC(i,j,bi,bj)
108       &               *_hFacS(I,J,K,bi,bj)             faceArea = _dxG(i,j,bi,bj)*drF(k)
109             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +       &               *_hFacS(i,j,k,bi,bj)
110       &      implicSurfPress*implicDiv2DFlow             aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
111       &           *faceArea*recip_dyC(I,J,bi,bj)       &              + implicSurfPress*implicDiv2DFlow
112         &               *faceArea*recip_dyC(i,j,bi,bj)
113            ENDDO            ENDDO
114           ENDDO           ENDDO
115          ENDDO          ENDDO
116            DO j=1,sNy
117             DO i=1,sNx
118  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
119          IF (useOBCS) THEN             aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
120           DO I=1,sNx       &                   *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj)
121            IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj),bi,bj)=0.             aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
122            IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj)+1,bi,bj)=0.       &                   *maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj)
123            IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj)+1,bi,bj)=0.  #endif /* ALLOW_OBCS */
124            IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj),bi,bj)=0.             myNorm = MAX(ABS(aW2d(i,j,bi,bj)),myNorm)
125           ENDDO             myNorm = MAX(ABS(aS2d(i,j,bi,bj)),myNorm)
          DO J=1,sNy  
           IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj),J,bi,bj)=0.  
           IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj)+1,J,bi,bj)=0.  
           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  
 #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        _GLOBAL_MAX_R4( myNorm, myThid )        _GLOBAL_MAX_RS( myNorm, myThid )
131        IF ( myNorm .NE. 0. _d 0 ) THEN        IF ( myNorm .NE. 0. _d 0 ) THEN
132         myNorm = 1. _d 0/myNorm         myNorm = 1. _d 0/myNorm
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  c     _EXCH_XY_R4(aW2d, myThid)        CALL EXCH_UV_XY_RS( aW2d, aS2d, .FALSE., myThid )
 c     _EXCH_XY_R4(aS2d, myThid)  
       CALL EXCH_UV_XY_RS(aW2d,aS2d,.FALSE.,myThid)  
153  CcnhDebugStarts  CcnhDebugStarts
154  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )
155  C     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  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 :  C-- Define the solver tolerance in the appropriate Unit :
162        cg2dNormaliseRHS = cg2dTargetResWunit.LE.0        cg2dNormaliseRHS = cg2dTargetResWunit.LE.0.
163        IF (cg2dNormaliseRHS) THEN        IF (cg2dNormaliseRHS) THEN
164  C-  when using a normalisation of RHS, tolerance has no unit => no conversion  C-  when using a normalisation of RHS, tolerance has no unit => no conversion
165          cg2dTolerance = cg2dTargetResidual          cg2dTolerance = cg2dTargetResidual
166        ELSE        ELSE
 C-  compute the total Area of the domain :  
         sumArea = 0.  
         DO bj=myByLo(myThid),myByHi(myThid)  
          DO bi=myBxLo(myThid),myBxHi(myThid)  
           DO j=1,sNy  
            DO i=1,sNx  
              IF (Ro_surf(i,j,bi,bj).GT.R_low(i,j,bi,bj)) THEN  
                sumArea = sumArea + rA(i,j,bi,bj)  
              ENDIF  
            ENDDO  
           ENDDO  
          ENDDO  
         ENDDO  
 c       WRITE(*,*) ' mythid, sumArea = ', mythid, sumArea  
         _GLOBAL_SUM_R8( sumArea, myThid )  
167  C-  convert Target-Residual (in W unit) to cg2d-solver residual unit [m^2/s^2]  C-  convert Target-Residual (in W unit) to cg2d-solver residual unit [m^2/s^2]
168          cg2dTolerance = cg2dNorm * cg2dTargetResWunit          cg2dTolerance = cg2dNorm * cg2dTargetResWunit
169       &                           * sumArea / deltaTmom       &                           * globalArea / deltaTmom
         _BEGIN_MASTER( myThid )  
         WRITE(standardMessageUnit,'(2A,1P2E22.14)') ' ini_cg2d: ',  
      &          'sumArea,cg2dTolerance =', sumArea,cg2dTolerance  
         _END_MASTER( myThid )  
170        ENDIF        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
187    
188  C--   Initialise preconditioner  C--   Initialise preconditioner
189  C     Note. 20th May 1998  C     Note. 20th May 1998
190  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 217  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*recip_Bo(I,J,bi,bj)*       &      +aS2d(i,j,bi,bj) + aS2d( i,j+1,bi,bj)
211       &     rA(I,J,bi,bj)/deltaTMom/deltaTfreesurf       &      +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*recip_Bo(I,J-1,bi,bj)*          DO j=1,sNy
217       &     rA(I,J-1,bi,bj)/deltaTMom/deltaTfreesurf           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*recip_Bo(I-1,J,bi,bj)*              pC(i,j,bi,bj) = 1. _d 0
223       &     rA(I-1,J,bi,bj)/deltaTMom/deltaTfreesurf             ELSE
224       &    )              pC(i,j,bi,bj) =  1. _d 0 / aC
225            IF ( aC .EQ. 0. ) THEN             ENDIF
226              pC(I,J,bi,bj) = 1. _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  c     _EXCH_XY_R4(pW, myThid)        CALL EXCH_UV_XY_RS( pW, pS, .FALSE., myThid )
 c     _EXCH_XY_R4(pS, myThid)  
       CALL EXCH_UV_XY_RS(pW,pS,.FALSE.,myThid)  
248  CcnhDebugStarts  CcnhDebugStarts
249  C     CALL PLOT_FIELD_XYRS( pC, 'pC   INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( pC, 'pC   INI_CG2D.2' , 1, myThid )
250  C     CALL PLOT_FIELD_XYRS( pW, 'pW   INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( pW, 'pW   INI_CG2D.2' , 1, myThid )
Legend:
Removed from v.1.37.6.2  
changed lines
  Added in v.1.51
  ViewVC Help
Powered by ViewVC 1.1.22