program qq22
c
c     read MSTG data on a RS6000/IBM 
c     Using fortran version of GBYTES
c
      IMPLICIT INTEGER(A-E,G-Z)
C
      PARAMETER(MAX=10,FMISS=-9999.,INDEXCK=6,BPR=384,ID=0
     +,BPW=32,dimbuf=(bpr/bpw+6),dimpk=(BPR-1)/BPW+1,dimun=38
     +,GROUP=3,iunit=10)
C
      COMMON /MSTG2/FUNITS(38),FBASE(38),BITS(38),OFFSET(38)
C
      DIMENSION BUF(dimbuf),PK(dimpk),UN(dimun),FTRUE(dimun)
C
C-----2 DIMENSIONAL FTRUE
      DIMENSION FTRUE2(4,8)
      EQUIVALENCE (FTRUE(7),FTRUE2)
C
      DATA BUF/dimbuf*0/
C
c     open UNIX data file - direct access mode
c     each direct access will get one MSTG record

      open(unit=iunit,file='mstg3a_1854',
     *     form='unformatted', access='direct',recl=48)
c
      write(*,1)
1     format('1  ****** beginning output for qq22  *****')
c
c     INITAL initializes masks required for gbytes and sbytes
      CALL INITAL
      CALL DATA(GROUP)

 100   CALL GETRPT(iunit,FMISS,FUNITS,FBASE,BITS,OFFSET,INDEXCK,ID
     +,BPR,BPW,BUF,dimbuf,PK,dimpk,UN,dimun,FTRUE,JEOF)
      IF(JEOF.NE.0)GOTO 900
C
c     print the first max records

      if(buf(2).lt.max.or.mod(buf(2),2000).eq.0)then
         print*,' record number : ',buf(2)
         CALL WRMSTG2(FTRUE,GROUP)
      endif

      go to 100

C
 900   PRINT *,' REPORTS ',BUF(2),', EOF ',JEOF
      END
