C $Header: /u/gcmpack/MITgcm/verification/global_ocean.cubed32x32x30/code/ini_curvilinear_grid.F,v 1.1 2002/12/12 13:03:34 cheisey Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

      SUBROUTINE INI_CURVILINEAR_GRID( myThid )
C     /==========================================================\
C     | SUBROUTINE INI_CURVILINEAR_GRID                          |
C     | o Initialise curvilinear coordinate system               |
C     |==========================================================|
C     \==========================================================/
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"

C     == Routine arguments ==
C     myThid -  Number of this instance of INI_CARTESIAN_GRID
      INTEGER myThid

C     == Local variables ==
      INTEGER bi,bj,I,J
      CHARACTER*(12) ff

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.
          cosFacU(J,bi,bj)=1.
          cosFacV(J,bi,bj)=1.
          sqcosFacU(J,bi,bj)=1.
          sqcosFacV(J,bi,bj)=1.
         ENDDO
        ENDDO

        CALL READSYMTILE_RS('DXF.bin',DXF,0,0,bi,bj,myThid)
        CALL READSYMTILE_RS('DYF.bin',DYF,0,0,bi,bj,myThid)
        CALL READSYMTILE_RS('RA.bin',RA,0,0,bi,bj,myThid)
        CALL READSYMTILE_RS('DXV.bin',DXV,1,1,bi,bj,myThid)
        CALL READSYMTILE_RS('DYU.bin',DYU,1,1,bi,bj,myThid)
        CALL READSYMTILE_RS('RAZ.bin',RAZ,1,1,bi,bj,myThid)
        CALL READSYMTILE_RS('DXC.bin',DXC,1,0,bi,bj,myThid)
        CALL READSYMTILE_RS('DYC.bin',DYC,0,1,bi,bj,myThid)
        CALL READSYMTILE_RS('RAW.bin',RAW,1,0,bi,bj,myThid)
        CALL READSYMTILE_RS('RAS.bin',RAS,0,1,bi,bj,myThid)
        CALL READSYMTILE_RS('DXG.bin',DXG,0,1,bi,bj,myThid)
        CALL READSYMTILE_RS('DYG.bin',DYG,1,0,bi,bj,myThid)

        write(ff(1:12),'(a,i3.3,a)') 'LONC.',bi+(bj-1)*nSx,'.bin'
        CALL READSYMTILE_RS(ff,XC,0,0,bi,bj,myThid)
        write(ff(1:12),'(a,i3.3,a)') 'LATC.',bi+(bj-1)*nSx,'.bin'
        CALL READSYMTILE_RS(ff,YC,0,0,bi,bj,myThid)
        write(ff(1:12),'(a,i3.3,a)') 'LONG.',bi+(bj-1)*nSx,'.bin'
        CALL READSYMTILE_RS(ff,XG,1,1,bi,bj,myThid)
        write(ff(1:12),'(a,i3.3,a)') 'LATG.',bi+(bj-1)*nSx,'.bin'
        CALL READSYMTILE_RS(ff,YG,1,1,bi,bj,myThid)
        write(0,*) 'icg: ',bi,bj

       ENDDO ! bi
      ENDDO ! bj

C     Here we make no assumptions about grid symmetry and simply
C     read the raw grid data from files

C-    Cell centered quantities
c     CALL MDSREADFIELD('LONC.bin',readBinaryPrec,'RL',1,XC,  1,myThid)
c     CALL MDSREADFIELD('LATC.bin',readBinaryPrec,'RL',1,YC,  1,myThid)
      _EXCH_XY_R4(XC,myThid)
      _EXCH_XY_R4(YC,myThid)
        write(0,*) 'icg: xc,yc'

c     CALL MDSREADFIELD('DXF.bin',readBinaryPrec,'RL',1,DXF,  1,myThid)
c     CALL MDSREADFIELD('DYF.bin',readBinaryPrec,'RL',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)

c     CALL MDSREADFIELD('RA.bin',readBinaryPrec,'RL',1,RA,  1,myThid)
      _EXCH_XY_R4(RA,myThid )
        write(0,*) 'icg: ra'

C-    Corner quantities
C       *********** this are not degbugged ************
c     CALL MDSREADFIELD('LONG.bin',readBinaryPrec,'RL',1,XG,  1,myThid)
c     CALL MDSREADFIELD('LATG.bin',readBinaryPrec,'RL',1,YG,  1,myThid)
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
      CALL EXCH_Z_XY_RS(XG,myThid)
      CALL EXCH_Z_XY_RS(YG,myThid)
        write(0,*) 'icg: xg,yg'

c     CALL MDSREADFIELD('DXV.bin',readBinaryPrec,'RL',1,DXV,  1,myThid)
c     CALL MDSREADFIELD('DYU.bin',readBinaryPrec,'RL',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
      CALL EXCH_Z_XY_RS(DXV,myThid)
      CALL EXCH_Z_XY_RS(DYU,myThid)
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DXV(sNx+1,1,bi,bj)=DXV(1,1,bi,bj)
        DXV(1,sNy+1,bi,bj)=DXV(1,1,bi,bj)
        DYU(sNx+1,1,bi,bj)=DYU(1,1,bi,bj)
        DYU(1,sNy+1,bi,bj)=DYU(1,1,bi,bj)
       ENDDO
      ENDDO
cs-   end block
C !!! _EXCH_ZUV_XY_R4(DXV, DYU, unSigned, myThid)
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)

