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

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

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

revision 1.2 by adcroft, Tue May 29 14:01:37 2001 UTC revision 1.32 by jmc, Sat Sep 2 21:27:58 2006 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    C- note: default is to use "new" grid files (OLD_GRID_IO undef)
7    C   but can still use (on 1 cpu, with MDSIO) OLD_GRID_IO and EXCH2 independently
8    #undef OLD_GRID_IO
9    
10    CBOP
11    C     !ROUTINE: INI_CURVILINEAR_GRID
12    C     !INTERFACE:
13        SUBROUTINE INI_CURVILINEAR_GRID( myThid )        SUBROUTINE INI_CURVILINEAR_GRID( myThid )
14  C     /==========================================================\  C     !DESCRIPTION: \bv
15  C     | SUBROUTINE INI_CURVILINEAR_GRID                          |  C     *==========================================================*
16  C     | o Initialise curvilinear coordinate system               |  C     | SUBROUTINE INI_CURVILINEAR_GRID
17  C     |==========================================================|  C     | o Initialise curvilinear coordinate system
18  C     \==========================================================/  C     *==========================================================*
19        IMPLICIT NONE  C     | Curvilinear grid settings are read from a file rather
20    C     | than coded in-line as for cartesian and spherical polar.
21    C     | This is more general but you have to create the grid
22    C     | yourself.
23    C     *==========================================================*
24    C     \ev
25    
26    C     !USES:
27          IMPLICIT NONE
28  C     === Global variables ===  C     === Global variables ===
29  #include "SIZE.h"  #include "SIZE.h"
30  #include "EEPARAMS.h"  #include "EEPARAMS.h"
31  #include "PARAMS.h"  #include "PARAMS.h"
32  #include "GRID.h"  #include "GRID.h"
33    #ifdef ALLOW_EXCH2
34    #include "W2_EXCH2_TOPOLOGY.h"
35    #include "W2_EXCH2_PARAMS.h"
36    #endif
37    #ifdef ALLOW_MNC
38    #include "MNC_PARAMS.h"
39    #endif
40    
41    C     !INPUT/OUTPUT PARAMETERS:
42  C     == Routine arguments ==  C     == Routine arguments ==
43  C     myThid -  Number of this instance of INI_CARTESIAN_GRID  C     myThid -  Number of this instance of INI_CURVILINEAR_GRID
44        INTEGER myThid        INTEGER myThid
45    
46    C     !LOCAL VARIABLES:
47  C     == Local variables ==  C     == Local variables ==
48        INTEGER bi,bj        INTEGER bi,bj, myIter
49        INTEGER I,J        INTEGER I,J
50          CHARACTER*(MAX_LEN_MBUF) msgBuf
51          LOGICAL anglesAreSet
52    #ifdef ALLOW_MNC
53          CHARACTER*(80) mncFn
54    #endif
55    #ifndef OLD_GRID_IO
56    # ifdef ALLOW_EXCH2
57          _RL buf(sNx*nSx*nPx+1)
58          INTEGER myTile
59    # else
60          _RL buf(sNx+1,sNy+1)
61    # endif
62          INTEGER iG, iL, iLen
63          CHARACTER*(MAX_LEN_FNAM) fName
64          CHARACTER*(MAX_LEN_MBUF) tmpBuf
65          INTEGER  ILNBLNK
66          EXTERNAL ILNBLNK
67    #endif
68    CEOP
69    
70  C--   Set everything to zero everywhere  C--   Set everything to zero everywhere
71        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
# Line 31  C--   Set everything to zero everywhere Line 73  C--   Set everything to zero everywhere
73    
74          DO J=1-Oly,sNy+Oly          DO J=1-Oly,sNy+Oly
75           DO I=1-Olx,sNx+Olx           DO I=1-Olx,sNx+Olx
76            XC(i,j,bi,bj)=0.            xC(i,j,bi,bj)=0.
77            YC(i,j,bi,bj)=0.            yC(i,j,bi,bj)=0.
78            XG(i,j,bi,bj)=0.            xG(i,j,bi,bj)=0.
79            YG(i,j,bi,bj)=0.            yG(i,j,bi,bj)=0.
80            DXC(i,j,bi,bj)=0.            dxC(i,j,bi,bj)=0.
81            DYC(i,j,bi,bj)=0.            dyC(i,j,bi,bj)=0.
82            DXG(i,j,bi,bj)=0.            dxG(i,j,bi,bj)=0.
83            DYG(i,j,bi,bj)=0.            dyG(i,j,bi,bj)=0.
84            DXF(i,j,bi,bj)=0.            dxF(i,j,bi,bj)=0.
85            DYF(i,j,bi,bj)=0.            dyF(i,j,bi,bj)=0.
86            DXV(i,j,bi,bj)=0.            dxV(i,j,bi,bj)=0.
87            DYU(i,j,bi,bj)=0.            dyU(i,j,bi,bj)=0.
88            RA(i,j,bi,bj)=0.            rA(i,j,bi,bj)=0.
89            RAZ(i,j,bi,bj)=0.            rAz(i,j,bi,bj)=0.
90            RAW(i,j,bi,bj)=0.            rAw(i,j,bi,bj)=0.
91            RAS(i,j,bi,bj)=0.            rAs(i,j,bi,bj)=0.
92            tanPhiAtU(i,j,bi,bj)=0.            tanPhiAtU(i,j,bi,bj)=0.
93            tanPhiAtV(i,j,bi,bj)=0.            tanPhiAtV(i,j,bi,bj)=0.
94              angleCosC(i,j,bi,bj)=1.
95              angleSinC(i,j,bi,bj)=0.
96            cosFacU(J,bi,bj)=1.            cosFacU(J,bi,bj)=1.
97            cosFacV(J,bi,bj)=1.            cosFacV(J,bi,bj)=1.
98            sqcosFacU(J,bi,bj)=1.            sqCosFacU(J,bi,bj)=1.
99            sqcosFacV(J,bi,bj)=1.            sqCosFacV(J,bi,bj)=1.
100           ENDDO           ENDDO
101          ENDDO          ENDDO
102    
103         ENDDO ! bi         ENDDO
104        ENDDO ! bj        ENDDO
105    
106    C--   Everyone must wait for the initialisation to be done
107          _BARRIER
108    
109  C     Here we make no assumptions about grid symmetry and simply  C     Here we make no assumptions about grid symmetry and simply
110  C     read the raw grid data from files  C     read the raw grid data from files
111    
112    #ifdef OLD_GRID_IO
113    
114  C-    Cell centered quantities  C-    Cell centered quantities
115        CALL MDSREADFIELD('LONC.bin',readBinaryPrec,'RL',1,XC,  1,myThid)        CALL MDSREADFIELD('LONC.bin',readBinaryPrec,'RS',1,xC,  1,myThid)
116        CALL MDSREADFIELD('LATC.bin',readBinaryPrec,'RL',1,YC,  1,myThid)        CALL MDSREADFIELD('LATC.bin',readBinaryPrec,'RS',1,yC,  1,myThid)
117        _EXCH_XY_R4(XC,myThid)        _EXCH_XY_R4(xC,myThid)
118        _EXCH_XY_R4(YC,myThid)        _EXCH_XY_R4(yC,myThid)
119    
120        CALL MDSREADFIELD('DXF.bin',readBinaryPrec,'RL',1,DXF,  1,myThid)        CALL MDSREADFIELD('DXF.bin',readBinaryPrec,'RS',1,dxF,  1,myThid)
121        CALL MDSREADFIELD('DYF.bin',readBinaryPrec,'RL',1,DYF,  1,myThid)        CALL MDSREADFIELD('DYF.bin',readBinaryPrec,'RS',1,dyF,  1,myThid)
122  C !!! _EXCH_OUV_XY_R4(DXF, DYF, unSigned, myThid )        CALL EXCH_UV_AGRID_3D_RS( dxF, dyF, .FALSE., 1, myThid )
 cs!   this is not correct! <= need paired exchange for DXF,DYF  
       _EXCH_XY_R4(DXF,myThid)  
       _EXCH_XY_R4(DYF,myThid)  
