/[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.9 by cnh, Wed May 27 05:57:02 1998 UTC revision 1.22 by adcroft, Tue Dec 8 19:44:28 1998 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    
3  #include "CPP_EEOPTIONS.h"  #include "CPP_OPTIONS.h"
4    
5  CStartOfInterface  CStartOfInterface
6        SUBROUTINE INI_CG2D( myThid )        SUBROUTINE INI_CG2D( myThid )
# Line 17  C     === Global variables === Line 17  C     === Global variables ===
17  #include "EEPARAMS.h"  #include "EEPARAMS.h"
18  #include "PARAMS.h"  #include "PARAMS.h"
19  #include "GRID.h"  #include "GRID.h"
20    #include "DYNVARS.h"
21  #include "CG2D.h"  #include "CG2D.h"
22    #include "OBCS.h"
23    
24  C     === Routine arguments ===  C     === Routine arguments ===
25  C     myThid - Thread no. that called this routine.  C     myThid - Thread no. that called this routine.
# Line 36  C     myNorm - Work variable used in clc Line 38  C     myNorm - Work variable used in clc
38        INTEGER bi, bj        INTEGER bi, bj
39        INTEGER I,  J, K        INTEGER I,  J, K
40        _RL     faceArea        _RL     faceArea
41        _RL     myNorm        _RS     myNorm
42        _RL     aC, aCw, aCs        _RL     aC, aCw, aCs
43    
44  C--   Initialise laplace operator  C--   Initialise laplace operator
# Line 51  C     aS2d: integral in Z Ay/dY Line 53  C     aS2d: integral in Z Ay/dY
53            aS2d(I,J,bi,bj) = 0. _d 0            aS2d(I,J,bi,bj) = 0. _d 0
54           ENDDO           ENDDO
55          ENDDO          ENDDO
56          DO K=1,Nz          DO K=1,Nr
57           DO J=1,sNy           DO J=1,sNy
58            DO I=1,sNx            DO I=1,sNx
59             faceArea = _dyG(I,J,bi,bj)*dzF(K)*HFacW(I,J,K,bi,bj)             faceArea = _dyG(I,J,bi,bj)*drF(K)
60         &               *_hFacW(I,J,K,bi,bj)
61             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +
62       &      gBaro*faceArea*rDxC(I,J,bi,bj)       &      gBaro*faceArea*recip_dxC(I,J,bi,bj)
63             faceArea = _dxG(I,J,bi,bj)*dzF(K)*HFacS(I,J,K,bi,bj)             faceArea = _dxG(I,J,bi,bj)*drF(K)
64         &               *_hFacS(I,J,K,bi,bj)
65             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +
66       &      gBaro*faceArea*rDyC(I,J,bi,bj)       &      gBaro*faceArea*recip_dyC(I,J,bi,bj)
67            ENDDO            ENDDO
68           ENDDO           ENDDO
69          ENDDO          ENDDO
70            IF (openBoundaries) THEN
71             DO I=1,sNx
72              IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj),bi,bj)=0.
73              IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj)+1,bi,bj)=0.
74              IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj)+1,bi,bj)=0.
75              IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj),bi,bj)=0.
76             ENDDO
77             DO J=1,sNy
78              IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj),J,bi,bj)=0.
79              IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj)+1,J,bi,bj)=0.
80              IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj)+1,J,bi,bj)=0.
81              IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj),J,bi,bj)=0.
82             ENDDO
83            ENDIF
84          DO J=1,sNy          DO J=1,sNy
85           DO I=1,sNx           DO I=1,sNx
86            myNorm = MAX(ABS(aW2d(I,J,bi,bj)),myNorm)            myNorm = MAX(ABS(aW2d(I,J,bi,bj)),myNorm)
# Line 72  C     aS2d: integral in Z Ay/dY Line 90  C     aS2d: integral in Z Ay/dY
90         ENDDO         ENDDO
91        ENDDO        ENDDO
92        cg2dNbuf(1,myThid) = myNorm        cg2dNbuf(1,myThid) = myNorm
93        _GLOBAL_MAX_R8( cg2dNbuf, myNorm, myThid )        _GLOBAL_MAX_R4( cg2dNbuf, myNorm, myThid )
94        IF ( cg2dNbuf(1,1) .NE. 0. _d 0 ) THEN        IF ( cg2dNbuf(1,1) .NE. 0. _d 0 ) THEN
95         myNorm = 1. _d 0/cg2dNbuf(1,1)         myNorm = 1. _d 0/cg2dNbuf(1,1)
96        ELSE        ELSE
# Line 81  C     aS2d: integral in Z Ay/dY Line 99  C     aS2d: integral in Z Ay/dY
99         cg2dNorm = myNorm         cg2dNorm = myNorm
100        _BEGIN_MASTER( myThid )        _BEGIN_MASTER( myThid )
101  CcnhDebugStarts  CcnhDebugStarts
102         WRITE(msgBuf,'(A,F)') '// CG2D normalisation factor = ', cg2dNorm         WRITE(msgBuf,'(A,E40.25)') '// CG2D normalisation factor = ',
103         &               cg2dNorm
104         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
105         WRITE(msgBuf,*) '                               '         WRITE(msgBuf,*) '                               '
106         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
# Line 100  CcnhDebugEnds Line 119  CcnhDebugEnds
119                
120  C--   Update overlap regions  C--   Update overlap regions
121  CcnhDebugStarts  CcnhDebugStarts
122  C     CALL PLOT_FIELD_XYR8( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )
123  C     CALL PLOT_FIELD_XYR8( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )
124  CcnhDebugEnds  CcnhDebugEnds
125        _EXCH_XY_R4(aW2d, myThid)        _EXCH_XY_R4(aW2d, myThid)
126        _EXCH_XY_R4(aS2d, myThid)        _EXCH_XY_R4(aS2d, myThid)
127  CcnhDebugStarts  CcnhDebugStarts
128  C     CALL PLOT_FIELD_XYR8( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )        CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )
129  C     CALL PLOT_FIELD_XYR8( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )        CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )
130  CcnhDebugEnds  CcnhDebugEnds
131    
132  C--   Initialise preconditioner  C--   Initialise preconditioner
# Line 132  C           defaults to 0.51 but can be Line 151  C           defaults to 0.51 but can be
151            aC = -(            aC = -(
152       &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)       &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)
153       &    +aS2d(I,J,bi,bj) + aS2D(I  ,J+1,bi,bj)       &    +aS2d(I,J,bi,bj) + aS2D(I  ,J+1,bi,bj)
154       &    +freeSurfFac*myNorm*       &    +freeSurfFac*myNorm* horiVertRatio*
155       &     zA(I,J,bi,bj)/deltaTMom/deltaTMom       &     rA(I,J,bi,bj)/deltaTMom/deltaTMom
156       &    )       &    )
157            aCs = -(            aCs = -(
158       &     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)
159       &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)       &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)
160       &    +freeSurfFac*myNorm*       &    +freeSurfFac*myNorm* horiVertRatio*
161       &     zA(I,J-1,bi,bj)/deltaTMom/deltaTMom       &     rA(I,J-1,bi,bj)/deltaTMom/deltaTMom
162       &    )       &    )
163            aCw = -(            aCw = -(
164       &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)       &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)
165       &    +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)
166       &    +freeSurfFac*myNorm*       &    +freeSurfFac*myNorm* horiVertRatio*
167       &     zA(I-1,J,bi,bj)/deltaTMom/deltaTMom       &     rA(I-1,J,bi,bj)/deltaTMom/deltaTMom
168       &    )       &    )
169            IF ( aC .EQ. 0. ) THEN            IF ( aC .EQ. 0. ) THEN
170              pC(I,J,bi,bj) = 0. _d 0              pC(I,J,bi,bj) = 0. _d 0
# Line 175  C--   Update overlap regions Line 194  C--   Update overlap regions
194        _EXCH_XY_R4(pC, myThid)        _EXCH_XY_R4(pC, myThid)
195        _EXCH_XY_R4(pW, myThid)        _EXCH_XY_R4(pW, myThid)
196        _EXCH_XY_R4(pS, myThid)        _EXCH_XY_R4(pS, myThid)
197    CcnhDebugStarts
198          CALL PLOT_FIELD_XYRS( pC, 'pC   INI_CG2D.2' , 1, myThid )
199          CALL PLOT_FIELD_XYRS( pW, 'pW   INI_CG2D.2' , 1, myThid )
200          CALL PLOT_FIELD_XYRS( pS, 'pS   INI_CG2D.2' , 1, myThid )
201    CcnhDebugEnds
202    
203  C--   Set default values for initial guess and RHS  C--   Set default values for initial guess and RHS
204        IF ( startTime .EQ. 0 ) THEN        IF ( startTime .EQ. 0 ) THEN
# Line 184  C--   Set default values for initial gue Line 208  C--   Set default values for initial gue
208            DO I=1,sNx            DO I=1,sNx
209             cg2d_x(I,J,bi,bj) = 0. _d 0             cg2d_x(I,J,bi,bj) = 0. _d 0
210             cg2d_b(I,J,bi,bj) = 0. _d 0             cg2d_b(I,J,bi,bj) = 0. _d 0
211    #ifdef INCLUDE_CD_CODE
212               cg2d_xNM1(I,J,bi,bj) = 0. _d 0
213    #endif
214            ENDDO            ENDDO
215           ENDDO           ENDDO
216          ENDDO          ENDDO
# Line 191  C--   Set default values for initial gue Line 218  C--   Set default values for initial gue
218  C--    Update overlap regions  C--    Update overlap regions
219         _EXCH_XY_R8(cg2d_x, myThid)         _EXCH_XY_R8(cg2d_x, myThid)
220         _EXCH_XY_R8(cg2d_b, myThid)         _EXCH_XY_R8(cg2d_b, myThid)
221    #ifdef INCLUDE_CD_CODE
222           _EXCH_XY_R8(cg2d_xNM1, myThid)
223    #endif
224        ENDIF        ENDIF
225    
226        RETURN        RETURN
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.22