c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
      subroutine emadata(ns,qs,nh,nz,qz,thresh,no,ql,qu,tl,tu)
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c    
c    This routine arranges a data set using the EMA algorithm 
c
c    This was prepared for Aquaterra by Tim Cohn, US Geological Survey
c       
c    Timothy A. Cohn        05 December 2003
c
c
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c       input variables:
c       ---------------------------------------------------------------------------
c            ns         i*4  number of systematic observations
c            qs(ns)     r*4  vector of systematic record floods
c            nh         i*4  number of years in historical period (censored or not)
c            nz         i*4  number of historical floods (q >= thresh)
c            qz(nz)     r*4  vector of historical floods
c            thresh     r*8  threshold for historical floods
c
c       output variables: (input to "emafit")
c       ---------------------------------------------------------------------------
c            no         i*4  total number of observations (equivalent to HISTPD)
c            ql(n)      r*8  vector of lower bounds on floods
c            qu(n)      r*8  vector of upper bounds on floods
c            tl(n)      r*8  vector of lower bounds on flood threshold
c            tu(n)      r*8  vector of upper bounds on flood threshold
c
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c       N.B.  The output arrays, {ql,qu} and {tl,tu}, represent two distinct 
c             concepts that need to be understood clearly:
c            a.  The observations, which are described by {ql,qu};
c                  specifically, the true value of q is assumed to lie
c                    wihin the interval ql(i) <= q(i) <= qu(i);
c                  if q(i) is known exactly, then ql(i) = q(i) = qu(i);
c
c                  {ql,qu} are required by EMA to fit P3 dist'n to data
c
c            b.  The censoring pattern, which is described by {tl,tu}
c                  {tl,tu} define interval on which data get observed:
c                   if q is not inside {tl,tu}, q is either left or right censored.
c                  Examples:
c                    {-inf,+inf}  => Systematic data (everything known exactly)
c                    {T,   +inf}  => Historical data (q>=t known exactly; q<t censored)
c                
c                  {tl,tu} are required by EMA to estimate confidence intervals
c
c        Also:
c            EMAFIT requires the logarithms of the floods; this routine will work 
c            with either the real-space flood data or the log-space data.
c
c
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      implicit none
      
      integer 
     1  ns,nh,nz,no,
     2  i
      
      real  qs(*),qz(*)
      double precision  thresh,ql(*),qu(*),tl(*),tu(*),
     2  qmin,qmax
     
      data
     1  qmin, qmax/-1.d099, 1.d099/

         no    = 0
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c  systematic data
c     
      do 10 i=1,ns
        if (qs(i).gt.0) then
          no = no + 1
          ql(no) = log(qs(i))
          qu(no) = ql(no)
          tl(no) = qmin
          tu(no) = qmax
        end if
10    continue
cprh        no    = no+ns

c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c  historical exceedance data
c
      do 20 i=1,nz
        ql(i+no) = log(qz(i))
        qu(i+no) = log(qz(i))
        tl(i+no) = thresh
        tu(i+no) = qmax
20    continue
        no    = no+nz
      
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c  historical less-than data
c     
      do 30 i=1,nh-nz
        ql(i+no) = qmin
        qu(i+no) = thresh
        tl(i+no) = thresh
        tu(i+no) = qmax
30    continue
      
      return
      end