/[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.26 - (hide annotations) (download)
Tue Jan 27 15:35:27 2009 UTC (15 years, 3 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62c, checkpoint62a, checkpoint62e, checkpoint62d, checkpoint62, checkpoint62b, checkpoint61n, checkpoint61q, checkpoint61o, checkpoint61l, checkpoint61m, checkpoint61j, checkpoint61k, checkpoint61i, checkpoint61v, checkpoint61w, checkpoint61t, checkpoint61u, checkpoint61r, checkpoint61s, checkpoint61p, checkpoint61z, checkpoint61x, checkpoint61y
Changes since 1.25: +8 -8 lines
rename thetaMin,phiMin -> xgOrigin,ygOrigin
 (temporary backward compatibility in ini_parms.F, until next checkpoint)

1 jmc 1.26 C $Header: /u/gcmpack/MITgcm/model/src/ini_spherical_polar_grid.F,v 1.25 2008/02/08 13:01:25 mlosch Exp $
2 cnh 1.17 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.19 CBOP
9     C !ROUTINE: INI_SPHERICAL_POLAR_GRID
10     C !INTERFACE:
11 cnh 1.1 SUBROUTINE INI_SPHERICAL_POLAR_GRID( myThid )
12 cnh 1.19 C !DESCRIPTION: \bv
13 cnh 1.1 C /==========================================================\
14     C | SUBROUTINE INI_SPHERICAL_POLAR_GRID |
15 cnh 1.19 C | o Initialise model coordinate system arrays |
16 cnh 1.1 C |==========================================================|
17     C | These arrays are used throughout the code in evaluating |
18     C | gradients, integrals and spatial avarages. This routine |
19     C | is called separately by each thread and initialise only |
20     C | the region of the domain it is "responsible" for. |
21     C | Under the spherical polar grid mode primitive distances |
22     C | in X and Y are in degrees. Distance in Z are in m or Pa |
23     C | depending on the vertical gridding mode. |
24     C \==========================================================/
25 cnh 1.19 C \ev
26    
27     C !USES:
28 adcroft 1.12 IMPLICIT NONE
29 cnh 1.1 C === Global variables ===
30     #include "SIZE.h"
31     #include "EEPARAMS.h"
32     #include "PARAMS.h"
33     #include "GRID.h"
34    
35 cnh 1.19 C !INPUT/OUTPUT PARAMETERS:
36 cnh 1.1 C == Routine arguments ==
37     C myThid - Number of this instance of INI_CARTESIAN_GRID
38     INTEGER myThid
39     CEndOfInterface
40    
41 cnh 1.19 C !LOCAL VARIABLES:
42 cnh 1.1 C == Local variables ==
43     C xG, yG - Global coordinate location.
44     C xBase - South-west corner location for process.
45     C yBase
46 jmc 1.26 C zUpper - Work arrays for upper and lower
47 adcroft 1.15 C zLower cell-face heights.
48     C phi - Temporary scalar
49 cnh 1.1 C iG, jG - Global coordinate index. Usually used to hold
50     C the south-west global coordinate of a tile.
51     C bi,bj - Loop counters
52     C zUpper - Temporary arrays holding z coordinates of
53     C zLower upper and lower faces.
54     C xBase - Lower coordinate for this threads cells
55     C yBase
56     C lat, latN, - Temporary variables used to hold latitude
57     C latS values.
58     C I,J,K
59     INTEGER iG, jG
60     INTEGER bi, bj
61 adcroft 1.14 INTEGER I, J
62 adcroft 1.15 _RL lat, dlat, dlon, xG0, yG0
63    
64    
65 cnh 1.19 C "Long" real for temporary coordinate calculation
66     C NOTICE the extended range of indices!!
67 adcroft 1.15 _RL xGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)
68     _RL yGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)
69    
70 cnh 1.19 C The functions iGl, jGl return the "global" index with valid values beyond
71     C halo regions
72     C cnh wrote:
73 adcroft 1.20 C > I dont understand why we would ever have to multiply the
74 cnh 1.19 C > overlap by the total domain size e.g
75     C > OLx*Nx, OLy*Ny.
76     C > Can anybody explain? Lines are in ini_spherical_polar_grid.F.
77 jmc 1.26 C > Surprised the code works if its wrong, so I am puzzled.
78 cnh 1.19 C jmc replied:
79     C Yes, I can explain this since I put this modification to work
80     C with small domain (where Oly > Ny, as for instance, zonal-average
81     C case):
82     C This has no effect on the acuracy of the evaluation of iGl(I,bi)
83     C and jGl(J,bj) since we take mod(a+OLx*Nx,Nx) and mod(b+OLy*Ny,Ny).
84     C But in case a or b is negative, then the FORTRAN function "mod"
85     C does not return the matematical value of the "modulus" function,
86     C and this is not good for your purpose.
87 jmc 1.26 C This is why I add +OLx*Nx and +OLy*Ny to be sure that the 1rst
88 cnh 1.19 C argument of the mod function is positive.
89 adcroft 1.15 INTEGER iGl,jGl
90     iGl(I,bi) = 1+mod(myXGlobalLo-1+(bi-1)*sNx+I+Olx*Nx-1,Nx)
91     jGl(J,bj) = 1+mod(myYGlobalLo-1+(bj-1)*sNy+J+Oly*Ny-1,Ny)
92 cnh 1.19 CEOP
93 cnh 1.1
94 adcroft 1.18
95 adcroft 1.15 C For each tile ...
96 cnh 1.1 DO bj = myByLo(myThid), myByHi(myThid)
97     DO bi = myBxLo(myThid), myBxHi(myThid)
98 adcroft 1.15
99     C-- "Global" index (place holder)
100     jG = myYGlobalLo + (bj-1)*sNy
101 cnh 1.1 iG = myXGlobalLo + (bi-1)*sNx
102    
103 adcroft 1.15 C-- First find coordinate of tile corner (meaning outer corner of halo)
104 jmc 1.26 xG0 = xgOrigin
105 adcroft 1.15 C Find the X-coordinate of the outer grid-line of the "real" tile
106     DO i=1, iG-1
107     xG0 = xG0 + delX(i)
108     ENDDO
109     C Back-step to the outer grid-line of the "halo" region
110     DO i=1, Olx
111     xG0 = xG0 - delX( 1+mod(Olx*Nx-1+iG-i,Nx) )
112     ENDDO
113     C Find the Y-coordinate of the outer grid-line of the "real" tile
114 jmc 1.26 yG0 = ygOrigin
115 adcroft 1.15 DO j=1, jG-1
116     yG0 = yG0 + delY(j)
117     ENDDO
118     C Back-step to the outer grid-line of the "halo" region
119     DO j=1, Oly
120     yG0 = yG0 - delY( 1+mod(Oly*Ny-1+jG-j,Ny) )
121     ENDDO
122    
123     C-- Calculate coordinates of cell corners for N+1 grid-lines
124     DO J=1-Oly,sNy+Oly +1
125     xGloc(1-Olx,J) = xG0
126     DO I=1-Olx,sNx+Olx
127     c xGloc(I+1,J) = xGloc(I,J) + delX(1+mod(Nx-1+iG-1+i,Nx))
128     xGloc(I+1,J) = xGloc(I,J) + delX( iGl(I,bi) )
129     ENDDO
130     ENDDO
131     DO I=1-Olx,sNx+Olx +1
132     yGloc(I,1-Oly) = yG0
133     DO J=1-Oly,sNy+Oly
134     c yGloc(I,J+1) = yGloc(I,J) + delY(1+mod(Ny-1+jG-1+j,Ny))
135     yGloc(I,J+1) = yGloc(I,J) + delY( jGl(J,bj) )
136     ENDDO
137     ENDDO
138    
139     C-- Make a permanent copy of [xGloc,yGloc] in [xG,yG]
140     DO J=1-Oly,sNy+Oly
141     DO I=1-Olx,sNx+Olx
142     xG(I,J,bi,bj) = xGloc(I,J)
143     yG(I,J,bi,bj) = yGloc(I,J)
144     ENDDO
145     ENDDO
146    
147     C-- Calculate [xC,yC], coordinates of cell centers
148     DO J=1-Oly,sNy+Oly
149     DO I=1-Olx,sNx+Olx
150     C by averaging
151 jmc 1.26 xC(I,J,bi,bj) = 0.25*(
152 adcroft 1.15 & xGloc(I,J)+xGloc(I+1,J)+xGloc(I,J+1)+xGloc(I+1,J+1) )
153 jmc 1.26 yC(I,J,bi,bj) = 0.25*(
154 adcroft 1.15 & yGloc(I,J)+yGloc(I+1,J)+yGloc(I,J+1)+yGloc(I+1,J+1) )
155     ENDDO
156     ENDDO
157    
158     C-- Calculate [dxF,dyF], lengths between cell faces (through center)
159     DO J=1-Oly,sNy+Oly
160     DO I=1-Olx,sNx+Olx
161     C by averaging
162     c dXF(I,J,bi,bj) = 0.5*(dXG(I,J,bi,bj)+dXG(I,J+1,bi,bj))
163     c dYF(I,J,bi,bj) = 0.5*(dYG(I,J,bi,bj)+dYG(I+1,J,bi,bj))
164     C by formula
165     lat = yC(I,J,bi,bj)
166     dlon = delX( iGl(I,bi) )
167     dlat = delY( jGl(J,bj) )
168     dXF(I,J,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
169     #ifdef USE_BACKWARD_COMPATIBLE_GRID
170     dXF(I,J,bi,bj) = delX(iGl(I,bi))*deg2rad*rSphere*
171     & COS(yc(I,J,bi,bj)*deg2rad)
172     #endif /* USE_BACKWARD_COMPATIBLE_GRID */
173     dYF(I,J,bi,bj) = rSphere*dlat*deg2rad
174     ENDDO
175     ENDDO
176    
177     C-- Calculate [dxG,dyG], lengths along cell boundaries
178     DO J=1-Oly,sNy+Oly
179     DO I=1-Olx,sNx+Olx
180     C by averaging
181     c dXG(I,J,bi,bj) = 0.5*(dXF(I,J,bi,bj)+dXF(I,J-1,bi,bj))
182     c dYG(I,J,bi,bj) = 0.5*(dYF(I,J,bi,bj)+dYF(I-1,J,bi,bj))
183     C by formula
184     lat = 0.5*(yGloc(I,J)+yGloc(I+1,J))
185     dlon = delX( iGl(I,bi) )
186     dlat = delY( jGl(J,bj) )
187     dXG(I,J,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
188 adcroft 1.18 if (dXG(I,J,bi,bj).LT.1.) dXG(I,J,bi,bj)=0.
189 adcroft 1.15 dYG(I,J,bi,bj) = rSphere*dlat*deg2rad
190     ENDDO
191     ENDDO
192    
193     C-- The following arrays are not defined in some parts of the halo
194     C region. We set them to zero here for safety. If they are ever
195     C referred to, especially in the denominator then it is a mistake!
196     DO J=1-Oly,sNy+Oly
197     DO I=1-Olx,sNx+Olx
198     dXC(I,J,bi,bj) = 0.
199     dYC(I,J,bi,bj) = 0.
200     dXV(I,J,bi,bj) = 0.
201     dYU(I,J,bi,bj) = 0.
202     rAw(I,J,bi,bj) = 0.
203     rAs(I,J,bi,bj) = 0.
204     ENDDO
205     ENDDO
206    
207     C-- Calculate [dxC], zonal length between cell centers
208     DO J=1-Oly,sNy+Oly
209     DO I=1-Olx+1,sNx+Olx ! NOTE range
210     C by averaging
211     dXC(I,J,bi,bj) = 0.5D0*(dXF(I,J,bi,bj)+dXF(I-1,J,bi,bj))
212     C by formula
213     c lat = 0.5*(yC(I,J,bi,bj)+yC(I-1,J,bi,bj))
214     c dlon = 0.5*(delX( iGl(I,bi) ) + delX( iGl(I-1,bi) ))
215     c dXC(I,J,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
216     C by difference
217     c lat = 0.5*(yC(I,J,bi,bj)+yC(I-1,J,bi,bj))
218     c dlon = (xC(I,J,bi,bj)+xC(I-1,J,bi,bj))
219     c dXC(I,J,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
220 cnh 1.1 ENDDO
221     ENDDO
222 adcroft 1.15
223     C-- Calculate [dyC], meridional length between cell centers
224     DO J=1-Oly+1,sNy+Oly ! NOTE range
225     DO I=1-Olx,sNx+Olx
226     C by averaging
227     dYC(I,J,bi,bj) = 0.5*(dYF(I,J,bi,bj)+dYF(I,J-1,bi,bj))
228     C by formula
229     c dlat = 0.5*(delY( jGl(J,bj) ) + delY( jGl(J-1,bj) ))
230     c dYC(I,J,bi,bj) = rSphere*dlat*deg2rad
231     C by difference
232     c dlat = (yC(I,J,bi,bj)+yC(I,J-1,bi,bj))
233     c dYC(I,J,bi,bj) = rSphere*dlat*deg2rad
234 cnh 1.1 ENDDO
235     ENDDO
236 adcroft 1.15
237     C-- Calculate [dxV,dyU], length between velocity points (through corners)
238     DO J=1-Oly+1,sNy+Oly ! NOTE range
239     DO I=1-Olx+1,sNx+Olx ! NOTE range
240     C by averaging (method I)
241     dXV(I,J,bi,bj) = 0.5*(dXG(I,J,bi,bj)+dXG(I-1,J,bi,bj))
242     dYU(I,J,bi,bj) = 0.5*(dYG(I,J,bi,bj)+dYG(I,J-1,bi,bj))
243     C by averaging (method II)
244     c dXV(I,J,bi,bj) = 0.5*(dXG(I,J,bi,bj)+dXG(I-1,J,bi,bj))
245     c dYU(I,J,bi,bj) = 0.5*(dYC(I,J,bi,bj)+dYC(I-1,J,bi,bj))
246 cnh 1.1 ENDDO
247     ENDDO
248 adcroft 1.15
249     C-- Calculate vertical face area (tracer cells)
250     DO J=1-Oly,sNy+Oly
251     DO I=1-Olx,sNx+Olx
252     lat=0.5*(yGloc(I,J)+yGloc(I+1,J))
253     dlon=delX( iGl(I,bi) )
254     dlat=delY( jGl(J,bj) )
255     rA(I,J,bi,bj) = rSphere*rSphere*dlon*deg2rad
256     & *abs( sin((lat+dlat)*deg2rad)-sin(lat*deg2rad) )
257     #ifdef USE_BACKWARD_COMPATIBLE_GRID
258     lat=yC(I,J,bi,bj)-delY( jGl(J,bj) )*0.5 _d 0
259     lon=yC(I,J,bi,bj)+delY( jGl(J,bj) )*0.5 _d 0
260 cnh 1.8 rA(I,J,bi,bj) = dyF(I,J,bi,bj)
261 adcroft 1.15 & *rSphere*(SIN(lon*deg2rad)-SIN(lat*deg2rad))
262     #endif /* USE_BACKWARD_COMPATIBLE_GRID */
263     ENDDO
264     ENDDO
265    
266     C-- Calculate vertical face area (u cells)
267     DO J=1-Oly,sNy+Oly
268     DO I=1-Olx+1,sNx+Olx ! NOTE range
269     C by averaging
270 adcroft 1.11 rAw(I,J,bi,bj) = 0.5*(rA(I,J,bi,bj)+rA(I-1,J,bi,bj))
271 adcroft 1.15 C by formula
272     c lat=yGloc(I,J)
273     c dlon=0.5*( delX( iGl(I,bi) ) + delX( iGl(I-1,bi) ) )
274     c dlat=delY( jGl(J,bj) )
275     c rAw(I,J,bi,bj) = rSphere*rSphere*dlon*deg2rad
276     c & *abs( sin((lat+dlat)*deg2rad)-sin(lat*deg2rad) )
277     ENDDO
278     ENDDO
279    
280     C-- Calculate vertical face area (v cells)
281     DO J=1-Oly,sNy+Oly
282     DO I=1-Olx,sNx+Olx
283     C by formula
284     lat=yC(I,J,bi,bj)
285     dlon=delX( iGl(I,bi) )
286     dlat=0.5*( delY( jGl(J,bj) ) + delY( jGl(J-1,bj) ) )
287     rAs(I,J,bi,bj) = rSphere*rSphere*dlon*deg2rad
288     & *abs( sin(lat*deg2rad)-sin((lat-dlat)*deg2rad) )
289     #ifdef USE_BACKWARD_COMPATIBLE_GRID
290     lon=yC(I,J,bi,bj)-delY( jGl(J,bj) )
291     lat=yC(I,J,bi,bj)
292     rAs(I,J,bi,bj) = rSphere*delX(iGl(I,bi))*deg2rad
293     & *rSphere*(SIN(lat*deg2rad)-SIN(lon*deg2rad))
294     #endif /* USE_BACKWARD_COMPATIBLE_GRID */
295     IF (abs(lat).GT.90..OR.abs(lat-dlat).GT.90.) rAs(I,J,bi,bj)=0.
296     ENDDO
297     ENDDO
298    
299     C-- Calculate vertical face area (vorticity points)
300     DO J=1-Oly,sNy+Oly
301     DO I=1-Olx,sNx+Olx
302     C by formula
303 jmc 1.21 lat =0.5 _d 0*(yGloc(I,J)+yGloc(I,J+1))
304     dlon=0.5 _d 0*( delX( iGl(I,bi) ) + delX( iGl(I-1,bi) ) )
305     dlat=0.5 _d 0*( delY( jGl(J,bj) ) + delY( jGl(J-1,bj) ) )
306 adcroft 1.15 rAz(I,J,bi,bj) = rSphere*rSphere*dlon*deg2rad
307     & *abs( sin(lat*deg2rad)-sin((lat-dlat)*deg2rad) )
308     IF (abs(lat).GT.90..OR.abs(lat-dlat).GT.90.) rAz(I,J,bi,bj)=0.
309     ENDDO
310     ENDDO
311    
312 jmc 1.23 C-- Calculate trigonometric terms & grid orientation:
313 adcroft 1.15 DO J=1-Oly,sNy+Oly
314     DO I=1-Olx,sNx+Olx
315     lat=0.5*(yGloc(I,J)+yGloc(I,J+1))
316 jmc 1.23 tanPhiAtU(I,J,bi,bj)=tan(lat*deg2rad)
317 adcroft 1.15 lat=0.5*(yGloc(I,J)+yGloc(I+1,J))
318 jmc 1.23 tanPhiAtV(I,J,bi,bj)=tan(lat*deg2rad)
319     angleCosC(I,J,bi,bj) = 1.
320     angleSinC(I,J,bi,bj) = 0.
321 adcroft 1.11 ENDDO
322 adcroft 1.18 ENDDO
323    
324     C-- Cosine(lat) scaling
325     DO J=1-OLy,sNy+OLy
326     jG = myYGlobalLo + (bj-1)*sNy + J-1
327     jG = min(max(1,jG),Ny)
328     IF (cosPower.NE.0.) THEN
329     cosFacU(J,bi,bj)=COS(yC(1,J,bi,bj)*deg2rad)
330     & **cosPower
331     cosFacV(J,bi,bj)=COS((yC(1,J,bi,bj)-0.5*delY(jG))*deg2rad)
332     & **cosPower
333 jmc 1.22 cosFacU(J,bi,bj)=ABS(cosFacU(J,bi,bj))
334     cosFacV(J,bi,bj)=ABS(cosFacV(J,bi,bj))
335 adcroft 1.18 sqcosFacU(J,bi,bj)=sqrt(cosFacU(J,bi,bj))
336     sqcosFacV(J,bi,bj)=sqrt(cosFacV(J,bi,bj))
337     ELSE
338     cosFacU(J,bi,bj)=1.
339     cosFacV(J,bi,bj)=1.
340     sqcosFacU(J,bi,bj)=1.
341     sqcosFacV(J,bi,bj)=1.
342     ENDIF
343 adcroft 1.11 ENDDO
344 adcroft 1.15
345     ENDDO ! bi
346     ENDDO ! bj
347    
348 mlosch 1.24 IF ( rotateGrid ) THEN
349     CALL ROTATE_SPHERICAL_POLAR_GRID( xC, yC, myThid )
350     CALL ROTATE_SPHERICAL_POLAR_GRID( xG, yG, myThid )
351 mlosch 1.25 CALL CALC_ANGLES( myThid )
352 mlosch 1.24 ENDIF
353    
354 cnh 1.1 RETURN
355     END

  ViewVC Help
Powered by ViewVC 1.1.22