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

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

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


Revision 1.10 - (show annotations) (download)
Mon Mar 8 21:20:02 2004 UTC (20 years, 2 months ago) by adcroft
Branch: MAIN
Changes since 1.9: +14 -14 lines
Renamed "USE_W2" to "ALLOW_EXCH2" so that it is no longer necessary to edit
CPP_EEOPTION.h as well as packages.conf to turn on/off exch2.
 - you can control the use of exch2 through packages.conf or -enable/disable.
 + need to add a run-time flag for this

1 C $Header: /u/gcmpack/MITgcm/model/src/ini_curvilinear_grid.F,v 1.9 2004/02/25 01:44:00 jmc Exp $
2 C $Name: $
3
4 #include "CPP_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: INI_CURVILINEAR_GRID
8 C !INTERFACE:
9 SUBROUTINE INI_CURVILINEAR_GRID( myThid )
10 C !DESCRIPTION: \bv
11 C *==========================================================*
12 C | SUBROUTINE INI_CURVILINEAR_GRID
13 C | o Initialise curvilinear coordinate system
14 C *==========================================================*
15 C | Curvilinear grid settings are read from a file rather
16 C | than coded in-line as for cartesian and spherical polar.
17 C | This is more general but you have to create the grid
18 C | yourself.
19 C *==========================================================*
20 C \ev
21
22 C !USES:
23 IMPLICIT NONE
24 C === Global variables ===
25 #include "SIZE.h"
26 #include "EEPARAMS.h"
27 #include "PARAMS.h"
28 #include "GRID.h"
29 #ifdef ALLOW_EXCH2
30 #include "W2_EXCH2_TOPOLOGY.h"
31 #include "W2_EXCH2_PARAMS.h"
32 #endif
33
34 #ifndef ALLOW_EXCH2
35 C- note: default is to use "new" grid files (OLD_GRID_IO undef) with EXCH2
36 C but can still use (on 1 cpu) OLD_GRID_IO and EXCH2 independently
37 #define OLD_GRID_IO
38 #endif /* ALLOW_EXCH2 */
39
40 C !INPUT/OUTPUT PARAMETERS:
41 C == Routine arguments ==
42 C myThid - Number of this instance of INI_CARTESIAN_GRID
43 INTEGER myThid
44
45 C !LOCAL VARIABLES:
46 C == Local variables ==
47 INTEGER bi,bj, myTile
48 INTEGER I,J
49 CHARACTER*(15) fName
50 _RL buf(sNx+1,sNy+1)
51 INTEGER iG, iL
52 CHARACTER*(MAX_LEN_MBUF) msgBuf
53 INTEGER ILNBLNK
54 EXTERNAL ILNBLNK
55 CEOP
56
57 C-- Set everything to zero everywhere
58 DO bj = myByLo(myThid), myByHi(myThid)
59 DO bi = myBxLo(myThid), myBxHi(myThid)
60
61 DO J=1-Oly,sNy+Oly
62 DO I=1-Olx,sNx+Olx
63 XC(i,j,bi,bj)=0.
64 YC(i,j,bi,bj)=0.
65 XG(i,j,bi,bj)=0.
66 YG(i,j,bi,bj)=0.
67 DXC(i,j,bi,bj)=0.
68 DYC(i,j,bi,bj)=0.
69 DXG(i,j,bi,bj)=0.
70 DYG(i,j,bi,bj)=0.
71 DXF(i,j,bi,bj)=0.
72 DYF(i,j,bi,bj)=0.
73 DXV(i,j,bi,bj)=0.
74 DYU(i,j,bi,bj)=0.
75 RA(i,j,bi,bj)=0.
76 RAZ(i,j,bi,bj)=0.
77 RAW(i,j,bi,bj)=0.
78 RAS(i,j,bi,bj)=0.
79 tanPhiAtU(i,j,bi,bj)=0.
80 tanPhiAtV(i,j,bi,bj)=0.
81 cosFacU(J,bi,bj)=1.
82 cosFacV(J,bi,bj)=1.
83 sqcosFacU(J,bi,bj)=1.
84 sqcosFacV(J,bi,bj)=1.
85 ENDDO
86 ENDDO
87
88 ENDDO
89 ENDDO
90
91 C Here we make no assumptions about grid symmetry and simply
92 C read the raw grid data from files
93
94 #ifdef OLD_GRID_IO
95
96 C- Cell centered quantities
97 CALL MDSREADFIELD('LONC.bin',readBinaryPrec,'RS',1,XC, 1,myThid)
98 CALL MDSREADFIELD('LATC.bin',readBinaryPrec,'RS',1,YC, 1,myThid)
99 _EXCH_XY_R4(XC,myThid)
100 _EXCH_XY_R4(YC,myThid)
101
102 CALL MDSREADFIELD('DXF.bin',readBinaryPrec,'RS',1,DXF, 1,myThid)
103 CALL MDSREADFIELD('DYF.bin',readBinaryPrec,'RS',1,DYF, 1,myThid)
104 C !!! _EXCH_OUV_XY_R4(DXF, DYF, unSigned, myThid )
105 cs! this is not correct! <= need paired exchange for DXF,DYF
106 _EXCH_XY_R4(DXF,myThid)
107 _EXCH_XY_R4(DYF,myThid)
108 cs! fix overlaps:
109 DO bj = myByLo(myThid), myByHi(myThid)
110 DO bi = myBxLo(myThid), myBxHi(myThid)
111 DO j=1,sNy
112 DO i=1,Olx
113 DXF(1-i,j,bi,bj)=DXF(i,j,bi,bj)
114 DXF(sNx+i,j,bi,bj)=DXF(sNx+1-i,j,bi,bj)
115 DYF(1-i,j,bi,bj)=DYF(i,j,bi,bj)
116 DYF(sNx+i,j,bi,bj)=DYF(sNx+1-i,j,bi,bj)
117 ENDDO
118 ENDDO
119 DO j=1,Oly
120 DO i=1,sNx
121 DXF(i,1-j,bi,bj)=DXF(i,j,bi,bj)
122 DXF(i,sNy+j,bi,bj)=DXF(i,sNy+1-j,bi,bj)
123 DYF(i,1-j,bi,bj)=DYF(i,j,bi,bj)
124 DYF(i,sNy+j,bi,bj)=DYF(i,sNy+1-j,bi,bj)
125 ENDDO
126 ENDDO
127 ENDDO
128 ENDDO
129 cs
130
131 CALL MDSREADFIELD('RA.bin',readBinaryPrec,'RS',1,RA, 1,myThid)
132 _EXCH_XY_R4(RA,myThid )
133
134 C- Corner quantities
135 C *********** this are not degbugged ************
136 CALL MDSREADFIELD('LONG.bin',readBinaryPrec,'RS',1,XG, 1,myThid)
137 CALL MDSREADFIELD('LATG.bin',readBinaryPrec,'RS',1,YG, 1,myThid)
138 cs- this block needed by cubed sphere until we write more useful I/O routines
139 bi=3
140 bj=1
141 YG(1,sNy+1,bj,1)=YG(1,1,bi,1)
142 bj=bj+2
143 YG(1,sNy+1,bj,1)=YG(1,1,bi,1)
144 bj=bj+2
145 YG(1,sNy+1,bj,1)=YG(1,1,bi,1)
146 bi=6
147 bj=2
148 YG(sNx+1,1,bj,1)=YG(1,1,bi,1)
149 bj=bj+2
150 YG(sNx+1,1,bj,1)=YG(1,1,bi,1)
151 bj=bj+2
152 YG(sNx+1,1,bj,1)=YG(1,1,bi,1)
153 cs- end block
154 CALL EXCH_Z_XY_RS(XG,myThid)
155 CALL EXCH_Z_XY_RS(YG,myThid)
156
157 CALL MDSREADFIELD('DXV.bin',readBinaryPrec,'RS',1,DXV, 1,myThid)
158 CALL MDSREADFIELD('DYU.bin',readBinaryPrec,'RS',1,DYU, 1,myThid)
159 cs- this block needed by cubed sphere until we write more useful I/O routines
160 C !!! _EXCH_ZUV_XY_R4(DXV, DYU, unSigned, myThid)
161 cs! this is not correct <= need paired exchange for dxv,dyu
162 cs CALL EXCH_Z_XY_RS(DXV,myThid)
163 cs CALL EXCH_Z_XY_RS(DYU,myThid)
164 DO bj = myByLo(myThid), myByHi(myThid)
165 DO bi = myBxLo(myThid), myBxHi(myThid)
166 cs! fix overlaps:
167 DO j=1,sNy
168 DO i=1,Olx
169 DXV(1-i,j,bi,bj)=DXV(1+i,j,bi,bj)
170 DXV(sNx+i,j,bi,bj)=DXV(i,j,bi,bj)
171 DYU(1-i,j,bi,bj)=DYU(1+i,j,bi,bj)
172 DYU(sNx+i,j,bi,bj)=DYU(i,j,bi,bj)
173 ENDDO
174 ENDDO
175 DO j=1,Oly
176 DO i=1-Olx,sNx+Olx
177 DXV(i,1-j,bi,bj)=DXV(i,1+j,bi,bj)
178 DXV(i,sNy+j,bi,bj)=DXV(i,j,bi,bj)
179 DYU(i,1-j,bi,bj)=DYU(i,1+j,bi,bj)
180 DYU(i,sNy+j,bi,bj)=DYU(i,j,bi,bj)
181 ENDDO
182 ENDDO
183 ENDDO
184 ENDDO
185 cs- end block
186
187 CALL MDSREADFIELD('RAZ.bin',readBinaryPrec,'RS',1,RAZ, 1,myThid)
188 cs- this block needed by cubed sphere until we write more useful I/O routines
189 CALL EXCH_Z_XY_RS(RAZ , myThid )
190 DO bj = myByLo(myThid), myByHi(myThid)
191 DO bi = myBxLo(myThid), myBxHi(myThid)
192 RAZ(sNx+1,1,bi,bj)=RAZ(1,1,bi,bj)
193 RAZ(1,sNy+1,bi,bj)=RAZ(1,1,bi,bj)
194 ENDDO
195 ENDDO
196 cs- end block
197 CALL EXCH_Z_XY_RS(RAZ,myThid)
198
199 C- Staggered (u,v pairs) quantities
200 CALL MDSREADFIELD('DXC.bin',readBinaryPrec,'RS',1,DXC, 1,myThid)
201 CALL MDSREADFIELD('DYC.bin',readBinaryPrec,'RS',1,DYC, 1,myThid)
202 CALL EXCH_UV_XY_RS(DXC,DYC,.FALSE.,myThid)
203
204 CALL MDSREADFIELD('RAW.bin',readBinaryPrec,'RS',1,RAW, 1,myThid)
205 CALL MDSREADFIELD('RAS.bin',readBinaryPrec,'RS',1,RAS, 1,myThid)
206 cs- this block needed by cubed sphere until we write more useful I/O routines
207 DO bj = myByLo(myThid), myByHi(myThid)
208 DO bi = myBxLo(myThid), myBxHi(myThid)
209 DO J = 1,sNy
210 c RAW(sNx+1,J,bi,bj)=RAW(1,J,bi,bj)
211 c RAS(J,sNy+1,bi,bj)=RAS(J,1,bi,bj)
212 ENDDO
213 ENDDO
214 ENDDO
215 cs- end block
216 CALL EXCH_UV_XY_RS(RAW,RAS,.FALSE.,myThid)
217
218 CALL MDSREADFIELD('DXG.bin',readBinaryPrec,'RS',1,DXG, 1,myThid)
219 CALL MDSREADFIELD('DYG.bin',readBinaryPrec,'RS',1,DYG, 1,myThid)
220 cs- this block needed by cubed sphere until we write more useful I/O routines
221 DO bj = myByLo(myThid), myByHi(myThid)
222 DO bi = myBxLo(myThid), myBxHi(myThid)
223 DO J = 1,sNy
224 c DYG(sNx+1,J,bi,bj)=DYG(1,J,bi,bj)
225 c DXG(J,sNy+1,bi,bj)=DXG(J,1,bi,bj)
226 ENDDO
227 ENDDO
228 ENDDO
229 cs- end block
230 CALL EXCH_UV_XY_RS(DYG,DXG,.FALSE.,myThid)
231
232 c write(10) XC
233 c write(10) YC
234 c write(10) DXF
235 c write(10) DYF
236 c write(10) RA
237 c write(10) XG
238 c write(10) YG
239 c write(10) DXV
240 c write(10) DYU
241 c write(10) RAZ
242 c write(10) DXC
243 c write(10) DYC
244 c write(10) RAW
245 c write(10) RAS
246 c write(10) DXG
247 c write(10) DYG
248
249 #else
250
251 C-- Only do I/O if I am the master thread
252 _BEGIN_MASTER(myThid)
253
254 DO bj = 1,nSy
255 DO bi = 1,nSx
256 iG=bi+(myXGlobalLo-1)/sNx
257 WRITE(fName(1:15),'("tile",I3.3,".mitgrid")') iG
258 WRITE(msgBuf,'(A,I4)') 'tile:',iG
259 #ifdef ALLOW_EXCH2
260 myTile = W2_myTileList(bi)
261 write(fName(1:15),'("tile",I3.3,".mitgrid")')
262 & exch2_myface(myTile)
263 WRITE(msgBuf,'(A,I4)') 'tile:',myTile
264 #endif
265 iL = ILNBLNK(msgBuf)
266 WRITE(msgBuf,'(3A)') msgBuf(1:iL),
267 & ' ; Read from file ',fName(1:15)
268 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
269 & SQUEEZE_RIGHT , myThid)
270 WRITE(msgBuf,'(A)') ' =>'
271
272 CALL READSYMTILE_RS(fName,1,XC,bi,bj,buf,myThid)
273 iL = ILNBLNK(msgBuf)
274 WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'XC'
275 CALL READSYMTILE_RS(fName,2,YC,bi,bj,buf,myThid)
276 iL = ILNBLNK(msgBuf)
277 WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'YC'
278 CALL READSYMTILE_RS(fName,3,DXF,bi,bj,buf,myThid)
279 iL = ILNBLNK(msgBuf)
280 WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'DXF'
281 CALL READSYMTILE_RS(fName,4,DYF,bi,bj,buf,myThid)
282 iL = ILNBLNK(msgBuf)
283 WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'DYF'
284 CALL READSYMTILE_RS(fName,5,RA,bi,bj,buf,myThid)
285 iL = ILNBLNK(msgBuf)
286 WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'RA'
287 CALL READSYMTILE_RS(fName,6,XG,bi,bj,buf,myThid)
288 iL = ILNBLNK(msgBuf)
289 WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'XG'
290 CALL READSYMTILE_RS(fName,7,YG,bi,bj,buf,myThid)
291 iL = ILNBLNK(msgBuf)
292 WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'YG'
293 CALL READSYMTILE_RS(fName,8,DXV,bi,bj,buf,myThid)
294 iL = ILNBLNK(msgBuf)
295 WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'DXV'
296 CALL READSYMTILE_RS(fName,9,DYU,bi,bj,buf,myThid)
297 iL = ILNBLNK(msgBuf)
298 WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'DYU'
299 CALL READSYMTILE_RS(fName,10,RAZ,bi,bj,buf,myThid)
300 iL = ILNBLNK(msgBuf)
301 WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'RAZ'
302 CALL READSYMTILE_RS(fName,11,DXC,bi,bj,buf,myThid)
303 iL = ILNBLNK(msgBuf)
304 WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'DXC'
305 CALL READSYMTILE_RS(fName,12,DYC,bi,bj,buf,myThid)
306 iL = ILNBLNK(msgBuf)
307 WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'DYC'
308 CALL READSYMTILE_RS(fName,13,RAW,bi,bj,buf,myThid)
309 iL = ILNBLNK(msgBuf)
310 WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'RAW'
311 CALL READSYMTILE_RS(fName,14,RAS,bi,bj,buf,myThid)
312 iL = ILNBLNK(msgBuf)
313 WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'RAS'
314 CALL READSYMTILE_RS(fName,15,DXG,bi,bj,buf,myThid)
315 iL = ILNBLNK(msgBuf)
316 WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'DXG'
317 CALL READSYMTILE_RS(fName,16,DYG,bi,bj,buf,myThid)
318 iL = ILNBLNK(msgBuf)
319 WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'DYG'
320
321 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
322 & SQUEEZE_RIGHT , myThid)
323
324 ENDDO
325 ENDDO
326 _END_MASTER(myThid)
327
328 CALL EXCH_XY_RS(XC,myThid)
329 CALL EXCH_XY_RS(YC,myThid)
330 C !!! _EXCH_OUV_XY_R4(DXF, DYF, unSigned, myThid )
331 c CALL EXCH_XY_RS(DXF,myThid)
332 c CALL EXCH_XY_RS(DYF,myThid)
333 CALL EXCH_XY_RS(RA,myThid )
334 #ifndef ALLOW_EXCH2
335 CALL EXCH_Z_XY_RS(XG,myThid)
336 CALL EXCH_Z_XY_RS(YG,myThid)
337 C !!! _EXCH_ZUV_XY_R4(DXV, DYU, unSigned, myThid)
338 c CALL EXCH_Z_XY_RS(DXV,myThid)
339 c CALL EXCH_Z_XY_RS(DYU,myThid)
340 CALL EXCH_Z_XY_RS(RAZ,myThid)
341 #endif /* ALLOW_EXCH2 */
342 CALL EXCH_UV_XY_RS(DXC,DYC,.FALSE.,myThid)
343 CALL EXCH_UV_XY_RS(RAW,RAS,.FALSE.,myThid)
344 CALL EXCH_UV_XY_RS(DYG,DXG,.FALSE.,myThid)
345
346 #endif
347
348 c CALL WRITE_FULLARRAY_RL('DXV',DXV,1,0,myThid)
349 c CALL WRITE_FULLARRAY_RL('DYU',DYU,1,0,myThid)
350 c CALL WRITE_FULLARRAY_RL('RAZ',RAZ,1,0,myThid)
351 c CALL WRITE_FULLARRAY_RL('XG',XG,1,0,myThid)
352 c CALL WRITE_FULLARRAY_RL('YG',YG,1,0,myThid)
353
354
355 RETURN
356 END
357
358 C --------------------------------------------------------------------------
359
360 SUBROUTINE READSYMTILE_RS(fName,irec,array,bi,bj,buf,myThid)
361 C /==========================================================\
362 C | SUBROUTINE READSYMTILE_RS |
363 C |==========================================================|
364 C \==========================================================/
365 IMPLICIT NONE
366
367 C === Global variables ===
368 #include "SIZE.h"
369 #include "EEPARAMS.h"
370 #ifdef ALLOW_EXCH2
371 #include "W2_EXCH2_TOPOLOGY.h"
372 #include "W2_EXCH2_PARAMS.h"
373 #endif /* ALLOW_EXCH2 */
374
375 C == Routine arguments ==
376 CHARACTER*(*) fName
377 INTEGER irec
378 _RS array(1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
379 INTEGER bi,bj,myThid
380 #ifdef ALLOW_EXCH2
381 _RL buf(1:sNx*nSx*nPx+1)
382 #else
383 _RL buf(1:sNx+1,1:sNy+1)
384 #endif /* ALLOW_EXCH2 */
385
386 C == Local variables ==
387 INTEGER I,J,dUnit
388 INTEGER length_of_rec
389 INTEGER MDS_RECLEN
390 INTEGER TN, DNX, DNY, TBX, TBY, TNX, TNY, II, iBase
391
392 #ifdef ALLOW_EXCH2
393 C Figure out offset of tile within face
394 TN = W2_myTileList(bi)
395 DNX = exch2_mydnx(TN)
396 DNY = exch2_mydny(TN)
397 TBX = exch2_tbasex(TN)
398 TBY = exch2_tbasey(TN)
399 TNX = exch2_tnx(TN)
400 TNY = exch2_tny(TN)
401
402 CALL MDSFINDUNIT( dUnit, mythid )
403 length_of_rec=MDS_RECLEN( 64, (dNx+1), myThid )
404 OPEN( dUnit, file=fName, status='old',
405 & access='direct', recl=length_of_rec )
406 J=0
407 iBase=(irec-1)*(dny+1)
408 DO I=1+TBY,SNY+1+TBY
409 READ(dUnit,rec=I+iBase)(buf(ii),ii=1,dnx+1)
410 #ifdef _BYTESWAPIO
411 #ifdef REAL4_IS_SLOW
412 CALL MDS_BYTESWAPR8((dNx+1), buf)
413 #else
414 CALL MDS_BYTESWAPR4((dNx+1), buf)
415 #endif
416 #endif
417 J=J+1
418 DO II=1,sNx+1
419 array(II,J,bi,bj)=buf(II+TBX)
420 ENDDO
421 ENDDO
422 CLOSE( dUnit )
423
424 #else /* ALLOW_EXCH2 */
425
426 CALL MDSFINDUNIT( dUnit, mythid )
427 length_of_rec=MDS_RECLEN( 64, (sNx+1)*(sNy+1), myThid )
428 OPEN( dUnit, file=fName, status='old',
429 & access='direct', recl=length_of_rec )
430 READ(dUnit,rec=irec) buf
431 CLOSE( dUnit )
432
433 #ifdef _BYTESWAPIO
434 #ifdef REAL4_IS_SLOW
435 CALL MDS_BYTESWAPR8((sNx+1)*(sNy+1), buf)
436 #else
437 CALL MDS_BYTESWAPR4((sNx+1)*(sNy+1), buf)
438 #endif
439 #endif
440
441 DO J=1,sNy+1
442 DO I=1,sNx+1
443 array(I,J,bi,bj)=buf(I,J)
444 ENDDO
445 ENDDO
446 c write(0,*) irec,buf(1,1),array(1,1,1,1)
447
448 #endif /* ALLOW_EXCH2 */
449
450 RETURN
451 END

  ViewVC Help
Powered by ViewVC 1.1.22