/[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.38 by jmc, Sun Feb 10 20:04:11 2002 UTC revision 1.47 by jmc, Sat Jul 11 22:00:40 2009 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"
31  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
# Line 39  C     myThid - Thread no. that called th Line 39  C     myThid - Thread no. that called th
39    
40  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
41  C     === Local variables ===  C     === Local variables ===
42  C     xG, yG - Global coordinate location.  C     bi,bj  :: tile indices
43  C     zG  C     I,J,K  :: Loop counters
 C     iG, jG - Global coordinate index  
 C     bi,bj  - Loop counters  
44  C     faceArea - Temporary used to hold cell face areas.  C     faceArea - Temporary used to hold cell face areas.
 C     I,J,K  
45  C     myNorm - Work variable used in calculating normalisation factor  C     myNorm - Work variable used in calculating normalisation factor
46  C     sumArea - Work variable used to compute the total Domain Area  C     sumArea - Work variable used to compute the total Domain Area
47        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
48        INTEGER bi, bj        INTEGER bi, bj
49        INTEGER I,  J, K        INTEGER I, J, K, ks
50        _RL     faceArea        _RL     faceArea
51        _RS     myNorm        _RS     myNorm
52        _RL     sumArea        _RS     aC, aCw, aCs
       _RL     aC, aCw, aCs  
53  CEOP  CEOP
54    
55  C--   Initialize arrays in common blocs (CG2D.h) ; not really necessary  C--   Initialize arrays in common blocs (CG2D.h) ; not really necessary
56  C     but safer when EXCH do not fill all the overlap regions.  C     but safer when EXCH do not fill all the overlap regions.
57        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
58         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
# Line 64  C     but safer when EXCH do not fill al Line 60  C     but safer when EXCH do not fill al
60           DO I=1-OLx,sNx+OLx           DO I=1-OLx,sNx+OLx
61            aW2d(I,J,bi,bj) = 0. _d 0            aW2d(I,J,bi,bj) = 0. _d 0
62            aS2d(I,J,bi,bj) = 0. _d 0            aS2d(I,J,bi,bj) = 0. _d 0
63              aC2d(I,J,bi,bj) = 0. _d 0
64            pW(I,J,bi,bj) = 0. _d 0            pW(I,J,bi,bj) = 0. _d 0
65            pS(I,J,bi,bj) = 0. _d 0            pS(I,J,bi,bj) = 0. _d 0
66            pC(I,J,bi,bj) = 0. _d 0            pC(I,J,bi,bj) = 0. _d 0
           cg2d_q(I,J,bi,bj) = 0. _d 0  
67           ENDDO           ENDDO
68          ENDDO          ENDDO
69          DO J=1-1,sNy+1          DO J=1-1,sNy+1
70           DO I=1-1,sNx+1           DO I=1-1,sNx+1
71              cg2d_q(I,J,bi,bj) = 0. _d 0
72            cg2d_r(I,J,bi,bj) = 0. _d 0            cg2d_r(I,J,bi,bj) = 0. _d 0
73            cg2d_s(I,J,bi,bj) = 0. _d 0            cg2d_s(I,J,bi,bj) = 0. _d 0
74    #ifdef ALLOW_CG2D_NSA
75              cg2d_z(I,J,bi,bj) = 0. _d 0
76    #endif /* ALLOW_CG2D_NSA */
77           ENDDO           ENDDO
78          ENDDO          ENDDO
79         ENDDO         ENDDO
# Line 94  C     aS2d: integral in Z Ay/dY Line 94  C     aS2d: integral in Z Ay/dY
94          DO K=1,Nr          DO K=1,Nr
95           DO J=1,sNy           DO J=1,sNy
96            DO I=1,sNx            DO I=1,sNx
97    C  deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect
98             faceArea = _dyG(I,J,bi,bj)*drF(K)             faceArea = _dyG(I,J,bi,bj)*drF(K)
99       &               *_hFacW(I,J,K,bi,bj)       &               *_hFacW(I,J,K,bi,bj)
100             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj)
101       &      implicSurfPress*implicDiv2DFlow       &              + implicSurfPress*implicDiv2DFlow
102       &           *faceArea*recip_dxC(I,J,bi,bj)       &               *faceArea*recip_dxC(I,J,bi,bj)
103             faceArea = _dxG(I,J,bi,bj)*drF(K)             faceArea = _dxG(I,J,bi,bj)*drF(K)
104       &               *_hFacS(I,J,K,bi,bj)       &               *_hFacS(I,J,K,bi,bj)
105             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj)
106       &      implicSurfPress*implicDiv2DFlow       &              + implicSurfPress*implicDiv2DFlow
107       &           *faceArea*recip_dyC(I,J,bi,bj)       &               *faceArea*recip_dyC(I,J,bi,bj)
108            ENDDO            ENDDO
109           ENDDO           ENDDO
110          ENDDO          ENDDO
# Line 131  C     aS2d: integral in Z Ay/dY Line 132  C     aS2d: integral in Z Ay/dY
132          ENDDO          ENDDO
133         ENDDO         ENDDO
134        ENDDO        ENDDO
135        _GLOBAL_MAX_R4( myNorm, myThid )        _GLOBAL_MAX_RS( myNorm, myThid )
136        IF ( myNorm .NE. 0. _d 0 ) THEN        IF ( myNorm .NE. 0. _d 0 ) THEN
137         myNorm = 1. _d 0/myNorm         myNorm = 1. _d 0/myNorm
138        ELSE        ELSE
139         myNorm = 1. _d 0         myNorm = 1. _d 0
140        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 )  
141        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
142         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
143          DO J=1,sNy          DO J=1,sNy
# Line 157  CcnhDebugEnds Line 148  CcnhDebugEnds
148          ENDDO          ENDDO
149         ENDDO         ENDDO
150        ENDDO        ENDDO
151          
152  C--   Update overlap regions  C--   Update overlap regions
153  CcnhDebugStarts  CcnhDebugStarts
154  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )
155  C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )
156  CcnhDebugEnds  CcnhDebugEnds
157  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)  
158  CcnhDebugStarts  CcnhDebugStarts
159  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )
160  C     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 :  C-- Define the solver tolerance in the appropriate Unit :
167        cg2dNormaliseRHS = cg2dTargetResWunit.LE.0        cg2dNormaliseRHS = cg2dTargetResWunit.LE.0.
168        IF (cg2dNormaliseRHS) THEN        IF (cg2dNormaliseRHS) THEN
169  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
170          cg2dTolerance = cg2dTargetResidual          cg2dTolerance = cg2dTargetResidual
171        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 )  
172  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]
173          cg2dTolerance = cg2dNorm * cg2dTargetResWunit          cg2dTolerance = cg2dNorm * cg2dTargetResWunit
174       &                           * sumArea / deltaTMom       &                           * globalArea / deltaTmom
         WRITE(*,'(2A,1P2E22.14)') ' ini_cg2d: ',  
      &          'sumArea,cg2dTolerance =', sumArea,cg2dTolerance  