123    
124        CALL MDSREADFIELD('RA.bin',readBinaryPrec,'RL',1,RA,  1,myThid)        CALL MDSREADFIELD('RA.bin',readBinaryPrec,'RS',1,rA,  1,myThid)
125        _EXCH_XY_R4(RA,myThid )        _EXCH_XY_R4(rA,myThid )
126    
127  C-    Corner quantities  C-    Corner quantities
128  C       *********** this are not degbugged ************        CALL MDSREADFIELD('LONG.bin',readBinaryPrec,'RS',1,xG,  1,myThid)
129        CALL MDSREADFIELD('LONG.bin',readBinaryPrec,'RL',1,XG,  1,myThid)        CALL MDSREADFIELD('LATG.bin',readBinaryPrec,'RS',1,yG,  1,myThid)
130        CALL MDSREADFIELD('LATG.bin',readBinaryPrec,'RL',1,YG,  1,myThid)        IF (useCubedSphereExchange) THEN
131  cs-   this block needed by cubed sphere until we write more useful I/O routines  cs-   this block needed by cubed sphere until we write more useful I/O routines
132        bi=3        bi=3
133        bj=1        bj=1
134        YG(1,sNy+1,bj,1)=YG(1,1,bi,1)        yG(1,sNy+1,bj,1)=yG(1,1,bi,1)
135        bj=bj+2        bj=bj+2
136        YG(1,sNy+1,bj,1)=YG(1,1,bi,1)        yG(1,sNy+1,bj,1)=yG(1,1,bi,1)
137        bj=bj+2        bj=bj+2
138        YG(1,sNy+1,bj,1)=YG(1,1,bi,1)        yG(1,sNy+1,bj,1)=yG(1,1,bi,1)
139        bi=6        bi=6
140        bj=2        bj=2
141        YG(sNx+1,1,bj,1)=YG(1,1,bi,1)        yG(sNx+1,1,bj,1)=yG(1,1,bi,1)
142        bj=bj+2        bj=bj+2
143        YG(sNx+1,1,bj,1)=YG(1,1,bi,1)        yG(sNx+1,1,bj,1)=yG(1,1,bi,1)
144        bj=bj+2        bj=bj+2
145        YG(sNx+1,1,bj,1)=YG(1,1,bi,1)        yG(sNx+1,1,bj,1)=yG(1,1,bi,1)
146  cs-   end block  cs-   end block
147        CALL EXCH_Z_XY_RS(XG,myThid)        ENDIF
148        CALL EXCH_Z_XY_RS(YG,myThid)        CALL EXCH_Z_3D_RS( xG, 1, myThid )
149          CALL EXCH_Z_3D_RS( yG, 1, myThid )
150    
151        CALL MDSREADFIELD('DXV.bin',readBinaryPrec,'RL',1,DXV,  1,myThid)        CALL MDSREADFIELD('DXV.bin',readBinaryPrec,'RS',1,dxV,  1,myThid)
152        CALL MDSREADFIELD('DYU.bin',readBinaryPrec,'RL',1,DYU,  1,myThid)        CALL MDSREADFIELD('DYU.bin',readBinaryPrec,'RS',1,dyU,  1,myThid)
153  cs-   this block needed by cubed sphere until we write more useful I/O routines  cs-   this block needed by cubed sphere until we write more useful I/O routines
154  C !!! _EXCH_ZUV_XY_R4(DXV, DYU, unSigned, myThid)  C !!! _EXCH_ZUV_XY_R4(dxV, dyU, unSigned, myThid)
155  cs!   this is not correct <= need paired exchange for dxv,dyu  cs!   this is not correct <= need paired exchange for dxv,dyu
156        CALL EXCH_Z_XY_RS(DXV,myThid)        IF (.NOT.useCubedSphereExchange) THEN
157        CALL EXCH_Z_XY_RS(DYU,myThid)        CALL EXCH_Z_3D_RS( dxV, 1, myThid )
158          CALL EXCH_Z_3D_RS( dyU, 1, myThid )
159          ELSE
160        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
161         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
162          DXV(sNx+1,1,bi,bj)=DXV(1,1,bi,bj)  cs! fix overlaps:
163          DXV(1,sNy+1,bi,bj)=DXV(1,1,bi,bj)          DO j=1,sNy
164          DYU(sNx+1,1,bi,bj)=DYU(1,1,bi,bj)           DO i=1,Olx
165          DYU(1,sNy+1,bi,bj)=DYU(1,1,bi,bj)            dxV(1-i,j,bi,bj)=dxV(1+i,j,bi,bj)
166              dxV(sNx+i,j,bi,bj)=dxV(i,j,bi,bj)
167              dyU(1-i,j,bi,bj)=dyU(1+i,j,bi,bj)
168              dyU(sNx+i,j,bi,bj)=dyU(i,j,bi,bj)
169             ENDDO
170            ENDDO
171            DO j=1,Oly
172             DO i=1-Olx,sNx+Olx
173              dxV(i,1-j,bi,bj)=dxV(i,1+j,bi,bj)
174              dxV(i,sNy+j,bi,bj)=dxV(i,j,bi,bj)
175              dyU(i,1-j,bi,bj)=dyU(i,1+j,bi,bj)
176              dyU(i,sNy+j,bi,bj)=dyU(i,j,bi,bj)
177             ENDDO
178            ENDDO
179         ENDDO         ENDDO
180        ENDDO        ENDDO
181  cs-   end block  cs-   end block
182  C !!! _EXCH_ZUV_XY_R4(DXV, DYU, unSigned, myThid)        ENDIF
 cs!   this is not correct <= need paired exchange for dxv,dyu  
       CALL EXCH_Z_XY_RS(DXV,myThid)  
       CALL EXCH_Z_XY_RS(DYU,myThid)  
