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

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

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

revision 1.14 by adcroft, Mon Mar 27 22:25:44 2000 UTC revision 1.15 by adcroft, Fri Feb 2 21:04:48 2001 UTC
# Line 2  C $Header$ Line 2  C $Header$
2    
3  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
4    
5    #undef  USE_BACKWARD_COMPATIBLE_GRID
6    
7  CStartOfInterface  CStartOfInterface
8        SUBROUTINE INI_SPHERICAL_POLAR_GRID( myThid )        SUBROUTINE INI_SPHERICAL_POLAR_GRID( myThid )
9  C     /==========================================================\  C     /==========================================================\
# Line 50  C     == Local variables == Line 52  C     == Local variables ==
52  C     xG, yG - Global coordinate location.  C     xG, yG - Global coordinate location.
53  C     xBase  - South-west corner location for process.  C     xBase  - South-west corner location for process.
54  C     yBase  C     yBase
55    C     zUpper - Work arrays for upper and lower
56    C     zLower   cell-face heights.
57    C     phi    - Temporary scalar
58  C     iG, jG - Global coordinate index. Usually used to hold  C     iG, jG - Global coordinate index. Usually used to hold
59  C              the south-west global coordinate of a tile.  C              the south-west global coordinate of a tile.
60  C     bi,bj  - Loop counters  C     bi,bj  - Loop counters
# Line 60  C     yBase Line 65  C     yBase
65  C     lat, latN, - Temporary variables used to hold latitude  C     lat, latN, - Temporary variables used to hold latitude
66  C     latS         values.  C     latS         values.
67  C     I,J,K  C     I,J,K
       _RL    xG, yG  
       _RL    xBase, yBase  
68        INTEGER iG, jG        INTEGER iG, jG
69        INTEGER bi, bj        INTEGER bi, bj
70        INTEGER  I,  J        INTEGER  I,  J
71        _RL lat, latS, latN        _RL lat, dlat, dlon, xG0, yG0
72    
73    
74    C "Long" real for temporary coordinate calculation
75    C  NOTICE the extended range of indices!!
76          _RL xGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)
77          _RL yGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)
78    
79    C These functions return the "global" index with valid values beyond
80    C halo regions
81          INTEGER iGl,jGl
82          iGl(I,bi) = 1+mod(myXGlobalLo-1+(bi-1)*sNx+I+Olx*Nx-1,Nx)
83          jGl(J,bj) = 1+mod(myYGlobalLo-1+(bj-1)*sNy+J+Oly*Ny-1,Ny)
84    
85  C--   Example of inialisation for spherical polar grid  C     For each tile ...
 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  
 C     Note: In the spherical polar case delX and delY are given in  
 C           degrees and are relative to some starting latitude and  
 C           longitude - phiMin and thetaMin.  
       xC0 = thetaMin  
       yC0 = phiMin  
86        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
        jG = myYGlobalLo + (bj-1)*sNy  
87         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
88    
89    C--     "Global" index (place holder)
90            jG = myYGlobalLo + (bj-1)*sNy
91          iG = myXGlobalLo + (bi-1)*sNx          iG = myXGlobalLo + (bi-1)*sNx
         yBase = phiMin  
         xBase = thetaMin  
         DO i=1,iG-1  
          xBase = xBase + delX(i)  
         ENDDO  
         DO j=1,jG-1  
          yBase = yBase + delY(j)  
         ENDDO  
         yG = yBase  
         DO J=1,sNy  
          xG = xBase  
          DO I=1,sNx  
           xc(I,J,bi,bj)  = xG + delX(iG+i-1)*0.5 _d 0  
           yc(I,J,bi,bj)  = yG + delY(jG+j-1)*0.5 _d 0  
           xG = xG + delX(iG+I-1)  
           dxF(I,J,bi,bj) = delX(iG+i-1)*deg2rad  
      &    *rSphere*COS(yc(I,J,bi,bj)*deg2rad)  
           dyF(I,J,bi,bj) = delY(jG+j-1)*deg2rad*rSphere  
          ENDDO  
          yG = yG + 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 )  