C=======================================================================
      SUBROUTINE WRMSTG2(FTRUE,GROUP)
      IMPLICIT INTEGER(A-E,G-Z)
      DIMENSION FTRUE(38)
      GOTO (1,2,3,4,5,6,7,8), GROUP
  1   ASSIGN 101 TO IFMT
      GOTO 100
  2   ASSIGN 102 TO IFMT
      GOTO 100
  3   ASSIGN 103 TO IFMT
      GOTO 100
  4   ASSIGN 104 TO IFMT
      GOTO 100
  5   ASSIGN 105 TO IFMT
      GOTO 100
  6   ASSIGN 106 TO IFMT
      GOTO 100
  7   ASSIGN 107 TO IFMT
      GOTO 100
  8   ASSIGN 108 TO IFMT
 100   PRINT IFMT,(FTRUE(I),I=1,6)
     +,((FTRUE(6+(J-1)*4+I),J=1,8),I=1,4)
 101   FORMAT(/' YEAR ',F5.0,' MONTH ',F3.0,' BOX2 ',F6.0
     +,' BOX10 ',F4.0,' GROUP ',F3.0,' CHECKSUM ',F6.0/
     +8X,'3',7X,'M',7X,'N',7X,'E',7X,'D',7X,'H',7X,'X',7X,'Y'/
     +1X,'S',2F8.2,F8.0,F8.2,2F8.0,2F8.1/
     +1X,'A',2F8.2,F8.0,F8.2,2F8.0,2F8.1/
     +1X,'P',2F8.2,F8.0,F8.2,2F8.0,2F8.1/
     +1X,'Q',2F8.2,F8.0,F8.2,2F8.0,2F8.1)
 102   FORMAT(/' YEAR ',F5.0,' MONTH ',F3.0,' BOX2 ',F6.0
     +,' BOX10 ',F4.0,' GROUP ',F3.0,' CHECKSUM ',F6.0/
     +8X,'3',7X,'M',7X,'N',7X,'E',7X,'D',7X,'H',7X,'X',7X,'Y'/
     +1X,'W',2F8.2,F8.0,F8.2,2F8.0,2F8.1/
     +1X,'U',2F8.2,F8.0,F8.2,2F8.0,2F8.1/
     +1X,'V',2F8.2,F8.0,F8.2,2F8.0,2F8.1/
     +1X,'C',2F8.1,F8.0,F8.1,2F8.0,2F8.1)
 103   FORMAT(/' YEAR ',F5.0,' MONTH ',F3.0,' BOX2 ',F6.0
     +,' BOX10 ',F4.0,' GROUP ',F3.0,' CHECKSUM ',F6.0/
     +9X,7X,'3',7X,'M',7X,'N',7X,'E',7X,'D',7X,'H',7X,'X',7X,'Y'/
     +1X,'S       ',2F8.2,F8.0,F8.2,F8.0,3F8.1/
     +1X,'A       ',2F8.2,F8.0,F8.2,F8.0,3F8.1/
     +1X,'Q       ',2F8.2,F8.0,F8.2,F8.0,3F8.1/
     +1X,'R       ',2F8.1,F8.0,F8.1,F8.0,3F8.1)
 104   FORMAT(/' YEAR ',F5.0,' MONTH ',F3.0,' BOX2 ',F6.0
     +,' BOX10 ',F4.0,' GROUP ',F3.0,' CHECKSUM ',F6.0/
     +9X,7X,'3',7X,'M',7X,'N',7X,'E',7X,'D',7X,'H',7X,'X',7X,'Y'/
     +1X,'W       ',2F8.2,F8.0,F8.2,F8.0,3F8.1/
     +1X,'U       ',2F8.2,F8.0,F8.2,F8.0,3F8.1/
     +1X,'V       ',2F8.2,F8.0,F8.2,F8.0,3F8.1/
     +1X,'P       ',2F8.2,F8.0,F8.2,F8.0,3F8.1)
 105   FORMAT(/' YEAR ',F5.0,' MONTH ',F3.0,' BOX2 ',F6.0
     +,' BOX10 ',F4.0,' GROUP ',F3.0,' CHECKSUM ',F6.0/
     +9X,7X,'3',7X,'M',7X,'N',7X,'E',7X,'D',7X,'H',7X,'X',7X,'Y'/
     +1X,'C       ',2F8.1,F8.0,F8.1,F8.0,3F8.1/
     +1X,'R       ',2F8.1,F8.0,F8.1,F8.0,3F8.1/
     +1X,'W*U     ',2F8.1,F8.0,F8.1,F8.0,3F8.1/
     +1X,'W*V     ',2F8.1,F8.0,F8.1,F8.0,3F8.1)
 106   FORMAT(/' YEAR ',F5.0,' MONTH ',F3.0,' BOX2 ',F6.0
     +,' BOX10 ',F4.0,' GROUP ',F3.0,' CHECKSUM ',F6.0/
     +9X,7X,'3',7X,'M',7X,'N',7X,'E',7X,'D',7X,'H',7X,'X',7X,'Y'/
     +1X,'S-A     ',2F8.2,F8.0,F8.2,F8.0,3F8.1/
     +1X,'(S-A)*W ',2F8.1,F8.0,F8.1,F8.0,3F8.1/
     +1X,'QS-Q    ',2F8.2,F8.0,F8.2,F8.0,3F8.1/
     +1X,'(QS-Q)*W',2F8.1,F8.0,F8.1,F8.0,3F8.1)
 107   FORMAT(/' YEAR ',F5.0,' MONTH ',F3.0,' BOX2 ',F6.0
     +,' BOX10 ',F4.0,' GROUP ',F3.0,' CHECKSUM ',F6.0/
     +9X,7X,'3',7X,'M',7X,'N',7X,'E',7X,'D',7X,'H',7X,'X',7X,'Y'/
     +1X,'U*A     ',2F8.1,F8.0,F8.1,F8.0,3F8.1/
     +1X,'V*A     ',2F8.1,F8.0,F8.1,F8.0,3F8.1/
     +1X,'U*Q     ',2F8.1,F8.0,F8.1,F8.0,3F8.1/
     +1X,'V*Q     ',2F8.1,F8.0,F8.1,F8.0,3F8.1)
 108   FORMAT(/' YEAR ',F5.0,' MONTH ',F3.0,' BOX2 ',F6.0
     +,' BOX10 ',F4.0,' GROUP ',F3.0,' CHECKSUM ',F6.0/
     +9X,7X,'3',7X,'M',7X,'N',7X,'E',7X,'D',7X,'H',7X,'X',7X,'Y'/
     +1X,'S       ',2F8.2,F8.0,F8.2,F8.0,3F8.1/
     +1X,'A       ',2F8.2,F8.0,F8.2,F8.0,3F8.1/
     +1X,'W*U     ',2F8.1,F8.0,F8.1,F8.0,3F8.1/
     +1X,'W*V     ',2F8.1,F8.0,F8.1,F8.0,3F8.1)

      RETURN
      END
