-
Notifications
You must be signed in to change notification settings - Fork 13
/
ini_spherical_polar_grid.F
284 lines (259 loc) · 9.86 KB
/
ini_spherical_polar_grid.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
#include "CPP_OPTIONS.h"
#undef USE_BACKWARD_COMPATIBLE_GRID
CBOP
C !ROUTINE: INI_SPHERICAL_POLAR_GRID
C !INTERFACE:
SUBROUTINE INI_SPHERICAL_POLAR_GRID( myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE INI_SPHERICAL_POLAR_GRID
C | o Initialise model coordinate system arrays
C *==========================================================*
C | These arrays are used throughout the code in evaluating
C | gradients, integrals and spatial avarages. This routine
C | is called separately by each thread and initialise only
C | the region of the domain it is "responsible" for.
C | Under the spherical polar grid mode primitive distances
C | in X and Y are in degrees. Distance in Z are in m or Pa
C | depending on the vertical gridding mode.
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
C myThid :: my Thread Id Number
INTEGER myThid
C !LOCAL VARIABLES:
C == Local variables ==
C bi,bj :: tile indices
C i, j :: loop counters
C lat :: Temporary variables used to hold latitude values.
C dlat :: Temporary variables used to hold latitudes increment.
C dlon :: Temporary variables used to hold longitude increment.
C delXloc :: mesh spacing in X direction
C delYloc :: mesh spacing in Y direction
C xGloc :: mesh corner-point location (local "Long" real array type)
C yGloc :: mesh corner-point location (local "Long" real array type)
LOGICAL skipCalcAngleC
INTEGER bi, bj
INTEGER i, j
INTEGER gridNx, gridNy
_RL lat, dlat, dlon
C NOTICE the extended range of indices!!
_RL delXloc(0-OLx:sNx+OLx)
_RL delYloc(0-OLy:sNy+OLy)
C NOTICE the extended range of indices!!
_RL xGloc(1-OLx:sNx+OLx+1,1-OLy:sNy+OLy+1)
_RL yGloc(1-OLx:sNx+OLx+1,1-OLy:sNy+OLy+1)
CEOP
C-- For each tile ...
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
C-- set tile local mesh (same units as delX,deY)
C corresponding to coordinates of cell corners for N+1 grid-lines
CALL INI_LOCAL_GRID(
O xGloc, yGloc,
O delXloc, delYloc,
O gridNx, gridNy,
I bi, bj, myThid )
C-- Make a permanent copy of [xGloc,yGloc] in [xG,yG]
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
xG(i,j,bi,bj) = xGloc(i,j)
yG(i,j,bi,bj) = yGloc(i,j)
ENDDO
ENDDO
C-- Calculate [xC,yC], coordinates of cell centers
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
C by averaging
xC(i,j,bi,bj) = 0.25 _d 0*(
& xGloc(i,j)+xGloc(i+1,j)+xGloc(i,j+1)+xGloc(i+1,j+1) )
yC(i,j,bi,bj) = 0.25 _d 0*(
& yGloc(i,j)+yGloc(i+1,j)+yGloc(i,j+1)+yGloc(i+1,j+1) )
ENDDO
ENDDO
C-- Calculate [dxF,dyF], lengths between cell faces (through center)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
C by averaging
c dxF(i,j,bi,bj) = 0.5*(dxG(i,j,bi,bj)+dxG(i,j+1,bi,bj))
c dyF(i,j,bi,bj) = 0.5*(dyG(i,j,bi,bj)+dyG(i+1,j,bi,bj))
C by formula
lat = yC(i,j,bi,bj)
dlon = delXloc(i)
dlat = delYloc(j)
dxF(i,j,bi,bj) = rSphere*COS(lat*deg2rad)*dlon*deg2rad
dyF(i,j,bi,bj) = rSphere*dlat*deg2rad
ENDDO
ENDDO
C-- Calculate [dxG,dyG], lengths along cell boundaries
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
C by averaging
c dxG(i,j,bi,bj) = 0.5*(dxF(i,j,bi,bj)+dxF(i,j-1,bi,bj))
c dyG(i,j,bi,bj) = 0.5*(dyF(i,j,bi,bj)+dyF(i-1,j,bi,bj))
C by formula
lat = 0.5 _d 0*(yGloc(i,j)+yGloc(i+1,j))
dlon = delXloc(i)
dlat = delYloc(j)
dxG(i,j,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
IF (dxG(i,j,bi,bj).LT.1.) dxG(i,j,bi,bj)=0.
dyG(i,j,bi,bj) = rSphere*dlat*deg2rad
ENDDO
ENDDO
C-- The following arrays are not defined in some parts of the halo
C region. We set them to zero here for safety. If they are ever
C referred to, especially in the denominator then it is a mistake!
C Note: this is now done earlier in main S/R INI_GRID
c DO j=1-OLy,sNy+OLy
c DO i=1-OLx,sNx+OLx
c dxC(i,j,bi,bj) = 0.
c dyC(i,j,bi,bj) = 0.
c dxV(i,j,bi,bj) = 0.
c dyU(i,j,bi,bj) = 0.
c rAw(i,j,bi,bj) = 0.
c rAs(i,j,bi,bj) = 0.
c ENDDO
c ENDDO
C-- Calculate [dxC], zonal length between cell centers
DO j=1-OLy,sNy+OLy
DO i=1-OLx+1,sNx+OLx ! NOTE range
C by averaging
dxC(i,j,bi,bj) = 0.5 _d 0*(dxF(i,j,bi,bj)+dxF(i-1,j,bi,bj))
C by formula
c lat = 0.5*(yC(i,j,bi,bj)+yC(i-1,j,bi,bj))
c dlon = 0.5*( delXloc(i) + delXloc(i-1) )
c dxC(i,j,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
C by difference
c lat = 0.5*(yC(i,j,bi,bj)+yC(i-1,j,bi,bj))
c dlon = (xC(i,j,bi,bj)-xC(i-1,j,bi,bj))
c dxC(i,j,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
ENDDO
ENDDO
C-- Calculate [dyC], meridional length between cell centers
DO j=1-OLy+1,sNy+OLy ! NOTE range
DO i=1-OLx,sNx+OLx
C by averaging
dyC(i,j,bi,bj) = 0.5 _d 0*(dyF(i,j,bi,bj)+dyF(i,j-1,bi,bj))
C by formula
c dlat = 0.5*( delYloc(j) + delYloc(j-1) )
c dyC(i,j,bi,bj) = rSphere*dlat*deg2rad
C by difference
c dlat = (yC(i,j,bi,bj)-yC(i,j-1,bi,bj))
c dyC(i,j,bi,bj) = rSphere*dlat*deg2rad
ENDDO
ENDDO
C-- Calculate [dxV,dyU], length between velocity points (through corners)
DO j=1-OLy+1,sNy+OLy ! NOTE range
DO i=1-OLx+1,sNx+OLx ! NOTE range
C by averaging (method I)
dxV(i,j,bi,bj) = 0.5 _d 0*(dxG(i,j,bi,bj)+dxG(i-1,j,bi,bj))
dyU(i,j,bi,bj) = 0.5 _d 0*(dyG(i,j,bi,bj)+dyG(i,j-1,bi,bj))
C by averaging (method II)
c dxV(i,j,bi,bj) = 0.5*(dxG(i,j,bi,bj)+dxG(i-1,j,bi,bj))
c dyU(i,j,bi,bj) = 0.5*(dyC(i,j,bi,bj)+dyC(i-1,j,bi,bj))
ENDDO
ENDDO
C-- Calculate vertical face area (tracer cells)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
lat=0.5 _d 0*(yGloc(i,j)+yGloc(i+1,j))
dlon = delXloc(i)
dlat = delYloc(j)
rA(i,j,bi,bj) = rSphere*rSphere*dlon*deg2rad
& *ABS( SIN((lat+dlat)*deg2rad)-SIN(lat*deg2rad) )
#ifdef USE_BACKWARD_COMPATIBLE_GRID
lat = yC(i,j,bi,bj) - delYloc(j)*0.5 _d 0
dlat= yC(i,j,bi,bj) + delYloc(j)*0.5 _d 0
rA(i,j,bi,bj) = dyF(i,j,bi,bj)
& *rSphere*(SIN(dlat*deg2rad)-SIN(lat*deg2rad))
#endif /* USE_BACKWARD_COMPATIBLE_GRID */
ENDDO
ENDDO
C-- Calculate vertical face area (u cells)
DO j=1-OLy,sNy+OLy
DO i=1-OLx+1,sNx+OLx ! NOTE range
C by averaging
rAw(i,j,bi,bj) = 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj))
C by formula
c lat=yGloc(i,j)
c dlon = 0.5*( delXloc(i) + delXloc(i-1) )
c dlat = delYloc(j)
c rAw(i,j,bi,bj) = rSphere*rSphere*dlon*deg2rad
c & *abs( sin((lat+dlat)*deg2rad)-sin(lat*deg2rad) )
ENDDO
ENDDO
C-- Calculate vertical face area (v cells)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
C by formula
lat=yC(i,j,bi,bj)
dlon = delXloc(i)
dlat = 0.5 _d 0*( delYloc(j) + delYloc(j-1) )
#ifdef USE_BACKWARD_COMPATIBLE_GRID
dlat= delYloc(j)
#endif /* USE_BACKWARD_COMPATIBLE_GRID */
rAs(i,j,bi,bj) = rSphere*rSphere*dlon*deg2rad
& *ABS( SIN(lat*deg2rad)-SIN((lat-dlat)*deg2rad) )
IF (ABS(lat).GT.90..OR.ABS(lat-dlat).GT.90.) rAs(i,j,bi,bj)=0.
ENDDO
ENDDO
C-- Calculate vertical face area (vorticity points)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
C by formula
lat = 0.5 _d 0*(yGloc(i,j)+yGloc(i,j+1))
dlon = 0.5 _d 0*( delXloc(i) + delXloc(i-1) )
dlat = 0.5 _d 0*( delYloc(j) + delYloc(j-1) )
rAz(i,j,bi,bj) = rSphere*rSphere*dlon*deg2rad
& *ABS( SIN(lat*deg2rad)-SIN((lat-dlat)*deg2rad) )
IF (ABS(lat).GT.90..OR.ABS(lat-dlat).GT.90.) rAz(i,j,bi,bj)=0.
ENDDO
ENDDO
C-- Calculate trigonometric terms & grid orientation:
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
lat=0.5 _d 0*(yGloc(i,j)+yGloc(i,j+1))
tanPhiAtU(i,j,bi,bj)=TAN(lat*deg2rad)
lat=0.5 _d 0*(yGloc(i,j)+yGloc(i+1,j))
tanPhiAtV(i,j,bi,bj)=TAN(lat*deg2rad)
C Note: this is now done earlier in main S/R INI_GRID
c angleCosC(i,j,bi,bj) = 1.
c angleSinC(i,j,bi,bj) = 0.
ENDDO
ENDDO
C-- Cosine(lat) scaling
DO j=1-OLy,sNy+OLy
i = 1
IF (cosPower.NE.0.) THEN
lat = 0.5 _d 0*(yGloc(i,j)+yGloc(i,j+1))
cosFacU(j,bi,bj) = ABS( COS(lat*deg2rad) )**cosPower
lat = 0.5 _d 0*(yGloc(i,j)+yGloc(i+1,j))
cosFacV(j,bi,bj) = ABS( COS(lat*deg2rad) )**cosPower
sqcosFacU(j,bi,bj) = SQRT(cosFacU(j,bi,bj))
sqcosFacV(j,bi,bj) = SQRT(cosFacV(j,bi,bj))
ELSE
cosFacU(j,bi,bj) = 1.
cosFacV(j,bi,bj) = 1.
sqcosFacU(j,bi,bj)=1.
sqcosFacV(j,bi,bj)=1.
ENDIF
ENDDO
C-- end bi,bj loops
ENDDO
ENDDO
IF ( rotateGrid ) THEN
CALL ROTATE_SPHERICAL_POLAR_GRID( xC, yC, myThid )
CALL ROTATE_SPHERICAL_POLAR_GRID( xG, yG, myThid )
skipCalcAngleC = .FALSE.
CALL CALC_GRID_ANGLES( skipCalcAngleC, myThid )
ENDIF
RETURN
END