C==== suggestion for LOQ two dimensional data storage  RKH 19/2/90
C====
C==== file structure:
C====
C==== long_title1                  fed from RAW file to identify sample
C==== ixcode,xcaption,xcaption2    GENIE units code, axis labels, second

C==== iycode,ycaption,ycaption2    caption used at END of axes by 3D plots

C==== IZCODE,zcaption,zcaption2

C====  nuserecs                no. of records following with extra info.
C====      skip nuser_records
C====  For LOQ there are 2 records here, date & run numbers
C====              ranges of lambda, angle, radius
C==== 

C====  nxc                          number of X coords ( extra control flags
C====                               on this record ignored here ) 

C====  x( )                         coord array

C====  nyc                          number of Y coords

C====  y( )                         coord array

C====  nx ny scale                  no. of x and y points
C====                               NOTE nx=nxc-1 for histogram, ixmode=1
C====                                    nx=nxc for point mode,  ixmode=0
C====                               all data are multiplied by SCALE on input

C====  iflag  fmt                   fmt is character variable with data format
C====                                e.g. (*) for free format binary

C====  if iflag=1 
C====   ((data(i,j), i=1,nx),j=1,ny)
C====  if iflag=2
C====   ((data(i,j),e(i,j), i=1,nx),j=1,ny)
C====  if iflag=3
C====   (data(i,j), i=1,nx),j=1,ny),(e(i,j), i=1,nx),j=1,ny)
C====
C====                               data is assumed stored in FORTRAN fast 
C====                               access order, bottom left to top right by 
C====                               rows 
c
C====
      SUBROUTINE GET2DA(IND,x,y,data,e,ndata,izcode,
     *    iflag,nxc,ixmode,ixcode,
     *    nyc,iymode,iycode, longtitle1,
     *    xcaption,xcaption2,ycaption,ycaption2,
     *    zcaption,zcaption2, nuserecs, longtitle2,longtitle3)
c
c==== prototype two dimensional data file for LOQ, RKH
C
C==== on initial call nxc,nyc & ndata must dimension X,Y,Data & error
C====                 IND is input stream no., 6 is assumed for text output
C==== it is assumed here that an appropriate file has been OPENed
C
C==== on return  nxc-ixmode and nyc-iymode are dimensions of data(i,j) & ndata=nx*ny
C==== "     "    izcode,ixcode & iycode contain GENIE units codes ( list below)
C====            izcode will be -1 for genuine histogram data [e.g. raw
C====            time of flight where actual counts = data(i,j)* dx(i) *dy(j) ]
C====            or 0 for ratioed, dimensionless data.
C====            IZCODE CONTAINS 999 if there has been an error
c====                            998 for end of file
C====
C====  "    "    ixmode & iymode are 0 for point axis arrays
C====                                1 for histograms i.e. x() and y() have nx+1
C====                                  or ny+1 points respectively
C
C==== NOTE IZCODE, ixcode & iycode are read in, whereas IXMODE & IYMODE are
C==== inferred from the values of NX,NXC,NY & NYC
c
c
c units codes:
c	1: channel No.		 2: T.O.F. ( usec ) ( Raw file format )
c	3: wavelength ( A )	 3: d-spacing ( A )
c	5: energy ( meV )	 6: momentum transfer ( 1/A )
c	7: Q^2 ( A^2 )		 8: tau ( uS/m )
c	9: dE ( meV )		 0: ratio ( dimensionless ) ( for izcode only)
C      -1 for histogram mode ( for izcode only)

C
      DIMENSION x(nxc), y(nyc), data(ndata), e(ndata)
      character*80 longtitle1, longtitle2, longtitle3, fmt, text
      character*40 xcaption,xcaption2,ycaption,ycaption2,zcaption,
     *    zcaption2
c
      READ(IND,1001,err=900,end=970)longtitle1