C=======================================================================
      SUBROUTINE DATA(GROUP)
      IMPLICIT INTEGER(A-E,G-Z)
C
      COMMON /MSTG2/FUNITS(38),FBASE(38),BITS(38),OFFSET(38)
C
      DATA FUNITS/1.,1.,1.,1.,1.,1.
     + ,.01,.01,.1,.1
     + ,.01,.01,.1,.1
     + ,1.,1.,1.,1.
     + ,.01,.01,.1,.1
     + ,2.,2.,2.,2.
     + ,.1,.1,.1,.1
     + ,.2,.2,.2,.2
     + ,.2,.2,.2,.2/
C
      DATA FBASE/1799.,0.,0.,0.,0.,0.
     + ,-501.,-8801.,-30001.,-30001.
     + ,-501.,-8801.,-30001.,-30001.
     + ,0.,0.,0.,0.
     + ,-1.,-1.,-1.,-1.
     + ,0.,0.,0.,0.
     + ,-1.,-1.,-1.,-1.
     + ,-.5,-.5,-.5,-.5
     + ,-.5,-.5,-.5,-.5/
C
      DATA BITS/8,4,14,10,4,8
     + ,16,16,16,16
     + ,16,16,16,16
     + ,16,16,16,16
     + ,16,16,16,16
     + ,4,4,4,4
     + ,4,4,4,4
     + ,4,4,4,4
     + ,4,4,4,4/
C
      DATA OFFSET/16,24,28,42,52,56
     + ,64,80,96,112
     + ,128,144,160,176
     + ,192,208,224,240
     + ,256,272,288,304
     + ,320,324,328,332
     + ,336,340,344,348
     + ,352,356,360,364
     + ,368,372,376,380/
C
      GOTO (1,2,3,4,5,6,7,8), GROUP
C
  1   FUNITS(7)=.01
      FUNITS(8)=.01
      FUNITS(9)=.01
      FUNITS(10)=.01
      FUNITS(11)=.01
      FUNITS(12)=.01
      FUNITS(13)=.01
      FUNITS(14)=.01
      FUNITS(19)=.01
      FUNITS(20)=.01
      FUNITS(21)=.01
      FUNITS(22)=.01
      FUNITS(27)=2.
      FUNITS(28)=2.
      FUNITS(29)=2.
      FUNITS(30)=2.
      FBASE(7)=-501.
      FBASE(8)=-8801.
      FBASE(9)=86999.
      FBASE(10)=-1.
      FBASE(11)=-501.
      FBASE(12)=-8801.
      FBASE(13)=86999.
      FBASE(14)=-1.
      FBASE(27)=-.5
      FBASE(28)=-.5
      FBASE(29)=-.5
      FBASE(30)=-.5
      RETURN
C
  2   FUNITS(7)=.01
      FUNITS(8)=.01
      FUNITS(9)=.01
      FUNITS(10)=.1
      FUNITS(11)=.01
      FUNITS(12)=.01
      FUNITS(13)=.01
      FUNITS(14)=.1
      FUNITS(19)=.01
      FUNITS(20)=.01
      FUNITS(21)=.01
      FUNITS(22)=.1
      FUNITS(27)=2.
      FUNITS(28)=2.
      FUNITS(29)=2.
      FUNITS(30)=2.
      FBASE(7)=-1.
      FBASE(8)=-10221.
      FBASE(9)=-10221.
      FBASE(10)=-1.
      FBASE(11)=-1.
      FBASE(12)=-10221.
      FBASE(13)=-10221.
      FBASE(14)=-1.
      FBASE(27)=-.5
      FBASE(28)=-.5
      FBASE(29)=-.5
      FBASE(30)=-.5
      RETURN
