/[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.44 by jmc, Tue Oct 17 18:52:34 2006 UTC revision 1.45 by jmc, Tue Dec 5 05:25:08 2006 UTC
# 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     aC, aCw, aCs        _RS     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 63  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
# Line 93  C     aS2d: integral in Z Ay/dY Line 91  C     aS2d: integral in Z Ay/dY
91          DO K=1,Nr          DO K=1,Nr
92           DO J=1,sNy           DO J=1,sNy
93            DO I=1,sNx            DO I=1,sNx
94    C  deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect
95             faceArea = _dyG(I,J,bi,bj)*drF(K)             faceArea = _dyG(I,J,bi,bj)*drF(K)
96       &               *_hFacW(I,J,K,bi,bj)       &               *_hFacW(I,J,K,bi,bj)
97             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj)
98       &      implicSurfPress*implicDiv2DFlow       &              + implicSurfPress*implicDiv2DFlow
99       &           *faceArea*recip_dxC(I,J,bi,bj)       &               *faceArea*recip_dxC(I,J,bi,bj)
100             faceArea = _dxG(I,J,bi,bj)*drF(K)             faceArea = _dxG(I,J,bi,bj)*drF(K)
101       &               *_hFacS(I,J,K,bi,bj)       &               *_hFacS(I,J,K,bi,bj)
102             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj)
103       &      implicSurfPress*implicDiv2DFlow       &              + implicSurfPress*implicDiv2DFlow
104       &           *faceArea*recip_dyC(I,J,bi,bj)       &               *faceArea*recip_dyC(I,J,bi,bj)
105            ENDDO            ENDDO
106           ENDDO           ENDDO
107          ENDDO          ENDDO
# Line 146  C     aS2d: integral in Z Ay/dY Line 145  C     aS2d: integral in Z Ay/dY
145          ENDDO          ENDDO
146         ENDDO         ENDDO
147        ENDDO        ENDDO
148          
149  C--   Update overlap regions  C--   Update overlap regions
150  CcnhDebugStarts  CcnhDebugStarts
151  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )
# Line 189  CcnhDebugStarts Line 188  CcnhDebugStarts
188         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
189        _END_MASTER( myThid )        _END_MASTER( myThid )
190  CcnhDebugEnds  CcnhDebugEnds
191        
192  C--   Initialise preconditioner  C--   Initialise preconditioner
193  C     Note. 20th May 1998  C     Note. 20th May 1998
194  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 208  C           defaults to 0.51 but can be Line 207  C           defaults to 0.51 but can be
207         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
208          DO J=1,sNy          DO J=1,sNy
209           DO I=1,sNx           DO I=1,sNx
210              ks = ksurfC(I,J,bi,bj)
211            pC(I,J,bi,bj) = 1. _d 0            pC(I,J,bi,bj) = 1. _d 0
212            aC = -(            aC = -(
213       &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)       &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)
214       &    +aS2d(I,J,bi,bj) + aS2d(I  ,J+1,bi,bj)       &    +aS2d(I,J,bi,bj) + aS2d(I  ,J+1,bi,bj)
215       &    +freeSurfFac*myNorm*recip_Bo(I,J,bi,bj)*       &    +freeSurfFac*myNorm*recip_Bo(I,J,bi,bj)*deepFac2F(ks)
216       &     rA(I,J,bi,bj)/deltaTMom/deltaTfreesurf       &                *rA(I,J,bi,bj)/deltaTMom/deltaTfreesurf
217       &    )       &    )
218            aCs = -(            aCs = -(
219       &     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)
220       &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)       &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)
221       &    +freeSurfFac*myNorm*recip_Bo(I,J-1,bi,bj)*       &    +freeSurfFac*myNorm*recip_Bo(I,J-1,bi,bj)*deepFac2F(ks)
222       &     rA(I,J-1,bi,bj)/deltaTMom/deltaTfreesurf       &                *rA(I,J-1,bi,bj)/deltaTMom/deltaTfreesurf
223       &    )       &    )
224            aCw = -(            aCw = -(
225       &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)       &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)
226       &    +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)
227       &    +freeSurfFac*myNorm*recip_Bo(I-1,J,bi,bj)*       &    +freeSurfFac*myNorm*recip_Bo(I-1,J,bi,bj)*deepFac2F(ks)
228       &     rA(I-1,J,bi,bj)/deltaTMom/deltaTfreesurf       &                *rA(I-1,J,bi,bj)/deltaTMom/deltaTfreesurf
229       &    )       &    )
230            IF ( aC .EQ. 0. ) THEN            IF ( aC .EQ. 0. ) THEN
231              pC(I,J,bi,bj) = 1. _d 0              pC(I,J,bi,bj) = 1. _d 0
# Line 235  C           defaults to 0.51 but can be Line 235  C           defaults to 0.51 but can be
235            IF ( aC + aCw .EQ. 0. ) THEN            IF ( aC + aCw .EQ. 0. ) THEN
236             pW(I,J,bi,bj) = 0.             pW(I,J,bi,bj) = 0.
237            ELSE            ELSE
238             pW(I,J,bi,bj) =             pW(I,J,bi,bj) =
239       &     -aW2d(I  ,J  ,bi,bj)/((cg2dpcOffDFac *(aCw+aC))**2 )       &     -aW2d(I  ,J  ,bi,bj)/((cg2dpcOffDFac *(aCw+aC))**2 )
240            ENDIF            ENDIF
241            IF ( aC + aCs .EQ. 0. ) THEN            IF ( aC + aCs .EQ. 0. ) THEN
# Line 244  C           defaults to 0.51 but can be Line 244  C           defaults to 0.51 but can be
244             pS(I,J,bi,bj) =             pS(I,J,bi,bj) =
245       &     -aS2d(I  ,J  ,bi,bj)/((cg2dpcOffDFac *(aCs+aC))**2 )       &     -aS2d(I  ,J  ,bi,bj)/((cg2dpcOffDFac *(aCs+aC))**2 )
246            ENDIF            ENDIF
247    C-    store solver main diagonal :
248              aC2d(I,J,bi,bj) = aC
249  C         pC(I,J,bi,bj) = 1.  C         pC(I,J,bi,bj) = 1.
250  C         pW(I,J,bi,bj) = 0.  C         pW(I,J,bi,bj) = 0.
251  C         pS(I,J,bi,bj) = 0.  C         pS(I,J,bi,bj) = 0.
Legend:
Removed from v.1.44  
changed lines
  Added in v.1.45
  ViewVC Help
Powered by ViewVC 1.1.22