/[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.15 by cnh, Sun Feb 4 14:38:47 2001 UTC revision 1.16 by adcroft, Tue May 29 14:01:37 2001 UTC
# Line 48  C     myThid -  Number of this instance Line 48  C     myThid -  Number of this instance
48  CEndOfInterface  CEndOfInterface
49    
50  C     == Local variables ==  C     == Local variables ==
51  C     xG, yG - Global coordinate location.        INTEGER iG, jG, bi, bj, I,  J
52  C     xBase  - South-west corner location for process.        _RL xG0, yG0
53  C     yBase  
54  C     zUpper - Work arrays for upper and lower  C "Long" real for temporary coordinate calculation
55  C     zLower   cell-face heights.  C  NOTICE the extended range of indices!!
56  C     phi    - Temporary scalar        _RL xGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)
57  C     xBase  - Temporaries for lower corner coordinate        _RL yGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)
58  C     yBase  
59  C     iG, jG - Global coordinate index. Usually used to hold  C These functions return the "global" index with valid values beyond
60  C              the south-west global coordinate of a tile.  C halo regions
61  C     bi,bj  - Loop counters        INTEGER iGl,jGl
62  C     zUpper - Temporary arrays holding z coordinates of        iGl(I,bi) = 1+mod(myXGlobalLo-1+(bi-1)*sNx+I+Olx*Nx-1,Nx)
63  C     zLower   upper and lower faces.        jGl(J,bj) = 1+mod(myYGlobalLo-1+(bj-1)*sNy+J+Oly*Ny-1,Ny)
64  C     I,J,K  
65        _RL    xGloc, yGloc  C     For each tile ...
       _RL    xBase, yBase  
       INTEGER iG, jG  
       INTEGER bi, bj  
       INTEGER  I,  J  
   
 C--   Simple example of inialisation on cartesian grid  
 C--   First set coordinates of cell centers  
 C     This operation is only performed at start up so for more  
 C     complex configurations it is usually OK to pass iG, jG to a custom  
 C     function and have it return xG and yG.  
 C     Set up my local grid first  
       xC0 = 0. _d 0  
       yC0 = 0. _d 0  
66        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
        jG = myYGlobalLo + (bj-1)*sNy  
67         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
68    
69    C--     "Global" index (place holder)
70            jG = myYGlobalLo + (bj-1)*sNy
71          iG = myXGlobalLo + (bi-1)*sNx          iG = myXGlobalLo + (bi-1)*sNx
         yBase = 0. _d 0  
         xBase = 0. _d 0  
         DO i=1,iG-1  
          xBase = xBase + delX(i)  
         ENDDO  
         DO j=1,jG-1  
          yBase = yBase + delY(j)  
         ENDDO  
         yGloc = yBase  
         DO J=1,sNy  
          xGloc = xBase  
          DO I=1,sNx  
           xG(I,J,bi,bj)  = xGloc  
           yG(I,J,bi,bj)  = yGloc  
           xc(I,J,bi,bj)  = xGloc + delX(iG+i-1)*0.5 _d 0  
           yc(I,J,bi,bj)  = yGloc + delY(jG+j-1)*0.5 _d 0  
           xGloc = xGloc + delX(iG+I-1)  
           dxF(I,J,bi,bj) = delX(iG+i-1)  
           dyF(I,J,bi,bj) = delY(jG+j-1)  
          ENDDO  
          yGloc = yGloc + delY(jG+J-1)  
         ENDDO  
        ENDDO  
       ENDDO  
 C     Now sync. and get edge regions from other threads and/or processes.  
 C     Note: We could just set the overlap regions ourselves here but  
 C           exchanging edges is safer and is good practice!  
       _EXCH_XY_R4( xc, myThid )  
       _EXCH_XY_R4( yc, myThid )  
       _EXCH_XY_R4(dxF, myThid )  
       _EXCH_XY_R4(dyF, myThid )  