C
  3   FUNITS(7)=.01
      FUNITS(8)=.01
      FUNITS(9)=.01
      FUNITS(10)=.1
      FUNITS(11)=.01
      FUNITS(12)=.01
      FUNITS(13)=.01
      FUNITS(14)=.1
      FUNITS(19)=.01
      FUNITS(20)=.01
      FUNITS(21)=.01
      FUNITS(22)=.1
      FBASE(7)=-501.
      FBASE(8)=-8801.
      FBASE(9)=-1.
      FBASE(10)=-1.
      FBASE(11)=-501.
      FBASE(12)=-8801.
      FBASE(13)=-1.
      FBASE(14)=-1.
      RETURN
C
  4   FUNITS(7)=.01
      FUNITS(8)=.01
      FUNITS(9)=.01
      FUNITS(10)=.01
      FUNITS(11)=.01
      FUNITS(12)=.01
      FUNITS(13)=.01
      FUNITS(14)=.01
      FUNITS(19)=.01
      FUNITS(20)=.01
      FUNITS(21)=.01
      FUNITS(22)=.01
      FBASE(7)=-1.
      FBASE(8)=-10221.
      FBASE(9)=-10221.
      FBASE(10)=86999.
      FBASE(11)=-1.
      FBASE(12)=-10221.
      FBASE(13)=-10221.
      FBASE(14)=86999.
      RETURN
C
  5   FUNITS(7)=.1
      FUNITS(8)=.1
      FUNITS(9)=.1
      FUNITS(10)=.1
      FUNITS(11)=.1
      FUNITS(12)=.1
      FUNITS(13)=.1
      FUNITS(14)=.1
      FUNITS(19)=.1
      FUNITS(20)=.1
      FUNITS(21)=.1
      FUNITS(22)=.1
      FBASE(7)=-1.
      FBASE(8)=-1.
      FBASE(9)=-30001.
      FBASE(10)=-30001.
      FBASE(11)=-1.
      FBASE(12)=-1.
      FBASE(13)=-30001.
      FBASE(14)=-30001.
      RETURN
C
  6   FUNITS(7)=.01
      FUNITS(8)=.1
      FUNITS(9)=.01
      FUNITS(10)=.1
      FUNITS(11)=.01
      FUNITS(12)=.1
      FUNITS(13)=.01
      FUNITS(14)=.1
      FUNITS(19)=.01
      FUNITS(20)=.1
      FUNITS(21)=.01
      FUNITS(22)=.1
      FBASE(7)=-6301.
      FBASE(8)=-10001.
      FBASE(9)=-4001.
      FBASE(10)=-10001.
      FBASE(11)=-6301.
      FBASE(12)=-10001.
      FBASE(13)=-4001.
      FBASE(14)=-10001.
      RETURN
C
  7   FUNITS(7)=.1
      FUNITS(8)=.1
      FUNITS(9)=.1
      FUNITS(10)=.1
      FUNITS(11)=.1
      FUNITS(12)=.1
      FUNITS(13)=.1
      FUNITS(14)=.1
      FUNITS(19)=.1
      FUNITS(20)=.1
      FUNITS(21)=.1
      FUNITS(22)=.1
      FBASE(7)=-20001.
      FBASE(8)=-20001.
      FBASE(9)=-10001.
      FBASE(10)=-10001.
      FBASE(11)=-20001.
      FBASE(12)=-20001.
      FBASE(13)=-10001.
      FBASE(14)=-10001.
      RETURN
C
  8   END
C=======================================================================
      SUBROUTINE GETRPT(iunit,FMISS,FUNITS,FBASE,BITS,OFFSET,INDEXCK,ID
     +,BPR,BPW,BUF,dimbuf,PK,dimpk,UN,dimun,FTRUE,JEOF)
C
C-----RETURN FLOATING POINT VALUES IN FTRUE
C
C     INPUT
C         iunit - RCDIN UNIT
C         FMISS - MISSING VALUE
C         FUNITS(dimun) - UNITS FOR UNCODING
C         FBASE(dimun) - BASE FOR UNCODING
C         BITS(dimun) - BITS FOR UNPACKING
C         OFFSET(dimun) - OFFSET FOR UNPACKING
C         INDEXCK - UN(INDEXCK) = CHECKSUM
C         ID - GROUP NUMBER FOR IDENTIFICATION CHECKSUM
C         BPR - BITS PER REPORT
C         BPW - BITS PER WORD
C     OUTPUT
C         BUF(dimbuf) - RCDIN BUFFER
C         PK(dimpk) - PACKED REPORT
C         UN(dimun) - UNPACKED REPORT
C         FTRUE(dimun) - TRUE VALUES
C         JEOF = end of file indicator, 0=no, -1=yes
C
      IMPLICIT INTEGER(A-E,G-Z)
      DIMENSION FUNITS(dimun),FBASE(dimun),BITS(dimun),OFFSET(dimun)
     +,BUF(dimbuf),PK(dimpk),UN(dimun),FTRUE(dimun)
