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