72    
73  C--   Calculate separation between other points  C--   First find coordinate of tile corner (meaning outer corner of halo)
74  C     dxG, dyG are separations between cell corners along cell faces.          xG0 = 0.
75        DO bj = myByLo(myThid), myByHi(myThid)  C       Find the X-coordinate of the outer grid-line of the "real" tile
76         DO bi = myBxLo(myThid), myBxHi(myThid)          DO i=1, iG-1
77          DO J=1,sNy           xG0 = xG0 + delX(i)
78           DO I=1,sNx          ENDDO
79            dxG(I,J,bi,bj) = (dxF(I,J,bi,bj)+dxF(I,J-1,bi,bj))*0.5 _d 0  C       Back-step to the outer grid-line of the "halo" region
80            dyG(I,J,bi,bj) = (dyF(I,J,bi,bj)+dyF(I-1,J,bi,bj))*0.5 _d 0          DO i=1, Olx
81             xG0 = xG0 - delX( 1+mod(Olx*Nx-1+iG-i,Nx) )
82            ENDDO
83    C       Find the Y-coordinate of the outer grid-line of the "real" tile
84            yG0 = 0.
85            DO j=1, jG-1
86             yG0 = yG0 + delY(j)
87            ENDDO
88    C       Back-step to the outer grid-line of the "halo" region
89            DO j=1, Oly
90             yG0 = yG0 - delY( 1+mod(Oly*Ny-1+jG-j,Ny) )
91            ENDDO
92    
93    C--     Calculate coordinates of cell corners for N+1 grid-lines
94            DO J=1-Oly,sNy+Oly +1
95             xGloc(1-Olx,J) = xG0
96             DO I=1-Olx,sNx+Olx
97    c         xGloc(I+1,J) = xGloc(I,J) + delX(1+mod(Nx-1+iG-1+i,Nx))
98              xGloc(I+1,J) = xGloc(I,J) + delX( iGl(I,bi) )
99           ENDDO           ENDDO
100          ENDDO          ENDDO
101         ENDDO          DO I=1-Olx,sNx+Olx +1
102        ENDDO           yGloc(I,1-Oly) = yG0
103        _EXCH_XY_R4(dxG, myThid )           DO J=1-Oly,sNy+Oly
104        _EXCH_XY_R4(dyG, myThid )  c         yGloc(I,J+1) = yGloc(I,J) + delY(1+mod(Ny-1+jG-1+j,Ny))
105  C     dxV, dyU are separations between velocity points along cell faces.            yGloc(I,J+1) = yGloc(I,J) + delY( jGl(J,bj) )
       DO bj = myByLo(myThid), myByHi(myThid)  
        DO bi = myBxLo(myThid), myBxHi(myThid)  
         DO J=1,sNy  
          DO I=1,sNx  
           dxV(I,J,bi,bj) = (dxG(I,J,bi,bj)+dxG(I-1,J,bi,bj))*0.5 _d 0  
           dyU(I,J,bi,bj) = (dyG(I,J,bi,bj)+dyG(I,J-1,bi,bj))*0.5 _d 0  
