/[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.16 - (hide annotations) (download)
Sun Feb 4 14:38:47 2001 UTC (23 years, 3 months ago) by cnh
Branch: MAIN
Changes since 1.15: +2 -1 lines
Made sure each .F and .h file had
the CVS keywords Header and Name at its start.
Most had header but very few currently have Name, so
lots of changes!

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

  ViewVC Help
Powered by ViewVC 1.1.22