/[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.13 - (show annotations) (download)
Fri Apr 2 00:53:18 2004 UTC (20 years, 1 month ago) by dimitri
Branch: MAIN
CVS Tags: checkpoint53, checkpoint53d_post, checkpoint52m_post, checkpoint53c_post, checkpoint53a_post, checkpoint52n_post, checkpoint53b_pre, checkpoint53b_post, checkpoint53d_pre
Changes since 1.12: +16 -3 lines
added IF(useCubedSphereExchange)THEN's for hacks that seem specific to cube

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

  ViewVC Help
Powered by ViewVC 1.1.22