92    
93  C--   Calculate separation between other points  C--   First find coordinate of tile corner (meaning outer corner of halo)
94  C     dxG, dyG are separations between cell corners along cell faces.          xG0 = thetaMin
95        DO bj = myByLo(myThid), myByHi(myThid)  C       Find the X-coordinate of the outer grid-line of the "real" tile
96         DO bi = myBxLo(myThid), myBxHi(myThid)          DO i=1, iG-1
97          DO J=1,sNy           xG0 = xG0 + delX(i)
98           DO I=1,sNx          ENDDO
99            jG = myYGlobalLo + (bj-1)*sNy + J-1  C       Back-step to the outer grid-line of the "halo" region
100            iG = myXGlobalLo + (bi-1)*sNx + I-1          DO i=1, Olx
101            lat = yc(I,J,bi,bj)-delY(jG) * 0.5 _d 0           xG0 = xG0 - delX( 1+mod(Olx*Nx-1+iG-i,Nx) )
102            dxG(I,J,bi,bj) = rSphere*COS(lat*deg2rad)*delX(iG)*deg2rad          ENDDO
103            dyG(I,J,bi,bj) = (dyF(I,J,bi,bj)+dyF(I-1,J,bi,bj))*0.5 _d 0  C       Find the Y-coordinate of the outer grid-line of the "real" tile
104            yG0 = phiMin
105            DO j=1, jG-1
106             yG0 = yG0 + delY(j)
107            ENDDO
108    C       Back-step to the outer grid-line of the "halo" region
109            DO j=1, Oly
110             yG0 = yG0 - delY( 1+mod(Oly*Ny-1+jG-j,Ny) )
111            ENDDO
112    
113    C--     Calculate coordinates of cell corners for N+1 grid-lines
114            DO J=1-Oly,sNy+Oly +1
115             xGloc(1-Olx,J) = xG0
116             DO I=1-Olx,sNx+Olx
117    c         xGloc(I+1,J) = xGloc(I,J) + delX(1+mod(Nx-1+iG-1+i,Nx))
118              xGloc(I+1,J) = xGloc(I,J) + delX( iGl(I,bi) )
119           ENDDO           ENDDO
120          ENDDO          ENDDO
121         ENDDO          DO I=1-Olx,sNx+Olx +1
122        ENDDO           yGloc(I,1-Oly) = yG0
123        _EXCH_XY_R4(dxG, myThid )           DO J=1-Oly,sNy+Oly
124        _EXCH_XY_R4(dyG, myThid )  c         yGloc(I,J+1) = yGloc(I,J) + delY(1+mod(Ny-1+jG-1+j,Ny))
125  C     dxC, dyC is separation between cell centers            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  
           dxC(I,J,bi,bj)    = (dxF(I,J,bi,bj)+dxF(I-1,J,bi,bj))*0.5 _d 0  
           dyC(I,J,bi,bj)    = (dyF(I,J,bi,bj)+dyF(I,J-1,bi,bj))*0.5 _d 0  
126           ENDDO           ENDDO
127          ENDDO          ENDDO
128         ENDDO  
129        ENDDO  C--     Make a permanent copy of [xGloc,yGloc] in [xG,yG]
130        _EXCH_XY_R4(dxC, myThid )          DO J=1-Oly,sNy+Oly
131        _EXCH_XY_R4(dyC, myThid )           DO I=1-Olx,sNx+Olx
132  C     dxV, dyU are separations between velocity points along cell faces.            xG(I,J,bi,bj) = xGloc(I,J)
133        DO bj = myByLo(myThid), myByHi(myThid)            yG(I,J,bi,bj) = yGloc(I,J)
        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  
 #ifdef OLD_UV_GEOMETRY  
           dyU(I,J,bi,bj) = (dyG(I,J,bi,bj)+dyG(I,J-1,bi,bj))*0.5 _d 0  
 #else  
           dyU(I,J,bi,bj) = (dyC(I,J,bi,bj)+dyC(I-1,J,bi,bj))*0.5 _d 0  
 #endif  