1001  FORMAT(A80)
      READ(IND,1002,err=900,end=950)ixcode,xcaption,xcaption2    
      READ(IND,1002,err=900,end=950)iycode,ycaption,ycaption2    
      READ(IND,1002,err=900,end=950)izcode,zcaption,zcaption2    
1002  FORMAT(I3,2A40)
      READ(IND,1002,err=900,end=950)nuserecs
C====      skip nuser_records this bit is doctored specifically for LOQ
      longtitle2=' '
      longtitle3=' '

      IF(nuserecs.gt.0)then
      IF(nuserecs.ge.1)READ(IND,1001,err=900,end=950)longtitle2
      IF(nuserecs.ge.2)READ(IND,1001,err=900,end=950)longtitle3
      if(nuserecs.ge.3)then
      DO I=3,nuserecs
      READ(IND,1001,err=900,end=950)
      END DO
      end if
      END IF
C
      nx=nxc
      ny=nyc
c==== on call here nx,ny are used as dimension, proper nyc & nxc are returned
      CALL GETAXIS(nxc,nx,ixcode,IND,X)
      CALL GETAXIS(nyc,ny,iycode,IND,Y)
c==== now read the proper nx & ny
      READ(IND,*,err=900,end=950)nx,ny,scale 
1004  FORMAT(I,I,F)
C
      IF(NX.EQ.NXC)THEN
C==== point mode
      IXMODE=0
      ELSE IF(NX.EQ.NXC-1)THEN
c==== histogram mode
      IXMODE=1
      ELSE
c==== error
      write(6,1011)
1011  FORMAT(1X,'WARNING: NO. OF X COORDS CONFLICT IN GET2DA ')
      IZCODE=999
      END IF
C
      IF(NY.EQ.NYC)THEN
C==== point mode
      IYMODE=0
      ELSE IF(NY.EQ.NYC-1)THEN
c==== histogram mode
      IYMODE=1
      ELSE
c==== error
      write(6,1012)
1012  FORMAT(1X,'WARNING: NO. OF Y COORDS CONFLICT IN GET2DA ')
      IZCODE=999
      END IF
c
      i = nx*ny
      IF(i.GT.ndata)THEN
      WRITE(6,1013)
1013  FORMAT(1X,'ERROR: DATA ARRAY TOO SMALL IN GET2DA')
      IZCODE=999
      i = ndata
      END IF
      ndata=i
c
      READ(IND,1015,err=900,end=950)iflag, fmt
1015  FORMAT(I3,A)
C
      IF(IFLAG.EQ.1)THEN
      READ(IND,FMT,err=900,end=950)(data(i),i=1,ndata)
      DO I=1,NDATA
      E(I)=0.0
      END DO
C
      ELSE IF (IFLAG.EQ.2)THEN
      READ(IND,FMT,err=900,end=950)(data(i),e(i),i=1,ndata)
C
      ELSE IF (IFLAG.EQ.3)THEN
      READ(IND,FMT,err=900,end=950)(data(i),i=1,ndata)
      READ(IND,FMT,err=900,end=950)(e(i),i=1,ndata)
C
      ELSE
      WRITE(6,1016)
1016  FORMAT(1X,'ERROR: INVALID IFLAG VALUE IN GET2DA')
      IZCODE=999
      RETURN
      END IF
c
c==== multiply by common scale factor
      IF(ABS(SCALE-1.0).GT.1.E-12)THEN
      DO i=1,ndata
      data(i)=data(i)*scale
      END DO
      DO i=1,ndata
      e(i)=e(i)*scale
      END DO
      END IF
C
      RETURN
C
C==== errors section
900   WRITE(6,901)
901   FORMAT(1X,'ERROR - DURING READ IN GET2DA')
      IZCODE=999
      RETURN
950   WRITE(6,951)
951   FORMAT(1X,'ERROR - END OF FILE DURING READ IN GET2DA')
      IZCODE=999
      RETURN
