/[MITgcm]/MITgcm_contrib/high_res_cube/code-mods/ini_curvilinear_grid.F
ViewVC logotype

Annotation of /MITgcm_contrib/high_res_cube/code-mods/ini_curvilinear_grid.F

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


Revision 1.4 - (hide annotations) (download)
Fri Mar 26 15:19:35 2004 UTC (21 years, 8 months ago) by edhill
Branch: MAIN
CVS Tags: HEAD
Changes since 1.3: +1 -1 lines
FILE REMOVED
 o removing files (as requested by CNH) to fix the recent CVS repository
   corruption (the dead-files-returning problem)

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

  ViewVC Help
Powered by ViewVC 1.1.22