program create_qcstats
c
c   program to process hires sounding data (in uniform height resolution) 
c   creates a table of QC stats to be used by assign_qcf.f

c
c   gfortran create_qcstats_hgt_L4.f -o create_qcstats_hgt_L4.x
c
c  1) read in data
c  2) compute stats 
c  3) input from upaqf_stnid
c  4) output to stats_stnid

      real                           xmis
      parameter (maxz=10000, ntmx=250, xmis=-999.)
      real pres(maxz), zs(ntmx)
      real z(maxz), p(ntmx,maxz), tc(ntmx,maxz), td(ntmx,maxz)
      real pp, zz, tt, dd, dir, spd, zinc
      real uc(ntmx,maxz), vc(ntmx,maxz)
      integer qp(ntmx,maxz), qh(ntmx,maxz), qt(ntmx,maxz) 
      integer qd(ntmx,maxz), qw(ntmx,maxz)
      integer qqp, qqh, qqt, qqd, qqw
      integer nt, nz, nu, nv, nd
      integer idate(ntmx), itime(ntmx)

      real pave(maxz),  psd(maxz)
      real tave(maxz),  tsd(maxz)
      real tdave(maxz), tdsd(maxz)
      real uave(maxz),  usd(maxz)
      real vave(maxz),  vsd(maxz)

   
      real      xlon, xlat, alt
      integer   iyr, mon, idy, ihr, min
      integer   imm, idd, ihh
     
      character fin*11, fstats*15, header(5)*64, resp*1

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

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

      print *, 'read in stn number (i5) ...'
      read(5,'(i5)') 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 *, 'read in height resolution (m) of data ...'
      read(5,*) zinc


c   initialize fields

      do i=1,ntmx
      do k=1,maxz
         p(i,k) = xmis 
         tc(i,k) = xmis 
         td(i,k) = xmis 
         uc(i,k) = xmis 
         vc(i,k) = xmis 
      enddo
      enddo

      print *, 'enter stn elev ...>'
      read(5,*) zmax  
      zstart = int(zmax)/zinc * zinc
      print *, zstart, zmax

      z(1) = zmax
      do k=2,maxz
         z(k) = ifix((zmax + (k-1)*zinc)/zinc)*zinc
c        print *, k, z(k)
      enddo

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

      write(fstats, 9) istnid
    9 format('stats_hgt_',i5)
      open (24, file=fstats, status='unknown', form='formatted')

c   read header information

      nt = 0
    5 nt = nt + 1

      do il=1,4
         read(20,'(a64)',end=50) header(il)
      enddo

      read(header(2)(9:14),*)   idate(nt) 
      read(header(2)(16:19),*)   itime(nt) 
      read(header(3)(7:10),'(i4)') nlv
c     print *, idate(nt), itime(nt), nlv

      do k=1,nlv
  
          read (20,200,end=50) pp, zz, tt, dd, dir, spd,
     2    qqp, qqh, qqt, qqd, qqw
c         print 200, pp, zz, tt, dd, dir, spd,
c    2    qqp, qqh, qqt, qqd, qqw
          if(k .eq. 1) then
             izindex = 1
             zs(nt) = zz
             print *, idate(nt), itime(nt), nlv, nt, tt
          else
             izindex = (zz - zstart)/zinc + 2
          endif
c         print *, k, zz, izindex
          p(nt,izindex) = pp
          tc(nt,izindex) = tt
          td(nt,izindex) = dd
          call spdr_uv (spd, dir, uc(nt,izindex), vc(nt,izindex), xmis)
          qp(nt,izindex) = qqp
          qh(nt,izindex) = qqh
          qt(nt,izindex) = qqt
          qd(nt,izindex) = qqd
          qw(nt,izindex) = qqw
      enddo

      go to 5

   50 continue

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

  999 continue

      nsd = nt-1
      print *, 'number of times read in ', nsd

      do k=1,maxz

         call mean(p(1,k), nsd, pave(k), xmis, nz)
         call varc(p(1,k), nsd, var, psd(k), pave(k), xmis)
         call mean(tc(1,k), nsd, tave(k), xmis, ntc)
         call varc(tc(1,k), nsd, var, tsd(k), tave(k), xmis)
         call mean(td(1,k), nsd, tdave(k), xmis, nd)
         call varc(td(1,k), nsd, var, tdsd(k), tdave(k), xmis)
         call mean(uc(1,k), nsd, uave(k), xmis, nu)
         call varc(uc(1,k), nsd, var, usd(k), uave(k), xmis)
         call mean(vc(1,k), nsd, vave(k), xmis, nv)
         call varc(vc(1,k), nsd, var, vsd(k), uave(k), xmis)
