program read_upaL4    
c
c   program to read upper-air Level 4 in CSU format 
c   level 4 data in GTS-res or 5-hPa interpolated with QC flags
c  
c   written by P. Ciesielski 
c
c   to create an executable on a LINUX box:
c   gfortran read_upaL4.f -o read_upaL4.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*20

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

c
c   read in station wmo number
c

      print *, 'enter the name of file to read (e.g.,upaqi_72797) ...>'
      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 

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

          nsnd = nsnd + 1

         write(6,30) stnid, iyr, mon, idy, ihr, min, nl

c        do l=1,nl
c           write(6,35) p(l), h(l), tc(l), td(l), dir(l), spd(l), 
c    2      xlat(l), xlon(l)
c        enddo

         if(iflg .eq. 1) exit

      enddo     

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

c   close files

      close(8)
   30 format('read in data for ',a,2x,3i2.2,2x,2i2.2, ' with ',
     2        i3,' lines')
   35 format(8f10.2)

      end

      subroutine readdt (iyr, mon, idy, ihr, min, p, h, t, td, 
     2   dir, spd, qp, qh, qt, qd, qw, xlat, xlon, sfcht, n, wmon, 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,'(a80)',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)
c            print *, line 
                 if(dir(n) .eq. 0. .and. spd(n) .gt. 0.) then
                    dir(n) = 360.
                 endif
            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