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

Contents of /MITgcm/model/src/ini_cartesian_grid.F

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


Revision 1.1 - (show annotations) (download)
Wed Apr 22 19:15:30 1998 UTC (26 years, 1 month ago) by cnh
Branch: MAIN
Branch point for: cnh
Initial revision

1 C $Id$
2
3 #include "CPP_EEOPTIONS.h"
4
5 CStartOfInterface
6 SUBROUTINE INI_CARTESIAN_GRID( myThid )
7 C /==========================================================\
8 C | SUBROUTINE INI_CARTESIAN_GRID |
9 C | o Initialise model coordinate system |
10 C |==========================================================|
11 C | These arrays are used throughout the code in evaluating |
12 C | gradients, integrals and spatial avarages. This routine |
13 C | is called separately by each thread and initialise only |
14 C | the region of the domain it is "responsible" for. |
15 C | Notes: |
16 C | Two examples are included. One illustrates the |
17 C | initialisation of a cartesian grid. The other shows the |
18 C | inialisation of a spherical polar grid. Other orthonormal|
19 C | grids can be fitted into this design. In this case |
20 C | custom metric terms also need adding to account for the |
21 C | projections of velocity vectors onto these grids. |
22 C | The structure used here also makes it possible to |
23 C | implement less regular grid mappings. In particular |
24 C | o Schemes which leave out blocks of the domain that are |
25 C | all land could be supported. |
26 C | o Multi-level schemes such as icosohedral or cubic |
27 C | grid projections onto a sphere can also be fitted |
28 C | within the strategy we use. |
29 C | Both of the above also require modifying the support |
30 C | routines that map computational blocks to simulation |
31 C | domain blocks. |
32 C | Under the cartesian grid mode primitive distances in X |
33 C | and Y are in metres. Disktance in Z are in m or Pa |
34 C | depending on the vertical gridding mode. |
35 C \==========================================================/
36
37 C === Global variables ===
38 #include "SIZE.h"
39 #include "EEPARAMS.h"
40 #include "PARAMS.h"
41 #include "GRID.h"
42
43 C == Routine arguments ==
44 C myThid - Number of this instance of INI_CARTESIAN_GRID
45 INTEGER myThid
46 CEndOfInterface
47
48 C == Local variables ==
49 C xG, yG - Global coordinate location.
50 C zG
51 C xBase - South-west corner location for process.
52 C yBase
53 C zUpper - Work arrays for upper and lower
54 C zLower cell-face heights.
55 C phi - Temporary scalar
56 C xBase - Temporaries for lower corner coordinate
57 C yBase
58 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 I,J,K
64 _RL xG, yG, zG
65 _RL phi
66 _RL zUpper(Nz), zLower(Nz)
67 _RL xBase, yBase
68 INTEGER iG, jG
69 INTEGER bi, bj
70 INTEGER I, J, K
71
72 C-- Simple example of inialisation on cartesian grid
73 C-- First set coordinates of cell centers
74 C This operation is only performed at start up so for more
75 C complex configurations it is usually OK to pass iG, jG to a custom
76 C function and have it return xG and yG.
77 C Set up my local grid first
78 DO bj = myByLo(myThid), myByHi(myThid)
79 jG = myYGlobalLo + (bj-1)*sNy
80 DO bi = myBxLo(myThid), myBxHi(myThid)
81 iG = myXGlobalLo + (bi-1)*sNx
82 yBase = 0. _d 0
83 xBase = 0. _d 0
84 DO i=1,iG-1
85 xBase = xBase + delX(i)
86 ENDDO
87 DO j=1,jG-1
88 yBase = yBase + delY(j)
89 ENDDO
90 yG = yBase
91 DO J=1,sNy
92 xG = xBase
93 DO I=1,sNx
94 xc(I,J,bi,bj) = xG + delX(iG+i-1)*0.5 _d 0
95 yc(I,J,bi,bj) = yG + delY(jG+j-1)*0.5 _d 0
96 xG = xG + delX(iG+I-1)
97 dxF(I,J,bi,bj) = delX(iG+i-1)
98 dyF(I,J,bi,bj) = delY(jG+j-1)
99 ENDDO
100 yG = yG + delY(jG+J-1)
101 ENDDO
102 ENDDO
103 ENDDO
104 C Now sync. and get edge regions from other threads and/or processes.
105 C Note: We could just set the overlap regions ourselves here but
106 C exchanging edges is safer and is good practice!
107 _EXCH_XY_R4( xc, myThid )
108 _EXCH_XY_R4( yc, myThid )
109 _EXCH_XY_R4(dxF, myThid )
110 _EXCH_XY_R4(dyF, myThid )
111
112 C-- Calculate separation between other points
113 C dxG, dyG are separations between cell corners along cell faces.
114 DO bj = myByLo(myThid), myByHi(myThid)
115 DO bi = myBxLo(myThid), myBxHi(myThid)
116 DO J=1,sNy
117 DO I=1,sNx
118 dxG(I,J,bi,bj) = (dxF(I,J,bi,bj)+dxF(I,J-1,bi,bj))*0.5 _d 0
119 dyG(I,J,bi,bj) = (dyF(I,J,bi,bj)+dyF(I-1,J,bi,bj))*0.5 _d 0
120 ENDDO
121 ENDDO
122 ENDDO
123 ENDDO
124 _EXCH_XY_R4(dxG, myThid )
125 _EXCH_XY_R4(dyG, myThid )
126 C dxV, dyU are separations between velocity points along cell faces.
127 DO bj = myByLo(myThid), myByHi(myThid)
128 DO bi = myBxLo(myThid), myBxHi(myThid)
129 DO J=1,sNy
130 DO I=1,sNx
131 dxV(I,J,bi,bj) = (dxG(I,J,bi,bj)+dxG(I-1,J,bi,bj))*0.5 _d 0
132 dyU(I,J,bi,bj) = (dyG(I,J,bi,bj)+dyG(I,J-1,bi,bj))*0.5 _d 0
133 ENDDO
134 ENDDO
135 ENDDO
136 ENDDO
137 _EXCH_XY_R4(dxV, myThid )
138 _EXCH_XY_R4(dyU, myThid )
139 C dxC, dyC is separation between cell centers
140 DO bj = myByLo(myThid), myByHi(myThid)
141 DO bi = myBxLo(myThid), myBxHi(myThid)
142 DO J=1,sNy
143 DO I=1,sNx
144 dxC(I,J,bi,bj) = (dxF(I,J,bi,bj)+dxF(I-1,J,bi,bj))*0.5 D0
145 dyC(I,J,bi,bj) = (dyF(I,J,bi,bj)+dyF(I,J-1,bi,bj))*0.5 D0
146 ENDDO
147 ENDDO
148 ENDDO
149 ENDDO
150 _EXCH_XY_R4(dxC, myThid )
151 _EXCH_XY_R4(dyC, myThid )
152 C Calculate recipricols
153 DO bj = myByLo(myThid), myByHi(myThid)
154 DO bi = myBxLo(myThid), myBxHi(myThid)
155 DO J=1,sNy
156 DO I=1,sNx
157 rDxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj)
158 rDyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj)
159 rDxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj)
160 rDyC(I,J,bi,bj)=1.d0/dyC(I,J,bi,bj)
161 rDxF(I,J,bi,bj)=1.d0/dxF(I,J,bi,bj)
162 rDyF(I,J,bi,bj)=1.d0/dyF(I,J,bi,bj)
163 rDxV(I,J,bi,bj)=1.d0/dxV(I,J,bi,bj)
164 rDyU(I,J,bi,bj)=1.d0/dyU(I,J,bi,bj)
165 ENDDO
166 ENDDO
167 ENDDO
168 ENDDO
169 _EXCH_XY_R4(rDxG, myThid )
170 _EXCH_XY_R4(rDyG, myThid )
171 _EXCH_XY_R4(rDxC, myThid )
172 _EXCH_XY_R4(rDyC, myThid )
173 _EXCH_XY_R4(rDxF, myThid )
174 _EXCH_XY_R4(rDyF, myThid )
175 _EXCH_XY_R4(rDxV, myThid )
176 _EXCH_XY_R4(rDyU, myThid )
177 C Calculate vertical face area
178 DO bj = myByLo(myThid), myByHi(myThid)
179 DO bi = myBxLo(myThid), myBxHi(myThid)
180 DO J=1,sNy
181 DO I=1,sNx
182 zA(I,J,bi,bj) = dxF(I,J,bi,bj)*dyF(I,J,bi,bj)
183 ENDDO
184 ENDDO
185 ENDDO
186 ENDDO
187
188 DO bj = myByLo(myThid), myByHi(myThid)
189 DO bi = myBxLo(myThid), myBxHi(myThid)
190 DO K=1,Nz
191 DO J=1,sNy
192 DO I=1,sNx
193 IF (HFacC(I,J,K,bi,bj) .NE. 0. D0 ) THEN
194 rHFacC(I,J,K,bi,bj) = 1. D0 / HFacC(I,J,K,bi,bj)
195 ELSE
196 rHFacC(I,J,K,bi,bj) = 0. D0
197 ENDIF
198 IF (HFacW(I,J,K,bi,bj) .NE. 0. D0 ) THEN
199 rHFacW(I,J,K,bi,bj) = 1. D0 / HFacW(I,J,K,bi,bj)
200 maskW(I,J,K,bi,bj) = 1. D0
201 ELSE
202 rHFacW(I,J,K,bi,bj) = 0. D0
203 maskW(I,J,K,bi,bj) = 0.0 D0
204 ENDIF
205 IF (HFacS(I,J,K,bi,bj) .NE. 0. D0 ) THEN
206 rHFacS(I,J,K,bi,bj) = 1. D0 / HFacS(I,J,K,bi,bj)
207 maskS(I,J,K,bi,bj) = 1. D0
208 ELSE
209 rHFacS(I,J,K,bi,bj) = 0. D0
210 maskS(I,J,K,bi,bj) = 0. D0
211 ENDIF
212 ENDDO
213 ENDDO
214 ENDDO
215 ENDDO
216 ENDDO
217 C Now sync. and get/send edge regions that are shared with
218 C other threads.
219 _EXCH_XYZ_R4(rHFacC , myThid )
220 _EXCH_XYZ_R4(rHFacW , myThid )
221 _EXCH_XYZ_R4(rHFacS , myThid )
222 _EXCH_XYZ_R4(maskW , myThid )
223 _EXCH_XYZ_R4(maskS , myThid )
224
225 C
226 RETURN
227 END
228
229 C $Id: $

  ViewVC Help
Powered by ViewVC 1.1.22