/[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.7 by adcroft, Thu Jul 2 14:16:24 1998 UTC revision 1.18 by adcroft, Tue May 29 14:01:37 2001 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_EEOPTIONS.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 )
# 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(Nz), zLower(Nz)  
       _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 The functions iGl, jGl return the "global" index with valid values beyond
81    C halo regions
82    C cnh wrote:
83    C >    I don't understand why we would ever have to multiply the
84    C >    overlap by the total domain size e.g
85    C >    OLx*Nx, OLy*Ny.
86    C >    Can anybody explain? Lines are in ini_spherical_polar_grid.F.
87    C >    Surprised the code works if its wrong, so I am puzzled.        
88    C jmc relied:
89    C Yes, I can explain this since I put this modification to work
90    C with small domain (where Oly > Ny, as for instance, zonal-average
91    C case):
92    C This has no effect on the acuracy of the evaluation of iGl(I,bi)
93    C and jGl(J,bj) since we take mod(a+OLx*Nx,Nx) and mod(b+OLy*Ny,Ny).
94    C But in case a or b is negative, then the FORTRAN function "mod"
95    C does not return the matematical value of the "modulus" function,
96    C and this is not good for your purpose.
97    C This is why I add +OLx*Nx and +OLy*Ny to be sure that the 1rst
98    C argument of the mod function is positive.
99          INTEGER iGl,jGl
100          iGl(I,bi) = 1+mod(myXGlobalLo-1+(bi-1)*sNx+I+Olx*Nx-1,Nx)
101          jGl(J,bj) = 1+mod(myYGlobalLo-1+(bj-1)*sNy+J+Oly*Ny-1,Ny)
102    
103  C--   Example of inialisation for spherical polar grid  
104  C--   First set coordinates of cell centers  C     For each tile ...
 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  
105        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
        jG = myYGlobalLo + (bj-1)*sNy  
106         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
107    
108    C--     "Global" index (place holder)
109            jG = myYGlobalLo + (bj-1)*sNy
110          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 )  
111    
112  C--   Calculate separation between other points  C--   First find coordinate of tile corner (meaning outer corner of halo)
113  C     dxG, dyG are separations between cell corners along cell faces.          xG0 = thetaMin
114        DO bj = myByLo(myThid), myByHi(myThid)  C       Find the X-coordinate of the outer grid-line of the "real" tile
115         DO bi = myBxLo(myThid), myBxHi(myThid)          DO i=1, iG-1
116          DO J=1,sNy           xG0 = xG0 + delX(i)
117           DO I=1,sNx          ENDDO
118            jG = myYGlobalLo + (bj-1)*sNy + J-1  C       Back-step to the outer grid-line of the "halo" region
119            iG = myXGlobalLo + (bi-1)*sNx + I-1          DO i=1, Olx
120            lat = yc(I,J,bi,bj)-delY(jG) * 0.5 _d 0           xG0 = xG0 - delX( 1+mod(Olx*Nx-1+iG-i,Nx) )
121            dxG(I,J,bi,bj) = rSphere*COS(lat*deg2rad)*delX(iG)*deg2rad          ENDDO
122            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
123            yG0 = phiMin
124            DO j=1, jG-1
125             yG0 = yG0 + delY(j)
126            ENDDO
127    C       Back-step to the outer grid-line of the "halo" region
128            DO j=1, Oly
129             yG0 = yG0 - delY( 1+mod(Oly*Ny-1+jG-j,Ny) )
130            ENDDO
131    
132    C--     Calculate coordinates of cell corners for N+1 grid-lines
133            DO J=1-Oly,sNy+Oly +1
134             xGloc(1-Olx,J) = xG0
135             DO I=1-Olx,sNx+Olx
136    c         xGloc(I+1,J) = xGloc(I,J) + delX(1+mod(Nx-1+iG-1+i,Nx))
137              xGloc(I+1,J) = xGloc(I,J) + delX( iGl(I,bi) )
138           ENDDO           ENDDO
139          ENDDO          ENDDO
140         ENDDO          DO I=1-Olx,sNx+Olx +1
141        ENDDO           yGloc(I,1-Oly) = yG0
142        _EXCH_XY_R4(dxG, myThid )           DO J=1-Oly,sNy+Oly
143        _EXCH_XY_R4(dyG, myThid )  c         yGloc(I,J+1) = yGloc(I,J) + delY(1+mod(Ny-1+jG-1+j,Ny))
144  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  
145           ENDDO           ENDDO
146          ENDDO          ENDDO
147         ENDDO  
148        ENDDO  C--     Make a permanent copy of [xGloc,yGloc] in [xG,yG]
149        _EXCH_XY_R4(dxV, myThid )          DO J=1-Oly,sNy+Oly
150        _EXCH_XY_R4(dyU, myThid )           DO I=1-Olx,sNx+Olx
151  C     dxC, dyC is separation between cell centers            xG(I,J,bi,bj) = xGloc(I,J)
152        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  
           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  
