CMAX: [Version 1.0 (6/22/93)]
CMAX: ----------------------------------------------------
CMAX: Transformation of RDMIX from rdmix2.f
CMAX: ----------------------------------------------------
c------------------------------------------------------------------c
c                                                                  c
c  Program Name: RDMIX2                                            c
c  Date:         December 10, 1993                                 c
c  Programmer:   Mike Seablom                                      c
c                                                                  c
c  Program RDMIX reads in relative humitdity and mixing ratio data c
c  from a non-repack dataset (flexible format).                    c
c                                                                  c
c------------------------------------------------------------------c
      SUBROUTINE rdmix(nymd,nhms,data,bq,nobs,pmand,psfc,im,jnp,mlev,
     &    mmix)
      
c
c     Maximal observation points
      INTEGER MaxObservPoints
c
      PARAMETER (MaxObservPoints = 150000)
c
c
      REAL data (MaxObservPoints * 10), quality_mixtime_displacement, ps
     &fc (im, jnp), pmand (mlev)
      LOGICAL end_of_file
      INTEGER error_code, undergnd_count, bq (MaxObservPoints)
      CHARACTER * 8 dataname, obsname
      CHARACTER * 60 message
!hpf$ distribute bq(block)
!hpf$ distribute data(block)
!hpf$ distribute pmand(block)
!hpf$ distribute psfc(block,block)
      DATA fill / 1.e15 /
      COMMON / dataptrs / ikx, ilat, ilon, ipmix, itype, imix, ifcer
C  Initialize data quality flags.
      INTEGER Ti1
      message(1:60) = ' '
      bq = 1
C  Define highest pressure level for mixing ratio analysis.
      lmix = mlev - mmix + 1
C  Initialize the starting addresses for array "DATA".
      ikx = 1
      ilat = 1 * maxobservpoints + 1
      ilon = 2 * maxobservpoints + 1
      ipmix = 3 * maxobservpoints + 1
      itype = 4 * maxobservpoints + 1
      imix = 5 * maxobservpoints + 1
      irh = 6 * maxobservpoints + 1
C  Check for date and time match and dataset positioning.
      REWIND( 78 ) 
      READ( 78,* ) date,time
      IF (date .NE. float(nymd) .OR. time .NE. float(nhms)) THEN
        WRITE( 6,'(/," Analysis time: ",2i9,/," Time on dataset: ",     
     &    2i9)' ) nymd,nhms,date,time
        message = ' RDMIX2: Dataset positioned at wrong ' 
     &      // 'date/time.'
        CALL flash(5,'ABORT in RDMIX2. Dataset positioned to wrong date/
     &time.     ',1)
      ENDIF
      nobs = 0
      end_of_file = .FALSE.
      i = 1
 7000 IF ((.NOT.end_of_file)) THEN
        READ( 78,*,IOSTAT=error_code ) time_displacement,
     &      data(ikx + i - 1),data(ilat + i - 1),data(ilon + i - 1),
     &      data(ipmix + i - 1),data(irh + i - 1),data(imix + i - 1),
     &      quality_mix
        IF (error_code .NE. 0) THEN
          end_of_file = .TRUE.
        ELSE
c  Store only reports passing the proximity check.
          IF (data(ikx + i - 1) .LT. 1.E7 .AND. quality_mix .LT. 1.99
     &          ) THEN
c  Interplolate model surface pressure to observation location.
            CALL terpv_v1(data,ilat + i - 1,maxobservpoints * 10,data,
     &          ilon + i - 1,maxobservpoints * 10,psobs,1,psfc)
c Trap for any mixing ratio observations below the interpolated surface
c pressure.  This check is quite crude and will sometimes remove and
c sometimes keep data near the surface in areas with rapidly varying
c topography.  Data well below the suface will be caught, however.
            IF (data(ipmix + i - 1) .LT. psobs .AND. data(ipmix + i - 1)
     &            .GE.pmand(lmix)) THEN
              i = i + 1
              nobs = nobs + 1
              IF (abs(quality_mix - 1.) .LT. .5) THEN
                bq(i - 1) = 0
              ENDIF
            ELSE
              undergnd_count = undergnd_count + 1
            ENDIF
          ENDIF
        ENDIF
        GOTO 7000
      ENDIF
 7001 CONTINUE
      WRITE( 6,'(/," Total number of mixing observations read:",        
     &i9)' ) nobs
      WRITE( 6,'(" Number of obs. below model surface:",i7)' ) 
     &    undergnd_count
C  Compress the data array in the surface  analysis.
      DO i = 1 , nobs
        data(i + nobs) = data(ilat + i - 1)
      ENDDO
      DO i = 1 , nobs
        data(i + 2 * nobs) = data(ilon + i - 1)
      ENDDO
      DO i = 1 , nobs
        data(i + 3 * nobs) = data(ipmix + i - 1)
      ENDDO
      DO i = 1 , nobs
        data(i + 5 * nobs) = data(imix + i - 1)
      ENDDO
      DO i = 1 , nobs
        data(i + 6 * nobs) = data(irh + i - 1)
      ENDDO
      ikx = 1
      ilat = ikx + nobs
      ilon = ilat + nobs
      ipmix = ilon + nobs
      itype = ipmix + nobs
      imix = itype + nobs
      RETURN
      END