c     CALL MDSREADFIELD('RAZ.bin',readBinaryPrec,'RL',1,RAZ,  1,myThid)
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
      CALL EXCH_Z_XY_RS(RAZ,myThid)

C-    Staggered (u,v pairs) quantities
c     CALL MDSREADFIELD('DXC.bin',readBinaryPrec,'RL',1,DXC,  1,myThid)
c     CALL MDSREADFIELD('DYC.bin',readBinaryPrec,'RL',1,DYC,  1,myThid)
      CALL EXCH_UV_XY_RS(DXC,DYC,.FALSE.,myThid)

c     CALL MDSREADFIELD('RAW.bin',readBinaryPrec,'RL',1,RAW,  1,myThid)
c     CALL MDSREADFIELD('RAS.bin',readBinaryPrec,'RL',1,RAS,  1,myThid)
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
      CALL EXCH_UV_XY_RS(RAW,RAS,.FALSE.,myThid)
        write(0,*) 'icg: raw,ras'

c     CALL MDSREADFIELD('DXG.bin',readBinaryPrec,'RL',1,DXG,  1,myThid)
c     CALL MDSREADFIELD('DYG.bin',readBinaryPrec,'RL',1,DYG,  1,myThid)
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
      CALL EXCH_UV_XY_RS(DYG,DXG,.FALSE.,myThid)

C Debugging:
      write(10) XC
      write(10) YC
      write(10) DXF
      write(10) DYF
      write(10) RA
      write(10) XG
      write(10) YG
      write(10) DXV
      write(10) DYU
      write(10) RAZ
      write(10) DXC
      write(10) DYC
      write(10) RAW
      write(10) RAS
      write(10) DXG
      write(10) DYG
      write(10) DYG
#ifdef X_COMMENTS_X
C Use the following M-script to read:
 fid=fopen('fort.10','r')
 n=24o=3
 xc=f77read(fid,n+2*o n+2*o 6,'real*8')
 yc=f77read(fid,n+2*o n+2*o 6,'real*8')
 dxf=f77read(fid,n+2*o n+2*o 6,'real*8')
 dyf=f77read(fid,n+2*o n+2*o 6,'real*8')
 ra=f77read(fid,n+2*o n+2*o 6,'real*8')
 xg=f77read(fid,n+2*o n+2*o 6,'real*8')
 yg=f77read(fid,n+2*o n+2*o 6,'real*8')
 dxv=f77read(fid,n+2*o n+2*o 6,'real*8')
 dyu=f77read(fid,n+2*o n+2*o 6,'real*8')
 raz=f77read(fid,n+2*o n+2*o 6,'real*8')
 dxc=f77read(fid,n+2*o n+2*o 6,'real*8')
 dyc=f77read(fid,n+2*o n+2*o 6,'real*8')
 raw=f77read(fid,n+2*o n+2*o 6,'real*8')
 ras=f77read(fid,n+2*o n+2*o 6,'real*8')
 dxg=f77read(fid,n+2*o n+2*o 6,'real*8')
 dyg=f77read(fid,n+2*o n+2*o 6,'real*8')
 fclose(fid)
#endif /* X_COMMENTS_X */

        write(0,*) 'icg: end'

      RETURN
      END


      SUBROUTINE READSYMTILE_RS(fileName,array,Xol,Yol,bi,bj,myThid)
C     /==========================================================\
C     | SUBROUTINE READSYMTILE_RS                                |
C     |==========================================================|
C     \==========================================================/
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"

C     == Routine arguments ==
      CHARACTER*(*) fileName
      _RS array(1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
      INTEGER Xol,Yol,bi,bj,myThid

C     == Local variables ==
      INTEGER I,J

      _BEGIN_MASTER(myThid)
      OPEN(36,FILE=fileName,STATUS='OLD',ACCESS='DIRECT',
#ifdef REAL4_IS_SLOW
     &     RECL=((sNx+Xol)*(sNy+Yol))*WORDLENGTH*2 )
#else
     &     RECL=((sNx+Xol)*(sNy+Yol))*WORDLENGTH )
#endif
      READ(36,REC=1) ((array(I,J,bi,bj),I=1,sNx+Xol),J=1,sNy+Yol)
      CLOSE(36)
#ifdef _BYTESWAPIO
#ifdef REAL4_IS_SLOW
      CALL MDS_BYTESWAPR8((sNx+2*Olx)*(sNy+2*Oly),
     &          array(1-Olx,1-Oly,bi,bj))
#else
      CALL MDS_BYTESWAPR4((sNx+2*Olx)*(sNy+2*Oly),
     &          array(1-Olx,1-Oly,bi,bj))
#endif
#endif

C Avoid broken exchanges
      DO J=1,Oly
       DO I=1,sNx+Xol
        array(I,1-J,bi,bj)=array(I,J+Yol,bi,bj)
       ENDDO
      ENDDO
      DO J=1+Yol,Oly
       DO I=1,sNx+Xol
        array(I,sNy+J,bi,bj)=array(I,sNy+1-J+Yol,bi,bj)
       ENDDO
      ENDDO
      DO J=1,sNy+Yol
c     DO J=1-Oly,sNy+Oly
       DO I=1,Olx
        array(1-I,J,bi,bj)=array(I+Xol,J,bi,bj)
       ENDDO
      ENDDO
      DO J=1,sNy+Yol
c     DO J=1-Oly,sNy+Oly
       DO I=1+Xol,Olx
        array(sNx+I,J,bi,bj)=array(sNy+1-I+Xol,J,bi,bj)
       ENDDO
      ENDDO
      _END_MASTER(myThid)

      RETURN
      END