      program separate
c**************************************************************************
c                                                                         *
c    This program performs hydrograph separation and supplies             *
c    streamflow and baseflow using methods described in Heath and Trainer *
c    1968 "introduction to groundwater hydrology".  Programmed by:        *
c    Cynthia E. Sellinger,Physical Scientist, Great Lakes Environmental   *
c    Research Laboratory,                                                 *
c    National Oceanic and Atmospheric Administration, Department of       *
c    Commerce.                                                            *
c**************************************************************************
      implicit none
      integer count4,nom1,nom2,count5,count1,count15
      real a,b,tmin,area,days
      character*40 infile1,infile2,outfile
c
      print *,' '
      print *,' '
      print *,' '
      print *,' '
      print *,' '
      print *,' '
      print *,' '
      print *,' '
      print *,' '
      print *,' '
      print *,' '
      print *,' '
      print *,' '
      print *,' '
      print *,' '
      print *,' '
      print *,' Note:  Precipitation values must be in units of'
      print *,' Millimeters per day, Flow values should be in units'
      print *,' of cubic meters per second, groundwater stage'
      print *,' files should be daily levels in meters above a '
      print *,' common datum, and the watershed area should be '
      print *,' entered in square kilometers '
      print *,' '
      print *, ' Enter the precipitation and streamflow file name   '
      read(5,2)infile1
      print *,' '
      print *,' '
      print *, ' Enter the stage file name   '
      read(5,2)infile2
      print *,' '
      print *,' '
      print *,' Enter the output file name   '
      read(5,2)outfile
      print *,' '
      print *,' '
      print *,' Enter the watershed area:    '
      read(5,*)area
      print *,' '
      print *,' '
      print *,' Enter the bottom elevation of the stream:    '
      read(5,*)tmin
c
c  find rainless periods
      area = area / 2.589998
      days = area**0.2
      nom1 = int(days)
      nom2 = nom1/2
      
c
      call rainless(count4,infile1,infile2,nom1,nom2,
     +count5,count1,count15)
c
c
c
c  Find the rating curve for the flow values
c
      call regress(count4,a,b,tmin)
c
c  compute component
c
c
c
      call component(a,b,tmin,outfile,count5,count15)
 
c
  2   format(a)
      close(11)
      close(12)
      close(15)
      close(13)
      stop
      end
c
c
      subroutine rainless(count4,
     +infile1,infile2,nom1,nom2,count5,count1,count15)
c
c  this program examines precipitation and streamflow data for rainless
c  periods then matches it up to stage data
c
      implicit none
      real precip,flow,stage
      integer day1,mon1,year1,i,j,day,
     +month,year,count4,count5,count1,count15,nom1,nom2
      character*40 infile1,infile2
c
c
      open(1,file=infile1,status='old')
      open(2,file=infile2,status='old')
      open(11,status='scratch',access='direct',recl=16,
     +form='unformatted')
      open(15,status='scratch',access='direct',recl=20,
     +form='unformatted')
      open(12,status='scratch',access='direct',recl=16,
     +form='unformatted')
c
c
c      write stage data to a scratch file
c
          count5 = 0
          do j = 1,50000
	  read(2,*,end=21)year,month,day,stage
	  count5 = count5 + 1
	  write(12,rec=count5)year,month,day,stage
          enddo
   21     continue
	  close(2)
c
c
cc  find rainless periods
c
         count1 = 0
	 count15 = 0
         do i = 1,50000
c
         read(1,*,end=20)year1,mon1,day1,precip,flow
	  if(flow.gt.0.00) then
	      count15 = count15 + 1
              write(15,rec=count15)year1,mon1,day1,flow
	      else
	   endif
c
c
          if(precip.le.0.1.and.flow.gt.0.00) then
              count1 = count1 + 1
              write(11,rec=count1)year1,mon1,day1,flow
	      print *, 'processing rainless period:   ',count1
c
          else
	  endif
         enddo
   20 continue
      close(1)
c
c
c
c  select number of days after the rainless period begins
c
       call select(count1,count4,nom1,nom2,count5)
c
c
       return
       end
c      
c
c
      subroutine select(count1,count4,nom1,nom2,count5)
c
c    this program examines the stage-flow-precipitation output file
c    and creates a similar file for consistant periods after three
c    days or less.
c
      implicit none
      integer i,j,mon1,day1,year1,count1,
     +day_diff,start,end,count2,count4,nom1,nom2,count5,
     +year2,mon2,day2,k,year1_2,mon1_2,day1_2,year1_3,
     +mon1_3,day1_3
      real flow1,flow1_2,flow1_3,
     +stage
c
c
      open(unit=13,status='scratch',access='direct',
     +recl=8,form='unformatted')
c

      count2 = 0
      count4 = 0
      do i = 1,count1
c
      if(i.lt.count1) then
      read(11,rec=i)year1,mon1,day1,flow1
      read(11,rec=i+1)year1_2,mon1_2,day1_2,flow1_2
      day_diff = day1_2 - day1
c
         if(year1.eq.year1_2.and.mon1.eq.mon1_2.and.
     +   day_diff.eq.1) then
         count2 = count2 + 1
	 else
           if(i.gt.1) then
           read(11,rec=i)year1,mon1,day1,flow1
	   read(11,rec=i-1)year1_3,mon1_3,day1_3,flow1_3
	   day_diff = day1 - day1_3
	   if(year1.eq.year1_3.and.mon1.eq.
     +     mon1_3.and.day_diff.eq.1) then
           count2 = count2 + 1       
           endif
	   endif
	    if(count2.ge.nom1) then
