C $Header: /u/gcmpack/MITgcm/model/src/ini_curvilinear_grid.F,v 1.20 2005/07/13 00:41:44 jmc Exp $
C $Name: $
#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"
CBOP
C !ROUTINE: INI_CURVILINEAR_GRID
C !INTERFACE:
SUBROUTINE INI_CURVILINEAR_GRID( myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE INI_CURVILINEAR_GRID
C | o Initialise curvilinear coordinate system
C *==========================================================*
C | Curvilinear grid settings are read from a file rather
C | than coded in-line as for cartesian and spherical polar.
C | This is more general but you have to create the grid
C | yourself.
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#ifdef ALLOW_EXCH2
#include "W2_EXCH2_TOPOLOGY.h"
#include "W2_EXCH2_PARAMS.h"
#endif
#ifdef ALLOW_MNC
#include "MNC_PARAMS.h"
#endif
#ifndef ALLOW_EXCH2
C- note: default is to use "new" grid files (OLD_GRID_IO undef) with EXCH2
C but can still use (on 1 cpu) OLD_GRID_IO and EXCH2 independently
#define OLD_GRID_IO
#endif /* ALLOW_EXCH2 */
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
C myThid - Number of this instance of INI_CARTESIAN_GRID
INTEGER myThid
C !LOCAL VARIABLES:
C == Local variables ==
INTEGER bi,bj, myTile, myiter
INTEGER I,J
CHARACTER*(MAX_LEN_FNAM) fName
#ifdef ALLOW_MNC
CHARACTER*(80) mncFn
#endif
#ifdef ALLOW_EXCH2
_RL buf(sNx*nSx*nPx+1)
#else
_RL buf(sNx+1,sNy+1)
#endif
INTEGER iG, iL, iLen
CHARACTER*(MAX_LEN_MBUF) msgBuf, tmpBuf
INTEGER ILNBLNK
EXTERNAL
CEOP
C-- Set everything to zero everywhere
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
DO J=1-Oly,sNy+Oly
DO I=1-Olx,sNx+Olx
XC(i,j,bi,bj)=0.
YC(i,j,bi,bj)=0.
XG(i,j,bi,bj)=0.
YG(i,j,bi,bj)=0.
DXC(i,j,bi,bj)=0.
DYC(i,j,bi,bj)=0.
DXG(i,j,bi,bj)=0.
DYG(i,j,bi,bj)=0.
DXF(i,j,bi,bj)=0.
DYF(i,j,bi,bj)=0.
DXV(i,j,bi,bj)=0.
DYU(i,j,bi,bj)=0.
RA(i,j,bi,bj)=0.
RAZ(i,j,bi,bj)=0.
RAW(i,j,bi,bj)=0.
RAS(i,j,bi,bj)=0.
tanPhiAtU(i,j,bi,bj)=0.
tanPhiAtV(i,j,bi,bj)=0.
angleCosC(i,j,bi,bj)=1.
angleSinC(i,j,bi,bj)=0.
cosFacU(J,bi,bj)=1.
cosFacV(J,bi,bj)=1.
sqcosFacU(J,bi,bj)=1.
sqcosFacV(J,bi,bj)=1.
ENDDO
ENDDO
ENDDO
ENDDO
#ifdef ALLOW_MNC
IF (useMNC .AND. readgrid_mnc) THEN
_BEGIN_MASTER(myThid)
DO i = 1,80
mncFn(i:i) = ' '
ENDDO
write(mncFn,'(a)') 'mitgrid'
DO i = 1,MAX_LEN_MBUF
msgBuf(i:i) = ' '
ENDDO
WRITE(msgBuf,'(2A)') msgBuf,' ; Reading grid info using MNC'
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT , myThid)
CALL MNC_FILE_CLOSE_ALL_MATCHING(mncFn, myThid)
CALL MNC_CW_SET_UDIM(mncFn, 1, myThid)
CALL MNC_CW_RS_R('R',mncFn,0,0,'XC', XC, myThid)
CALL MNC_CW_RS_R('R',mncFn,0,0,'XG', XG, myThid)
CALL MNC_CW_RS_R('R',mncFn,0,0,'YC', YC, myThid)
CALL MNC_CW_RS_R('R',mncFn,0,0,'YG', YG, myThid)
CALL MNC_CW_RS_R('R',mncFn,0,0,'dxC',DXC, myThid)
CALL MNC_CW_RS_R('R',mncFn,0,0,'dyC',DYC, myThid)
CALL MNC_CW_RS_R('R',mncFn,0,0,'dxF',DXF, myThid)
CALL MNC_CW_RS_R('R',mncFn,0,0,'dyF',DYF, myThid)
CALL MNC_CW_RS_R('R',mncFn,0,0,'dxG',DXG, myThid)
CALL MNC_CW_RS_R('R',mncFn,0,0,'dyG',DYG, myThid)
CALL MNC_CW_RS_R('R',mncFn,0,0,'dxV',DXV, myThid)
CALL MNC_CW_RS_R('R',mncFn,0,0,'dyU',DYU, myThid)
CALL MNC_CW_RS_R('R',mncFn,0,0,'rA', RA, myThid)
CALL MNC_CW_RS_R('R',mncFn,0,0,'rAz',RAZ, myThid)
CALL MNC_CW_RS_R('R',mncFn,0,0,'rAw',RAW, myThid)
CALL MNC_CW_RS_R('R',mncFn,0,0,'rAs',RAS, myThid)
_END_MASTER(myThid)
CALL EXCH_XY_RS(XC,myThid)
CALL EXCH_XY_RS(YC,myThid)
#ifdef HRCUBE
CALL EXCH_XY_RS(DXF,myThid)
CALL EXCH_XY_RS(DYF,myThid)
#endif
CALL EXCH_XY_RS(RA,myThid )
CALL EXCH_Z_XY_RS(XG,myThid)
CALL EXCH_Z_XY_RS(YG,myThid)
CALL EXCH_Z_XY_RS(RAZ,myThid)
CALL EXCH_UV_XY_RS(DXC,DYC,.FALSE.,myThid)
CALL EXCH_UV_XY_RS(RAW,RAS,.FALSE.,myThid)
CALL EXCH_UV_XY_RS(DYG,DXG,.FALSE.,myThid)
CALL EXCH_UV_AGRID_XY_RS(angleSinC,angleCosC,.TRUE.,myThid)
ELSE
#endif
C Here we make no assumptions about grid symmetry and simply
C read the raw grid data from files
#ifdef OLD_GRID_IO
C- Cell centered quantities
CALL MDSREADFIELD('LONC.bin',readBinaryPrec,'RS',1,XC, 1,myThid)
CALL MDSREADFIELD('LATC.bin',readBinaryPrec,'RS',1,YC, 1,myThid)
_EXCH_XY_R4(XC,myThid)
_EXCH_XY_R4(YC,myThid)
CALL MDSREADFIELD('DXF.bin',readBinaryPrec,'RS',1,DXF, 1,myThid)
CALL MDSREADFIELD('DYF.bin',readBinaryPrec,'RS',1,DYF, 1,myThid)
C !!! _EXCH_OUV_XY_R4(DXF, DYF, unSigned, myThid )
cs! this is not correct! <= need paired exchange for DXF,DYF
_EXCH_XY_R4(DXF,myThid)
_EXCH_XY_R4(DYF,myThid)
IF (useCubedSphereExchange) THEN
cs! fix overlaps:
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
DO j=1,sNy
DO i=1,Olx
DXF(1-i,j,bi,bj)=DXF(i,j,bi,bj)
DXF(sNx+i,j,bi,bj)=DXF(sNx+1-i,j,bi,bj)
DYF(1-i,j,bi,bj)=DYF(i,j,bi,bj)
DYF(sNx+i,j,bi,bj)=DYF(sNx+1-i,j,bi,bj)
ENDDO
ENDDO
DO j=1,Oly
DO i=1,sNx
DXF(i,1-j,bi,bj)=DXF(i,j,bi,bj)
DXF(i,sNy+j,bi,bj)=DXF(i,sNy+1-j,bi,bj)
DYF(i,1-j,bi,bj)=DYF(i,j,bi,bj)
DYF(i,sNy+j,bi,bj)=DYF(i,sNy+1-j,bi,bj)
ENDDO
ENDDO
ENDDO
ENDDO
ENDIF
cs
CALL MDSREADFIELD('RA.bin',readBinaryPrec,'RS',1,RA, 1,myThid)
_EXCH_XY_R4(RA,myThid )
C- Corner quantities
C *********** this are not degbugged ************
CALL MDSREADFIELD('LONG.bin',readBinaryPrec,'RS',1,XG, 1,myThid)
CALL MDSREADFIELD('LATG.bin',readBinaryPrec,'RS',1,YG, 1,myThid)
IF (useCubedSphereExchange) THEN
cs- this block needed by cubed sphere until we write more useful I/O routines
bi=3
bj=1
YG(1,sNy+1,bj,1)=YG(1,1,bi,1)
bj=bj+2
YG(1,sNy+1,bj,1)=YG(1,1,bi,1)
bj=bj+2
YG(1,sNy+1,bj,1)=YG(1,1,bi,1)
bi=6
bj=2
YG(sNx+1,1,bj,1)=YG(1,1,bi,1)
bj=bj+2
YG(sNx+1,1,bj,1)=YG(1,1,bi,1)
bj=bj+2
YG(sNx+1,1,bj,1)=YG(1,1,bi,1)
cs- end block
ENDIF
CALL EXCH_Z_XY_RS(XG,myThid)
CALL EXCH_Z_XY_RS(YG,myThid)
CALL MDSREADFIELD('DXV.bin',readBinaryPrec,'RS',1,DXV, 1,myThid)
CALL MDSREADFIELD('DYU.bin',readBinaryPrec,'RS',1,DYU, 1,myThid)
cs- this block needed by cubed sphere until we write more useful I/O routines
C !!! _EXCH_ZUV_XY_R4(DXV, DYU, unSigned, myThid)
cs! this is not correct <= need paired exchange for dxv,dyu
IF (.NOT.useCubedSphereExchange) THEN
CALL EXCH_Z_XY_RS(DXV,myThid)
CALL EXCH_Z_XY_RS(DYU,myThid)
ELSE
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
cs! fix overlaps:
DO j=1,sNy
DO i=1,Olx
DXV(1-i,j,bi,bj)=DXV(1+i,j,bi,bj)
DXV(sNx+i,j,bi,bj)=DXV(i,j,bi,bj)
DYU(1-i,j,bi,bj)=DYU(1+i,j,bi,bj)
DYU(sNx+i,j,bi,bj)=DYU(i,j,bi,bj)
ENDDO
ENDDO
DO j=1,Oly
DO i=1-Olx,sNx+Olx
DXV(i,1-j,bi,bj)=DXV(i,1+j,bi,bj)
DXV(i,sNy+j,bi,bj)=DXV(i,j,bi,bj)
DYU(i,1-j,bi,bj)=DYU(i,1+j,bi,bj)
DYU(i,sNy+j,bi,bj)=DYU(i,j,bi,bj)
ENDDO
ENDDO
ENDDO
ENDDO
cs- end block
ENDIF
CALL MDSREADFIELD('RAZ.bin',readBinaryPrec,'RS',1,RAZ, 1,myThid)
IF (useCubedSphereExchange) THEN
cs- this block needed by cubed sphere until we write more useful I/O routines
CALL EXCH_Z_XY_RS(RAZ , myThid )
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
RAZ(sNx+1,1,bi,bj)=RAZ(1,1,bi,bj)
RAZ(1,sNy+1,bi,bj)=RAZ(1,1,bi,bj)
ENDDO
ENDDO
cs- end block
ENDIF
CALL EXCH_Z_XY_RS(RAZ,myThid)
C- Staggered (u,v pairs) quantities
CALL MDSREADFIELD('DXC.bin',readBinaryPrec,'RS',1,DXC, 1,myThid)
CALL MDSREADFIELD('DYC.bin',readBinaryPrec,'RS',1,DYC, 1,myThid)
CALL EXCH_UV_XY_RS(DXC,DYC,.FALSE.,myThid)
CALL MDSREADFIELD('RAW.bin',readBinaryPrec,'RS',1,RAW, 1,myThid)
CALL MDSREADFIELD('RAS.bin',readBinaryPrec,'RS',1,RAS, 1,myThid)
IF (useCubedSphereExchange) THEN
cs- this block needed by cubed sphere until we write more useful I/O routines
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
DO J = 1,sNy
c RAW(sNx+1,J,bi,bj)=RAW(1,J,bi,bj)
c RAS(J,sNy+1,bi,bj)=RAS(J,1,bi,bj)
ENDDO
ENDDO
ENDDO
cs- end block
ENDIF
CALL EXCH_UV_XY_RS(RAW,RAS,.FALSE.,myThid)
CALL MDSREADFIELD('DXG.bin',readBinaryPrec,'RS',1,DXG, 1,myThid)
CALL MDSREADFIELD('DYG.bin',readBinaryPrec,'RS',1,DYG, 1,myThid)
IF (useCubedSphereExchange) THEN
cs- this block needed by cubed sphere until we write more useful I/O routines
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
DO J = 1,sNy
c DYG(sNx+1,J,bi,bj)=DYG(1,J,bi,bj)
c DXG(J,sNy+1,bi,bj)=DXG(J,1,bi,bj)
ENDDO
ENDDO
ENDDO
cs- end block
ENDIF
CALL EXCH_UV_XY_RS(DYG,DXG,.FALSE.,myThid)
CALL EXCH_UV_AGRID_XY_RS(angleSinC,angleCosC,.TRUE.,myThid)
c write(10) XC
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
#else /* ifndef OLD_GRID_IO */
C-- Only do I/O if I am the master thread
_BEGIN_MASTER(myThid)
DO bj = 1,nSy
DO bi = 1,nSx
iG=bi+(myXGlobalLo-1)/sNx
WRITE(tmpBuf,'(A,I4)') 'tile:',iG
#ifdef ALLOW_EXCH2
myTile = W2_myTileList(bi)
WRITE(tmpBuf,'(A,I4)') 'tile:',myTile
iG = exch2_myface(myTile)
#endif
iLen = ILNBLNK(horizGridFile)
IF ( iLen .EQ. 0 ) THEN
WRITE(fName,'("tile",I3.3,".mitgrid")') iG
ELSE
WRITE(fName,'(2A,I3.3,A)') horizGridFile(1:iLen),
& '.face',iG,'.bin'
ENDIF
iLen = ILNBLNK(fName)
iL = ILNBLNK(tmpBuf)
WRITE(msgBuf,'(3A)') tmpBuf(1:iL),
& ' ; Read from file ',fName(1:iLen)
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT , myThid)
WRITE(msgBuf,'(A)') ' =>'
CALL READSYMTILE_RS(fName,1,XC,bi,bj,buf,myThid)
iL = ILNBLNK(msgBuf)
WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'XC'
CALL READSYMTILE_RS(fName,2,YC,bi,bj,buf,myThid)
iL = ILNBLNK(tmpBuf)
WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'YC'
CALL READSYMTILE_RS(fName,3,DXF,bi,bj,buf,myThid)
iL = ILNBLNK(msgBuf)
WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'DXF'
CALL READSYMTILE_RS(fName,4,DYF,bi,bj,buf,myThid)
iL = ILNBLNK(tmpBuf)
WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'DYF'
CALL READSYMTILE_RS(fName,5,RA,bi,bj,buf,myThid)
iL = ILNBLNK(msgBuf)
WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'RA'
CALL READSYMTILE_RS(fName,6,XG,bi,bj,buf,myThid)
iL = ILNBLNK(tmpBuf)
WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'XG'
CALL READSYMTILE_RS(fName,7,YG,bi,bj,buf,myThid)
iL = ILNBLNK(msgBuf)
WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'YG'
CALL READSYMTILE_RS(fName,8,DXV,bi,bj,buf,myThid)
iL = ILNBLNK(tmpBuf)
WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'DXV'
CALL READSYMTILE_RS(fName,9,DYU,bi,bj,buf,myThid)
iL = ILNBLNK(msgBuf)
WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'DYU'
CALL READSYMTILE_RS(fName,10,RAZ,bi,bj,buf,myThid)
iL = ILNBLNK(tmpBuf)
WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'RAZ'
CALL READSYMTILE_RS(fName,11,DXC,bi,bj,buf,myThid)
iL = ILNBLNK(msgBuf)
WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'DXC'
CALL READSYMTILE_RS(fName,12,DYC,bi,bj,buf,myThid)
iL = ILNBLNK(tmpBuf)
WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'DYC'
CALL READSYMTILE_RS(fName,13,RAW,bi,bj,buf,myThid)
iL = ILNBLNK(msgBuf)
WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'RAW'
CALL READSYMTILE_RS(fName,14,RAS,bi,bj,buf,myThid)
iL = ILNBLNK(tmpBuf)
WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'RAS'
CALL READSYMTILE_RS(fName,15,DXG,bi,bj,buf,myThid)
iL = ILNBLNK(msgBuf)
WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'DXG'
CALL READSYMTILE_RS(fName,16,DYG,bi,bj,buf,myThid)
iL = ILNBLNK(tmpBuf)
WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'DYG'
iLen = ILNBLNK(horizGridFile)
IF ( iLen.GT.0 ) THEN
CALL READSYMTILE_RS(fName,17,angleCosC,bi,bj,buf,myThid)
iL = ILNBLNK(msgBuf)
WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'AngleCS'
CALL READSYMTILE_RS(fName,18,angleSinC,bi,bj,buf,myThid)
iL = ILNBLNK(tmpBuf)
WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'AngleSN'
ENDIF
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT , myThid)
ENDDO
ENDDO
_END_MASTER(myThid)
CALL EXCH_XY_RS(XC,myThid)
CALL EXCH_XY_RS(YC,myThid)
C !!! _EXCH_OUV_XY_R4(DXF, DYF, unSigned, myThid )
#ifdef HRCUBE
CALL EXCH_XY_RS(DXF,myThid)
CALL EXCH_XY_RS(DYF,myThid)
#endif
CALL EXCH_XY_RS(RA,myThid )
CALL EXCH_Z_XY_RS(XG,myThid)
CALL EXCH_Z_XY_RS(YG,myThid)
C !!! _EXCH_ZUV_XY_R4(DXV, DYU, unSigned, myThid)
c CALL EXCH_Z_XY_RS(DXV,myThid)
c CALL EXCH_Z_XY_RS(DYU,myThid)
CALL EXCH_Z_XY_RS(RAZ,myThid)
CALL EXCH_UV_XY_RS(DXC,DYC,.FALSE.,myThid)
CALL EXCH_UV_XY_RS(RAW,RAS,.FALSE.,myThid)
CALL EXCH_UV_XY_RS(DYG,DXG,.FALSE.,myThid)
CALL EXCH_UV_AGRID_XY_RS(angleSinC,angleCosC,.TRUE.,myThid)
#endif /* OLD_GRID_IO */
#ifdef ALLOW_MNC
ENDIF
#endif /* ALLOW_MNC */
c CALL WRITE_FULLARRAY_RL('DXV',DXV,1,0,0,0,myThid)
c CALL WRITE_FULLARRAY_RL('DYU',DYU,1,0,0,0,myThid)
c CALL WRITE_FULLARRAY_RL('RAZ',RAZ,1,0,0,0,myThid)
c CALL WRITE_FULLARRAY_RL('XG',XG,1,0,0,0,myThid)
c CALL WRITE_FULLARRAY_RL('YG',YG,1,0,0,0,myThid)
C-- Require that 0 <= longitude < 360 if using exf package
#ifdef ALLOW_EXF
DO bj = 1,nSy
DO bi = 1,nSx
DO J=1-Oly,sNy+Oly
DO I=1-Olx,sNx+Olx
IF (XC(i,j,bi,bj).lt.0.) XC(i,j,bi,bj)=XC(i,j,bi,bj)+360.
IF (XG(i,j,bi,bj).lt.0.) XG(i,j,bi,bj)=XG(i,j,bi,bj)+360.
ENDDO
ENDDO
ENDDO
ENDDO
#endif /* ALLOW_EXF */
C-- Now let's look at all these beasts
IF ( debugLevel .GE. debLevB ) THEN
myiter = 1
CALL PLOT_FIELD_XYRL( XC , 'Current XC ' ,
& myIter, myThid )
CALL PLOT_FIELD_XYRL( YC , 'Current YC ' ,
& myIter, myThid )
CALL PLOT_FIELD_XYRL( DXF , 'Current DXF ' ,
& myIter, myThid )
CALL PLOT_FIELD_XYRL( XC , 'Current XC ' ,
& myIter, myThid )
CALL PLOT_FIELD_XYRL( DYF , 'Current DYF ' ,
& myIter, myThid )
CALL PLOT_FIELD_XYRL( RA , 'Current RA ' ,
& myIter, myThid )
CALL PLOT_FIELD_XYRL( XG , 'Current XG ' ,
& myIter, myThid )
CALL PLOT_FIELD_XYRL( YG , 'Current YG ' ,
& myIter, myThid )
CALL PLOT_FIELD_XYRL( DXV , 'Current DXV ' ,
& myIter, myThid )
CALL PLOT_FIELD_XYRL( DYU , 'Current DYU ' ,
& myIter, myThid )
CALL PLOT_FIELD_XYRL( RAZ , 'Current RAZ ' ,
& myIter, myThid )
CALL PLOT_FIELD_XYRL( DXC , 'Current DXC ' ,
& myIter, myThid )
CALL PLOT_FIELD_XYRL( DYC , 'Current DYC ' ,
& myIter, myThid )
CALL PLOT_FIELD_XYRL( RAW , 'Current RAW ' ,
& myIter, myThid )
CALL PLOT_FIELD_XYRL( RAS , 'Current RAS ' ,
& myIter, myThid )
CALL PLOT_FIELD_XYRL( DXG , 'Current DXG ' ,
& myIter, myThid )
CALL PLOT_FIELD_XYRL( DYG , 'Current DYG ' ,
& myIter, myThid )
CALL PLOT_FIELD_XYRL(angleCosC, 'Current AngleCS ' ,
& myIter, myThid )
CALL PLOT_FIELD_XYRL(angleSinC, 'Current AngleSN ' ,
& myIter, myThid )
ENDIF
RETURN
END
C --------------------------------------------------------------------------
SUBROUTINE READSYMTILE_RS(fName,irec,array,bi,bj,buf,myThid)
C /==========================================================\
C | SUBROUTINE READSYMTILE_RS |
C |==========================================================|
C \==========================================================/
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#ifdef ALLOW_EXCH2
#include "W2_EXCH2_TOPOLOGY.h"
#include "W2_EXCH2_PARAMS.h"
#endif /* ALLOW_EXCH2 */
C == Routine arguments ==
CHARACTER*(*) fName
INTEGER irec
_RS array(1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
INTEGER bi,bj,myThid
#ifdef ALLOW_EXCH2
_RL buf(1:sNx*nSx*nPx+1)
#else
_RL buf(1:sNx+1,1:sNy+1)
#endif /* ALLOW_EXCH2 */
C == Local variables ==
INTEGER I,J,dUnit, iLen
INTEGER length_of_rec
INTEGER MDS_RECLEN
INTEGER TN, dNx, dNy, TBX, TBY, TNX, TNY, II, iBase
INTEGER ILNBLNK
EXTERNAL
iLen = ILNBLNK(fName)
#ifdef ALLOW_EXCH2
C Figure out offset of tile within face
TN = W2_myTileList(bi)
dNx = exch2_mydnx(TN)
dNy = exch2_mydny(TN)
TBX = exch2_tbasex(TN)
TBY = exch2_tbasey(TN)
TNX = exch2_tnx(TN)
TNY = exch2_tny(TN)
CALL MDSFINDUNIT( dUnit, myThid )
length_of_rec=MDS_RECLEN( 64, (dNx+1), myThid )
OPEN( dUnit, file=fName(1:iLen), status='old',
& access='direct', recl=length_of_rec )
J=0
iBase=(irec-1)*(dny+1)
DO I=1+TBY,sNy+1+TBY
READ(dUnit,rec=I+iBase)(buf(ii),ii=1,dNx+1)
#ifdef _BYTESWAPIO
#ifdef REAL4_IS_SLOW
CALL MDS_BYTESWAPR8((dNx+1), buf)
#else
CALL MDS_BYTESWAPR4((dNx+1), buf)
#endif
#endif
J=J+1
DO II=1,sNx+1
array(II,J,bi,bj)=buf(II+TBX)
ENDDO
ENDDO
CLOSE( dUnit )
#else /* ALLOW_EXCH2 */
CALL MDSFINDUNIT( dUnit, myThid )
length_of_rec=MDS_RECLEN( 64, (sNx+1)*(sNy+1), myThid )
OPEN( dUnit, file=fName(1:iLen), status='old',
& access='direct', recl=length_of_rec )
READ(dUnit,rec=irec) buf
CLOSE( dUnit )
#ifdef _BYTESWAPIO
#ifdef REAL4_IS_SLOW
CALL MDS_BYTESWAPR8((sNx+1)*(sNy+1), buf)
#else
CALL MDS_BYTESWAPR4((sNx+1)*(sNy+1), buf)
#endif
#endif
DO J=1,sNy+1
DO I=1,sNx+1
array(I,J,bi,bj)=buf(I,J)
ENDDO
ENDDO
c write(0,*) irec,buf(1,1),array(1,1,1,1)
#endif /* ALLOW_EXCH2 */
RETURN
END