134           ENDDO           ENDDO
135          ENDDO          ENDDO
136         ENDDO  
137        ENDDO  C--     Calculate [xC,yC], coordinates of cell centers
138        _EXCH_XY_R4(dxV, myThid )          DO J=1-Oly,sNy+Oly
139        _EXCH_XY_R4(dyU, myThid )           DO I=1-Olx,sNx+Olx
140  C     Calculate vertical face area and trigonometric terms  C         by averaging
141        DO bj = myByLo(myThid), myByHi(myThid)            xC(I,J,bi,bj) = 0.25*(
142         DO bi = myBxLo(myThid), myBxHi(myThid)       &     xGloc(I,J)+xGloc(I+1,J)+xGloc(I,J+1)+xGloc(I+1,J+1) )
143          DO J=1,sNy            yC(I,J,bi,bj) = 0.25*(
144           DO I=1,sNx       &     yGloc(I,J)+yGloc(I+1,J)+yGloc(I,J+1)+yGloc(I+1,J+1) )
145            jG = myYGlobalLo + (bj-1)*sNy + J-1           ENDDO
146            iG = myXGlobalLo + (bi-1)*sNx + I-1          ENDDO
147            latS = yc(i,j,bi,bj)-delY(jG)*0.5 _d 0  
148            latN = yc(i,j,bi,bj)+delY(jG)*0.5 _d 0  C--     Calculate [dxF,dyF], lengths between cell faces (through center)
149  #ifdef OLD_UV_GEOMETRY          DO J=1-Oly,sNy+Oly
150             DO I=1-Olx,sNx+Olx
151    C         by averaging
152    c         dXF(I,J,bi,bj) = 0.5*(dXG(I,J,bi,bj)+dXG(I,J+1,bi,bj))
153    c         dYF(I,J,bi,bj) = 0.5*(dYG(I,J,bi,bj)+dYG(I+1,J,bi,bj))
154    C         by formula
155              lat = yC(I,J,bi,bj)
156              dlon = delX( iGl(I,bi) )
157              dlat = delY( jGl(J,bj) )
158              dXF(I,J,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
159    #ifdef    USE_BACKWARD_COMPATIBLE_GRID
160              dXF(I,J,bi,bj) = delX(iGl(I,bi))*deg2rad*rSphere*
161         &                     COS(yc(I,J,bi,bj)*deg2rad)
162    #endif    /* USE_BACKWARD_COMPATIBLE_GRID */
163              dYF(I,J,bi,bj) = rSphere*dlat*deg2rad
164             ENDDO
165            ENDDO
166    
167    C--     Calculate [dxG,dyG], lengths along cell boundaries
168            DO J=1-Oly,sNy+Oly
169             DO I=1-Olx,sNx+Olx
170    C         by averaging
171    c         dXG(I,J,bi,bj) = 0.5*(dXF(I,J,bi,bj)+dXF(I,J-1,bi,bj))
172    c         dYG(I,J,bi,bj) = 0.5*(dYF(I,J,bi,bj)+dYF(I-1,J,bi,bj))
173    C         by formula
174              lat = 0.5*(yGloc(I,J)+yGloc(I+1,J))
175              dlon = delX( iGl(I,bi) )
176              dlat = delY( jGl(J,bj) )
177              dXG(I,J,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
178              dYG(I,J,bi,bj) = rSphere*dlat*deg2rad
179             ENDDO
180            ENDDO
181    
182    C--     The following arrays are not defined in some parts of the halo
183    C       region. We set them to zero here for safety. If they are ever
184    C       referred to, especially in the denominator then it is a mistake!
185            DO J=1-Oly,sNy+Oly
186             DO I=1-Olx,sNx+Olx
187              dXC(I,J,bi,bj) = 0.
188              dYC(I,J,bi,bj) = 0.
189              dXV(I,J,bi,bj) = 0.
190              dYU(I,J,bi,bj) = 0.
191              rAw(I,J,bi,bj) = 0.
192              rAs(I,J,bi,bj) = 0.
193             ENDDO
194            ENDDO
195    
196    C--     Calculate [dxC], zonal length between cell centers
197            DO J=1-Oly,sNy+Oly
198             DO I=1-Olx+1,sNx+Olx ! NOTE range
199    C         by averaging
200              dXC(I,J,bi,bj) = 0.5D0*(dXF(I,J,bi,bj)+dXF(I-1,J,bi,bj))
201    C         by formula
202    c         lat = 0.5*(yC(I,J,bi,bj)+yC(I-1,J,bi,bj))
203    c         dlon = 0.5*(delX( iGl(I,bi) ) + delX( iGl(I-1,bi) ))
204    c         dXC(I,J,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
205    C         by difference
206    c         lat = 0.5*(yC(I,J,bi,bj)+yC(I-1,J,bi,bj))
207    c         dlon = (xC(I,J,bi,bj)+xC(I-1,J,bi,bj))
208    c         dXC(I,J,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
209             ENDDO
210            ENDDO
211    
212    C--     Calculate [dyC], meridional length between cell centers
213            DO J=1-Oly+1,sNy+Oly ! NOTE range
214             DO I=1-Olx,sNx+Olx
215    C         by averaging
216              dYC(I,J,bi,bj) = 0.5*(dYF(I,J,bi,bj)+dYF(I,J-1,bi,bj))
217    C         by formula
218    c         dlat = 0.5*(delY( jGl(J,bj) ) + delY( jGl(J-1,bj) ))
219    c         dYC(I,J,bi,bj) = rSphere*dlat*deg2rad
220    C         by difference
221    c         dlat = (yC(I,J,bi,bj)+yC(I,J-1,bi,bj))
222    c         dYC(I,J,bi,bj) = rSphere*dlat*deg2rad
223             ENDDO
224            ENDDO
225    
226    C--     Calculate [dxV,dyU], length between velocity points (through corners)
227            DO J=1-Oly+1,sNy+Oly ! NOTE range
228             DO I=1-Olx+1,sNx+Olx ! NOTE range
229    C         by averaging (method I)
230              dXV(I,J,bi,bj) = 0.5*(dXG(I,J,bi,bj)+dXG(I-1,J,bi,bj))
231              dYU(I,J,bi,bj) = 0.5*(dYG(I,J,bi,bj)+dYG(I,J-1,bi,bj))
232    C         by averaging (method II)
233    c         dXV(I,J,bi,bj) = 0.5*(dXG(I,J,bi,bj)+dXG(I-1,J,bi,bj))
234    c         dYU(I,J,bi,bj) = 0.5*(dYC(I,J,bi,bj)+dYC(I-1,J,bi,bj))
235             ENDDO
236            ENDDO
237    
238    C--     Calculate vertical face area (tracer cells)
239            DO J=1-Oly,sNy+Oly
240             DO I=1-Olx,sNx+Olx
241              lat=0.5*(yGloc(I,J)+yGloc(I+1,J))
242              dlon=delX( iGl(I,bi) )
243              dlat=delY( jGl(J,bj) )
244              rA(I,J,bi,bj) = rSphere*rSphere*dlon*deg2rad
245         &        *abs( sin((lat+dlat)*deg2rad)-sin(lat*deg2rad) )
246    #ifdef    USE_BACKWARD_COMPATIBLE_GRID
247              lat=yC(I,J,bi,bj)-delY( jGl(J,bj) )*0.5 _d 0
248              lon=yC(I,J,bi,bj)+delY( jGl(J,bj) )*0.5 _d 0
249            rA(I,J,bi,bj) = dyF(I,J,bi,bj)            rA(I,J,bi,bj) = dyF(I,J,bi,bj)
250       &    *rSphere*(SIN(latN*deg2rad)-SIN(latS*deg2rad))       &    *rSphere*(SIN(lon*deg2rad)-SIN(lat*deg2rad))
251  #else  #endif    /* USE_BACKWARD_COMPATIBLE_GRID */
252            rA(I,J,bi,bj) = rSphere*delX(iG)*deg2rad           ENDDO
253       &    *rSphere*(SIN(latN*deg2rad)-SIN(latS*deg2rad))          ENDDO
254  #endif  
255  C         Area cannot be zero but sin can be if lat if < -90.  C--     Calculate vertical face area (u cells)
256            IF ( rA(I,J,bi,bj) .LT. 0. ) rA(I,J,bi,bj) = -rA(I,J,bi,bj)          DO J=1-Oly,sNy+Oly
257            tanPhiAtU(i,j,bi,bj)=tan(_yC(i,j,bi,bj)*deg2rad)           DO I=1-Olx+1,sNx+Olx ! NOTE range
258            tanPhiAtV(i,j,bi,bj)=tan(latS*deg2rad)  C         by averaging
          ENDDO  
         ENDDO  
        ENDDO  
       ENDDO  
       _EXCH_XY_R4 (rA       , myThid )  
       _EXCH_XY_R4 (tanPhiAtU , myThid )  
       _EXCH_XY_R4 (tanPhiAtV , myThid )  
       DO bj = myByLo(myThid), myByHi(myThid)  
        DO bi = myBxLo(myThid), myBxHi(myThid)  
         DO J=1,sNy  
          DO I=1,sNx  
           iG = myXGlobalLo + (bi-1)*sNx + I-1  
           jG = myYGlobalLo + (bj-1)*sNy + J-1  
           latS = yc(i,j-1,bi,bj)  
           latN = yc(i,j,bi,bj)  
 #ifdef OLD_UV_GEOMETRY  
           rAw(I,J,bi,bj) = 0.5*(rA(I,J,bi,bj)+rA(I-1,J,bi,bj))  
           rAs(I,J,bi,bj) = 0.5*(rA(I,J,bi,bj)+rA(I,J-1,bi,bj))  
 #else  
259            rAw(I,J,bi,bj) = 0.5*(rA(I,J,bi,bj)+rA(I-1,J,bi,bj))            rAw(I,J,bi,bj) = 0.5*(rA(I,J,bi,bj)+rA(I-1,J,bi,bj))
260            rAs(I,J,bi,bj) = rSphere*delX(iG)*deg2rad  C         by formula
261       &    *rSphere*(SIN(latN*deg2rad)-SIN(latS*deg2rad))  c         lat=yGloc(I,J)
262  #endif  c         dlon=0.5*( delX( iGl(I,bi) ) + delX( iGl(I-1,bi) ) )
263    c         dlat=delY( jGl(J,bj) )
264    c         rAw(I,J,bi,bj) = rSphere*rSphere*dlon*deg2rad
265    c    &        *abs( sin((lat+dlat)*deg2rad)-sin(lat*deg2rad) )
266             ENDDO
267            ENDDO
268    
269    C--     Calculate vertical face area (v cells)
270            DO J=1-Oly,sNy+Oly
271             DO I=1-Olx,sNx+Olx
272    C         by formula
273              lat=yC(I,J,bi,bj)
274              dlon=delX( iGl(I,bi) )
275              dlat=0.5*( delY( jGl(J,bj) ) + delY( jGl(J-1,bj) ) )
276              rAs(I,J,bi,bj) = rSphere*rSphere*dlon*deg2rad
277         &        *abs( sin(lat*deg2rad)-sin((lat-dlat)*deg2rad) )
278    #ifdef    USE_BACKWARD_COMPATIBLE_GRID
279              lon=yC(I,J,bi,bj)-delY( jGl(J,bj) )
280              lat=yC(I,J,bi,bj)
281              rAs(I,J,bi,bj) = rSphere*delX(iGl(I,bi))*deg2rad
282         &    *rSphere*(SIN(lat*deg2rad)-SIN(lon*deg2rad))
283    #endif    /* USE_BACKWARD_COMPATIBLE_GRID */
284              IF (abs(lat).GT.90..OR.abs(lat-dlat).GT.90.) rAs(I,J,bi,bj)=0.
285             ENDDO
286            ENDDO
287    
288    C--     Calculate vertical face area (vorticity points)
289            DO J=1-Oly,sNy+Oly
290             DO I=1-Olx,sNx+Olx
291    C         by formula
292              lat=yC(I,J,bi,bj)
293              dlon=delX( iGl(I,bi) )
294              dlat=0.5*( delY( jGl(J,bj) ) + delY( jGl(J-1,bj) ) )
295              rAz(I,J,bi,bj) = rSphere*rSphere*dlon*deg2rad
296         &     *abs( sin(lat*deg2rad)-sin((lat-dlat)*deg2rad) )
297              IF (abs(lat).GT.90..OR.abs(lat-dlat).GT.90.) rAz(I,J,bi,bj)=0.
298             ENDDO
299            ENDDO
300    
301    C--     Calculate trigonometric terms
302            DO J=1-Oly,sNy+Oly
303             DO I=1-Olx,sNx+Olx
304              lat=0.5*(yGloc(I,J)+yGloc(I,J+1))
305              tanPhiAtU(i,j,bi,bj)=tan(lat*deg2rad)
306              lat=0.5*(yGloc(I,J)+yGloc(I+1,J))
307              tanPhiAtV(i,j,bi,bj)=tan(lat*deg2rad)
308           ENDDO           ENDDO
309          ENDDO          ENDDO
310         ENDDO  
311        ENDDO         ENDDO ! bi
312        _EXCH_XY_R4 (rAw      , myThid )        ENDDO ! bj
313        _EXCH_XY_R4 (rAs      , myThid )  
314  C        write(0,*) ' yC=', (yC(1,j,1,1),j=1,sNy)
315          write(0,*) 'dxF=', (dXF(1,j,1,1),j=1,sNy)
316          write(0,*) 'dyF=', (dYF(1,j,1,1),j=1,sNy)
317          write(0,*) 'dxG=', (dXG(1,j,1,1),j=1,sNy)
318          write(0,*) 'dyG=', (dYG(1,j,1,1),j=1,sNy)
319          write(0,*) 'dxC=', (dXC(1,j,1,1),j=1,sNy)
320          write(0,*) 'dyC=', (dYC(1,j,1,1),j=1,sNy)
321          write(0,*) 'dxV=', (dXV(1,j,1,1),j=1,sNy)
322          write(0,*) 'dyU=', (dYU(1,j,1,1),j=1,sNy)
323          write(0,*) ' rA=', (rA(1,j,1,1),j=1,sNy)
324          write(0,*) 'rAw=', (rAw(1,j,1,1),j=1,sNy)
325          write(0,*) 'rAs=', (rAs(1,j,1,1),j=1,sNy)
326    
327        RETURN        RETURN
328        END        END
Legend:
Removed from v.1.14  
changed lines
  Added in v.1.15
  ViewVC Help
Powered by ViewVC 1.1.22