program process_snd
c
c   program to process hires sounding data obtained from ASPEN in 
c   EOL or CLS format
c   creates a uniform height dataset ... upazi_stnid
c
c   gfortran create_upa_hgt_L4.f -o create_upa_hgt_L4.x
c
c  1) read in data
c  2) interpolate data to specific uniform vertical height levels 
c  3) put qc flags on data (all flags are set as good)

      real               xmis,      zmis,        ixmis
      parameter (ip=10000, xmis=-999., zmis=-999., ixmis=-999.)
      real p(ip), z(ip), tc(ip), td(ip), uc(ip), vc(ip)
      real fuc(ip), fvc(ip), xlni(ip), xlti(ip)
      real zsfc(1000)
      integer qp(ip), qh(ip), qt(ip), qd(ip), qu(ip)
      integer istnid, nf
   
      parameter (m=10000)

      real      xlon, xlat, alt
      integer   iyr, mon, idy, ihr, min
      integer   imm, idd, ihh
     
      common /header/ xlon, xlat, alt, iyr, mon, idy, ihr, min

      common /data/ p0(m), z0(m), t0(m), d0(m), u0(m), v0(m), 
     2        qpp(m), qtt(m), qhh(m), quu(m), qvv(m), xln(m), xlt(m)

      character fname*33, cid*5, fout*12, resp*1

      integer igood, iques, idoub, ibad, iint, iest, imis, iunc
      data igood, iques,  ibad, iint, iest, iunc, imis
     2    /    1,     2,   4,    6,    7,    8,    9/

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

c   specify constants
                             
      cid = '     '          
      print *, 'read in stnid (e.g., 46810) ...>'
      read(5,*) istnid 
      print *, 'istnid = ', istnid
      print *, 'Is this correct (y/n) ...'
      read(5,*) resp
      if(resp .eq. 'n' .or. resp .eq. 'N') then
         print *, 'incorrect cid and istnid'
         stop
      endif
      print *, 'if humidity has been corrected type 1, else 0 ...>'
      read(5,*) ihumcor

      write(fout, 8) istnid
    8 format('upazi_',i5)
      open (20, file=fout, status='unknown', form='formatted')

c
c  specify resolution of dataset in meter 
c  note: rise rate of the balloon is typically 5 m/s and will 1 to 2 sec data
c        5 meters represents a reasonable limit for highest resolution 
c
      print *, 'Specify desired uniform resolution in meters ' 
      print *, '(5 m greater, where 50 meter is generally adequate '
      print *, 'for most research applications ...'
      read (5,*) zinc

      r = 287.054
      g = 9.806
      pi = acos(-1.)
      degran = pi / 180.

c  contains list of files to be processed

      open(2, file='files_hires', form='formatted', status='old')

      nf = 0

    5 continue

c  set quality control flags to good at all levels
      do l=1,ip
         qp(l) = igood
         qh(l) = igood
         qt(l) = igood
         qd(l) = igood
         qu(l) = igood
      enddo

c     read in file 

      read(2,'(a)',end=999) fname

      nf = nf + 1
      print *, 'processing data for ', nf, fname

c   open data files

c   file 8 contain input hi-res data 

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

c
c   read in one data for one time from file 8
c   

c   select format(cls or eol) that input data is in by 
c   uncommenting appropriate line
c     call readcls(8, nl)
      call readeol(8, nl)

c
c  get height level to interpolate the data to
c
      zsfc(nf) = z0(1)
      z(1) = z0(1)
      z(2) = ifix((z0(1)+zinc)/zinc) * zinc
      print *, z(1), z(2) 
      if(z(1) .eq. z(2)) then
         print 233, mon, idy, ihr
 233     format(40x, 'check ', 3i3)
         z(2) = z(1) - zinc
      endif
      n = 2
   60 n = n + 1
      z(n) = z(n-1) + zinc
      if(z(n) .lt. z0(nl)) go to 60
      nz = n - 1
      if(nz .gt. ip) then
         print *, 'ip too small ...', nz, ip
         stop
      endif

cc     do i=1,10
cc        print *, i, z(i)
cc     enddo
cc     do i=nz-10,nz
cc        print *, i, z(i)
cc     enddo