183    
184        CALL MDSREADFIELD('RAZ.bin',readBinaryPrec,'RL',1,RAZ,  1,myThid)        CALL MDSREADFIELD('RAZ.bin',readBinaryPrec,'RS',1,rAz,  1,myThid)
185          IF (useCubedSphereExchange) THEN
186  cs-   this block needed by cubed sphere until we write more useful I/O routines  cs-   this block needed by cubed sphere until we write more useful I/O routines
187        CALL EXCH_Z_XY_RS(RAZ , myThid )        CALL EXCH_Z_3D_RS( rAz, 1, myThid )
188        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
189         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
190          RAZ(sNx+1,1,bi,bj)=RAZ(1,1,bi,bj)          rAz(sNx+1,1,bi,bj)=rAz(1,1,bi,bj)
191          RAZ(1,sNy+1,bi,bj)=RAZ(1,1,bi,bj)          rAz(1,sNy+1,bi,bj)=rAz(1,1,bi,bj)
192         ENDDO         ENDDO
193        ENDDO        ENDDO
194  cs-   end block  cs-   end block
195        CALL EXCH_Z_XY_RS(RAZ,myThid)        ENDIF
196          CALL EXCH_Z_3D_RS( rAz, 1, myThid )
197    
198  C-    Staggered (u,v pairs) quantities  C-    Staggered (u,v pairs) quantities
199        CALL MDSREADFIELD('DXC.bin',readBinaryPrec,'RL',1,DXC,  1,myThid)        CALL MDSREADFIELD('DXC.bin',readBinaryPrec,'RS',1,dxC,  1,myThid)
200        CALL MDSREADFIELD('DYC.bin',readBinaryPrec,'RL',1,DYC,  1,myThid)        CALL MDSREADFIELD('DYC.bin',readBinaryPrec,'RS',1,dyC,  1,myThid)
201        CALL EXCH_UV_XY_RS(DXC,DYC,.FALSE.,myThid)        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          CALL EXCH_UV_XY_RS(rAw,rAs,.FALSE.,myThid)
206    
207          CALL MDSREADFIELD('DXG.bin',readBinaryPrec,'RS',1,dxG,  1,myThid)
208          CALL MDSREADFIELD('DYG.bin',readBinaryPrec,'RS',1,dyG,  1,myThid)
209          CALL EXCH_UV_XY_RS(dyG,dxG,.FALSE.,myThid)
210          CALL EXCH_UV_AGRID_3D_RS(angleSinC,angleCosC,.TRUE., 1, myThid)
211          anglesAreSet = .FALSE.
212    
213    c     write(10) xC
214    c     write(10) yC
215    c     write(10) dxF
216    c     write(10) dyF
217    c     write(10) rA
218    c     write(10) xG
219    c     write(10) yG
220    c     write(10) dxV
221    c     write(10) dyU
222    c     write(10) rAz
223    c     write(10) dxC
224    c     write(10) dyC
225    c     write(10) rAw
226    c     write(10) rAs
227    c     write(10) dxG
228    c     write(10) dyG
229    
230    #else /* ifndef OLD_GRID_IO */
231    
232    C--   Only do I/O if I am the master thread
233          _BEGIN_MASTER(myThid)
234    
235    #ifdef ALLOW_MNC
236          IF (useMNC .AND. readgrid_mnc) THEN
237    C--   read NetCDF files:
238    
239        CALL MDSREADFIELD('RAW.bin',readBinaryPrec,'RL',1,RAW,  1,myThid)          DO i = 1,80
240        CALL MDSREADFIELD('RAS.bin',readBinaryPrec,'RL',1,RAS,  1,myThid)            mncFn(i:i) = ' '
241  cs-   this block needed by cubed sphere until we write more useful I/O routines          ENDDO
242        DO bj = myByLo(myThid), myByHi(myThid)          write(mncFn,'(a)') 'mitgrid'
243         DO bi = myBxLo(myThid), myBxHi(myThid)          DO i = 1,MAX_LEN_MBUF
244          DO J = 1,sNy            msgBuf(i:i) = ' '
 c        RAW(sNx+1,J,bi,bj)=RAW(1,J,bi,bj)  
 c        RAS(J,sNy+1,bi,bj)=RAS(J,1,bi,bj)  