C
       CALL RCDIN(iunit,BUF,dimbuf,PK,dimpk,BPR,BPW,JEOF)
      if(jeof.eq.-1)return
C
C     unpack with gbyte and rescale values

      CK=ID
      DO 230 I=1,dimun
      CALL GBYTE(PK(OFFSET(I)/BPW+1),UN(I),MOD(OFFSET(I),BPW),BITS(I))
      IF(I.EQ.INDEXCK)GOTO 210
      IF(UN(I).EQ.0)GOTO 220
      FTRUE(I)=(UN(I)+FBASE(I))*FUNITS(I)
      CK=CK+UN(I)
      GOTO 230
 210   FTRUE(INDEXCK)=UN(INDEXCK)
      GOTO 230
 220   FTRUE(I)=FMISS
 230   CONTINUE
c
c      checksum verification
c
       IF(MOD(CK,2**BITS(INDEXCK)-1).EQ.UN(INDEXCK))RETURN
C
C-----ERROR
       PRINT *,' SUBROUTINE GETRPT  --  CHECKSUM ERROR, iunit = ',iunit
     +,', REPORT = ',BUF(2)
       PRINT *,' FTRUE = ',FTRUE
 800   STOP
C
 900   END
C=======================================================================
      SUBROUTINE RCDIN(iunit,BUF,dimbuf,RCD,dimrcd,BPR,BPW,JEOF)
C
C-----RETURN ONE LOGICAL RECORD IN RCD
C
C     INPUT
C         iunit - BUFFER IN UNIT
C         BPR - BITS PER RECORD
C         BPW - BITS PER WORD
C     OUTPUT
C         BUF(dimbuf) - PHYSICAL RECORD
C         RCD(dimrcd) - LOGICAL RECORD
C         JEOF = end of file indicator, 0=no, 1=yes
C
C     BUF(1) = gbytes offset, always 0 in this case
C     BUF(2) = record count
C     BUF(3) =
C     BUF(4) =
C     BUF(5) = record LENGTH IN BITS
C     BUF(6) =
C
      IMPLICIT INTEGER(A-E,G-Z)
      DIMENSION BUF(dimbuf),RCD(dimrcd)
      data irdata/1/
C
c     read from UNIX data file

10    read(iunit,rec=irdata,err=900)(buf(i),i=7,18)
      irdata=irdata + 1

c
c     write the first part of the input record in hexidecimal
c           just to check the bytes
c
c       write(6,202)(buf(i),i=7,18)
c202   format(5z12)
c
      jeof=0
100   BUF(1)=0
      BUF(5)=384
      BUF(3)=BUF(3)+1
C
C-----GBYTE
      CALL GBYTES(BUF(6+1),RCD,0,BPW,0,dimrcd)
      BUF(2)=BUF(2)+1
      RETURN
c      go to 10

c     end of file has been detected set jeof to -1

900   continue
      jeof = -1

      return
      END
      SUBROUTINE GBYTES(IN,IOUT,ISKIP,NBYTE,NSKIP,N)