175        ENDIF        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 217  C           defaults to 0.51 but can be Line 208  C           defaults to 0.51 but can be
208         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
209          DO J=1,sNy          DO J=1,sNy
210           DO I=1,sNx           DO I=1,sNx
211              ks = ksurfC(I,J,bi,bj)
212            pC(I,J,bi,bj) = 1. _d 0            pC(I,J,bi,bj) = 1. _d 0
213            aC = -(            aC = -(
214       &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)       &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)
215       &    +aS2d(I,J,bi,bj) + aS2d(I  ,J+1,bi,bj)       &    +aS2d(I,J,bi,bj) + aS2d(I  ,J+1,bi,bj)
216       &    +freeSurfFac*myNorm*recip_Bo(I,J,bi,bj)*       &    +freeSurfFac*myNorm*recip_Bo(I,J,bi,bj)*deepFac2F(ks)
217       &     rA(I,J,bi,bj)/deltaTMom/deltaTMom       &                *rA(I,J,bi,bj)/deltaTMom/deltaTfreesurf
218       &    )       &    )
219            aCs = -(            aCs = -(
220       &     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)
221       &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)       &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)
222       &    +freeSurfFac*myNorm*recip_Bo(I,J-1,bi,bj)*       &    +freeSurfFac*myNorm*recip_Bo(I,J-1,bi,bj)*deepFac2F(ks)
223       &     rA(I,J-1,bi,bj)/deltaTMom/deltaTMom       &                *rA(I,J-1,bi,bj)/deltaTMom/deltaTfreesurf
224       &    )       &    )
225            aCw = -(            aCw = -(
226       &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)       &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)
227       &    +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)
228       &    +freeSurfFac*myNorm*recip_Bo(I-1,J,bi,bj)*       &    +freeSurfFac*myNorm*recip_Bo(I-1,J,bi,bj)*deepFac2F(ks)
229       &     rA(I-1,J,bi,bj)/deltaTMom/deltaTMom       &                *rA(I-1,J,bi,bj)/deltaTMom/deltaTfreesurf
230       &    )       &    )
231            IF ( aC .EQ. 0. ) THEN            IF ( aC .EQ. 0. ) THEN
232              pC(I,J,bi,bj) = 1. _d 0              pC(I,J,bi,bj) = 1. _d 0
# Line 244  C           defaults to 0.51 but can be Line 236  C           defaults to 0.51 but can be
236            IF ( aC + aCw .EQ. 0. ) THEN            IF ( aC + aCw .EQ. 0. ) THEN
237             pW(I,J,bi,bj) = 0.             pW(I,J,bi,bj) = 0.
238            ELSE            ELSE
239             pW(I,J,bi,bj) =             pW(I,J,bi,bj) =
240       &     -aW2d(I  ,J  ,bi,bj)/((cg2dpcOffDFac *(aCw+aC))**2 )       &     -aW2d(I  ,J  ,bi,bj)/((cg2dpcOffDFac *(aCw+aC))**2 )
241            ENDIF            ENDIF
242            IF ( aC + aCs .EQ. 0. ) THEN            IF ( aC + aCs .EQ. 0. ) THEN
# Line 253  C           defaults to 0.51 but can be Line 245  C           defaults to 0.51 but can be
245             pS(I,J,bi,bj) =             pS(I,J,bi,bj) =
246       &     -aS2d(I  ,J  ,bi,bj)/((cg2dpcOffDFac *(aCs+aC))**2 )       &     -aS2d(I  ,J  ,bi,bj)/((cg2dpcOffDFac *(aCs+aC))**2 )
247            ENDIF            ENDIF
248    C-    store solver main diagonal :
249              aC2d(I,J,bi,bj) = aC
250  C         pC(I,J,bi,bj) = 1.  C         pC(I,J,bi,bj) = 1.
251  C         pW(I,J,bi,bj) = 0.  C         pW(I,J,bi,bj) = 0.
252  C         pS(I,J,bi,bj) = 0.  C         pS(I,J,bi,bj) = 0.
# Line 261  C         pS(I,J,bi,bj) = 0. Line 255  C         pS(I,J,bi,bj) = 0.
255         ENDDO         ENDDO
256        ENDDO        ENDDO
257  C--   Update overlap regions  C--   Update overlap regions
258        _EXCH_XY_R4(pC, myThid)        CALL EXCH_XY_RS( pC, myThid )
259  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)  
260  CcnhDebugStarts  CcnhDebugStarts
261  C     CALL PLOT_FIELD_XYRS( pC, 'pC   INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( pC, 'pC   INI_CG2D.2' , 1, myThid )
262  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.38  
changed lines
  Added in v.1.47
  ViewVC Help
Powered by ViewVC 1.1.22