c
c
	        start = i - count2 + nom2
		end = i
		  do j = start,end
		   read(11,rec=j)year1,mon1,day1,flow1
                    do k = 1,count5
		    read(12,rec=k)year2,mon2,day2,stage
                    if(year1.eq.year2.and.mon1.eq.mon2.
     +              and.day1.eq.day2) then		      
		    count4 = count4 + 1
                    write(13,rec=count4)stage,flow1
       print *, 'matched flow data to rainless period number:',
     +j, '  of    ',count1
                    else
		    endif
		    enddo
                   enddo
	     endif
	     count2 = 0
	     endif
      else
	   read(11,rec=i-1)year1_3,mon1_3,day1_3,flow1_3
           read(11,rec=i)year1,mon1,day1,flow1
	   day_diff = day1 - day1_3
             if(year1.eq.year1_3.and.mon1.eq.
     +       mon1_3.and.day_diff.eq.1) then
             count2 = count2 + 1
	           if(count2.ge.nom1) then
c
c  skip the first two days to allow overland runoff to occur
c
 	             start = i - count2 + nom2
		     end = i
		     do j = start,end
        	     read(11,rec=j)year1,mon1,day1,flow1
        	    do k = 1,count5
		    read(12,rec=k)year2,mon2,day2,stage
                    if(year1.eq.year2.and.mon1.eq.mon2.
     +              and.day1.eq.day2) then		      
		    count4 = count4 + 1
                    write(13,rec=count4)stage,flow1
		    print *, 'match stage number ',i
		    else
		    endif
		    enddo
                     enddo
	          endif
              endif
      endif
      enddo
  222 format(3i7,2x,2f10.2)
      return
      end
		  
c
c
      subroutine regress(n,a,b,tmin)
      implicit none
      real x,y,meany,meanx,a,b,denomenator,
     +sumx,sumy,numerator,
     +x_trans,y_trans,tmin,
     +sum_xbar_squared,sumx_ybar
      integer i,n,n2
c
c
c compute the regression
c
      sumx_ybar = 0.0
      sumx = 0.0
      sumy = 0.0
      sum_xbar_squared = 0.0
      n2 = 0
      do i = 1,n
c      
c  must first transform the data
c
      read(13,rec=i)x,y
c
      if(x.gt.tmin) then
      n2 = n2 + 1
      x_trans = log(x - tmin)
      y_trans = log(y)
c
      sumx = sumx + x_trans
      sumy = sumy + y_trans
      else
      endif
      enddo
c
c  if the stream bottom is higher in elevation than the water table
c  then no baseflows enters the stream 
c
c
      if(n2.eq.0) then
      print *, ' '
      print *, ' '
      print *, ' '
      print *, ' '
      print *, ' '
      print *, ' '
      print *, ' '
      print *, ' '
      print *, ' ***********  Attention::!!!!!!!!!!!!!!!!!!!!!!'
      print *, ' '
      print *, ' '
      print *, ' '
      print *, ' '
      print *, '    The stream bottom has a higher elevation than'
      print *, '    the groundwater table, thus baseflow separation'
      print *, '    did not occur!!!!!!!!!!!'
      print *, ' '
      print *, ' '
      print *, ' '
      print *, ' '
      print *, ' '
      print *, ' '
      stop
      else
      endif
      meanx = sumx/float(n2)
      meany = sumy/float(n2)
c
c
      do i = 1,n
      read(13,rec=i)x,y
      if(x.gt.tmin) then
      x_trans = log(x - tmin)
      y_trans = log(y)
      sumx_ybar = sumx_ybar + ((x_trans - meanx)*(y_trans-meany))
      sum_xbar_squared = sum_xbar_squared + ((x_trans - meanx)**2)
      else
      endif
      enddo
      numerator = sumx_ybar
      denomenator = sum_xbar_squared
c
      a = numerator/denomenator
      b = meany - (a*meanx)
            b = exp(b)

c
c  compute discharges for the available stages, and change b to
c  normal form--take the exponential
c
c
c
      return
      end            
c      
c
      subroutine component(a,b,tmin,outfile,count5,count15)
c
c  this program examines precipitation and streamflow data for rainless
c  periods then matches it up to stage data
c
      implicit none
      real flow,stage,a,b,tmin,rating,time,mon(12)
      integer day1,mon1,mon2,year1,day2,year2,i,j,time2,
     +count5,count11,count15
      character*40 outfile
c
      data mon/0.,31.,59.,90.,120.,151.,181.,212.,243.,273.,304.
     +,335./
c
      open(3,file=outfile,status='new')
c
c
c  read begining dates
c
      count11 = 0
      time2 = 0.
      do i = 1,count15
      read(15,rec=i)year1,mon1,day1,flow
         do j = 1,count5
	 read(12,rec=j)year2,mon2,day2,stage
c
	 if(year1.eq.year2.and.mon1.eq.mon2.
     +   and.day1.eq.day2) then
          if(stage.gt.tmin) then
	  count11 = count11 + 1
	  time = float(year1) + ((mon(mon1) + float(day1))/365.)
c	  time = time - 1900.
	   rating = b*((stage - tmin)**a)
	   time2 = time2 + 1.
	   write(3,4)year1,mon1,day1,TIME,stage,flow,rating
      print *, 'computing baseflow #',count11     ,' of a possible',
     +count5
          goto 21
          else
	  endif
	  else
          endif
         enddo
   21  continue
      enddo
c
   4  format(2x,3i5,4f10.3)
      return
      end