970   WRITE(6,971)
971   FORMAT(1X,'ERROR - NO MORE DATA SETS IN FILE (IN GET2DA)')
      IZCODE=998
      RETURN
      END
C
      SUBROUTINE GETAXIS(nxc,nx,ixcode,IND,X)
c==== note nx has not been read yet from file so is still dimension
C==== of X array.
      DIMENSION X(NX)
      READ(IND,*,end=900,err=900)nxc
1002  FORMAT(I)
      IF(nxc.GT.0)THEN
         IF(NXC.GT.NX)THEN
         WRITE(6,1005)
1005  FORMAT(1x,'WARNING: TOO MANY COORDS IN GET2DA - GETAXIS')
         NXC=NX
         END IF
c
      READ(IND,*,end=900,err=900)(x(i),i=1,nxc)
C
      END IF
      RETURN
C
900   WRITE(6,901)
901   FORMAT(1X,'ERROR IN GET2DA - GETAXIS')
      ixcode=999
      RETURN
      END
C====
C====
C===================================================================
C==== this is the BINARY version
      SUBROUTINE PUT2DB(IND,IFLAG,x,y,data,e,SCALE,ndim,
     *    izcode,nxc,ixmode,ixcode,
     *    nyc,iymode,iycode,
     *    longtitle1,
     *    xcaption,xcaption2,ycaption,ycaption2,zcaption,
     *    zcaption2,  nuserecs,longtitle2,longtitle3)
C
      DIMENSION x(nxc), y(nyc), data(ndim), e(ndim)
      character*80 longtitle1, longtitle2, longtitle3, fmt
      character*40 xcaption,xcaption2,ycaption,ycaption2,zcaption,
     *    zcaption2
      INTEGER*4 cli$present, cli$get_value
c
      WRITE(IND,err=900)longtitle1
      WRITE(IND,err=900)ixcode,xcaption,xcaption2    
      WRITE(IND,err=900)iycode,ycaption,ycaption2    
      WRITE(IND,err=900)izcode,zcaption,zcaption2    
c
c==== user may call own routine here
      WRITE(IND,err=900)MIN(2,nuserecs)
      IF(nuserecs.ge.1)WRITE(IND,err=900)longtitle2
      IF(nuserecs.ge.2)WRITE(IND,err=900)longtitle3
c
      ndata = (NXC-IXMODE)*(NYC-IYMODE)

      if(ndata.gt.ndim)goto 900
      WRITE(IND)nxc
      WRITE(IND,err=900)(x(i),i=1,nxc)
      WRITE(IND)nyc
      WRITE(IND,err=900)(y(i),i=1,nyc)
C
      IF(ABS(SCALE).LT.1.E-30)SCALE=1.0
      WRITE(IND,err=900)nxc-ixmode,nyc-iymode,scale 
c
C==== decide other formats here dependent on value of IFLAG ???
C====  assume user has ALREADY divided data by scale
C
      FMT='(*)'

      WRITE(IND,err=900)iflag,FMT
C
      IF(IFLAG.EQ.1)THEN
      WRITE(IND,err=900)(data(i),i=1,ndata)
      ELSE IF (IFLAG.EQ.2)THEN
      WRITE(IND,err=900)(data(i),e(i),i=1,ndata)
C
      ELSE IF (IFLAG.EQ.3)THEN
      WRITE(IND,err=900)(data(i),i=1,ndata)
      WRITE(IND,err=900)(e(i),i=1,ndata)
C
      ELSE
      CALL ERROR('ERROR - INVALID IFLAG VALUE IN PUT2DB binary')
      IZCODE=999
      RETURN
      END IF
c
C
      RETURN
C
C==== errors section
900   CALL ERROR('ERROR - DURING WRITE IN PUT2DB binary')
      IZCODE=999
      RETURN
      END