106           ENDDO           ENDDO
107          ENDDO          ENDDO
108         ENDDO  
109        ENDDO  C--     Make a permanent copy of [xGloc,yGloc] in [xG,yG]
110        _EXCH_XY_R4(dxV, myThid )          DO J=1-Oly,sNy+Oly
111        _EXCH_XY_R4(dyU, myThid )           DO I=1-Olx,sNx+Olx
112  C     dxC, dyC is separation between cell centers            xG(I,J,bi,bj) = xGloc(I,J)
113        DO bj = myByLo(myThid), myByHi(myThid)            yG(I,J,bi,bj) = yGloc(I,J)
114         DO bi = myBxLo(myThid), myBxHi(myThid)           ENDDO
115          DO J=1,sNy          ENDDO
116           DO I=1,sNx  
117            dxC(I,J,bi,bj)    = (dxF(I,J,bi,bj)+dxF(I-1,J,bi,bj))*0.5 _d 0  C--     Calculate [xC,yC], coordinates of cell centers
118            dyC(I,J,bi,bj)    = (dyF(I,J,bi,bj)+dyF(I,J-1,bi,bj))*0.5 _d 0          DO J=1-Oly,sNy+Oly
119             DO I=1-Olx,sNx+Olx
120    C         by averaging
121              xC(I,J,bi,bj) = 0.25*(
122         &     xGloc(I,J)+xGloc(I+1,J)+xGloc(I,J+1)+xGloc(I+1,J+1) )
123              yC(I,J,bi,bj) = 0.25*(
124         &     yGloc(I,J)+yGloc(I+1,J)+yGloc(I,J+1)+yGloc(I+1,J+1) )
125             ENDDO
126            ENDDO
127    
128    C--     Calculate [dxF,dyF], lengths between cell faces (through center)
129            DO J=1-Oly,sNy+Oly
130             DO I=1-Olx,sNx+Olx
131              dXF(I,J,bi,bj) = delX( iGl(I,bi) )
132              dYF(I,J,bi,bj) = delY( jGl(J,bj) )
133             ENDDO
134            ENDDO
135    
136    C--     Calculate [dxG,dyG], lengths along cell boundaries
137            DO J=1-Oly,sNy+Oly
138             DO I=1-Olx,sNx+Olx
139              dXG(I,J,bi,bj) = delX( iGl(I,bi) )
140              dYG(I,J,bi,bj) = delY( jGl(J,bj) )
141             ENDDO
142            ENDDO
143    
144    C--     The following arrays are not defined in some parts of the halo
145    C       region. We set them to zero here for safety. If they are ever
146    C       referred to, especially in the denominator then it is a mistake!
147            DO J=1-Oly,sNy+Oly
148             DO I=1-Olx,sNx+Olx
149              dXC(I,J,bi,bj) = 0.
150              dYC(I,J,bi,bj) = 0.
151              dXV(I,J,bi,bj) = 0.
152              dYU(I,J,bi,bj) = 0.
153              rAw(I,J,bi,bj) = 0.
154              rAs(I,J,bi,bj) = 0.
155             ENDDO
156            ENDDO
157    
158    C--     Calculate [dxC], zonal length between cell centers
159            DO J=1-Oly,sNy+Oly
160             DO I=1-Olx+1,sNx+Olx ! NOTE range
161              dXC(I,J,bi,bj) = 0.5D0*(dXF(I,J,bi,bj)+dXF(I-1,J,bi,bj))
162             ENDDO
163            ENDDO
164    
165    C--     Calculate [dyC], meridional length between cell centers
166            DO J=1-Oly+1,sNy+Oly ! NOTE range
167             DO I=1-Olx,sNx+Olx
168              dYC(I,J,bi,bj) = 0.5*(dYF(I,J,bi,bj)+dYF(I,J-1,bi,bj))
169           ENDDO           ENDDO
170          ENDDO          ENDDO
171         ENDDO  
172        ENDDO  C--     Calculate [dxV,dyU], length between velocity points (through corners)
173        _EXCH_XY_R4(dxC, myThid )          DO J=1-Oly+1,sNy+Oly ! NOTE range
174        _EXCH_XY_R4(dyC, myThid )           DO I=1-Olx+1,sNx+Olx ! NOTE range
175    C         by averaging (method I)
176              dXV(I,J,bi,bj) = 0.5*(dXG(I,J,bi,bj)+dXG(I-1,J,bi,bj))
177              dYU(I,J,bi,bj) = 0.5*(dYG(I,J,bi,bj)+dYG(I,J-1,bi,bj))
178    C         by averaging (method II)
179    c         dXV(I,J,bi,bj) = 0.5*(dXG(I,J,bi,bj)+dXG(I-1,J,bi,bj))
180    c         dYU(I,J,bi,bj) = 0.5*(dYC(I,J,bi,bj)+dYC(I-1,J,bi,bj))
181             ENDDO
182            ENDDO
183    
184  C     Calculate vertical face area  C     Calculate vertical face area
185        DO bj = myByLo(myThid), myByHi(myThid)          DO J=1-Oly,sNy+Oly
186         DO bi = myBxLo(myThid), myBxHi(myThid)           DO I=1-Olx,sNx+Olx
         DO J=1,sNy  
          DO I=1,sNx  
187            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)
188            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)
189            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)
190            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)
191            tanPhiAtU(I,J,bi,bj) = 0. _d 0            tanPhiAtU(I,J,bi,bj) = 0.
192            tanPhiAtV(I,J,bi,bj) = 0. _d 0            tanPhiAtV(I,J,bi,bj) = 0.
193           ENDDO           ENDDO
194          ENDDO          ENDDO
        ENDDO  
       ENDDO  
       _EXCH_XY_R4 (rA       , myThid )  
       _EXCH_XY_R4 (rAw      , myThid )  
       _EXCH_XY_R4 (rAs      , myThid )  
       _EXCH_XY_R4 (tanPhiAtU , myThid )  
       _EXCH_XY_R4 (tanPhiAtV , myThid )  
195    
196  C  C--     Cosine(lat) scaling
197            DO J=1-OLy,sNy+OLy
198             cosFacU(J,bi,bj)=1.
199             cosFacV(J,bi,bj)=1.
200             sqcosFacU(J,bi,bj)=1.
201             sqcosFacV(J,bi,bj)=1.
202            ENDDO
203    
204           ENDDO ! bi
205          ENDDO ! bj
206    
207        RETURN        RETURN
208        END        END
Legend:
Removed from v.1.15  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.22