/[MITgcm]/MITgcm/model/src/ini_cartesian_grid.F
ViewVC logotype

Diff of /MITgcm/model/src/ini_cartesian_grid.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.21 by jmc, Tue Jan 27 15:35:27 2009 UTC revision 1.22 by jmc, Sat Apr 17 18:25:12 2010 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4    c#include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7  CBOP  CBOP
8  C     !ROUTINE: INI_CARTESIAN_GRID  C     !ROUTINE: INI_CARTESIAN_GRID
9  C     !INTERFACE:  C     !INTERFACE:
10        SUBROUTINE INI_CARTESIAN_GRID( myThid )        SUBROUTINE INI_CARTESIAN_GRID( myThid )
11    
12  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
13  C     *==========================================================*  C     *==========================================================*
14  C     | SUBROUTINE INI_CARTESIAN_GRID  C     | SUBROUTINE INI_CARTESIAN_GRID
# Line 48  C     === Global variables === Line 50  C     === Global variables ===
50  #include "EEPARAMS.h"  #include "EEPARAMS.h"
51  #include "PARAMS.h"  #include "PARAMS.h"
52  #include "GRID.h"  #include "GRID.h"
53    c#ifdef ALLOW_EXCH2
54    c#include "W2_EXCH2_SIZE.h"
55    c#include "W2_EXCH2_TOPOLOGY.h"
56    c#include "W2_EXCH2_PARAMS.h"
57    c#endif /* ALLOW_EXCH2 */
58    
59  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
60  C     == Routine arguments ==  C     == Routine arguments ==
# Line 56  C     myThid  ::  Number of this instanc Line 63  C     myThid  ::  Number of this instanc
63    
64  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
65  C     == Local variables ==  C     == Local variables ==
66        INTEGER iG, jG, bi, bj, I,  J        INTEGER iG, jG, bi, bj, i, j
67        _RL xG0, yG0        _RL xG0, yG0
68  C     "Long" real for temporary coordinate calculation  C     "Long" real for temporary coordinate calculation
69  C      NOTICE the extended range of indices!!  C      NOTICE the extended range of indices!!
# Line 65  C      NOTICE the extended range of indi Line 72  C      NOTICE the extended range of indi
72  C     These functions return the "global" index with valid values beyond  C     These functions return the "global" index with valid values beyond
73  C     halo regions  C     halo regions
74        INTEGER iGl,jGl        INTEGER iGl,jGl
75        iGl(I,bi) = 1+mod(myXGlobalLo-1+(bi-1)*sNx+I+Olx*Nx-1,Nx)        iGl(i,bi) = 1+MOD(myXGlobalLo-1+(bi-1)*sNx+i+Olx*Nx-1,Nx)
76        jGl(J,bj) = 1+mod(myYGlobalLo-1+(bj-1)*sNy+J+Oly*Ny-1,Ny)        jGl(j,bj) = 1+MOD(myYGlobalLo-1+(bj-1)*sNy+j+Oly*Ny-1,Ny)
77    c#ifdef ALLOW_EXCH2
78    c      INTEGER tN
79    c#endif /* ALLOW_EXCH2 */
80  CEOP  CEOP
81    
82  C     For each tile ...  C     For each tile ...
# Line 76  C     For each tile ... Line 86  C     For each tile ...
86  C--     "Global" index (place holder)  C--     "Global" index (place holder)
87          jG = myYGlobalLo + (bj-1)*sNy          jG = myYGlobalLo + (bj-1)*sNy
88          iG = myXGlobalLo + (bi-1)*sNx          iG = myXGlobalLo + (bi-1)*sNx
89    c#ifdef ALLOW_EXCH2
90    c        IF ( W2_useE2ioLayOut ) THEN
91    cC- note: does not work for non-uniform delX or delY
92    c          tN = W2_myTileList(bi,bj)
93    c          iG = exch2_txGlobalo(tN)
94    c          jG = exch2_tyGlobalo(tN)
95    c        ENDIF
96    c#endif /* ALLOW_EXCH2 */
97    
98  C--   First find coordinate of tile corner (meaning outer corner of halo)  C--   First find coordinate of tile corner (meaning outer corner of halo)
99          xG0 = xgOrigin          xG0 = xgOrigin
# Line 85  C       Find the X-coordinate of the out Line 103  C       Find the X-coordinate of the out
103          ENDDO          ENDDO
104  C       Back-step to the outer grid-line of the "halo" region  C       Back-step to the outer grid-line of the "halo" region
105          DO i=1, Olx          DO i=1, Olx
106           xG0 = xG0 - delX( 1+mod(Olx*Nx-1+iG-i,Nx) )           xG0 = xG0 - delX( 1+MOD(Olx*Nx-1+iG-i,Nx) )
107          ENDDO          ENDDO
108  C       Find the Y-coordinate of the outer grid-line of the "real" tile  C       Find the Y-coordinate of the outer grid-line of the "real" tile
109          yG0 = ygOrigin          yG0 = ygOrigin
# Line 94  C       Find the Y-coordinate of the out Line 112  C       Find the Y-coordinate of the out
112          ENDDO          ENDDO
113  C       Back-step to the outer grid-line of the "halo" region  C       Back-step to the outer grid-line of the "halo" region
114          DO j=1, Oly          DO j=1, Oly
115           yG0 = yG0 - delY( 1+mod(Oly*Ny-1+jG-j,Ny) )           yG0 = yG0 - delY( 1+MOD(Oly*Ny-1+jG-j,Ny) )
116          ENDDO          ENDDO
117    
118  C--     Calculate coordinates of cell corners for N+1 grid-lines  C--     Calculate coordinates of cell corners for N+1 grid-lines
119          DO J=1-Oly,sNy+Oly +1          DO j=1-Oly,sNy+Oly +1
120           xGloc(1-Olx,J) = xG0           xGloc(1-Olx,j) = xG0
121           DO I=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
122  c         xGloc(I+1,J) = xGloc(I,J) + delX(1+mod(Nx-1+iG-1+i,Nx))  c         xGloc(i+1,j) = xGloc(i,j) + delX(1+mod(Nx-1+iG-1+i,Nx))
123            xGloc(I+1,J) = xGloc(I,J) + delX( iGl(I,bi) )            xGloc(i+1,j) = xGloc(i,j) + delX( iGl(i,bi) )
124           ENDDO           ENDDO
125          ENDDO          ENDDO
126          DO I=1-Olx,sNx+Olx +1          DO i=1-Olx,sNx+Olx +1
127           yGloc(I,1-Oly) = yG0           yGloc(i,1-Oly) = yG0
128           DO J=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
129  c         yGloc(I,J+1) = yGloc(I,J) + delY(1+mod(Ny-1+jG-1+j,Ny))  c         yGloc(i,j+1) = yGloc(i,j) + delY(1+mod(Ny-1+jG-1+j,Ny))
130            yGloc(I,J+1) = yGloc(I,J) + delY( jGl(J,bj) )            yGloc(i,j+1) = yGloc(i,j) + delY( jGl(j,bj) )
131           ENDDO           ENDDO
132          ENDDO          ENDDO
133    
134  C--     Make a permanent copy of [xGloc,yGloc] in [xG,yG]  C--     Make a permanent copy of [xGloc,yGloc] in [xG,yG]
135          DO J=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
136           DO I=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
137            xG(I,J,bi,bj) = xGloc(I,J)            xG(i,j,bi,bj) = xGloc(i,j)
138            yG(I,J,bi,bj) = yGloc(I,J)            yG(i,j,bi,bj) = yGloc(i,j)
139           ENDDO           ENDDO
140          ENDDO          ENDDO
141    
142  C--     Calculate [xC,yC], coordinates of cell centers  C--     Calculate [xC,yC], coordinates of cell centers
143          DO J=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
144           DO I=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
145  C         by averaging  C         by averaging
146            xC(I,J,bi,bj) = 0.25*(            xC(i,j,bi,bj) = 0.25 _d 0*(
147       &     xGloc(I,J)+xGloc(I+1,J)+xGloc(I,J+1)+xGloc(I+1,J+1) )       &     xGloc(i,j)+xGloc(i+1,j)+xGloc(i,j+1)+xGloc(i+1,j+1) )
148            yC(I,J,bi,bj) = 0.25*(            yC(i,j,bi,bj) = 0.25 _d 0*(
149       &     yGloc(I,J)+yGloc(I+1,J)+yGloc(I,J+1)+yGloc(I+1,J+1) )       &     yGloc(i,j)+yGloc(i+1,j)+yGloc(i,j+1)+yGloc(i+1,j+1) )
150           ENDDO           ENDDO
151          ENDDO          ENDDO
152    
153  C--     Calculate [dxF,dyF], lengths between cell faces (through center)  C--     Calculate [dxF,dyF], lengths between cell faces (through center)
154          DO J=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
155           DO I=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
156            dXF(I,J,bi,bj) = delX( iGl(I,bi) )            dxF(i,j,bi,bj) = delX( iGl(i,bi) )
157            dYF(I,J,bi,bj) = delY( jGl(J,bj) )            dyF(i,j,bi,bj) = delY( jGl(j,bj) )
158           ENDDO           ENDDO
159          ENDDO          ENDDO
160    
161  C--     Calculate [dxG,dyG], lengths along cell boundaries  C--     Calculate [dxG,dyG], lengths along cell boundaries
162          DO J=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
163           DO I=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
164            dXG(I,J,bi,bj) = delX( iGl(I,bi) )            dxG(i,j,bi,bj) = delX( iGl(i,bi) )
165            dYG(I,J,bi,bj) = delY( jGl(J,bj) )            dyG(i,j,bi,bj) = delY( jGl(j,bj) )
166           ENDDO           ENDDO
167          ENDDO          ENDDO
168    
169  C--     The following arrays are not defined in some parts of the halo  C--     The following arrays are not defined in some parts of the halo
170  C       region. We set them to zero here for safety. If they are ever  C       region. We set them to zero here for safety. If they are ever
171  C       referred to, especially in the denominator then it is a mistake!  C       referred to, especially in the denominator then it is a mistake!
172          DO J=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
173           DO I=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
174            dXC(I,J,bi,bj) = 0.            dxC(i,j,bi,bj) = 0.
175            dYC(I,J,bi,bj) = 0.            dyC(i,j,bi,bj) = 0.
176            dXV(I,J,bi,bj) = 0.            dxV(i,j,bi,bj) = 0.
177            dYU(I,J,bi,bj) = 0.            dyU(i,j,bi,bj) = 0.
178            rAw(I,J,bi,bj) = 0.            rAw(i,j,bi,bj) = 0.
179            rAs(I,J,bi,bj) = 0.            rAs(i,j,bi,bj) = 0.
180           ENDDO           ENDDO
181          ENDDO          ENDDO
182    
183  C--     Calculate [dxC], zonal length between cell centers  C--     Calculate [dxC], zonal length between cell centers
184          DO J=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
185           DO I=1-Olx+1,sNx+Olx ! NOTE range           DO i=1-Olx+1,sNx+Olx ! NOTE range
186            dXC(I,J,bi,bj) = 0.5D0*(dXF(I,J,bi,bj)+dXF(I-1,J,bi,bj))            dxC(i,j,bi,bj) = 0.5 _d 0*(dxF(i,j,bi,bj)+dxF(i-1,j,bi,bj))
187           ENDDO           ENDDO
188          ENDDO          ENDDO
189    
190  C--     Calculate [dyC], meridional length between cell centers  C--     Calculate [dyC], meridional length between cell centers
191          DO J=1-Oly+1,sNy+Oly ! NOTE range          DO j=1-Oly+1,sNy+Oly ! NOTE range
192           DO I=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
193            dYC(I,J,bi,bj) = 0.5*(dYF(I,J,bi,bj)+dYF(I,J-1,bi,bj))            dyC(i,j,bi,bj) = 0.5 _d 0*(dyF(i,j,bi,bj)+dyF(i,j-1,bi,bj))
194           ENDDO           ENDDO
195          ENDDO          ENDDO
196    
197  C--     Calculate [dxV,dyU], length between velocity points (through corners)  C--     Calculate [dxV,dyU], length between velocity points (through corners)
198          DO J=1-Oly+1,sNy+Oly ! NOTE range          DO j=1-Oly+1,sNy+Oly ! NOTE range
199           DO I=1-Olx+1,sNx+Olx ! NOTE range           DO i=1-Olx+1,sNx+Olx ! NOTE range
200  C         by averaging (method I)  C         by averaging (method I)
201            dXV(I,J,bi,bj) = 0.5*(dXG(I,J,bi,bj)+dXG(I-1,J,bi,bj))            dxV(i,j,bi,bj) = 0.5 _d 0*(dxG(i,j,bi,bj)+dxG(i-1,j,bi,bj))
202            dYU(I,J,bi,bj) = 0.5*(dYG(I,J,bi,bj)+dYG(I,J-1,bi,bj))            dyU(i,j,bi,bj) = 0.5 _d 0*(dyG(i,j,bi,bj)+dyG(i,j-1,bi,bj))
203  C         by averaging (method II)  C         by averaging (method II)
204  c         dXV(I,J,bi,bj) = 0.5*(dXG(I,J,bi,bj)+dXG(I-1,J,bi,bj))  c         dxV(i,j,bi,bj) = 0.5*(dxG(i,j,bi,bj)+dxG(i-1,j,bi,bj))
205  c         dYU(I,J,bi,bj) = 0.5*(dYC(I,J,bi,bj)+dYC(I-1,J,bi,bj))  c         dyU(i,j,bi,bj) = 0.5*(dyC(i,j,bi,bj)+dyC(i-1,j,bi,bj))
206           ENDDO           ENDDO
207          ENDDO          ENDDO
208    
209  C--     Calculate vertical face area  C--     Calculate vertical face area
210          DO J=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
211           DO I=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
212            rA (I,J,bi,bj) = dxF(I,J,bi,bj)*dyF(I,J,bi,bj)            rA (i,j,bi,bj) = dxF(i,j,bi,bj)*dyF(i,j,bi,bj)
213            rAw(I,J,bi,bj) = dxC(I,J,bi,bj)*dyG(I,J,bi,bj)            rAw(i,j,bi,bj) = dxC(i,j,bi,bj)*dyG(i,j,bi,bj)
214            rAs(I,J,bi,bj) = dxG(I,J,bi,bj)*dyC(I,J,bi,bj)            rAs(i,j,bi,bj) = dxG(i,j,bi,bj)*dyC(i,j,bi,bj)
215            rAz(I,J,bi,bj) = dxV(I,J,bi,bj)*dyU(I,J,bi,bj)            rAz(i,j,bi,bj) = dxV(i,j,bi,bj)*dyU(i,j,bi,bj)
216  C--     Set trigonometric terms & grid orientation:  C--     Set trigonometric terms & grid orientation:
217            tanPhiAtU(I,J,bi,bj) = 0.            tanPhiAtU(i,j,bi,bj) = 0.
218            tanPhiAtV(I,J,bi,bj) = 0.            tanPhiAtV(i,j,bi,bj) = 0.
219            angleCosC(I,J,bi,bj) = 1.            angleCosC(i,j,bi,bj) = 1.
220            angleSinC(I,J,bi,bj) = 0.            angleSinC(i,j,bi,bj) = 0.
221           ENDDO           ENDDO
222          ENDDO          ENDDO
223    
224  C--     Cosine(lat) scaling  C--     Cosine(lat) scaling
225          DO J=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
226           cosFacU(J,bi,bj)=1.           cosFacU(j,bi,bj)=1.
227           cosFacV(J,bi,bj)=1.           cosFacV(j,bi,bj)=1.
228           sqcosFacU(J,bi,bj)=1.           sqcosFacU(j,bi,bj)=1.
229           sqcosFacV(J,bi,bj)=1.           sqcosFacV(j,bi,bj)=1.
230          ENDDO          ENDDO
231    
232  C--   end bi,bj loops  C--   end bi,bj loops
Legend:
Removed from v.1.21  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.22