c  get data at interpolated levels

      call interz(z0, p0, nl, z, p, nz, zmis)
      
      call interz(z0, t0, nl, z, tc, nz, xmis)
      call interz(z0, d0, nl, z, td, nz, xmis)
      call layave(z0, zinc, u0, nl, z, uc, nz, xmis)
      call layave(z0, zinc, v0, nl, z, vc, nz, xmis)
c     call interz(z0, u0, nl, z, uc, nz, xmis)
c     call interz(z0, v0, nl, z, vc, nz, xmis)
      call fil121(uc, nz, fuc, ixmis)
      call fil121(vc, nz, fvc, ixmis)
      call interz(z0, xln, nl, z, xlni, nz, xmis)
      call interz(z0, xlt, nl, z, xlti, nz, xmis)

c  compute date/time tag

      call dattim(iyr, mon, idy, ihr, min, imm, idd, ihh)  
      write(11,222) fname, mon, idy, ihr*100+min
  222 format(a19,3i5)

c   print out header information

      write(20,700)
  700 format(' STN     DATE   GMT   HTS     LAT     LONG    ID')
c     write(20,701) cid, iyr,imm,idd,ihh*100, nint(alt), xlat,
c    2              xlon, istnid
      write(20,701) cid, iyr,mon,idy,ihr*100+min, nint(alt), xlat,
     2              xlon, istnid
  701 format(a4, 4x, 3i2.2, 1x, i4.4, 1x, i5, 2(1x,f8.2), i7)

      if (ihumcor .eq. 0) then
         write(20,702) nz
      else
         write(20,704) nz
      endif
  702 format('NLVL =',i4,' COMMENTS:')
  704 format('NLVL =',i4,' COMMENTS: humidity corrected')
      write(20,703)
  703 format('     P      HT      TC       TD     DIR      SPD  QP QH', 
     2       ' QT QD QW   LON      LAT')

      do it=1,nz
  
          if(p(it) .eq. xmis) qp(it) = imis
          if(z(it) .eq. xmis) qh(it) = imis
          if(tc(it) .eq. xmis) qt(it) = imis
          if(td(it) .eq. xmis) qd(it) = imis
          if(fuc(it) .ne. xmis .and. fvc(it) .ne. xmis) then
             call spd_dir(fuc(it), fvc(it), spd, dir, ixmis)
             qu(it) = igood
          else
             spd = ixmis
             dir = ixmis
             qu(it) = imis
          endif
          if(xlni(it) .eq. ixmis) then
             xlni(it) = xlni(it-1)
          endif
          if(xlti(it) .eq. ixmis) then
             xlti(it) = xlti(it-1)
          endif

          write (20,200) p(it),z(it),tc(it),td(it), dir, spd,
     2      qp(it), qh(it), qt(it), qd(it), qu(it), xlni(it), xlti(it)

      enddo

  200 format(6f8.1,1x,5i3,2f8.2)

      go to 5

  999 continue

      ntot = nf - 1

      print *, 'number of sondes processed ', ntot

      call mean (zsfc, ntot, azsfc, xmis, nn)

      print *, 'mean surface height  = ', azsfc
      print *, 'number of obs  = ', nn


      end

      subroutine fil121(fld, n, ffld, xmis)
c
c   filter data with 1-2-1 filter
c
      real fld(n), ffld(n), xmis

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

      ffld(1) = fld(1)
      ffld(n) = fld(n)
      do i=2,n-1
 
         ip = i + 1
         im = i - 1
         if(fld(ip).ne.xmis .and. fld(i).ne.xmis 
     2                      .and. fld(im).ne.xmis) then 
            ffld(i) = 0.25*(fld(im)+fld(ip)) + 0.5*fld(i)
         else
            ffld(i) = fld(i)
         endif

      enddo

      return
      end 
      subroutine dattim(iyr, mon, idy, ihr, min, imm, idd, ihh)

      integer iyr, mon, idy, ihr, min, imm, idd, ihh

c   add one hour to sonde time (sondes are typically launched 1 hr
c   prior to time that obs are needed)

      ihh = ihr + nint(float(min)/60.) + 1 
      imm = mon
      idd = idy 

      if(ihh .eq. 24) then
         ihh = 1
         idd = idy + 1
         if(mon .eq. 6 .and. idd .gt. 30) then
            imm = 7 
            idd = 1
         elseif(mon .eq. 7 .and. idd .gt. 31) then
            imm = 8 
            idd = 1
         elseif(mon .eq. 8 .and. idd .gt. 31) then
            imm = 9 
            idd = 1
         endif
      elseif(ihh .gt. 25) then
         print *, 'error in dattim', ihh, imm, idd
         stop
      endif

      return
      end

      subroutine spd_dir(u, v, spd, dir, xmis)