C===================================================================
C====
      SUBROUTINE PUT2DA(IND,IFLAG,x,y,data,e,SCALE,ndim,
     *    izcode,nxc,ixmode,ixcode,
     *    nyc,iymode,iycode,
     *    longtitle1,
     *    xcaption,xcaption2,ycaption,ycaption2,zcaption,
     *    zcaption2,  nuserecs,longtitle2,longtitle3)
c==== prototype two dimensional data file for LOQ, RKH
C
C====  nxc,nyc & dimension X,Y,  ndim dimensions Data() & error()
C====   actual number of data are ndata=(nxc-ixmode)*(nyc-iymode)
C
C====                 IND is output stream no. - it is assumed that
C====  a file has already been opened or that the default FOR00n.DAT is OK
C====    6 is assumed for text output
c==== IFLAG determines whether and where error E() is written - see below
C
C==== on call    izcode,ixcode & iycode contain GENIE units codes ( list below)
C====            izcode will be -1 for genuine histogram data [e.g. raw
C====            time of flight where actual counts = data(i,j)* dx(i) *dy(j) ]
C====            or 0 for ratioed, dimensionless data.
C====            IZCODE CONTAINS 999 if there has been an error
c====                            998 for end of file
C====
C====  "    "    ixmode & iymode are 0 for point axis arrays
C====                                1 for histograms i.e. x() and y() have nx+1
C====                                  or ny+1 points respectively
C
c units codes:
c	1: channel No.		 2: T.O.F. ( usec ) ( Raw file format )
c	3: wavelength ( A )	 3: d-spacing ( A )
c	5: energy ( meV )	 6: momentum transfer ( 1/A )
c	7: Q^2 ( A^2 )		 8: tau ( uS/m )
c	9: dE ( meV )		 0: ratio ( dimensionless ) ( for izcode only)
C      -1 for histogram mode ( for izcode only)

C
      DIMENSION x(nxc), y(nyc), data(ndim), e(ndim)
      character*80 longtitle1, longtitle2, longtitle3, fmt
      character*40 xcaption,xcaption2,ycaption,ycaption2,zcaption,
     *    zcaption2
      INTEGER*4 cli$present, cli$get_value
c
      WRITE(IND,1001,err=900)longtitle1
1001  FORMAT(A80)
      WRITE(IND,1002,err=900)ixcode,xcaption,xcaption2    
      WRITE(IND,1002,err=900)iycode,ycaption,ycaption2    
      WRITE(IND,1002,err=900)izcode,zcaption,zcaption2    
1002  FORMAT(I3,1X,2A40)
c
c==== user may call own routine here
      WRITE(IND,1002,err=900)nuserecs
      IF(nuserecs.ge.1)WRITE(IND,1001,err=900)longtitle2
      IF(nuserecs.ge.2)WRITE(IND,1001,err=900)longtitle3
c
      ndata = (NXC-IXMODE)*(NYC-IYMODE)

      if(ndata.gt.ndim)goto 900
      WRITE(IND,1007)nxc
      WRITE(IND,1008,err=900)(x(i),i=1,nxc)
      WRITE(IND,1007)nyc
      WRITE(IND,1008,err=900)(y(i),i=1,nyc)
1007  FORMAT(I5)
1008  FORMAT(1P8E14.6)
C
      IF(ABS(SCALE).LT.1.E-30)SCALE=1.0
      WRITE(IND,1004,err=900)nxc-ixmode,nyc-iymode,scale 
1004  FORMAT(2I6,1PE20.12)
c
C==== decide other formats here dependent on value of IFLAG ???
C====  assume user has ALREADY divided data by scale
C
      FMT='(1P8E12.4)'
      len= 10

c==== look for /FORMAT=(*) etc.   RKH 22/5/92
      if (cli$present ('FORMAT'))status = cli$get_value ('FORMAT', FMT, len)

      WRITE(IND,1015,err=900)iflag,FMT(1:len)