c        print *, p(k), nsd, nz, ntc, nd, nu, nv
c        write(22,222) k, pres(k), zave(k), zsd(k), nz, tave(k), tsd(k), 
c    2                    ntc,tdave(k), tdsd(k), nd, uave(k), usd(k), 
c    3                    nu, vave(k), vsd(k), nv
      enddo

c remove outliers in data fields (> fac*std) and recompute stats

      fac = 4.

      do k=1,maxz
         do i=1,nt
  
            if(p(i,k) .ne. xmis) then
               dif = abs(p(i,k) - pave(k))
               if(dif .gt. fac*psd(k)) then
                  p(i,k) = xmis 
               endif
            endif
            if(tc(i,k) .ne. xmis) then
               dif = abs(tc(i,k) - tave(k))
               if(dif .gt. fac*tsd(k)) then
                  tc(i,k) = xmis 
               endif
            endif
            if(td(i,k) .ne. xmis) then
               dif = abs(td(i,k) - tdave(k))
               if(dif .gt. fac*tdsd(k)) then
                  td(i,k) = xmis 
               endif
            endif
            if(uc(i,k) .ne. xmis) then
               dif = abs(uc(i,k) - uave(k))
               if(dif .gt. fac*usd(k)) then
                  uc(i,k) = xmis 
               endif
            endif
            if(vc(i,k) .ne. xmis) then
               dif = abs(vc(i,k) - vave(k))
               if(dif .gt. fac*vsd(k)) then
                  vc(i,k) = xmis 
               endif
            endif

         enddo
          
         call mean(p(1,k), nsd, pave(k), xmis, nz)
         call varc(p(1,k), nsd, var, psd(k), pave(k), xmis)

         call mean(tc(1,k), nsd, tave(k), xmis, ntc)
         call varc(tc(1,k), nsd, var, tsd(k), tave(k), xmis)

         call mean(td(1,k), nsd, tdave(k), xmis, nd)
         call varc(td(1,k), nsd, var, tdsd(k), tdave(k), xmis)

         call mean(uc(1,k), nsd, uave(k), xmis, nu)
         call varc(uc(1,k), nsd, var, usd(k), uave(k), xmis)

         call mean(vc(1,k), nsd, vave(k), xmis, nv)
         call varc(vc(1,k), nsd, var, vsd(k), uave(k), xmis)

         if(nz.eq.0 .and. ntc.eq.0 .and. nd.eq.0 .and. nu.eq.0) then
            go to 9989
         endif

         write(24,222) k, z(k), pave(k), psd(k), nz, tave(k), tsd(k), 
     2                    ntc,tdave(k), tdsd(k), nd, uave(k), usd(k), 
     3                    nu, vave(k), vsd(k), nv
     
      enddo

 9989 continue
      close(24)

  222 format(i5,f7.0,1x,3(f6.1,f6.2,i4),2(f7.1,f6.1,i4))
  223 format(3i6,f8.1)

      end

      SUBROUTINE VARC(DATA, N, VAR, SDEV, XMEAN, XMIS)
c
c  routine to compute the standard deviation and variance
c
      REAL DATA(1), VAR, XMEAN
      M = 0
      SUM = 0.
      VAR = -99. 
      SDEV = -99. 
      DO 20  I=1,N
         IF(DATA(I) .EQ. XMIS) GO TO 20
         M = M + 1
         SUM = SUM + (DATA(I) - XMEAN)**2
   20 CONTINUE
      if(m .gt. 1) then
         VAR =  SUM / FLOAT(M-1)
         SDEV = SQRT(VAR)
      endif
                                                                                      
      RETURN
      END
                                                                                      
      SUBROUTINE MEAN(DATA, N, EZMEAN, XMIS, nn)
C
C   THIS CALCULATES THE MEAN OF A SET OF DATA
C
      REAL DATA(1)
c   *********************************************
      SUM = 0.
      nn = 0
      ezmean = xmis
      DO 10  I=1,N
         if(data(i) .ne. xmis) then
            SUM = SUM + DATA(I)
            nn = nn + 1
         endif
   10 CONTINUE
      if(nn .gt. 1) then
         EZMEAN = SUM/FLOAT(NN)
      endif
                                                                                      
      RETURN
      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 spdr_uv (spd, dir, u, v, xmis)