c
c   convert u and v wind components into speed and direction
c
      pi = acos(-1.)
      ran2deg = 180. / pi

      if(u.ne.xmis .and. v.ne.xmis) then
         spd = sqrt (u*u + v*v)
         if (spd .eq. 0.0) then
            dir = 0.
         else
            if(u .eq. 0.) then
               if (v .ge. 0.) then
                  dir = 180.
               else
                  dir = 360.
               endif
            elseif(v .eq. 0.) then
               if (u .ge. 0.) then
                  dir = 270.
               else
                  dir = 90.
               endif
            else
               dir = acos(-v/spd) * ran2deg
               if(u .gt. 0) dir = 360. - dir
            endif
         endif
      else
         spd = xmis
         dir = xmis
      endif

      return
      end

      subroutine layave(z, zinc, fld, nl, zi, fldi, nz, xmis) 

c   this routine computes a smooth profile by taking layer averages
c
c   maximum interval over which to do linear interpolation
c
      real       zmax,       imis
      parameter (zmax = 500., imis=-999.)
      
c   raw fields
      real z(nl), fld(nl)

c  interpolated fields
      real zi(nz), fldi(nz)

c  missing value flag
      real xmis

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

c     set surface value to first raw value

      if(fld(1) .ne. xmis) then
         fldi(1) = fld(1)
      else
         fldi(1) = imis
      endif

      do 70  il=2,nz

c  set field to missing but fill in if possible

         fldi(il) = imis

         za = zi(il) + zinc/2.
         if(za .gt. z(nl)) then
            za = z(nl)
         endif

         zb = zi(il) - zinc/2.
         if(zb .lt. z(1)) then
            zb = z(1)
         endif

c        print *, il, zb, za

         sum = 0.
         nobs  = 0
         do iraw = 1,nl

            if(z(iraw) .gt. zb .and. z(iraw).lt.za) then
               if(fld(iraw) .ne. xmis) then
                  sum = sum + fld(iraw)
                  nobs  = nobs  + 1
               endif
            endif

         enddo

         if(nobs .gt. 0) then
            fldi(il) = sum / float(nobs)
         endif 

   70 continue

      return
      end

      subroutine interz(z, fld, nl, zi, fldi, nz, xmis) 

c   this routine linearly interpolates a field in height
c
c   maximum interval over which to do linear interpolation
c
      real       zmax,       imis
      parameter (zmax = 500., imis=-999.)
      
c   raw fields
      real z(nl), fld(nl)

c  interpolated fields
      real zi(nz), fldi(nz)

c  missing value flag
      real xmis

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

c     set surface value to first raw value

      fldi(1) = fld(1)

      do 70  il=2,nz

c  set field to missing but fill in if possible

         fldi(il) = imis

c   this loop gets fld if interpolated height level is a raw data level
         do iraw = 1,nl

           if(zi(il) .eq. z(iraw)) then
               if(fld(iraw) .eq. xmis) go to 10
               fldi(il) = fld(iraw)
c              print *, z(il), fldi(il)
               go to 70
           endif
         enddo

c   this loop get fld if interpolated pressure level is not a raw data level
c   then interpolated value in linearly interpolated in log of pressure

c   find izlv - zlevel in raw data closest to z(il)

   10    continue
         zdif = 1.e4
         do iraw = 1,nl
            dz = abs(z(iraw) - zi(il))
c           if(iraw .lt. 20 ) then
c             print *, iraw, z(iraw), zi(il), dz, zdif, izlv
c           endif
            if(dz .lt. zdif) then
               zdif = dz
               izlv = iraw
            endif
         enddo
c        print *, il, izlv, zi(il), z(izlv)


c  find first level below zi(il) in which data is not missing
         do 20 iz=izlv,1,-1
c           print *, z(iz), zi(il), fld(iz)
            if(z(iz) .gt. zi(il)) go to 20
            if(fld(iz) .ne. xmis) then
               izb = iz
               go to 30
            endif
   20    continue
         go to 70 

   30    continue
c  find first level above zi(il) in which data is not missing
         do 40 iz=izlv,nl