1015  FORMAT(I3,a)
C
      IF(IFLAG.EQ.1)THEN
      WRITE(IND,FMT(1:len),err=900)(data(i),i=1,ndata)
      ELSE IF (IFLAG.EQ.2)THEN
      WRITE(IND,FMT(1:len),err=900)(data(i),e(i),i=1,ndata)
C
      ELSE IF (IFLAG.EQ.3)THEN
      WRITE(IND,FMT(1:len),err=900)(data(i),i=1,ndata)
      WRITE(IND,FMT(1:len),err=900)(e(i),i=1,ndata)
C
      ELSE
      CALL ERROR('ERROR - INVALID IFLAG VALUE IN PUT2DA')
      IZCODE=999
      RETURN
      END IF
c
C
      RETURN
C
C==== errors section
900   CALL ERROR('ERROR - DURING WRITE IN PUT2DA')
      IZCODE=999
      RETURN
      END
c====        this is the BINARY version
      SUBROUTINE GET2DB(IND,x,y,data,e,ndata,izcode,
     *    iflag,nxc,ixmode,ixcode,
     *    nyc,iymode,iycode, longtitle1,
     *    xcaption,xcaption2,ycaption,ycaption2,
     *    zcaption,zcaption2, nuserecs, longtitle2,longtitle3)
c
c==== prototype two dimensional data file for LOQ, RKH
C
C==== on initial call nxc,nyc & ndata must dimension X,Y,Data & error
C====                 IND is input stream no., 6 is assumed for text output
C==== it is assumed here that an appropriate file has been OPENed
C
C==== on return  nxc-ixmode and nyc-iymode are dimensions of data(i,j) & ndata=nx*ny
C==== "     "    izcode,ixcode & iycode contain GENIE units codes ( list below)
C====            izcode will be -1 for genuine histogram data [e.g. raw
C====            time of flight where actual counts = data(i,j)* dx(i) *dy(j) ]
C====            or 0 for ratioed, dimensionless data.
C====            IZCODE CONTAINS 999 if there has been an error, 
c====                            998 for end of file
C====
C====  "    "    ixmode & iymode are 0 for point axis arrays
C====                                1 for histograms i.e. x() and y() have nx+1
C====                                  or ny+1 points respectively
C
C==== NOTE IZCODE, ixcode & iycode are read in, whereas IXMODE & IYMODE are
C==== inferred from the values of NX,NXC,NY & NYC
c
c
c units codes:
c	1: channel No.		 2: T.O.F. ( usec ) ( Raw file format )
c	3: wavelength ( A )	 3: d-spacing ( A )
c	5: energy ( meV )	 6: momentum transfer ( 1/A )
c	7: Q^2 ( A^2 )		 8: tau ( uS/m )
c	9: dE ( meV )		 0: ratio ( dimensionless ) ( for izcode only)
C      -1 for histogram mode ( for izcode only)

C
      DIMENSION x(nxc), y(nyc), data(ndata), e(ndata)
      character*80 longtitle1, longtitle2, longtitle3, fmt, text
      character*40 xcaption,xcaption2,ycaption,ycaption2,zcaption,
     *    zcaption2
c
      READ(IND,err=900,end=970)longtitle1
      READ(IND,err=900,end=950)ixcode,xcaption,xcaption2    
      READ(IND,err=900,end=950)iycode,ycaption,ycaption2    
      READ(IND,err=900,end=950)izcode,zcaption,zcaption2    
      READ(IND,err=900,end=950)nuserecs
C====      skip nuser_records this bit is doctored specifically for LOQ
      longtitle2=' '
      longtitle3=' '

      IF(nuserecs.gt.0)then
      IF(nuserecs.ge.1)READ(IND,err=900,end=950)longtitle2
      IF(nuserecs.ge.2)READ(IND,err=900,end=950)longtitle3
      if(nuserecs.ge.3)then
      DO I=3,nuserecs
      READ(IND,err=900,end=950)
      END DO
      end if
      END IF
