/[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.14 by cnh, Mon Jun 15 05:13:56 1998 UTC revision 1.23 by adcroft, Wed Dec 9 16:11:52 1998 UTC
# 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"
# Line 19  C     === Global variables === Line 20  C     === Global variables ===
20  #include "GRID.h"  #include "GRID.h"
21  #include "DYNVARS.h"  #include "DYNVARS.h"
22  #include "CG2D.h"  #include "CG2D.h"
23    #include "OBCS.h"
24    
25  C     === Routine arguments ===  C     === Routine arguments ===
26  C     myThid - Thread no. that called this routine.  C     myThid - Thread no. that called this routine.
# Line 37  C     myNorm - Work variable used in clc Line 39  C     myNorm - Work variable used in clc
39        INTEGER bi, bj        INTEGER bi, bj
40        INTEGER I,  J, K        INTEGER I,  J, K
41        _RL     faceArea        _RL     faceArea
42        _RL     myNorm        _RS     myNorm
43        _RL     aC, aCw, aCs        _RL     aC, aCw, aCs
44    
45  C--   Initialise laplace operator  C--   Initialise laplace operator
# Line 52  C     aS2d: integral in Z Ay/dY Line 54  C     aS2d: integral in Z Ay/dY
54            aS2d(I,J,bi,bj) = 0. _d 0            aS2d(I,J,bi,bj) = 0. _d 0
55           ENDDO           ENDDO
56          ENDDO          ENDDO
57          DO K=1,Nz          DO K=1,Nr
58           DO J=1,sNy           DO J=1,sNy
59            DO I=1,sNx            DO I=1,sNx
60             faceArea = _dyG(I,J,bi,bj)*dzF(K)*_hFacW(I,J,K,bi,bj)             faceArea = _dyG(I,J,bi,bj)*drF(K)
61         &               *_hFacW(I,J,K,bi,bj)
62             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +
63       &      gBaro*faceArea*_rdxC(I,J,bi,bj)       &      gBaro*faceArea*recip_dxC(I,J,bi,bj)
64             faceArea = _dxG(I,J,bi,bj)*dzF(K)*_hFacS(I,J,K,bi,bj)             faceArea = _dxG(I,J,bi,bj)*drF(K)
65         &               *_hFacS(I,J,K,bi,bj)
66             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +
67       &      gBaro*faceArea*_rdyC(I,J,bi,bj)       &      gBaro*faceArea*recip_dyC(I,J,bi,bj)
68            ENDDO            ENDDO
69           ENDDO           ENDDO
70          ENDDO          ENDDO
71            IF (openBoundaries) THEN
72             DO I=1,sNx
73              IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj),bi,bj)=0.
74              IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj)+1,bi,bj)=0.
75              IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj)+1,bi,bj)=0.
76              IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj),bi,bj)=0.
77             ENDDO
78             DO J=1,sNy
79              IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj),J,bi,bj)=0.
80              IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj)+1,J,bi,bj)=0.
81              IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj)+1,J,bi,bj)=0.
82              IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj),J,bi,bj)=0.
83             ENDDO
84            ENDIF
85          DO J=1,sNy          DO J=1,sNy
86           DO I=1,sNx           DO I=1,sNx
87            myNorm = MAX(ABS(aW2d(I,J,bi,bj)),myNorm)            myNorm = MAX(ABS(aW2d(I,J,bi,bj)),myNorm)
# Line 73  C     aS2d: integral in Z Ay/dY Line 91  C     aS2d: integral in Z Ay/dY
91         ENDDO         ENDDO
92        ENDDO        ENDDO
93        cg2dNbuf(1,myThid) = myNorm        cg2dNbuf(1,myThid) = myNorm
94        _GLOBAL_MAX_R8( cg2dNbuf, myNorm, myThid )        _GLOBAL_MAX_R4( cg2dNbuf, myNorm, myThid )
95        IF ( cg2dNbuf(1,1) .NE. 0. _d 0 ) THEN        IF ( cg2dNbuf(1,1) .NE. 0. _d 0 ) THEN
96         myNorm = 1. _d 0/cg2dNbuf(1,1)         myNorm = 1. _d 0/cg2dNbuf(1,1)
97        ELSE        ELSE
# Line 82  C     aS2d: integral in Z Ay/dY Line 100  C     aS2d: integral in Z Ay/dY
100         cg2dNorm = myNorm         cg2dNorm = myNorm
101        _BEGIN_MASTER( myThid )        _BEGIN_MASTER( myThid )
102  CcnhDebugStarts  CcnhDebugStarts
103         WRITE(msgBuf,'(A,F)') '// CG2D normalisation factor = ', cg2dNorm         WRITE(msgBuf,'(A,E40.25)') '// CG2D normalisation factor = ',
104         &               cg2dNorm
105         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
106         WRITE(msgBuf,*) '                               '         WRITE(msgBuf,*) '                               '
107         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
# Line 107  CcnhDebugEnds Line 126  CcnhDebugEnds
126        _EXCH_XY_R4(aW2d, myThid)        _EXCH_XY_R4(aW2d, myThid)
127        _EXCH_XY_R4(aS2d, myThid)        _EXCH_XY_R4(aS2d, myThid)
128  CcnhDebugStarts  CcnhDebugStarts
129  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )        CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )
130  C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )        CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )
131  CcnhDebugEnds  CcnhDebugEnds
132    
133  C--   Initialise preconditioner  C--   Initialise preconditioner
# Line 133  C           defaults to 0.51 but can be Line 152  C           defaults to 0.51 but can be
152            aC = -(            aC = -(
153       &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)       &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)
154       &    +aS2d(I,J,bi,bj) + aS2D(I  ,J+1,bi,bj)       &    +aS2d(I,J,bi,bj) + aS2D(I  ,J+1,bi,bj)
155       &    +freeSurfFac*myNorm*       &    +freeSurfFac*myNorm* horiVertRatio*
156       &     zA(I,J,bi,bj)/deltaTMom/deltaTMom       &     rA(I,J,bi,bj)/deltaTMom/deltaTMom
157       &    )       &    )
158            aCs = -(            aCs = -(
159       &     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)
160       &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)       &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)
161       &    +freeSurfFac*myNorm*       &    +freeSurfFac*myNorm* horiVertRatio*
162       &     zA(I,J-1,bi,bj)/deltaTMom/deltaTMom       &     rA(I,J-1,bi,bj)/deltaTMom/deltaTMom
163       &    )       &    )
164            aCw = -(            aCw = -(
165       &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)       &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)
166       &    +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)
167       &    +freeSurfFac*myNorm*       &    +freeSurfFac*myNorm* horiVertRatio*
168       &     zA(I-1,J,bi,bj)/deltaTMom/deltaTMom       &     rA(I-1,J,bi,bj)/deltaTMom/deltaTMom
169       &    )       &    )
170            IF ( aC .EQ. 0. ) THEN            IF ( aC .EQ. 0. ) THEN
171              pC(I,J,bi,bj) = 0. _d 0              pC(I,J,bi,bj) = 0. _d 0
# Line 176  C--   Update overlap regions Line 195  C--   Update overlap regions
195        _EXCH_XY_R4(pC, myThid)        _EXCH_XY_R4(pC, myThid)
196        _EXCH_XY_R4(pW, myThid)        _EXCH_XY_R4(pW, myThid)
197        _EXCH_XY_R4(pS, myThid)        _EXCH_XY_R4(pS, myThid)
198    CcnhDebugStarts
199          CALL PLOT_FIELD_XYRS( pC, 'pC   INI_CG2D.2' , 1, myThid )
200          CALL PLOT_FIELD_XYRS( pW, 'pW   INI_CG2D.2' , 1, myThid )
201          CALL PLOT_FIELD_XYRS( pS, 'pS   INI_CG2D.2' , 1, myThid )
202    CcnhDebugEnds
203    
204  C--   Set default values for initial guess and RHS  C--   Set default values for initial guess and RHS
205        IF ( startTime .EQ. 0 ) THEN        IF ( startTime .EQ. 0 ) THEN
# Line 185  C--   Set default values for initial gue Line 209  C--   Set default values for initial gue
209            DO I=1,sNx            DO I=1,sNx
210             cg2d_x(I,J,bi,bj) = 0. _d 0             cg2d_x(I,J,bi,bj) = 0. _d 0
211             cg2d_b(I,J,bi,bj) = 0. _d 0             cg2d_b(I,J,bi,bj) = 0. _d 0
212  #ifdef ALLOW_CD  #ifdef INCLUDE_CD_CODE
213             cg2d_xNM1(I,J,bi,bj) = 0. _d 0             cg2d_xNM1(I,J,bi,bj) = 0. _d 0
214  #endif  #endif
215            ENDDO            ENDDO
# Line 195  C--   Set default values for initial gue Line 219  C--   Set default values for initial gue
219  C--    Update overlap regions  C--    Update overlap regions
220         _EXCH_XY_R8(cg2d_x, myThid)         _EXCH_XY_R8(cg2d_x, myThid)
221         _EXCH_XY_R8(cg2d_b, myThid)         _EXCH_XY_R8(cg2d_b, myThid)
222  #ifdef ALLOW_CD  #ifdef INCLUDE_CD_CODE
223         _EXCH_XY_R8(cg2d_xNM1, myThid)         _EXCH_XY_R8(cg2d_xNM1, myThid)
224  #endif  #endif
225        ENDIF        ENDIF
Legend:
Removed from v.1.14  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.22