c           print *, z(iz), zi(il), fld(iz)
            if(z(iz) .lt. zi(il)) go to 40
            if(fld(iz) .ne. xmis) then
               iza = iz
               go to 50
            endif
   40    continue
         go to 70 

   50    continue
c        print *, iza, izb

         if(abs(z(izb)-z(iza)) .gt. zmax) go to 70

          az  = zi(il)
          azb = z(izb)
          aza = z(iza)
c         print *, il, izlv, izb, iza
c         print *, zi(il), z(izb), z(iza)
c         print *, fld(izb), fld(iza)

          sl = (az - azb) / (aza - azb)
          fldi(il) = fld(izb) + sl*(fld(iza) - fld(izb))

   70 continue

      return
      end

      subroutine readeol(iun, lvl)
c
c   read in data for one time
c   assuming data is in NCAR CLass EOL

      parameter (m=10000, pmis=-999.)

      real    time, rh, wc, wspd, wdir, dz, xlnt, xltt
      integer iyr, mon, idy, ihr, min, ih, im
      real    sec, gpsalt

      common /header/ xlon, xlat, alt, iyr, mon, idy, ihr, min

      common /data/ p0(m), z0(m), t0(m), d0(m), u0(m), v0(m), 
     2        qp(m), qt(m), qh(m), qu(m), qv(m), xln(m), xlt(m)

      character hline*100, dline*135
                     
c   *************************************************************
c
c routine to read in uncorrected souding data
c
      iflag = 0
 10   read (iun,'(a100)',end=100) hline
      if (hline(1:4) .eq. 'Data') then
c  skip over a few lines
         read(iun, '(a100)',end=100) hline
         read(iun, '(a100)',end=100) hline
         read(iun, '(a100)',end=100) hline
c  read in lon, lat, alt info
         read(iun, '(a100)',end=100) hline
         read (hline(56:65),*) xlon
         read (hline(79:87),*) xlat
         read (hline(91:97),*) alt
c  read in date / time information
         read (iun,'(a100)') hline
         read (hline(46:47),*) iyr
         read (hline(50:51),*) mon
         read (hline(54:55),*) idy
         read (hline(58:59),*) ihr
         read (hline(61:62),*) min
c  read remainder of header lines
         do i=1,8
            read (iun,'(a82)') hline
         enddo
      else
         go to 10
      endif
           
C------- read in data     ---------------------------------------------------
      lvl = 0
      nl  = 0
   20 read (iun,'(a133)',end=100) dline
        nl = nl + 1
        read (dline,*) time, ih, im, sec, pres
        if(pres .ne. pmis) then
           lvl = lvl + 1
           if(lvl .gt. m) go to 100
           read (dline,*) time, ih, im, sec, p0(lvl), t0(lvl), d0(lvl), 
     2        rh, u0(lvl), v0(lvl), wspd, wdir, dz, z0(lvl), xlnt, 
     3        xltt, gpsalt 
           if(xlnt .eq. -999.) then
              if(lvl .eq. 1) then
                 xln(lvl)  = xlon 
              else
                 xln(lvl)  =  xln(lvl-1)
              endif
           else
              xln(lvl)  =  xlnt
           endif
           if(xlt(lvl) .eq. -999.) then
              if(lvl .eq. 1) then
                 xlt(lvl)  = xlat 
              else
                 xlt(lvl)  =  xlt(lvl-1)
              endif
           else
              xlt(lvl)  =  xltt
           endif
c          print 702,  p0(lvl), z0(lvl), t0(lvl), d0(lvl), 
c    2               u0(lvl), v0(lvl), qp(lvl), qt(lvl), qu(lvl)
c          print 702,  p0(lvl), xlnt, xln(lvl), xltt, xlt(lvl)
           write(49,*)  lvl, p0(lvl), z0(lvl)
        endif
        go to 20

  100 continue
      print 222, xlon, xlat, alt, 
     2           iyr, mon, idy, ihr, min,nl,lvl

  222 format (2f7.2,f4.1,7i5,f7.2)
  702 format (9(1x,f7.1))

      return
      end

      subroutine readcls(iun, lvl)
