program read_upa
c
c   program to read upper-air DYNAMO data in CSU format 
c
c   NOTE: this is preliminary GTS-resolution data with limited QC 
c  
c   written by P. Ciesielski (on 09/23/09) from Richard Johnson's 
c   group at CSU 
c
c   to create an executable on a LINUX box:
c   gfortran read_upa.f -o read_upa.x 
c
c   this program assumes that the "upa_file" file resides in 
c   same directory as the executable program
c
c   ml - max number of levels
c   xmis - missing value flag

      parameter (ml=220, xmis=-999.0)

      integer iyr, mon, idy, ihr, min, sfcht
      integer qp(ml), qh(ml), qt(ml), qd(ml), qw(ml)

      real    p(ml), h(ml), tc(ml), td(ml), dir(ml), spd(ml) 
      real    xlat(ml), xlon(ml)

      character stnid*5, upa_file*14

c   *******************************************************************

c
c   read in station wmo number
c

      print *, 'name of upa data file to read (e.g., upaqf.yymmddhh...>'
      read(5,'(a)') upa_file

c  open data file for specified station 

      open(8, file=upa_file, form='formatted', status='old')

      nsnd  = 0

c   read in data, iflg = 0 is normal return, iflg=1 if at EOF 

   25 call readdt (iyr, mon, idy, ihr, min, p, h, tc, td, dir, spd,
     2             qp, qh, qt, qd, qw,
     3             xlat, xlon, sfcht, nl, iflg)

      nsnd = nsnd + 1

c     do l=1,nl
c        if(p(l) .eq. 1000.) then
c           write(33,233) mon, idy, ihr, min, tc(l)
c        endif
c     enddo

  233 format(4i4,f7.2)

      write(6,30) iyr, mon, idy, ihr, min, nl
   30 format('read in data for ',5i2.2, ' with ',i3,' lines')

      if(iflg .eq. 1) go to 99 
      go to 25

   99 continue

      write(6,*) 'number of sounding read in  = ', nsnd

c   close files

      close(8)

      end

      subroutine readdt (iyr, mon, idy, ihr, min, p, h, t, td, 
     2   dir, spd, qp, qh, qt, qd, qw, xlat, xlon, sfcht, n, iflg)

c
c   routine to read data in one line at a time
c   data is in CSU format
c
c   Description of variables:
c
c   iyr - year
c   imn - month
c   idy - day
c   ihr - hour
c   sfcht - surface height
c   xlt  - latitude of site
c   xln  - longitude of site
c   wmon  - station wmo number
c           
c   p - pressure (mb)
c   h - height or altitude (m)
c   t - temperature (C)
c   td - dew point temperature (C)
c   dir - wind direction (degrees)
c   spd - wind speed (m/s)
c   xlat  - latitudinal position of balloon 
c   xlon  - longitude position of balloon
c
c   qp - quality flag on pressure
c   qh - quality flag on height
c   qt - quality flag on temperature
c   qd - quality flag on dew point
c   qw - quality flag on winds
c
c         Flag    Meaning
c        Value
c          1     parameter good
c          2     parameter questionable
c          3     parameter "visually" questionable
c          4     parameter bad
c          5     parameter "visually" bad
c          6     parameter interpolated
c          7     parameter estimated
c          8     parameter unchecked
c          9     parameter missing

c   ml - max number of data levels
c   n - actual number of data levels

      parameter (ml=220)
      real    p(ml), h(ml), t(ml), td(ml), dir(ml), spd(ml)
      real    xlat(ml), xlon(ml)
      integer qp(ml), qh(ml), qt(ml), qd(ml), qw(ml)
      integer iyr, mon, idy, ihr, min, n, iflg, sfcht
      character line*80, wmon*5

c  *****************************************************************

      n = 0
      iflg = 0

   10 read(8,'(A64)',end=21) line
      if (line(2:4) .eq. 'STN') then
         read (8,'(a64)',end=21) line
         read (line,702) iyr, mon, idy, ihr, min, sfcht, xlt, xln, 
     2                   wmon 
         read (8,'(a64)',end=21) line
         read (line(7:10),'(i4)') nlvl 
         read (8,'(a64)',end=21) line
   20    read (8,'(a64)',end=21) line
         if (line(2:4) .eq. 'STN') then
            backspace(8)
            if(n .ne. abs(nlvl)) then
               print *, mon, idy, ihr, min
               print *, 'n and nlvl not equal ', n, nlvl
               stop
            endif
            return
         else
            n = n + 1
            read (line,701) p(n), h(n), t(n), td(n), dir(n), spd(n),
     2           qp(n), qh(n), qt(n), qd(n), qw(n), xlon(n), xlat(n)
            go to 20
         endif
      else
         go to 10
      endif
   21 continue
      iflg = 1

c   format statements

  701 format(6(1x,f7.1),1x,5i3,2f8.2)
  702 format(8x,3i2,1x,2i2,1x,i5,2(1x,f8.2),2x,a5)

      return
      end