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

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

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


Revision 1.15 - (hide annotations) (download)
Fri Feb 2 21:04:48 2001 UTC (23 years, 3 months ago) by adcroft
Branch: MAIN
Changes since 1.14: +247 -133 lines
Merged changes from branch "branch-atmos-merge" into MAIN (checkpoint34)
 - substantial modifications to algorithm sequence (dynamics.F)
 - packaged OBCS, Shapiro filter, Zonal filter, Atmospheric Physics

1 adcroft 1.15 C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_spherical_polar_grid.F,v 1.14.2.4 2001/01/30 19:12:15 jmc Exp $
2 cnh 1.1
3 cnh 1.10 #include "CPP_OPTIONS.h"
4 cnh 1.1
5 adcroft 1.15 #undef USE_BACKWARD_COMPATIBLE_GRID
6    
7 cnh 1.1 CStartOfInterface
8     SUBROUTINE INI_SPHERICAL_POLAR_GRID( myThid )
9     C /==========================================================\
10     C | SUBROUTINE INI_SPHERICAL_POLAR_GRID |
11     C | o Initialise model coordinate system |
12     C |==========================================================|
13     C | These arrays are used throughout the code in evaluating |
14     C | gradients, integrals and spatial avarages. This routine |
15     C | is called separately by each thread and initialise only |
16     C | the region of the domain it is "responsible" for. |
17     C | Notes: |
18     C | Two examples are included. One illustrates the |
19     C | initialisation of a cartesian grid. The other shows the |
20     C | inialisation of a spherical polar grid. Other orthonormal|
21     C | grids can be fitted into this design. In this case |
22     C | custom metric terms also need adding to account for the |
23     C | projections of velocity vectors onto these grids. |
24     C | The structure used here also makes it possible to |
25     C | implement less regular grid mappings. In particular |
26     C | o Schemes which leave out blocks of the domain that are |
27     C | all land could be supported. |
28     C | o Multi-level schemes such as icosohedral or cubic |
29     C | grid projections onto a sphere can also be fitted |
30     C | within the strategy we use. |
31     C | Both of the above also require modifying the support |
32     C | routines that map computational blocks to simulation |
33     C | domain blocks. |
34     C | Under the spherical polar grid mode primitive distances |
35     C | in X and Y are in degrees. Distance in Z are in m or Pa |
36     C | depending on the vertical gridding mode. |
37     C \==========================================================/
38 adcroft 1.12 IMPLICIT NONE
39 cnh 1.1
40     C === Global variables ===
41     #include "SIZE.h"
42     #include "EEPARAMS.h"
43     #include "PARAMS.h"
44     #include "GRID.h"
45    
46     C == Routine arguments ==
47     C myThid - Number of this instance of INI_CARTESIAN_GRID
48     INTEGER myThid
49     CEndOfInterface
50    
51     C == Local variables ==
52     C xG, yG - Global coordinate location.
53     C xBase - South-west corner location for process.
54     C yBase
55 adcroft 1.15 C zUpper - Work arrays for upper and lower
56     C zLower cell-face heights.
57     C phi - Temporary scalar
58 cnh 1.1 C iG, jG - Global coordinate index. Usually used to hold
59     C the south-west global coordinate of a tile.
60     C bi,bj - Loop counters
61     C zUpper - Temporary arrays holding z coordinates of
62     C zLower upper and lower faces.
63     C xBase - Lower coordinate for this threads cells
64     C yBase
65     C lat, latN, - Temporary variables used to hold latitude
66     C latS values.
67     C I,J,K
68     INTEGER iG, jG
69     INTEGER bi, bj
70 adcroft 1.14 INTEGER I, J
71 adcroft 1.15 _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 cnh 1.1
85 adcroft 1.15 C For each tile ...
86 cnh 1.1 DO bj = myByLo(myThid), myByHi(myThid)
87     DO bi = myBxLo(myThid), myBxHi(myThid)
88 adcroft 1.15
89     C-- "Global" index (place holder)
90     jG = myYGlobalLo + (bj-1)*sNy
91 cnh 1.1 iG = myXGlobalLo + (bi-1)*sNx
92    
93 adcroft 1.15 C-- First find coordinate of tile corner (meaning outer corner of halo)
94     xG0 = thetaMin
95     C Find the X-coordinate of the outer grid-line of the "real" tile
96     DO i=1, iG-1
97     xG0 = xG0 + delX(i)
98     ENDDO
99     C Back-step to the outer grid-line of the "halo" region
100     DO i=1, Olx
101     xG0 = xG0 - delX( 1+mod(Olx*Nx-1+iG-i,Nx) )
102     ENDDO
103     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
120     ENDDO
121     DO I=1-Olx,sNx+Olx +1
122     yGloc(I,1-Oly) = yG0
123     DO J=1-Oly,sNy+Oly
124     c yGloc(I,J+1) = yGloc(I,J) + delY(1+mod(Ny-1+jG-1+j,Ny))
125     yGloc(I,J+1) = yGloc(I,J) + delY( jGl(J,bj) )
126     ENDDO
127     ENDDO
128    
129     C-- Make a permanent copy of [xGloc,yGloc] in [xG,yG]
130     DO J=1-Oly,sNy+Oly
131     DO I=1-Olx,sNx+Olx
132     xG(I,J,bi,bj) = xGloc(I,J)
133     yG(I,J,bi,bj) = yGloc(I,J)
134     ENDDO
135     ENDDO
136    
137     C-- Calculate [xC,yC], coordinates of cell centers
138     DO J=1-Oly,sNy+Oly
139     DO I=1-Olx,sNx+Olx
140     C by averaging
141     xC(I,J,bi,bj) = 0.25*(
142     & xGloc(I,J)+xGloc(I+1,J)+xGloc(I,J+1)+xGloc(I+1,J+1) )
143     yC(I,J,bi,bj) = 0.25*(
144     & yGloc(I,J)+yGloc(I+1,J)+yGloc(I,J+1)+yGloc(I+1,J+1) )
145     ENDDO
146     ENDDO
147    
148     C-- Calculate [dxF,dyF], lengths between cell faces (through center)
149     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 cnh 1.1 ENDDO
210     ENDDO
211 adcroft 1.15
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 cnh 1.1 ENDDO
224     ENDDO
225 adcroft 1.15
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 cnh 1.1 ENDDO
236     ENDDO
237 adcroft 1.15
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 cnh 1.8 rA(I,J,bi,bj) = dyF(I,J,bi,bj)
250 adcroft 1.15 & *rSphere*(SIN(lon*deg2rad)-SIN(lat*deg2rad))
251     #endif /* USE_BACKWARD_COMPATIBLE_GRID */
252     ENDDO
253     ENDDO
254    
255     C-- Calculate vertical face area (u cells)
256     DO J=1-Oly,sNy+Oly
257     DO I=1-Olx+1,sNx+Olx ! NOTE range
258     C by averaging
259 adcroft 1.11 rAw(I,J,bi,bj) = 0.5*(rA(I,J,bi,bj)+rA(I-1,J,bi,bj))
260 adcroft 1.15 C by formula
261     c lat=yGloc(I,J)
262     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 adcroft 1.11 ENDDO
309     ENDDO
310 adcroft 1.15
311     ENDDO ! bi
312     ENDDO ! bj
313    
314     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 cnh 1.1 RETURN
328     END

  ViewVC Help
Powered by ViewVC 1.1.22