245          ENDDO          ENDDO
246            WRITE(msgBuf,'(2A)') msgBuf,' ; Reading grid info using MNC'
247            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
248         &       SQUEEZE_RIGHT , myThid)
249            CALL MNC_FILE_CLOSE_ALL_MATCHING(mncFn, myThid)
250            CALL MNC_CW_SET_UDIM(mncFn, 1, myThid)
251            CALL MNC_CW_SET_CITER(mncFn, 2, -1, -1, -1, myThid)
252            CALL MNC_CW_SET_UDIM(mncFn, 1, myThid)
253            CALL MNC_CW_RS_R('D',mncFn,0,0,'XC', xC,  myThid)
254            CALL MNC_CW_RS_R('D',mncFn,0,0,'XG', xG,  myThid)
255            CALL MNC_CW_RS_R('D',mncFn,0,0,'YC', yC,  myThid)
256            CALL MNC_CW_RS_R('D',mncFn,0,0,'YG', yG,  myThid)
257            CALL MNC_CW_RS_R('D',mncFn,0,0,'dxC',dxC, myThid)
258            CALL MNC_CW_RS_R('D',mncFn,0,0,'dyC',dyC, myThid)
259            CALL MNC_CW_RS_R('D',mncFn,0,0,'dxF',dxF, myThid)
260            CALL MNC_CW_RS_R('D',mncFn,0,0,'dyF',dyF, myThid)
261            CALL MNC_CW_RS_R('D',mncFn,0,0,'dxG',dxG, myThid)
262            CALL MNC_CW_RS_R('D',mncFn,0,0,'dyG',dyG, myThid)
263            CALL MNC_CW_RS_R('D',mncFn,0,0,'dxV',dxV, myThid)
264            CALL MNC_CW_RS_R('D',mncFn,0,0,'dyU',dyU, myThid)
265            CALL MNC_CW_RS_R('D',mncFn,0,0,'rA', rA,  myThid)
266            CALL MNC_CW_RS_R('D',mncFn,0,0,'rAz',rAz, myThid)
267            CALL MNC_CW_RS_R('D',mncFn,0,0,'rAw',rAw, myThid)
268            CALL MNC_CW_RS_R('D',mncFn,0,0,'rAs',rAs, myThid)
269            CALL MNC_CW_RS_R('D',mncFn,0,0,'AngleCS',angleCosC,myThid)
270            CALL MNC_CW_RS_R('D',mncFn,0,0,'AngleSN',angleSinC,myThid)
271            anglesAreSet = .TRUE.
272    
273          ELSE
274    C--   read Binary files:
275    #endif /* ALLOW_MNC */
276    
277          DO bj = 1,nSy
278           DO bi = 1,nSx
279            iG=bi+(myXGlobalLo-1)/sNx
280            WRITE(tmpBuf,'(A,I4)') 'tile:',iG
281    #ifdef ALLOW_EXCH2
282            myTile = W2_myTileList(bi)
283            WRITE(tmpBuf,'(A,I4)') 'tile:',myTile
284            iG = exch2_myface(myTile)
285    #endif
286            iLen = ILNBLNK(horizGridFile)
287            IF ( iLen .EQ. 0 ) THEN
288              WRITE(fName,'("tile",I3.3,".mitgrid")') iG
289            ELSE
290              WRITE(fName,'(2A,I3.3,A)') horizGridFile(1:iLen),
291         &                              '.face',iG,'.bin'
292            ENDIF
293            iLen = ILNBLNK(fName)
294            iL = ILNBLNK(tmpBuf)
295            WRITE(msgBuf,'(3A)') tmpBuf(1:iL),
296         &                   ' ; Read from file ',fName(1:iLen)
297            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
298         &                      SQUEEZE_RIGHT , myThid)
299            WRITE(msgBuf,'(A)') '  =>'
300    
301            CALL READSYMTILE_RS(fName,1,xC,buf,bi,bj,myThid)
302            iL = ILNBLNK(msgBuf)
303            WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'xC'
304            CALL READSYMTILE_RS(fName,2,yC,buf,bi,bj,myThid)
305            iL = ILNBLNK(tmpBuf)
306            WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'yC'
307            CALL READSYMTILE_RS(fName,3,dxF,buf,bi,bj,myThid)
308            iL = ILNBLNK(msgBuf)
309            WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'dxF'
310            CALL READSYMTILE_RS(fName,4,dyF,buf,bi,bj,myThid)
311            iL = ILNBLNK(tmpBuf)
312            WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'dyF'
313            CALL READSYMTILE_RS(fName,5,rA,buf,bi,bj,myThid)
314            iL = ILNBLNK(msgBuf)
315            WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'rA'
316            CALL READSYMTILE_RS(fName,6,xG,buf,bi,bj,myThid)
317            iL = ILNBLNK(tmpBuf)
318            WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'xG'
319            CALL READSYMTILE_RS(fName,7,yG,buf,bi,bj,myThid)
320            iL = ILNBLNK(msgBuf)
321            WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'yG'
322            CALL READSYMTILE_RS(fName,8,dxV,buf,bi,bj,myThid)
323            iL = ILNBLNK(tmpBuf)
324            WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'dxV'
325            CALL READSYMTILE_RS(fName,9,dyU,buf,bi,bj,myThid)
326            iL = ILNBLNK(msgBuf)
327            WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'dyU'
328            CALL READSYMTILE_RS(fName,10,rAz,buf,bi,bj,myThid)
329            iL = ILNBLNK(tmpBuf)
330            WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'rAz'
331            CALL READSYMTILE_RS(fName,11,dxC,buf,bi,bj,myThid)
332            iL = ILNBLNK(msgBuf)
333            WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'dxC'
334            CALL READSYMTILE_RS(fName,12,dyC,buf,bi,bj,myThid)
335            iL = ILNBLNK(tmpBuf)
336            WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'dyC'
337            CALL READSYMTILE_RS(fName,13,rAw,buf,bi,bj,myThid)
338            iL = ILNBLNK(msgBuf)
339            WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'rAw'
340            CALL READSYMTILE_RS(fName,14,rAs,buf,bi,bj,myThid)
341            iL = ILNBLNK(tmpBuf)
342            WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'rAs'
343            CALL READSYMTILE_RS(fName,15,dxG,buf,bi,bj,myThid)
344            iL = ILNBLNK(msgBuf)
345            WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'dxG'
346            CALL READSYMTILE_RS(fName,16,dyG,buf,bi,bj,myThid)
347            iL = ILNBLNK(tmpBuf)
348            WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'dyG'
349    
350            iLen = ILNBLNK(horizGridFile)
351            IF ( iLen.GT.0 ) THEN
352             CALL READSYMTILE_RS(fName,17,angleCosC,buf,bi,bj,myThid)
353             iL = ILNBLNK(msgBuf)
354             WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'AngleCS'
355             CALL READSYMTILE_RS(fName,18,angleSinC,buf,bi,bj,myThid)
356             iL = ILNBLNK(tmpBuf)
357             WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'AngleSN'
358             anglesAreSet = .TRUE.
359            ELSE
360             anglesAreSet = .FALSE.
361            ENDIF
362    
363            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
364         &                      SQUEEZE_RIGHT , myThid)
365    
366         ENDDO         ENDDO
367        ENDDO        ENDDO
 cs-   end block  
       CALL EXCH_UV_XY_RS(RAW,RAS,.FALSE.,myThid)  