c
c   convert u and v wind components into speed and direction
c

      pi = acos(-1.0)
      degran = pi / 180.

      u = xmis
      v = xmis
      if (dir.gt.360. .or. dir.lt.0. .or. spd.lt.0.) then
         return
      else
         u = -1.* spd * sin(dir*degran)
         v = -1.* spd * cos(dir*degran)
      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(p, pinc, fld, nl, pi, fldi, np, xmis) 

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

c  interpolated fields
      real pi(np), fldi(np)

c  missing value flag
      real xmis

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

c     set surface value to first raw value

      fldi(1) = fld(1)

      do 70  il=2,np

c  set field to missing but fill in if possible

         fldi(il) = imis

         pa = pi(il) - pinc/2.
         if(pa .lt. p(nl)) then
            pa = p(nl)
         endif
         pb = pi(il) + pinc/2.
         if(pb .gt. p(1)) then
            pb = p(1)
         endif

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

            if(p(iraw) .lt. pb .and. p(iraw).gt.pa) 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 interp(p, fld, nl, pi, fldi, np, xmis) 

c   this routine linearly interpolates a field in log pressure 
c
c   maximum interval over which to do linear interpolation
c
      real       pmax,       imis
      parameter (pmax = 50., imis=-999.)
      
c   raw fields
      real p(nl), fld(nl)

c  interpolated fields
      real pi(np), fldi(np)

c  missing value flag
      real xmis

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

c     set surface value to first raw value

      fldi(1) = fld(1)

      do 70  il=2,np

c  set field to missing but fill in if possible

         fldi(il) = imis

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

           if(pi(il) .eq. p(iraw)) then
               if(fld(iraw) .eq. xmis) go to 10
               fldi(il) = fld(iraw)
c              print *, p(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 iplv - pressure level in raw data closest to p(il)

   10    continue
         pdif = 1.e4
         do iraw = 1,nl
            dp = abs(p(iraw) - pi(il))
c           if(iraw .lt. 20 )print *, iraw, p(iraw), pi(il), dp, pdif, iplv
            if(dp .lt. pdif) then
               pdif = dp
               iplv = iraw
            endif
         enddo
c        print *, il, pi(il), p(iplv)

c  find first level below pi(il) in which data is not missing
         do 20 ip=iplv,1,-1
            if(p(ip) .lt. pi(il)) go to 20
            if(fld(ip) .ne. xmis) then
               ipb = ip
               go to 30
            endif
   20    continue
         go to 70 

   30    continue
c  find first level above pi(il) in which data is not missing
         do 40 ip=iplv,nl
            if(p(ip) .gt. pi(il)) go to 40
            if(fld(ip) .ne. xmis) then
               ipa = ip
               go to 50
            endif
   40    continue
         go to 70 

   50    continue

         if(abs(p(ipb)-p(ipa)) .gt. pmax) go to 70

          ap  = alog(pi(il))
          apb = alog(p(ipb))
          apa = alog(p(ipa))
c         print *, il, iplv, ipb, ipa
c         print *, pi(il), p(ipb), p(ipa)
c         print *, fld(ipb), fld(ipa)

          sl = (ap - apb) / (apa - apb)
          fldi(il) = fld(ipb) + sl*(fld(ipa) - fld(ipb))
c         print *, fld(ipb), fldi(il)

   70 continue

      return
      end

      subroutine readupa(iun, lvl)
c
c   read in data for one time

      parameter (m=10000, pmis=9999.)

      real    time, rh, wc, wspd, wdir, dz, xln, xlt, rng
      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)

      character hline*130, dline*130
                                                                                     
c   *************************************************************
c
c routine to read in uncorrected souding data
c
      iflag = 0
 10   read (iun,'(a82)',end=100) hline
      if (hline(1:10) .eq. 'Launch Loc') then
c  read in lon, lat, alt info
         read (hline(61:70),*) xlon
         read (hline(72:78),*) xlat
         read (hline(81:81),*) alt 
c  read in date / time information
         read (iun,'(a82)') 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,'(a130)',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), wc, wspd, wdir, dz, xln, xlt, rng, 
     3        z0(lvl), qp(lvl), qt(lvl), qh(lvl), qu(lvl), qv(lvl)
c          print 702,  p0(lvl), z0(lvl), t0(lvl), d0(lvl), 
c    2               u0(lvl), v0(lvl), qp(lvl), qt(lvl), qu(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