c----program hazallXL modified on 12/20/01 to fix alog(zero) problem
c---- log - log interpolation
c--- sums annual prob. of exceed.
c--- from output of hazgridX and hazFX
c--- finds Poisson prob. of exceed. for time t or
c--- finds ground motions for specified annual probability of exceedance
c--- also can output file with hazard curves
c--- to run: hazallX 
c--- this version gets parameter file from header of 1st input file
c--- input files have hazard curves for each site concatenated
c---parfile= parameter file as used in hazgrid, etc. 
      type header_type
        character*30 :: name(6)
        real*4 :: period
        integer*4 :: nlev
        real*4 :: xlev(20)
        real*4 :: extra(10)
      end type header_type
      type(header_type) headr, head
      integer*4 readn
      dimension prob(6000000),out(600000),xlev(20)
      dimension prob0(6000000)
      character infile*256,name*256,namout*256,namepar*256
c
      do 9 i=1,6000000
 9      prob0(i)= 0.
      write(6,*) "enter an input file"
      read(5,900) infile
      write(6,*) "enter number of input files"
      read(5,*) nfile
      do 11 ifil=1,nfile
        write(6,*) ifil
        write(6,*) "enter name of prob. file"
        read(5,900) name
        write(6,900) name
 900    format(a30)
        call openr(name)
        ndata= 308
        call gethead(head,ndata,readn)
c------ get map parameters
        if(ifil.eq.1) then
          write(6,900) head%name(1)
c         open(unit=1,file=head%name(1))
          open(unit=1,file=infile)
c         headr%name(1)=head%name(1)
          headr%name(1)=infile
          read(1,*) ymin,ymax,dy
          read(1,*) xmin,xmax,dx
          write(6,*) "Y Values ",ymin,ymax,dy
          write(6,*) "X Values ",xmin,xmax,dx
          close(1)
          nx= ceiling((xmax-xmin)/dx+1)
          ny= ceiling((ymax-ymin)/dy+1)
          nrec= nx*ny
        endif
cccc
        write(6,900) head%name(1)
        write(6,900) head%name(2)
        write(6,*) head%period
        write(6,*) head%nlev
        write(6,*) (head%xlev(i),i=1,head%nlev)
        write(6,*) (head%extra(i),i=1,10)
        nlev= head%nlev
        ndata= nrec*nlev
        do 1 i=1,nlev
  1       xlev(i)= head%xlev(i)
        call getbuf2(prob,ndata,readn)
        write(6,*) ndata,readn
        write(6,*) "enter weight for probs"
        read(5,*) factor
        write(6,*) factor
        probm = 1.E4
        do 60 i=1,ndata
          if(prob(i).lt.probm) then
            probm = prob(i)
          endif
 60       prob(i)= prob0(i)+prob(i)*factor
        do 61 i=1,ndata
 61       prob0(i)= prob(i)
        headr%nlev= head%nlev
        write(6,*) "Minimum probability: ",probm
        do 95 i=1,nlev
 95       headr%xlev(i)= head%xlev(i)
        if(ifil.le.5) headr%name(ifil+1)= name
        if(ifil.le.10) headr%extra(ifil)= factor
        headr%period=head%period
 11   continue
      write(6,*) "enter 1 to just output summed prob file"
      read(5,*) iout
      if(iout.eq.1) then
        write(6,*) "enter name of output file"
        read(5,900) namout
        call openwx(ifp,namout)
        ndata=308
        call puthead(ifp,headr,ndata,readn)
        ndata= nrec*nlev
        call putbufx(ifp,prob,ndata,readn)
        write(6,*) ndata,readn
        go to 99
      endif
      write(6,*) "enter name of output file"
      read(5,900) namout
      call openwx(ifp,namout)
      write(6,*) "enter annual probability of exceedance"
      read(5,*) pex
      pex= alog(pex)
c      write(6,*) "enter rec number"
c      read(5,*) irec
c      write(6,*) prob(irec)
c      do 1 i=1,ndata
c 1    prob(i)= 1.-exp(-p(i)*t)
c
c
      write(6,*) "enter 0 to output prob gnd motion; 1 for Poisson prob"
      read(5,*) ip
      if(ip.eq.1) then
         write(6,*) "enter gnd motion level index,t"
         read(5,*) ilev,t
         do 100 i=1,nrec
         ind= (i-1)*nlev + ilev
  100    out(i)= 1.- exp(-prob(ind)*t)
         call openwx(ifp,namout)
         ndata=308
         headr%extra(9)= ilev
         headr%extra(10)= t
         call puthead(ifp,headr,ndata,ndata)
         call putbufx(ifp,out,nrec,readn)
         write(6,*) nrec,readn
         go to 99
         endif
      nlev2= nlev-1
c--- interpolate to find ground motions with specified pex
      do 10 i=1,nrec
c      if(prob(i).eq.0.) then
c        out(i)= alog(xlev(1))
c        go to 10
c        endif
      do 20 j=1,nlev2
      ind= (i-1)*nlev +j
      if(prob(ind).eq.0.) then
        out(i)= 0.
        go to 20
        endif
      y1= alog(prob(ind))
      y2= alog(prob(ind+1))
      x1= alog(xlev(j))
      x2= alog(xlev(j+1))
      if((j.eq.1).and.(pex.gt.y1)) then
        out(i)= x1
        go to 10
        endif
      if((j.eq.nlev2).and.(pex.lt.y2)) then
        out(i)= x2
        go to 10
        endif
      if((y1.ge.pex).and.(y2.le.pex)) then
        tmp= y2-y1
        if(tmp.eq.0.) write(6,*) i,j,y2,y1
        out(i)= x1+ (x2-x1)*(pex-y1)/tmp
        go to 10
        endif
 20   continue
 10   continue
      write(6,*) "enter scale factor"
      read(5,*) scale
      ymax= 0.
      do 90 i=1,nrec
      if(out(i).ne.0.) out(i)= exp(out(i))
      if(out(i).eq.0.) out(i)=  xlev(1)
c      write(6,*) out(i)
      out(i)= scale*(out(i))
 90   if(out(i).gt.ymax) ymax= out(i)
      write(6,800) ymax
 800  format(' ymax=',f12.1)
      headr%extra(8)= pex
      headr%extra(9)= scale
      headr%extra(10)= ymax
c---- output probabilistic ground motions
      ndata=308
      call puthead(ifp,headr,ndata,readn)
      call putbufx(ifp,out,nrec,readn)
      write(6,*) nrec,readn
 99   continue
      end