368    
369        CALL MDSREADFIELD('DXG.bin',readBinaryPrec,'RL',1,DXG,  1,myThid)  #ifdef ALLOW_MNC
370        CALL MDSREADFIELD('DYG.bin',readBinaryPrec,'RL',1,DYG,  1,myThid)        ENDIF
371  cs-   this block needed by cubed sphere until we write more useful I/O routines  #endif /* ALLOW_MNC */
372        DO bj = myByLo(myThid), myByHi(myThid)  
373         DO bi = myBxLo(myThid), myBxHi(myThid)        _END_MASTER(myThid)
374          DO J = 1,sNy  
375  c        DYG(sNx+1,J,bi,bj)=DYG(1,J,bi,bj)        CALL EXCH_XY_RS(xC,myThid)
376  c        DXG(J,sNy+1,bi,bj)=DXG(J,1,bi,bj)        CALL EXCH_XY_RS(yC,myThid)
377          ENDDO        CALL EXCH_UV_AGRID_3D_RS( dxF, dyF, .FALSE., 1, myThid )
378          CALL EXCH_XY_RS(rA,myThid )
379          CALL EXCH_Z_3D_RS( xG, 1, myThid )
380          CALL EXCH_Z_3D_RS( yG, 1, myThid )
381    C !!! _EXCH_ZUV_XY_R4(dxV, dyU, unSigned, myThid)
382    c     CALL EXCH_Z_3D_RS( dxV, 1, myThid )
383    c     CALL EXCH_Z_3D_RS( dyU, 1, myThid )
384          CALL EXCH_Z_3D_RS( rAz, 1, myThid )
385          CALL EXCH_UV_XY_RS(dxC,dyC,.FALSE.,myThid)
386          CALL EXCH_UV_XY_RS(rAw,rAs,.FALSE.,myThid)
387          CALL EXCH_UV_XY_RS(dyG,dxG,.FALSE.,myThid)
388          CALL EXCH_UV_AGRID_3D_RS(angleSinC,angleCosC,.TRUE., 1, myThid)
389    
390    #endif /* OLD_GRID_IO */
391    
392    C--   Stop if Angle have not been loaded but are needed :
393          _BEGIN_MASTER(myThid)
394          IF ( .NOT.anglesAreSet .AND. use3dCoriolis ) THEN
395            WRITE(msgBuf,'(2A)')
396         &   'INI_CURVILINEAR_GRID: Angle of CurvilinearGrid not set',
397         &   ' but needed for 3-D Coriolis'
398            CALL PRINT_ERROR( msgBuf , myThid)
399            STOP 'ABNORMAL END: S/R INI_CURVILINEAR_GRID'
400          ENDIF
401          _END_MASTER(myThid)
402    
403    c     CALL WRITE_FULLARRAY_RL('dxV',dxV,1,0,0,0,myThid)
404    c     CALL WRITE_FULLARRAY_RL('dyU',dyU,1,0,0,0,myThid)
405    c     CALL WRITE_FULLARRAY_RL('rAz',rAz,1,0,0,0,myThid)
406    c     CALL WRITE_FULLARRAY_RL('xG',xG,1,0,0,0,myThid)
407    c     CALL WRITE_FULLARRAY_RL('yG',yG,1,0,0,0,myThid)
408    
409    C--   Now let's look at all these beasts
410          IF ( debugLevel .GE. debLevB ) THEN
411            CALL PLOT_FIELD_XYRS( xC      , 'Current xC      ', 0, myThid )
412            CALL PLOT_FIELD_XYRS( yC      , 'Current yC      ', 0, myThid )
413            CALL PLOT_FIELD_XYRS( dxF     , 'Current dxF     ', 0, myThid )
414            CALL PLOT_FIELD_XYRS( dyF     , 'Current dyF     ', 0, myThid )
415            CALL PLOT_FIELD_XYRS( rA      , 'Current rA      ', 0, myThid )
416            CALL PLOT_FIELD_XYRS( xG      , 'Current xG      ', 0, myThid )
417            CALL PLOT_FIELD_XYRS( yG      , 'Current yG      ', 0, myThid )
418            CALL PLOT_FIELD_XYRS( dxV     , 'Current dxV     ', 0, myThid )
419            CALL PLOT_FIELD_XYRS( dyU     , 'Current dyU     ', 0, myThid )
420            CALL PLOT_FIELD_XYRS( rAz     , 'Current rAz     ', 0, myThid )
421            CALL PLOT_FIELD_XYRS( dxC     , 'Current dxC     ', 0, myThid )
422            CALL PLOT_FIELD_XYRS( dyC     , 'Current dyC     ', 0, myThid )
423            CALL PLOT_FIELD_XYRS( rAw     , 'Current rAw     ', 0, myThid )
424            CALL PLOT_FIELD_XYRS( rAs     , 'Current rAs     ', 0, myThid )
425            CALL PLOT_FIELD_XYRS( dxG     , 'Current dxG     ', 0, myThid )
426            CALL PLOT_FIELD_XYRS( dyG     , 'Current dyG     ', 0, myThid )
427            CALL PLOT_FIELD_XYRS(angleCosC, 'Current AngleCS ', 0, myThid )
428            CALL PLOT_FIELD_XYRS(angleSinC, 'Current AngleSN ', 0, myThid )
429          ENDIF
430    
431          RETURN
432          END
433    
434    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
435    
436    CBOP
437    C     !ROUTINE: READSYMTILE_RS
438    C     !INTERFACE:
439          SUBROUTINE READSYMTILE_RS(
440         I                           fName, irec,
441         U                           array, buf,
442         I                           bi,bj, myThid )
443    C     !DESCRIPTION: \bv
444    C     *==========================================================*
445    C     | SUBROUTINE READSYMTILE_RS
446    C     *==========================================================*
447    C     *==========================================================*
448    C     \ev
449    
450    C     !USES:
451          IMPLICIT NONE
452    C     === Global variables ===
453    #include "SIZE.h"
454    #include "EEPARAMS.h"
455    #ifdef ALLOW_EXCH2
456    #include "W2_EXCH2_TOPOLOGY.h"
457    #include "W2_EXCH2_PARAMS.h"
458    #endif /* ALLOW_EXCH2 */
459    
460    C     !INPUT/OUTPUT PARAMETERS:
461    C     == Routine arguments ==
462          CHARACTER*(*) fName
463          INTEGER irec
464          _RS array(1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
465    #ifdef ALLOW_EXCH2
466          _RL buf(1:sNx*nSx*nPx+1)
467    #else
468          _RL buf(1:sNx+1,1:sNy+1)
469    #endif /* ALLOW_EXCH2 */
470          INTEGER bi,bj, myThid
471    CEOP
472    
473    C     !LOCAL VARIABLES:
474    C     == Local variables ==
475          INTEGER I,J,dUnit, iLen
476          INTEGER length_of_rec
477          INTEGER MDS_RECLEN
478    #ifdef ALLOW_EXCH2
479          INTEGER TN, dNx, dNy, TBX, TBY, TNX, TNY, II, iBase
480    #endif
481          INTEGER  ILNBLNK
482          EXTERNAL ILNBLNK
483    
484          iLen = ILNBLNK(fName)
485    #ifdef ALLOW_EXCH2
486    C     Figure out offset of tile within face
487          TN  = W2_myTileList(bi)
488          dNx = exch2_mydnx(TN)
489          dNy = exch2_mydny(TN)
490          TBX = exch2_tbasex(TN)
491          TBY = exch2_tbasey(TN)
492          TNX = exch2_tnx(TN)
493          TNY = exch2_tny(TN)
494    
495          CALL MDSFINDUNIT( dUnit, myThid )
496          length_of_rec=MDS_RECLEN( 64, (dNx+1), myThid )
497          OPEN( dUnit, file=fName(1:iLen), status='old',
498         &             access='direct', recl=length_of_rec )
499          J=0
500          iBase=(irec-1)*(dny+1)
501          DO I=1+TBY,sNy+1+TBY
502           READ(dUnit,rec=I+iBase)(buf(ii),ii=1,dNx+1)
503    #ifdef _BYTESWAPIO
504    #ifdef REAL4_IS_SLOW
505           CALL MDS_BYTESWAPR8((dNx+1), buf)
506    #else
507           CALL MDS_BYTESWAPR4((dNx+1), buf)
508    #endif
509    #endif
510           J=J+1
511           DO II=1,sNx+1
512            array(II,J,bi,bj)=buf(II+TBX)
513         ENDDO         ENDDO
514        ENDDO        ENDDO
515  cs-   end block        CLOSE( dUnit )
516        CALL EXCH_UV_XY_RS(DYG,DXG,.FALSE.,myThid)  
517    #else /* ALLOW_EXCH2 */
518    
519          CALL MDSFINDUNIT( dUnit, myThid )
520          length_of_rec=MDS_RECLEN( 64, (sNx+1)*(sNy+1), myThid )
521          OPEN( dUnit, file=fName(1:iLen), status='old',
522         &             access='direct', recl=length_of_rec )
523          READ(dUnit,rec=irec) buf
524          CLOSE( dUnit )
525    
526    #ifdef _BYTESWAPIO
527    #ifdef REAL4_IS_SLOW
528          CALL MDS_BYTESWAPR8((sNx+1)*(sNy+1), buf)
529    #else
530          CALL MDS_BYTESWAPR4((sNx+1)*(sNy+1), buf)
531    #endif
532    #endif
533    
534          DO J=1,sNy+1
535           DO I=1,sNx+1
536            array(I,J,bi,bj)=buf(I,J)
537           ENDDO
538          ENDDO
539    c       write(0,*) irec,buf(1,1),array(1,1,1,1)
540    
541  c     write(10) XC  #endif /* ALLOW_EXCH2 */
 c     write(10) YC  
 c     write(10) DXF  
 c     write(10) DYF  
 c     write(10) RA  
 c     write(10) XG  
 c     write(10) YG  
 c     write(10) DXV  
 c     write(10) DYU  
 c     write(10) RAZ  
 c     write(10) DXC  
 c     write(10) DYC  
 c     write(10) RAW  
 c     write(10) RAS  
 c     write(10) DXG  
 c     write(10) DYG  
542    
543        RETURN        RETURN
544        END        END
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.32
  ViewVC Help
Powered by ViewVC 1.1.22