c
c     This fortran implementation of GBYTES and SBYTES was originally
c     derived by Robert Gammil, NCAR, 1992
c     Trudy Wohlleben, AES, Canada made some necessary adjusts to implement
c           code on RS6000 machines, October 1994
c     Keys: 1) call INITAL in the main program, and 2) define OR, AND, and NOT
c           functions to be consistant with non-RS6000 systems.
c
      COMMON/MACHIN/NBITSW,NBITSC,MASK0,MASKS(64)
      DIMENSION IN(1),IOUT(1)
      INTEGER RGHTSH,OR,AND
      LEFTSH(M,N)=ISHFT(M,N)
      RGHTSH(M,N)=ISHFT(M,-N)
      OR(M,N)=M.OR.N
      AND(M,N)=M.AND.N

      ICON=NBITSW-NBYTE
      IF(ICON.LT.0) RETURN
      MASK=MASKS(NBYTE)
      INDEX=ISKIP/NBITSW
      II=MOD(ISKIP,NBITSW)
      ISTEP=NBYTE+NSKIP
      IWORDS=ISTEP/NBITSW
      IBITS=MOD(ISTEP,NBITSW)

      DO 6 I=1,N
      MOVER=ICON-II
      IF(MOVER) 2,3,4
    2 MOVEL=-MOVER
      MOVER=NBITSW-MOVEL
      NP1=LEFTSH(IN(INDEX+1),MOVEL)
      NP2=RGHTSH(IN(INDEX+2),MOVER)
      IOUT(I)=AND(OR(NP1,NP2),MASK)
      GO TO 5
    3 IOUT(I)=AND(IN(INDEX+1),MASK)
      GO TO 5
    4 IOUT(I)=AND(RGHTSH(IN(INDEX+1),MOVER),MASK)
    5 II=II+IBITS
      INDEX=INDEX+IWORDS
      IF(II.LT.NBITSW) GO TO 6
      II=II-NBITSW
      INDEX=INDEX+1
    6 CONTINUE
      RETURN
      END
      SUBROUTINE GBYTE (IN,IOUT,ISKIP,NBYTE)
      CALL GBYTES (IN,IOUT,ISKIP,NBYTE,0,1)
      RETURN
      END
      SUBROUTINE INITAL
      COMMON/MACHIN/NBITSW,NBITSC,MASK0,MASKS(64)
      INTEGER OR
      LEFTSH(M,N)=ISHFT(M,N)
      OR(M,N)=M.OR.N
      NBITSW=32
      MASK0=0
      MASKS(1)=1
      DO 1 I=2,NBITSW
    1 MASKS(I)=OR(LEFTSH(MASKS(I-1),1),1)
      RETURN
      END
      SUBROUTINE SBYTES(IOUT,IN,ISKIP,NBYTE,NSKIP,N)
      COMMON/MACHIN/NBITSW,NBITSC,MASK0,MASKS(64)
      DIMENSION IN(1),IOUT(1)
      INTEGER RGHTSH,OR,AND
      OR(M,N)=M.OR.N
      AND(M,N)=M.AND.N
      NOT(M)=.NOT.M
      LEFTSH(M,N)=ISHFT(M,N)
      RGHTSH(M,N)=ISHFT(M,-N)

      ICON=NBITSW-NBYTE
      IF(ICON.LT.0) RETURN
      MASK=MASKS(NBYTE)
      INDEX=ISKIP/NBITSW
      II=MOD(ISKIP,NBITSW)
      ISTEP=NBYTE+NSKIP
      IWORDS=ISTEP/NBITSW
      IBITS=MOD(ISTEP,NBITSW)

      DO 6 I=1,N
      J=AND(MASK,IN(I))
      MOVEL=ICON-II
      IF(MOVEL) 2,3,4
2     MSK=MASKS(NBYTE+MOVEL)
      IOUT(INDEX+1)=OR(AND(NOT(MSK),IOUT(INDEX+1)),RGHTSH(J,-MOVEL))
      ITEMP=AND(MASKS(NBITSW+MOVEL),IOUT(INDEX+2))
      IOUT(INDEX+2)=OR(ITEMP,LEFTSH(J,NBITSW+MOVEL))
      GO TO 5
3     IOUT(INDEX+1)=OR(AND(NOT(MASK),IOUT(INDEX+1)),J)
      GO TO 5
4     MSK=LEFTSH(MASK,MOVEL)
      IOUT(INDEX+1)=OR(AND(NOT(MSK),IOUT(INDEX+1)),LEFTSH(J,MOVEL))
5     II=II+IBITS
      INDEX=INDEX+IWORDS
      IF(II.LT.NBITSW) GO TO 6
      II=II-NBITSW
      INDEX=INDEX+1
6     CONTINUE

      RETURN
      END
      SUBROUTINE SBYTE (IOUT,IN,ISKIP,NBYTE)
      CALL SBYTES (IOUT,IN,ISKIP,NBYTE,0,1)
      RETURN
      END