153           ENDDO           ENDDO
154          ENDDO          ENDDO
        ENDDO  
       ENDDO  
       _EXCH_XY_R4(dxC, myThid )  
       _EXCH_XY_R4(dyC, myThid )  
 C     Calculate vertical face area and trigonometric terms  
       DO bj = myByLo(myThid), myByHi(myThid)  
        DO bi = myBxLo(myThid), myBxHi(myThid)  
         DO J=1,sNy  
          DO I=1,sNx  
           jG = myYGlobalLo + (bj-1)*sNy + J-1  
           latS = yc(i,j,bi,bj)-delY(jG)*0.5 _d 0  
           latN = yc(i,j,bi,bj)+delY(jG)*0.5 _d 0  
           zA(I,J,bi,bj) = dyF(I,J,bi,bj)  
      &    *rSphere*(SIN(latN*deg2rad)-SIN(latS*deg2rad))  
           tanPhiAtU(i,j,bi,bj)=tan(_yC(i,j,bi,bj)*deg2rad)  
           tanPhiAtV(i,j,bi,bj)=tan(latS*deg2rad)  
          ENDDO  
         ENDDO  
        ENDDO  
       ENDDO  
       _EXCH_XY_R4 (zA       , myThid )  
       _EXCH_XY_R4 (tanPhiAtU , myThid )  
       _EXCH_XY_R4 (tanPhiAtV , myThid )  