c
c   read in data for one time
c   assuming data is in NCAR CLass format

      parameter (m=10000, pmis=9999.)

      real    time, rh, wc, wspd, wdir, dz, rng, xlnt, xltt
      integer iyr, mon, idy, ihr, min

      common /header/ xlon, xlat, alt, iyr, mon, idy, ihr, min

      common /data/ p0(m), z0(m), t0(m), d0(m), u0(m), v0(m), 
     2        qp(m), qt(m), qh(m), qu(m), qv(m), xln(m), xlt(m)

      character hline*100, dline*100
                     
c   *************************************************************
c
c routine to read in uncorrected souding data
c
      iflag = 0
 10   read (iun,'(a100)',end=100) hline
      if (hline(1:4) .eq. 'Data') then
c  skip over a few lines
         read(iun, '(a100)',end=100) hline
         read(iun, '(a100)',end=100) hline
c  read in lon, lat, alt info
         read(iun, '(a100)',end=100) hline
c   read for 46734
         read (hline(61:69),*) xlon
         read (hline(71:78),*) xlat
         read (hline(80:82),*) alt
c        read (hline(63:73),*) xlon
c        read (hline(74:82),*) xlat
c        read (hline(84:85),*) alt
c  read in date / time information
         read (iun,'(a100)') hline
         read (hline(38:39),*) iyr
         read (hline(42:43),*) mon
         read (hline(46:47),*) idy
         read (hline(50:51),*) ihr
         read (hline(53:54),*) min
c  read remainder of header lines
         do i=1,10
            read (iun,'(a82)') hline
         enddo
      else
         go to 10
      endif
           
C------- read in data     ---------------------------------------------------
      lvl = 0
      nl  = 0
   20 read (iun,'(a100)',end=100) dline
        nl = nl + 1
        read (dline,*) time, pres
        if(pres .ne. pmis) then
           lvl = lvl + 1
           if(lvl .gt. m) go to 100
           read (dline,*) time, p0(lvl), t0(lvl), d0(lvl), rh, 
     2        u0(lvl), v0(lvl), wspd, wdir, dz, xlnt, 
     3        xltt, rng, az, z0(lvl)
           if(xlnt .ge. 999.) then
              if(lvl .eq. 1) then
                 xln(lvl)  = xlon 
              else
                 xln(lvl)  =  xln(lvl-1)
              endif
           else
              xln(lvl)  =  xlnt
           endif
           if(xlt(lvl) .ge. 999.) then
              if(lvl .eq. 1) then
                 xlt(lvl)  = xlat 
              else
                 xlt(lvl)  =  xlt(lvl-1)
              endif
           else
              xlt(lvl)  =  xltt
           endif
c          print 702,  p0(lvl), z0(lvl), t0(lvl), d0(lvl), 
c    2               u0(lvl), v0(lvl), qp(lvl), qt(lvl), qu(lvl)
c          print 702,  p0(lvl), xlnt, xln(lvl), xltt, xlt(lvl)
c          write(49,702)  p0(lvl), xlnt, xln(lvl), xltt, xlt(lvl)
        endif
        go to 20

  100 continue
      print 222, xlon, xlat, alt, 
     2           iyr, mon, idy, ihr, min,nl,lvl,p0(lvl)

  222 format (2f7.2,f4.1,7i5,f7.2)
  702 format (9(1x,f7.1))

      return
      end

      SUBROUTINE MEAN(DATA, N, EZMEAN, XMIS, nn)
C
C   THIS CALCULATES THE MEAN OF A SET OF DATA
C
      REAL DATA(1)
      SUM = 0.
      nn = 0
      DO 10  I=1,N
         if(data(i) .ne. xmis) then
            SUM = SUM + DATA(I)
            nn = nn + 1
         endif
   10 CONTINUE
      if(nn .gt. 0) then
         EZMEAN = SUM/FLOAT(NN)
      else
         EZMEAN = xmis
      endif
      RETURN
      END

      SUBROUTINE PFMIMA_mis( N, X, XMIN, XMAX, xmis )
      INTEGER            N
      REAL               X(N), XMIN, XMAX
C
C   Finds the minimum and maximum values in the vector  X  of length  N
C
      INTEGER    I
      REAL       BIG
      PARAMETER  ( BIG = 1.0E29 )
C
      XMIN = +BIG
      XMAX = -BIG
      DO 10 I=1,N
          if(x(i) .eq. xmis) go to 10
          XMIN = MIN( XMIN, X(I) )
          XMAX = MAX( XMAX, X(I) )
   10 CONTINUE
      RETURN
      END