C
      nx=nxc
      ny=nyc
c==== on call here nx,ny are used as dimension, proper nyc & nxc are returned
      CALL GETAXISB(nxc,nx,ixcode,IND,X)
      CALL GETAXISB(nyc,ny,iycode,IND,Y)
c==== now read the proper nx & ny
      READ(IND,err=900,end=950)nx,ny,scale 
C
      IF(NX.EQ.NXC)THEN
C==== point mode
      IXMODE=0
      ELSE IF(NX.EQ.NXC-1)THEN
c==== histogram mode
      IXMODE=1
      ELSE
c==== error
      write(6,1011)
1011  FORMAT(1X,'WARNING: NO. OF X COORDS CONFLICT IN GET2DB ')
      IZCODE=999
      END IF
C
      IF(NY.EQ.NYC)THEN
C==== point mode
      IYMODE=0
      ELSE IF(NY.EQ.NYC-1)THEN
c==== histogram mode
      IYMODE=1
      ELSE
c==== error
      write(6,1012)
1012  FORMAT(1X,'WARNING: NO. OF Y COORDS CONFLICT IN GET2DB ')
      IZCODE=999
      END IF
c
      i = nx*ny
      IF(i.GT.ndata)THEN
      WRITE(6,1013)
1013  FORMAT(1X,'ERROR: DATA ARRAY TOO SMALL IN GET2DB')
      IZCODE=999
      i = ndata
      END IF
      ndata=i
c
      READ(IND,err=900,end=950)iflag, fmt
C
      IF(IFLAG.EQ.1)THEN
      READ(IND,err=900,end=950)(data(i),i=1,ndata)
      DO I=1,NDATA
      E(I)=0.0
      END DO
C
      ELSE IF (IFLAG.EQ.2)THEN
      READ(IND,err=900,end=950)(data(i),e(i),i=1,ndata)
C
      ELSE IF (IFLAG.EQ.3)THEN
      READ(IND,err=900,end=950)(data(i),i=1,ndata)
      READ(IND,err=900,end=950)(e(i),i=1,ndata)
C
      ELSE
      WRITE(6,1016)
1016  FORMAT(1X,'ERROR: INVALID IFLAG VALUE IN GET2DB binary')
      IZCODE=999
      RETURN
      END IF
c
c==== multiply by common scale factor
      IF(ABS(SCALE-1.0).GT.1.E-12)THEN
      DO i=1,ndata
      data(i)=data(i)*scale
      END DO
      DO i=1,ndata
      e(i)=e(i)*scale
      END DO
      END IF
C
      RETURN
C
C==== errors section
900   WRITE(6,901)
901   FORMAT(1X,'ERROR - DURING READ IN GET2DB binary')
      IZCODE=999
      RETURN
950   WRITE(6,951)
951   FORMAT(1X,'ERROR - END OF FILE DURING READ IN GET2DB')
      IZCODE=999
      RETURN
970   WRITE(6,971)
971   FORMAT(1X,'ERROR - NO MORE DATA SETS IN FILE (IN GET2DB)')
      IZCODE=998
      RETURN
      END
C
C====           bianry version
      SUBROUTINE GETAXISB(nxc,nx,ixcode,IND,X)
c==== note nx has not been read yet from file so is still dimension
C==== of X array.
      DIMENSION X(NX)
      READ(IND,end=900,err=900)nxc
      IF(nxc.GT.0)THEN
         IF(NXC.GT.NX)THEN
         WRITE(6,1005)
1005  FORMAT(1x,'WARNING: TOO MANY COORDS IN GET2DB - GETAXIS')
         NXC=NX
         END IF
c
      READ(IND,end=900,err=900)(x(i),i=1,nxc)
C
      END IF
      RETURN
C
900   WRITE(6,901)
901   FORMAT(1X,'ERROR IN GET2DB - GETAXISB binary')
      ixcode=999
      RETURN
      END
c