155    
156  C  C--     Calculate [xC,yC], coordinates of cell centers
157            DO J=1-Oly,sNy+Oly
158             DO I=1-Olx,sNx+Olx
159    C         by averaging
160              xC(I,J,bi,bj) = 0.25*(
161         &     xGloc(I,J)+xGloc(I+1,J)+xGloc(I,J+1)+xGloc(I+1,J+1) )
162              yC(I,J,bi,bj) = 0.25*(
163         &     yGloc(I,J)+yGloc(I+1,J)+yGloc(I,J+1)+yGloc(I+1,J+1) )
164             ENDDO
165            ENDDO
166    
167    C--     Calculate [dxF,dyF], lengths between cell faces (through center)
168            DO J=1-Oly,sNy+Oly
169             DO I=1-Olx,sNx+Olx
170    C         by averaging
171    c         dXF(I,J,bi,bj) = 0.5*(dXG(I,J,bi,bj)+dXG(I,J+1,bi,bj))
172    c         dYF(I,J,bi,bj) = 0.5*(dYG(I,J,bi,bj)+dYG(I+1,J,bi,bj))
173    C         by formula
174              lat = yC(I,J,bi,bj)
175              dlon = delX( iGl(I,bi) )
176              dlat = delY( jGl(J,bj) )
177              dXF(I,J,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
178    #ifdef    USE_BACKWARD_COMPATIBLE_GRID
179              dXF(I,J,bi,bj) = delX(iGl(I,bi))*deg2rad*rSphere*
180         &                     COS(yc(I,J,bi,bj)*deg2rad)
181    #endif    /* USE_BACKWARD_COMPATIBLE_GRID */
182              dYF(I,J,bi,bj) = rSphere*dlat*deg2rad
183             ENDDO
184            ENDDO
185    
186    C--     Calculate [dxG,dyG], lengths along cell boundaries
187            DO J=1-Oly,sNy+Oly
188             DO I=1-Olx,sNx+Olx
189    C         by averaging
190    c         dXG(I,J,bi,bj) = 0.5*(dXF(I,J,bi,bj)+dXF(I,J-1,bi,bj))
191    c         dYG(I,J,bi,bj) = 0.5*(dYF(I,J,bi,bj)+dYF(I-1,J,bi,bj))
192    C         by formula
193              lat = 0.5*(yGloc(I,J)+yGloc(I+1,J))
194              dlon = delX( iGl(I,bi) )
195              dlat = delY( jGl(J,bj) )
196              dXG(I,J,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
197              if (dXG(I,J,bi,bj).LT.1.) dXG(I,J,bi,bj)=0.
198              dYG(I,J,bi,bj) = rSphere*dlat*deg2rad
199             ENDDO
200            ENDDO
201    
202    C--     The following arrays are not defined in some parts of the halo
203    C       region. We set them to zero here for safety. If they are ever
204    C       referred to, especially in the denominator then it is a mistake!
205            DO J=1-Oly,sNy+Oly
206             DO I=1-Olx,sNx+Olx
207              dXC(I,J,bi,bj) = 0.
208              dYC(I,J,bi,bj) = 0.
209              dXV(I,J,bi,bj) = 0.
210              dYU(I,J,bi,bj) = 0.
211              rAw(I,J,bi,bj) = 0.
212              rAs(I,J,bi,bj) = 0.
213             ENDDO
214            ENDDO
215    
216    C--     Calculate [dxC], zonal length between cell centers
217            DO J=1-Oly,sNy+Oly
218             DO I=1-Olx+1,sNx+Olx ! NOTE range
219    C         by averaging
220              dXC(I,J,bi,bj) = 0.5D0*(dXF(I,J,bi,bj)+dXF(I-1,J,bi,bj))
221    C         by formula
222    c         lat = 0.5*(yC(I,J,bi,bj)+yC(I-1,J,bi,bj))
223    c         dlon = 0.5*(delX( iGl(I,bi) ) + delX( iGl(I-1,bi) ))
224    c         dXC(I,J,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
225    C         by difference
226    c         lat = 0.5*(yC(I,J,bi,bj)+yC(I-1,J,bi,bj))
227    c         dlon = (xC(I,J,bi,bj)+xC(I-1,J,bi,bj))
228    c         dXC(I,J,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
229             ENDDO
230            ENDDO
231    
232    C--     Calculate [dyC], meridional length between cell centers
233            DO J=1-Oly+1,sNy+Oly ! NOTE range
234             DO I=1-Olx,sNx+Olx
235    C         by averaging
236              dYC(I,J,bi,bj) = 0.5*(dYF(I,J,bi,bj)+dYF(I,J-1,bi,bj))
237    C         by formula
238    c         dlat = 0.5*(delY( jGl(J,bj) ) + delY( jGl(J-1,bj) ))
239    c         dYC(I,J,bi,bj) = rSphere*dlat*deg2rad
240    C         by difference
241    c         dlat = (yC(I,J,bi,bj)+yC(I,J-1,bi,bj))
242    c         dYC(I,J,bi,bj) = rSphere*dlat*deg2rad
243             ENDDO
244            ENDDO
245    
246    C--     Calculate [dxV,dyU], length between velocity points (through corners)
247            DO J=1-Oly+1,sNy+Oly ! NOTE range
248             DO I=1-Olx+1,sNx+Olx ! NOTE range
249    C         by averaging (method I)
250              dXV(I,J,bi,bj) = 0.5*(dXG(I,J,bi,bj)+dXG(I-1,J,bi,bj))
251              dYU(I,J,bi,bj) = 0.5*(dYG(I,J,bi,bj)+dYG(I,J-1,bi,bj))
252    C         by averaging (method II)
253    c         dXV(I,J,bi,bj) = 0.5*(dXG(I,J,bi,bj)+dXG(I-1,J,bi,bj))
254    c         dYU(I,J,bi,bj) = 0.5*(dYC(I,J,bi,bj)+dYC(I-1,J,bi,bj))
255             ENDDO
256            ENDDO
257    
258    C--     Calculate vertical face area (tracer cells)
259            DO J=1-Oly,sNy+Oly
260             DO I=1-Olx,sNx+Olx
261              lat=0.5*(yGloc(I,J)+yGloc(I+1,J))
262              dlon=delX( iGl(I,bi) )
263              dlat=delY( jGl(J,bj) )
264              rA(I,J,bi,bj) = rSphere*rSphere*dlon*deg2rad
265         &        *abs( sin((lat+dlat)*deg2rad)-sin(lat*deg2rad) )
266    #ifdef    USE_BACKWARD_COMPATIBLE_GRID
267              lat=yC(I,J,bi,bj)-delY( jGl(J,bj) )*0.5 _d 0
268              lon=yC(I,J,bi,bj)+delY( jGl(J,bj) )*0.5 _d 0
269              rA(I,J,bi,bj) = dyF(I,J,bi,bj)
270         &    *rSphere*(SIN(lon*deg2rad)-SIN(lat*deg2rad))
271    #endif    /* USE_BACKWARD_COMPATIBLE_GRID */
272             ENDDO
273            ENDDO
274    
275    C--     Calculate vertical face area (u cells)
276            DO J=1-Oly,sNy+Oly
277             DO I=1-Olx+1,sNx+Olx ! NOTE range
278    C         by averaging
279              rAw(I,J,bi,bj) = 0.5*(rA(I,J,bi,bj)+rA(I-1,J,bi,bj))
280    C         by formula
281    c         lat=yGloc(I,J)
282    c         dlon=0.5*( delX( iGl(I,bi) ) + delX( iGl(I-1,bi) ) )
283    c         dlat=delY( jGl(J,bj) )
284    c         rAw(I,J,bi,bj) = rSphere*rSphere*dlon*deg2rad
285    c    &        *abs( sin((lat+dlat)*deg2rad)-sin(lat*deg2rad) )
286             ENDDO
287            ENDDO
288    
289    C--     Calculate vertical face area (v cells)
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              rAs(I,J,bi,bj) = rSphere*rSphere*dlon*deg2rad
297         &        *abs( sin(lat*deg2rad)-sin((lat-dlat)*deg2rad) )
298    #ifdef    USE_BACKWARD_COMPATIBLE_GRID
299              lon=yC(I,J,bi,bj)-delY( jGl(J,bj) )
300              lat=yC(I,J,bi,bj)
301              rAs(I,J,bi,bj) = rSphere*delX(iGl(I,bi))*deg2rad
302         &    *rSphere*(SIN(lat*deg2rad)-SIN(lon*deg2rad))
303    #endif    /* USE_BACKWARD_COMPATIBLE_GRID */
304              IF (abs(lat).GT.90..OR.abs(lat-dlat).GT.90.) rAs(I,J,bi,bj)=0.
305             ENDDO
306            ENDDO
307    
308    C--     Calculate vertical face area (vorticity points)
309            DO J=1-Oly,sNy+Oly
310             DO I=1-Olx,sNx+Olx
311    C         by formula
312              lat=yC(I,J,bi,bj)
313              dlon=delX( iGl(I,bi) )
314              dlat=0.5*( delY( jGl(J,bj) ) + delY( jGl(J-1,bj) ) )
315              rAz(I,J,bi,bj) = rSphere*rSphere*dlon*deg2rad
316         &     *abs( sin(lat*deg2rad)-sin((lat-dlat)*deg2rad) )
317              IF (abs(lat).GT.90..OR.abs(lat-dlat).GT.90.) rAz(I,J,bi,bj)=0.
318             ENDDO
319            ENDDO
320    
321    C--     Calculate trigonometric terms
322            DO J=1-Oly,sNy+Oly
323             DO I=1-Olx,sNx+Olx
324              lat=0.5*(yGloc(I,J)+yGloc(I,J+1))
325              tanPhiAtU(i,j,bi,bj)=tan(lat*deg2rad)
326              lat=0.5*(yGloc(I,J)+yGloc(I+1,J))
327              tanPhiAtV(i,j,bi,bj)=tan(lat*deg2rad)
328             ENDDO
329            ENDDO
330    
331    C--     Cosine(lat) scaling
332            DO J=1-OLy,sNy+OLy
333             jG = myYGlobalLo + (bj-1)*sNy + J-1
334             jG = min(max(1,jG),Ny)
335             IF (cosPower.NE.0.) THEN
336              cosFacU(J,bi,bj)=COS(yC(1,J,bi,bj)*deg2rad)
337         &                    **cosPower
338              cosFacV(J,bi,bj)=COS((yC(1,J,bi,bj)-0.5*delY(jG))*deg2rad)
339         &                    **cosPower
340              sqcosFacU(J,bi,bj)=sqrt(cosFacU(J,bi,bj))
341              sqcosFacV(J,bi,bj)=sqrt(cosFacV(J,bi,bj))
342             ELSE
343              cosFacU(J,bi,bj)=1.
344              cosFacV(J,bi,bj)=1.
345              sqcosFacU(J,bi,bj)=1.
346              sqcosFacV(J,bi,bj)=1.
347             ENDIF
348            ENDDO
349    
350           ENDDO ! bi
351          ENDDO ! bj
352    
353        RETURN        RETURN
354        END        END
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.18
  ViewVC Help
Powered by ViewVC 1.1.22