/[MITgcm]/MITgcm_contrib/ifenty/Curvi/code/ini_curvilinear_grid.F
ViewVC logotype

Contents of /MITgcm_contrib/ifenty/Curvi/code/ini_curvilinear_grid.F

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


Revision 1.1 - (show annotations) (download)
Wed Jul 5 16:36:02 2006 UTC (16 years, 10 months ago) by ifenty
Branch: MAIN
CVS Tags: HEAD
Initial checkin of curvilinear configuration of Labrador Sea region

1 C $Header: /u/gcmpack/MITgcm/model/src/ini_curvilinear_grid.F,v 1.24 2005/09/28 01:43:55 edhill 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 #ifdef ALLOW_MNC
35 #include "MNC_PARAMS.h"
36 #endif
37
38 #ifndef ALLOW_EXCH2
39 C- note: default is to use "new" grid files (OLD_GRID_IO undef) with EXCH2
40 C but can still use (on 1 cpu) OLD_GRID_IO and EXCH2 independently
41 #ifdef ALLOW_MDSIO
42 #define OLD_GRID_IO
43 #endif
44 #endif /* ALLOW_EXCH2 */
45
46 C !INPUT/OUTPUT PARAMETERS:
47 C == Routine arguments ==
48 C myThid - Number of this instance of INI_CURVILINEAR_GRID
49 INTEGER myThid
50
51 C !LOCAL VARIABLES:
52 C == Local variables ==
53 INTEGER bi,bj, myIter
54 INTEGER I,J
55 CHARACTER*(MAX_LEN_FNAM) fName
56 #ifdef ALLOW_MNC
57 CHARACTER*(80) mncFn
58 #endif
59 #ifdef ALLOW_EXCH2
60 _RL buf(sNx*nSx*nPx+1)
61 INTEGER myTile
62 #else
63 _RL buf(sNx+1,sNy+1)
64 #endif
65 INTEGER iG, iL, iLen
66 CHARACTER*(MAX_LEN_MBUF) msgBuf, tmpBuf
67 INTEGER ILNBLNK
68 EXTERNAL ILNBLNK
69 CEOP
70
71 C-- Set everything to zero everywhere
72 DO bj = myByLo(myThid), myByHi(myThid)
73 DO bi = myBxLo(myThid), myBxHi(myThid)
74
75 DO J=1-Oly,sNy+Oly
76 DO I=1-Olx,sNx+Olx
77 XC(i,j,bi,bj)=0.
78 YC(i,j,bi,bj)=0.
79 XG(i,j,bi,bj)=0.
80 YG(i,j,bi,bj)=0.
81 DXC(i,j,bi,bj)=0.
82 DYC(i,j,bi,bj)=0.
83 DXG(i,j,bi,bj)=0.
84 DYG(i,j,bi,bj)=0.
85 DXF(i,j,bi,bj)=0.
86 DYF(i,j,bi,bj)=0.
87 DXV(i,j,bi,bj)=0.
88 DYU(i,j,bi,bj)=0.
89 RA(i,j,bi,bj)=0.
90 RAZ(i,j,bi,bj)=0.
91 RAW(i,j,bi,bj)=0.
92 RAS(i,j,bi,bj)=0.
93 tanPhiAtU(i,j,bi,bj)=0.
94 tanPhiAtV(i,j,bi,bj)=0.
95 angleCosC(i,j,bi,bj)=1.
96 angleSinC(i,j,bi,bj)=0.
97 cosFacU(J,bi,bj)=1.
98 cosFacV(J,bi,bj)=1.
99 sqcosFacU(J,bi,bj)=1.
100 sqcosFacV(J,bi,bj)=1.
101 ENDDO
102 ENDDO
103
104 ENDDO
105 ENDDO
106
107
108 #ifdef ALLOW_MNC
109 IF (useMNC .AND. readgrid_mnc) THEN
110
111 _BEGIN_MASTER(myThid)
112 DO i = 1,80
113 mncFn(i:i) = ' '
114 ENDDO
115 write(mncFn,'(a)') 'mitgrid'
116 DO i = 1,MAX_LEN_MBUF
117 msgBuf(i:i) = ' '
118 ENDDO
119 WRITE(msgBuf,'(2A)') msgBuf,' ; Reading grid info using MNC'
120 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
121 & SQUEEZE_RIGHT , myThid)
122 CALL MNC_FILE_CLOSE_ALL_MATCHING(mncFn, myThid)
123 CALL MNC_CW_SET_UDIM(mncFn, 1, myThid)
124 CALL MNC_CW_SET_CITER(mncFn, 2, -1, -1, -1, myThid)
125 CALL MNC_CW_SET_UDIM(mncFn, 1, myThid)
126 CALL MNC_CW_RS_R('D',mncFn,0,0,'XC', XC, myThid)
127 CALL MNC_CW_RS_R('D',mncFn,0,0,'XG', XG, myThid)
128 CALL MNC_CW_RS_R('D',mncFn,0,0,'YC', YC, myThid)
129 CALL MNC_CW_RS_R('D',mncFn,0,0,'YG', YG, myThid)
130 CALL MNC_CW_RS_R('D',mncFn,0,0,'dxC',DXC, myThid)
131 CALL MNC_CW_RS_R('D',mncFn,0,0,'dyC',DYC, myThid)
132 CALL MNC_CW_RS_R('D',mncFn,0,0,'dxF',DXF, myThid)
133 CALL MNC_CW_RS_R('D',mncFn,0,0,'dyF',DYF, myThid)
134 CALL MNC_CW_RS_R('D',mncFn,0,0,'dxG',DXG, myThid)
135 CALL MNC_CW_RS_R('D',mncFn,0,0,'dyG',DYG, myThid)
136 CALL MNC_CW_RS_R('D',mncFn,0,0,'dxV',DXV, myThid)
137 CALL MNC_CW_RS_R('D',mncFn,0,0,'dyU',DYU, myThid)
138 CALL MNC_CW_RS_R('D',mncFn,0,0,'rA', RA, myThid)
139 CALL MNC_CW_RS_R('D',mncFn,0,0,'rAz',RAZ, myThid)
140 CALL MNC_CW_RS_R('D',mncFn,0,0,'rAw',RAW, myThid)
141 CALL MNC_CW_RS_R('D',mncFn,0,0,'rAs',RAS, myThid)
142 CALL MNC_CW_RS_R('D',mncFn,0,0,'AngleCS',angleCosC,myThid)
143 CALL MNC_CW_RS_R('D',mncFn,0,0,'AngleSN',angleSinC,myThid)
144
145 _END_MASTER(myThid)
146
147 CALL EXCH_XY_RS(XC,myThid)
148 CALL EXCH_XY_RS(YC,myThid)
149 #ifdef HRCUBE
150 CALL EXCH_XY_RS(DXF,myThid)
151 CALL EXCH_XY_RS(DYF,myThid)
152 #endif
153 CALL EXCH_XY_RS(RA,myThid )
154 CALL EXCH_Z_XY_RS(XG,myThid)
155 CALL EXCH_Z_XY_RS(YG,myThid)
156 CALL EXCH_Z_XY_RS(RAZ,myThid)
157 CALL EXCH_UV_XY_RS(DXC,DYC,.FALSE.,myThid)
158 CALL EXCH_UV_XY_RS(RAW,RAS,.FALSE.,myThid)
159 CALL EXCH_UV_XY_RS(DYG,DXG,.FALSE.,myThid)
160 CALL EXCH_UV_AGRID_XY_RS(angleSinC,angleCosC,.TRUE.,myThid)
161
162 ELSE
163 #endif
164
165 C Here we make no assumptions about grid symmetry and simply
166 C read the raw grid data from files
167
168 #ifdef MONKEY_BUTT_OLD_GRID_IO
169
170 C- Cell centered quantities
171 CALL MDSREADFIELD('LONC.bin',readBinaryPrec,'RS',1,XC, 1,myThid)
172 CALL MDSREADFIELD('LATC.bin',readBinaryPrec,'RS',1,YC, 1,myThid)
173 _EXCH_XY_R4(XC,myThid)
174 _EXCH_XY_R4(YC,myThid)
175
176 CALL MDSREADFIELD('DXF.bin',readBinaryPrec,'RS',1,DXF, 1,myThid)
177 CALL MDSREADFIELD('DYF.bin',readBinaryPrec,'RS',1,DYF, 1,myThid)
178 C !!! _EXCH_OUV_XY_R4(DXF, DYF, unSigned, myThid )
179 cs! this is not correct! <= need paired exchange for DXF,DYF
180 _EXCH_XY_R4(DXF,myThid)
181 _EXCH_XY_R4(DYF,myThid)
182 IF (useCubedSphereExchange) THEN
183 cs! fix overlaps:
184 DO bj = myByLo(myThid), myByHi(myThid)
185 DO bi = myBxLo(myThid), myBxHi(myThid)
186 DO j=1,sNy
187 DO i=1,Olx
188 DXF(1-i,j,bi,bj)=DXF(i,j,bi,bj)
189 DXF(sNx+i,j,bi,bj)=DXF(sNx+1-i,j,bi,bj)
190 DYF(1-i,j,bi,bj)=DYF(i,j,bi,bj)
191 DYF(sNx+i,j,bi,bj)=DYF(sNx+1-i,j,bi,bj)
192 ENDDO
193 ENDDO
194 DO j=1,Oly
195 DO i=1,sNx
196 DXF(i,1-j,bi,bj)=DXF(i,j,bi,bj)
197 DXF(i,sNy+j,bi,bj)=DXF(i,sNy+1-j,bi,bj)
198 DYF(i,1-j,bi,bj)=DYF(i,j,bi,bj)
199 DYF(i,sNy+j,bi,bj)=DYF(i,sNy+1-j,bi,bj)
200 ENDDO
201 ENDDO
202 ENDDO
203 ENDDO
204 ENDIF
205 cs
206
207 CALL MDSREADFIELD('RA.bin',readBinaryPrec,'RS',1,RA, 1,myThid)
208 _EXCH_XY_R4(RA,myThid )
209
210 C- Corner quantities
211 C *********** this are not degbugged ************
212 CALL MDSREADFIELD('LONG.bin',readBinaryPrec,'RS',1,XG, 1,myThid)
213 CALL MDSREADFIELD('LATG.bin',readBinaryPrec,'RS',1,YG, 1,myThid)
214 IF (useCubedSphereExchange) THEN
215 cs- this block needed by cubed sphere until we write more useful I/O routines
216 bi=3
217 bj=1
218 YG(1,sNy+1,bj,1)=YG(1,1,bi,1)
219 bj=bj+2
220 YG(1,sNy+1,bj,1)=YG(1,1,bi,1)
221 bj=bj+2
222 YG(1,sNy+1,bj,1)=YG(1,1,bi,1)
223 bi=6
224 bj=2
225 YG(sNx+1,1,bj,1)=YG(1,1,bi,1)
226 bj=bj+2
227 YG(sNx+1,1,bj,1)=YG(1,1,bi,1)
228 bj=bj+2
229 YG(sNx+1,1,bj,1)=YG(1,1,bi,1)
230 cs- end block
231 ENDIF
232 CALL EXCH_Z_XY_RS(XG,myThid)
233 CALL EXCH_Z_XY_RS(YG,myThid)
234
235 CALL MDSREADFIELD('DXV.bin',readBinaryPrec,'RS',1,DXV, 1,myThid)
236 CALL MDSREADFIELD('DYU.bin',readBinaryPrec,'RS',1,DYU, 1,myThid)
237 cs- this block needed by cubed sphere until we write more useful I/O routines
238 C !!! _EXCH_ZUV_XY_R4(DXV, DYU, unSigned, myThid)
239 cs! this is not correct <= need paired exchange for dxv,dyu
240 IF (.NOT.useCubedSphereExchange) THEN
241 CALL EXCH_Z_XY_RS(DXV,myThid)
242 CALL EXCH_Z_XY_RS(DYU,myThid)
243 ELSE
244 DO bj = myByLo(myThid), myByHi(myThid)
245 DO bi = myBxLo(myThid), myBxHi(myThid)
246 cs! fix overlaps:
247 DO j=1,sNy
248 DO i=1,Olx
249 DXV(1-i,j,bi,bj)=DXV(1+i,j,bi,bj)
250 DXV(sNx+i,j,bi,bj)=DXV(i,j,bi,bj)
251 DYU(1-i,j,bi,bj)=DYU(1+i,j,bi,bj)
252 DYU(sNx+i,j,bi,bj)=DYU(i,j,bi,bj)
253 ENDDO
254 ENDDO
255 DO j=1,Oly
256 DO i=1-Olx,sNx+Olx
257 DXV(i,1-j,bi,bj)=DXV(i,1+j,bi,bj)
258 DXV(i,sNy+j,bi,bj)=DXV(i,j,bi,bj)
259 DYU(i,1-j,bi,bj)=DYU(i,1+j,bi,bj)
260 DYU(i,sNy+j,bi,bj)=DYU(i,j,bi,bj)
261 ENDDO
262 ENDDO
263 ENDDO
264 ENDDO
265 cs- end block
266 ENDIF
267
268 CALL MDSREADFIELD('RAZ.bin',readBinaryPrec,'RS',1,RAZ, 1,myThid)
269 IF (useCubedSphereExchange) THEN
270 cs- this block needed by cubed sphere until we write more useful I/O routines
271 CALL EXCH_Z_XY_RS(RAZ , myThid )
272 DO bj = myByLo(myThid), myByHi(myThid)
273 DO bi = myBxLo(myThid), myBxHi(myThid)
274 RAZ(sNx+1,1,bi,bj)=RAZ(1,1,bi,bj)
275 RAZ(1,sNy+1,bi,bj)=RAZ(1,1,bi,bj)
276 ENDDO
277 ENDDO
278 cs- end block
279 ENDIF
280 CALL EXCH_Z_XY_RS(RAZ,myThid)
281
282 C- Staggered (u,v pairs) quantities
283 CALL MDSREADFIELD('DXC.bin',readBinaryPrec,'RS',1,DXC, 1,myThid)
284 CALL MDSREADFIELD('DYC.bin',readBinaryPrec,'RS',1,DYC, 1,myThid)
285 CALL EXCH_UV_XY_RS(DXC,DYC,.FALSE.,myThid)
286
287 CALL MDSREADFIELD('RAW.bin',readBinaryPrec,'RS',1,RAW, 1,myThid)
288 CALL MDSREADFIELD('RAS.bin',readBinaryPrec,'RS',1,RAS, 1,myThid)
289 IF (useCubedSphereExchange) THEN
290 cs- this block needed by cubed sphere until we write more useful I/O routines
291 DO bj = myByLo(myThid), myByHi(myThid)
292 DO bi = myBxLo(myThid), myBxHi(myThid)
293 DO J = 1,sNy
294 c RAW(sNx+1,J,bi,bj)=RAW(1,J,bi,bj)
295 c RAS(J,sNy+1,bi,bj)=RAS(J,1,bi,bj)
296 ENDDO
297 ENDDO
298 ENDDO
299 cs- end block
300 ENDIF
301 CALL EXCH_UV_XY_RS(RAW,RAS,.FALSE.,myThid)
302
303 CALL MDSREADFIELD('DXG.bin',readBinaryPrec,'RS',1,DXG, 1,myThid)
304 CALL MDSREADFIELD('DYG.bin',readBinaryPrec,'RS',1,DYG, 1,myThid)
305 IF (useCubedSphereExchange) THEN
306 cs- this block needed by cubed sphere until we write more useful I/O routines
307 DO bj = myByLo(myThid), myByHi(myThid)
308 DO bi = myBxLo(myThid), myBxHi(myThid)
309 DO J = 1,sNy
310 c DYG(sNx+1,J,bi,bj)=DYG(1,J,bi,bj)
311 c DXG(J,sNy+1,bi,bj)=DXG(J,1,bi,bj)
312 ENDDO
313 ENDDO
314 ENDDO
315 cs- end block
316 ENDIF
317 CALL EXCH_UV_XY_RS(DYG,DXG,.FALSE.,myThid)
318 CALL EXCH_UV_AGRID_XY_RS(angleSinC,angleCosC,.TRUE.,myThid)
319
320 c write(10) XC
321 c write(10) YC
322 c write(10) DXF
323 c write(10) DYF
324 c write(10) RA
325 c write(10) XG
326 c write(10) YG
327 c write(10) DXV
328 c write(10) DYU
329 c write(10) RAZ
330 c write(10) DXC
331 c write(10) DYC
332 c write(10) RAW
333 c write(10) RAS
334 c write(10) DXG
335 c write(10) DYG
336
337 #else /* ifndef OLD_GRID_IO */
338
339 C-- Only do I/O if I am the master thread
340 _BEGIN_MASTER(myThid)
341
342 write(msgbuf,'(3a)') 'Open Bondary File'
343 call print_message( msgbuf, standardmessageunit,
344 & SQUEEZE_RIGHT , mythid)
345
346 DO bj = 1,nSy
347 DO bi = 1,nSx
348 iG=bi+(myXGlobalLo-1)/sNx
349 WRITE(tmpBuf,'(A,I4)') 'tile:',iG
350 #ifdef ALLOW_EXCH2
351 myTile = W2_myTileList(bi)
352 WRITE(tmpBuf,'(A,I4)') 'tile:',myTile
353 iG = exch2_myface(myTile)
354 #endif
355 iLen = ILNBLNK(horizGridFile)
356 IF ( iLen .EQ. 0 ) THEN
357 WRITE(fName,'("tile",I3.3,".mitgrid")') iG
358 ELSE
359 WRITE(fName,'(2A,I3.3,A)') horizGridFile(1:iLen),
360 & '.face',iG,'.bin'
361 ENDIF
362 iLen = ILNBLNK(fName)
363 iL = ILNBLNK(tmpBuf)
364 WRITE(msgBuf,'(3A)') tmpBuf(1:iL),
365 & ' ; Read from file ',fName(1:iLen)
366 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
367 & SQUEEZE_RIGHT , myThid)
368 WRITE(msgBuf,'(A)') ' =>'
369 CALL READSYMTILE_RS(fName,1,XC,bi,bj,buf,myThid)
370 iL = ILNBLNK(msgBuf)
371 WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'XC'
372 CALL READSYMTILE_RS(fName,2,YC,bi,bj,buf,myThid)
373 iL = ILNBLNK(tmpBuf)
374 WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'YC'
375 CALL READSYMTILE_RS(fName,3,DXF,bi,bj,buf,myThid)
376 iL = ILNBLNK(msgBuf)
377 WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'DXF'
378 CALL READSYMTILE_RS(fName,4,DYF,bi,bj,buf,myThid)
379 iL = ILNBLNK(tmpBuf)
380 WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'DYF'
381 CALL READSYMTILE_RS(fName,5,RA,bi,bj,buf,myThid)
382 iL = ILNBLNK(msgBuf)
383 WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'RA'
384 CALL READSYMTILE_RS(fName,6,XG,bi,bj,buf,myThid)
385 iL = ILNBLNK(tmpBuf)
386 WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'XG'
387 CALL READSYMTILE_RS(fName,7,YG,bi,bj,buf,myThid)
388 iL = ILNBLNK(msgBuf)
389 WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'YG'
390 CALL READSYMTILE_RS(fName,8,DXV,bi,bj,buf,myThid)
391 iL = ILNBLNK(tmpBuf)
392 WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'DXV'
393 CALL READSYMTILE_RS(fName,9,DYU,bi,bj,buf,myThid)
394 iL = ILNBLNK(msgBuf)
395 WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'DYU'
396 CALL READSYMTILE_RS(fName,10,RAZ,bi,bj,buf,myThid)
397 iL = ILNBLNK(tmpBuf)
398 WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'RAZ'
399 CALL READSYMTILE_RS(fName,11,DXC,bi,bj,buf,myThid)
400 iL = ILNBLNK(msgBuf)
401 WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'DXC'
402 CALL READSYMTILE_RS(fName,12,DYC,bi,bj,buf,myThid)
403 iL = ILNBLNK(tmpBuf)
404 WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'DYC'
405 CALL READSYMTILE_RS(fName,13,RAW,bi,bj,buf,myThid)
406 iL = ILNBLNK(msgBuf)
407 WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'RAW'
408 CALL READSYMTILE_RS(fName,14,RAS,bi,bj,buf,myThid)
409 iL = ILNBLNK(tmpBuf)
410 WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'RAS'
411 CALL READSYMTILE_RS(fName,15,DXG,bi,bj,buf,myThid)
412 iL = ILNBLNK(msgBuf)
413 WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'DXG'
414 CALL READSYMTILE_RS(fName,16,DYG,bi,bj,buf,myThid)
415 iL = ILNBLNK(tmpBuf)
416 WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'DYG'
417 c iLen = ILNBLNK(horizGridFile)
418 c IF ( iLen.GT.0 ) THEN
419 c CALL READSYMTILE_RS(fName,17,angleCosC,bi,bj,buf,myThid)
420 c iL = ILNBLNK(msgBuf)
421 c WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'AngleCS'
422 c CALL READSYMTILE_RS(fName,18,angleSinC,bi,bj,buf,myThid)
423 c iL = ILNBLNK(tmpBuf)
424 c WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'AngleSN'
425 c ENDIF
426
427 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
428 & SQUEEZE_RIGHT , myThid)
429
430 ENDDO
431 ENDDO
432
433 _END_MASTER(myThid)
434
435 CALL EXCH_XY_RS(XC,myThid)
436 CALL EXCH_XY_RS(YC,myThid)
437 C !!! _EXCH_OUV_XY_R4(DXF, DYF, unSigned, myThid )
438 #ifdef HRCUBE
439 CALL EXCH_XY_RS(DXF,myThid)
440 CALL EXCH_XY_RS(DYF,myThid)
441 #endif
442 CALL EXCH_XY_RS(RA,myThid )
443 CALL EXCH_Z_XY_RS(XG,myThid)
444 CALL EXCH_Z_XY_RS(YG,myThid)
445 C !!! _EXCH_ZUV_XY_R4(DXV, DYU, unSigned, myThid)
446 c CALL EXCH_Z_XY_RS(DXV,myThid)
447 c CALL EXCH_Z_XY_RS(DYU,myThid)
448 CALL EXCH_Z_XY_RS(RAZ,myThid)
449 CALL EXCH_UV_XY_RS(DXC,DYC,.FALSE.,myThid)
450 CALL EXCH_UV_XY_RS(RAW,RAS,.FALSE.,myThid)
451 CALL EXCH_UV_XY_RS(DYG,DXG,.FALSE.,myThid)
452 CALL EXCH_UV_AGRID_XY_RS(angleSinC,angleCosC,.TRUE.,myThid)
453
454 #endif /* OLD_GRID_IO */
455
456 #ifdef ALLOW_MNC
457 ENDIF
458 #endif /* ALLOW_MNC */
459
460 c CALL WRITE_FULLARRAY_RL('DXV',DXV,1,0,0,0,myThid)
461 c CALL WRITE_FULLARRAY_RL('DYU',DYU,1,0,0,0,myThid)
462 c CALL WRITE_FULLARRAY_RL('RAZ',RAZ,1,0,0,0,myThid)
463 c CALL WRITE_FULLARRAY_RL('XG',XG,1,0,0,0,myThid)
464 c CALL WRITE_FULLARRAY_RL('YG',YG,1,0,0,0,myThid)
465
466 C-- Require that 0 <= longitude < 360 if using exf package
467 #ifdef ALLOW_EXF
468 DO bj = 1,nSy
469 DO bi = 1,nSx
470 DO J=1-Oly,sNy+Oly
471 DO I=1-Olx,sNx+Olx
472 IF (XC(i,j,bi,bj).lt.0.) XC(i,j,bi,bj)=XC(i,j,bi,bj)+360.
473 IF (XG(i,j,bi,bj).lt.0.) XG(i,j,bi,bj)=XG(i,j,bi,bj)+360.
474 ENDDO
475 ENDDO
476 ENDDO
477 ENDDO
478 #endif /* ALLOW_EXF */
479
480 C-- Now let's look at all these beasts
481 IF ( debugLevel .GE. debLevB ) THEN
482 myIter = 1
483 CALL PLOT_FIELD_XYRL( XC , 'Current XC ' ,
484 & myIter, myThid )
485 CALL PLOT_FIELD_XYRL( YC , 'Current YC ' ,
486 & myIter, myThid )
487 CALL PLOT_FIELD_XYRL( DXF , 'Current DXF ' ,
488 & myIter, myThid )
489 CALL PLOT_FIELD_XYRL( XC , 'Current XC ' ,
490 & myIter, myThid )
491 CALL PLOT_FIELD_XYRL( DYF , 'Current DYF ' ,
492 & myIter, myThid )
493 CALL PLOT_FIELD_XYRL( RA , 'Current RA ' ,
494 & myIter, myThid )
495 CALL PLOT_FIELD_XYRL( XG , 'Current XG ' ,
496 & myIter, myThid )
497 CALL PLOT_FIELD_XYRL( YG , 'Current YG ' ,
498 & myIter, myThid )
499 CALL PLOT_FIELD_XYRL( DXV , 'Current DXV ' ,
500 & myIter, myThid )
501 CALL PLOT_FIELD_XYRL( DYU , 'Current DYU ' ,
502 & myIter, myThid )
503 CALL PLOT_FIELD_XYRL( RAZ , 'Current RAZ ' ,
504 & myIter, myThid )
505 CALL PLOT_FIELD_XYRL( DXC , 'Current DXC ' ,
506 & myIter, myThid )
507 CALL PLOT_FIELD_XYRL( DYC , 'Current DYC ' ,
508 & myIter, myThid )
509 CALL PLOT_FIELD_XYRL( RAW , 'Current RAW ' ,
510 & myIter, myThid )
511 CALL PLOT_FIELD_XYRL( RAS , 'Current RAS ' ,
512 & myIter, myThid )
513 CALL PLOT_FIELD_XYRL( DXG , 'Current DXG ' ,
514 & myIter, myThid )
515 CALL PLOT_FIELD_XYRL( DYG , 'Current DYG ' ,
516 & myIter, myThid )
517 CALL PLOT_FIELD_XYRL(angleCosC, 'Current AngleCS ' ,
518 & myIter, myThid )
519 CALL PLOT_FIELD_XYRL(angleSinC, 'Current AngleSN ' ,
520 & myIter, myThid )
521 ENDIF
522
523 RETURN
524 END
525
526 C --------------------------------------------------------------------------
527
528 SUBROUTINE READSYMTILE_RS(fName,irec,array,bi,bj,buf,myThid)
529 C /==========================================================\
530 C | SUBROUTINE READSYMTILE_RS |
531 C |==========================================================|
532 C \==========================================================/
533 IMPLICIT NONE
534
535 C === Global variables ===
536 #include "SIZE.h"
537 #include "EEPARAMS.h"
538 #ifdef ALLOW_EXCH2
539 #include "W2_EXCH2_TOPOLOGY.h"
540 #include "W2_EXCH2_PARAMS.h"
541 #endif /* ALLOW_EXCH2 */
542
543 C == Routine arguments ==
544 CHARACTER*(*) fName
545 INTEGER irec
546 _RS array(1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
547 INTEGER bi,bj,myThid
548 #ifdef ALLOW_EXCH2
549 _RL buf(1:sNx*nSx*nPx+1)
550 #else
551 _RL buf(1:sNx+1,1:sNy+1)
552 #endif /* ALLOW_EXCH2 */
553
554 C == Local variables ==
555 INTEGER I,J,dUnit, iLen
556 INTEGER length_of_rec
557 INTEGER MDS_RECLEN
558 #ifdef ALLOW_EXCH2
559 INTEGER TN, dNx, dNy, TBX, TBY, TNX, TNY, II, iBase
560 #endif
561 INTEGER ILNBLNK
562 EXTERNAL ILNBLNK
563
564 iLen = ILNBLNK(fName)
565 #ifdef ALLOW_EXCH2
566 C Figure out offset of tile within face
567 TN = W2_myTileList(bi)
568 dNx = exch2_mydnx(TN)
569 dNy = exch2_mydny(TN)
570 TBX = exch2_tbasex(TN)
571 TBY = exch2_tbasey(TN)
572 TNX = exch2_tnx(TN)
573 TNY = exch2_tny(TN)
574
575 CALL MDSFINDUNIT( dUnit, myThid )
576 length_of_rec=MDS_RECLEN( 64, (dNx+1), myThid )
577 OPEN( dUnit, file=fName(1:iLen), status='old',
578 & access='direct', recl=length_of_rec )
579 J=0
580 iBase=(irec-1)*(dny+1)
581 DO I=1+TBY,sNy+1+TBY
582 READ(dUnit,rec=I+iBase)(buf(ii),ii=1,dNx+1)
583 #ifdef _BYTESWAPIO
584 #ifdef REAL4_IS_SLOW
585 CALL MDS_BYTESWAPR8((dNx+1), buf)
586 #else
587 CALL MDS_BYTESWAPR4((dNx+1), buf)
588 #endif
589 #endif
590 J=J+1
591 DO II=1,sNx+1
592 array(II,J,bi,bj)=buf(II+TBX)
593 ENDDO
594 ENDDO
595 CLOSE( dUnit )
596
597 #else /* ALLOW_EXCH2 */
598
599 CALL MDSFINDUNIT( dUnit, myThid )
600 length_of_rec=MDS_RECLEN( 64, (sNx+1)*(sNy+1), myThid )
601 OPEN( dUnit, file=fName(1:iLen), status='old',
602 & access='direct', recl=length_of_rec )
603 READ(dUnit,rec=irec) buf
604 CLOSE( dUnit )
605
606 #ifdef _BYTESWAPIO
607 #ifdef REAL4_IS_SLOW
608 CALL MDS_BYTESWAPR8((sNx+1)*(sNy+1), buf)
609 #else
610 CALL MDS_BYTESWAPR4((sNx+1)*(sNy+1), buf)
611 #endif
612 #endif
613
614 DO J=1,sNy+1
615 DO I=1,sNx+1
616 array(I,J,bi,bj)=buf(I,J)
617 ENDDO
618 ENDDO
619 c write(0,*) irec,buf(1,1),array(1,1,1,1)
620
621 #endif /* ALLOW_EXCH2 */
622
623 RETURN
624 END

  ViewVC Help
Powered by ViewVC 1.1.22