c  Uses the Gamma probability density function as a transfer function
c  to describe a time-series response function produced by a specified
c  time-series driving function.
c
c  e.g. Route effective infiltration (i.e. flux leaving the root zone --
c  below this point flux is assumed to always be downward) through the
c  unsaturated zone from the bottom of the root zone to the water
c  table.
c
c  NOTE: Gamma function overflows for arguments > 171 on flalsws026.
c        To avoid overflow, natural logs are used for all gamma
c        function and gamma pdf calculations -- gamma pdf untransformed
c        before used in other calculations.
c
c  Andrew M. O'Reilly
c
c  Beta Versions:
c  ==============
c  v1.0 11/30/2001 - Original code
c  v1.1 03/05/2002 - Corrected calculation of area under driving and
c    response function curves to use correct time period:
c      drf: (initial time) thru (final time)
c      rsf: (initial time) thru (final time + gamma pdf "memory")
c    For purposes of calculating rsf when memory extends to before
c    initial time for drf, drf assumed = 0 for times < initial time
c       03/08/2002 - Corrected error that can occur when subtracting
c    two "identical" floating point numbers.  This can occur when
c    calculating (t-tau) if the fractional part of each value is equal,
c      e.g. t=64.1 and tau=3.1
c           (t-tau)=61, which truncates to 61
c           however...if t=64.09999999999999
c           then (t-tau)=60.99999999999999, which truncates to 60
c  v2.0 03/11/2002 - Incorporated initial time lag (tlagi) in response
c    function, i.e. rsf does not start to respond to drf until t>tlagi
c       03/19/2002 - Corrected error in calculation of gamma cdf.
c    Revised calculation of 2nd term of gpdfln in gampdf to avoid
c    overflow by converting log(x**a) to a*log(x)
c       03/20/2002 - Corrected convol so drf = 0 when t>final time as
c    as well as when t<initial time
c       04/04/2002 - Corrected calculation of area under driving
c    function curve for drf(itf)
c  v3.0 07/15/2002 - converted to subroutine for MAIN water-balance/
c    transfer-function program; moved date and jdaycm subroutines 
c    into separate file
c       07/16/2002 - removed interactive input--moved to input file
c    read by MAIN program
c       07/19/2002 - added new output file rffil2 which contains 
c    response function for times concurrent with driving function
c       10/11/2002 - changed gcdfm from 0.99 to 0.95; changed mem
c    (memory of gamma pdf) to 2*X, where X is defined by 
c    P(tau <= X) = gcdfm; this is a crude fix so that extremely
c    small values of dtau (<0.01) are not required for small values
c    of xn and xk	(resulting in a gamma pdf with a sharp peak)
c  v4.0 10/23/2002 - modified convol to use dtau/2 as initial tau
c    (therefore, if dtau=0.1, tau=0.05, 0.15, 0.25, etc) this is more
c    accurate because the gamma pdf value at center of dtau interval is
c    used rather than at end of dtau interval
c    - appropriate dtau is now determined based on RMSE between actual
c    and estimated Gamma pdf slope; dtau initially = dt, but is 
c    reduced by a factor (dtaurf, currently 0.5) until RMSE as a percent
c    of average absolute value of actual slope is less than a specified
c    amount (praer, currently 5%); supercedes changes made 10/11/2002
c    - tlagi does NOT have to be an even multiple of dt or dtau
c       10/28/2002 - revised algorithm to select appropriate dtau
c    - modified convol to use estimated Gamma pdf at first tau value
c    (tau=dtau/2) based on prescribed Gamma pdf at tau=0;  this is
c    necessary when xn<=1 because Gamma pdf --> infinity as tau --> 0
c    therefore, dtau would have to be infintesimally small in order
c    to obtain a gamma pdf with an area of ~1
c       11/07/2002 - modified convol to use tau values such that
c    taumtl (tau-tlagi) values are the same as the tau values used to
c    compute Gamma pdf memory; this previously would not be the case if
c    tlagi was not an even multiple of dtau; now the mass balance check
c    on rsf should be exactly the same percentage as the Gamma cdf at
c    the selected memory
c  v4.1 11/08/2002 - response function averaged over driving function
c    time step is now computed and written to rffil2
c       11/15/2002 - corrected calculation of nttli to account for three
c    cases (1) tlagi < dtau/2, (2) tlagi is an uneven multiple of dtau,
c    and (3) tlagi is an even multiple of dtau
c       12/10/2002 - dtdrf now specified in input file
c       12/17/2002 - change gcdfer from 0.001 to 0.01 (for xn=0.6 and
c    xk=100: this only made ~0.04% difference in gamma cdf values;
c    changing gcdfer from 0.01 to 0.1 makes ~0.4% difference in gamma
c    cdf values)
c       12/31/2002 - modified to write gamma pdf to gam.csv for final
c    dtau value ONLY
c       03/11/2003 - modified to write gamma pdf to gam.csv ONLY if
c    dtau gt 0.01 to avoid creating giant file
c  v5.0 03/17/2003 - modified to work with v5.0 of wbal; modified for
c    drf to be "backward discretized," ie, data for unit event time
c    t(n) applies for t(n-1)<times<=t(n), previously applied for
c    t(n)<=times<t(n+1)
c
c                        |----------|----------|
c    unit event time:   n-1         n         n+1
c							  
c                        \----+----/\-----+----/
c                             | 	        |
c    discretization:     "backward"    "forward" 
c
c       - the length of the unit driving event should be small relative   
c    to the length of the corresponding response event so that it       
c    approximates an instantaneous pulse, e.g. a unit event length of   
c    1 day is sufficient for a unit effective infiltration event that   
c    produces a recharge event that spans 30 days, in contrast, a       
c    recharge event that spans 3 days probably should be simulated by   
c    a unit effective infiltration event length of 1 hr.; based on this,
c    dtdrf will always be 1, e.g. representing 7 days, 1 day,           
c    1 1/10th day, 1 hr, etc.; the real time unit conversion factor,    
c    truc, and offset, tri, are used to convert unit event times to     
c    a real time of the user's choosing
c       03/21/2003 - modified to write response function averaged over
c    user-specified time step (dtravg) to file rffil2, the averaged
c    response function has the same time units as dtravg
c       03/24/2003 - corrected error in calc dgpdfe and gamcdf when
c    writting gam.csv (only affected gam.csv, calcs correct elsewhere)
c       04/02/2003 - change praer to 500 (was 5), REMINDER: pra=RMSE of
c    actual to estimated gamma pdf slope values as a percentage of 
c    average absolute slope -- for convolution, use largest dtau for 
c    which pra<praer; this is OK based on: for n=0.2035, k=55.04 hrs, &
c    tlagi=18.37 hrs, average diff between rsf with pra=318 vs pra=1.4
c    was 0.006% (max diff was 0.36%)
c       06/24/2003 - corrected error in ittau function that incorrectly
c    computed ittau when tau>t, now ittau=-99 when tau>t and these 
c    values ignored by convol since this is before the start of drf
c       01/05/2004 - added computation of memory+tlagi
c  v6.0 01/08/2004 - minor changes to screen output since dtdrf changed;
c    modified ittau function to calc (t-tau)/dtdrf -- this is required
c    so that mtmtau in convol is the correct index for drf; corrected
c    calculation of averaged rsf -- divide by ntsra (not ntsrd)
c       01/21/2004 - change equation to calc tr for writing to files --
c    old equation provided corrected value of tr but was cumbersome;
c    changed what is written to rsffil and rffil2; added error check for
c    rsfa/drfa mass-balance discrepancy
c
c  Public Release Versions:
c  ========================
c  v1.0 05/28/2004 - minor changes to format of Main Output File
c        
c
c
      subroutine trfn(rsffil,rffil2,xn,xk,tlagi,dtdrf,dt,truc,tri,
     &dtravg,ntsd,drf)
 	double precision dtau,dtaurf,dt,dtdrf,dterr,tau,t,truc,tri,dtravg,
     &tlagi,gammln,xn,xk,tgtf
	character*50 rsffil,rffil2
      character*25 version
	dimension drf(ntsd)
	allocatable rsf(:),t(:)
c
	version='Version 1.0 -- 05/28/2004'
	write(*,*)
	write(*,*) '**************************************************'
      write(*,*) '***** Transfer-Function module of WBTF model *****'
	write(*,*) '*********** ',version,' ************'
	write(*,*) '**************************************************'
	write(*,*)
c
c  Open files
c
      open(unit=10,file='gam.csv',status='unknown')
	open(unit=30,file=rsffil,status='unknown')
	open(unit=40,file=rffil2,status='unknown')
c
      write(*,*) 'Gamma pdf shape parameter (N) =',xn
	write(*,*) 'Gamma pdf initial time lag (TAUI) =',tlagi
      write(*,*) 'Gamma pdf scale parameter (K) =',xk
	write(*,*) 'Time step for driving function EI (DTPE) =',dtdrf
	write(*,*) 'Time step for response function RCHIN (DTU) =',dt
	write(*,*) 'Time step for averaging period (DTRAVG) =',dtravg
c
c  Check time step
c
	dterr=dtdrf/dt-dint(dtdrf/dt)
	if (dterr.gt.1.d-14) then
	   write(*,*) 'ERROR--DTU must be an even multiple of',dtdrf
	   write(30,*) 'ERROR--DTU must be an even multiple of',dtdrf
	   write(40,*) 'ERROR--DTU must be an even multiple of',dtdrf
	   stop
	endif
	write(*,*)
c
c  Allocate arrays and fill time array
c
      ntsra=idnint(dtravg/truc/dt)  !# rsf time steps per averging period
	write(*,*) 'Number of unit-event time steps (DTU) per averaging pe
     &riod =',ntsra
	ntsrd=idnint(dtdrf/dt)  !# time steps in rsf per unit drf event
	write(*,*) 'Number of unit-event time steps (DTU) per driving-func
     &tion EI time step (DTPE/DTU) =',ntsrd
	ntsr=ntsd*ntsrd   !total # rsf time steps spanned by drf
	write(*,*) 'Number of unit-even time steps (DTU) spanned by drivin
     &g function EI =',ntsr
	write(*,*)
	allocate (rsf(ntsr),t(ntsr))
	do 10 m=1,ntsr
	   t(m)=dble(float(m))*dt	  !times to calc rsf
	   rsf(m)=-999.
   10 continue
c
c  Calculate Gamma pdf for specified parameters; determine the
c  memory, mem, of the pdf (i.e. probability(tau<=mem) > gcdfm) by
c  calculating the gamma cdf; and determine the appropriate dtau
c						  
	gcdfm=0.99    !value of gamcdf used to calc memory
	gcdff=0.999   !final value of gamcdf, used to stop calc gampdf
	if (xn.gt.1.d0) then
	   gcdfer=0.01  !xn>1 calc slope error stats when gamcdf>gcdfer
	   dger=1.e30    !xn>1 calc slope error stats when dgpdf<dger
	else
	   gcdfer=0.01  !xn<=1 calc slope error stats when gamcdf>gcdfer
	   dger=100.     !xn<=1 calc slope error stats when dgpdf<dger
	endif
	praer=500.     !RMSE as % avg abs slope, use dtau where pra<praer
	dtaurf=0.5d0   !reduction factor to calc smaller values of dtau
	dtauap=1.d-2   !approximate area 0 to dtau when dtau<dtauap
	dtaumi=1.d-4   !minimum dtau
	dtau=dt
	gampza=0.  !add to gammpz to approx Gamma pdf at zero for xn<=1
	gpdfav=xn*xk   !mean of Gamma pdf
c
	write(*,*) '+++++++++++++++++++++++++++++++++++++++++++++++++++++'
	write(*,*) 'The following information concerns the selection of th
     &e appropriate value for dtau...'
	write(*,*) 'Too large a dtau will result in discretization error.' 
   90	write(*,'(18h *****************,/,15h Trying dtau = ,g10.4)') dtau
      i=0
	j=0
	taum=-1.
	tauslp=0.   !calc slope error stats for only for tau>tauslp
	gammpd=999.
c
	tau=dtau  !est. gampdf at tau=0 by extrapolating slope at tau=dtau
	dgpdf=dexp(-tau/xk-xn*dlog(xk)-gammln(xn))	     !actual slope
     &      *((xn-1.d0)*tau**(xn-2.d0)-tau**(xn-1.d0)/xk)
	if (xn.gt.1.d0) then
	   gammpz=gampdf(xn,xk,tau)-dgpdf*dtau   !~gampdf at tau=0
	   if (gammpz.lt.0) gammpz=0.
	   gammpo=gammpz
	else
	   gammpz=abs(dgpdf)*dtau+gampdf(xn,xk,tau)+gampza
	   gammpo=gammpz
	endif
	write(*,*) 'Slope of Gamma pdf at tau =',tau,' is',dgpdf
	write(*,*) 'Prescribed Gamma pdf at tau = 0 is',gammpz
c
	gamcdf=0.
	resmax=0.
	sres=0.
	ssres=0.
	sadg=0.
  100	if (gamcdf.gt.gcdff) then
	   goto 110
	elseif (gamcdf.gt.gcdfm.and.gammpd.lt.1.e-6) then
	   write(*,*) 'Gamma cdf > ',gcdfm,' but dtau probably too large t
     &o calculate Gamma cdf > ',gcdff
	   goto 110
	elseif (tau.gt.gpdfav.and.gammpd.lt.1.e-6.and.gamcdf.le.gcdfm)
     &then   !exceed mean pdf, cdf inc. slowly, still < memory criterion
	   write(*,'(44h Cannot reach specified memory criterion of ,f4.2,
     &         33h for Gamma cdf using current dtau)') gcdfm
	   write(*,'(9h P(tau <=,f9.3,4h) = ,f7.5)') tau,gamcdf
	   if (dtau.lt.dtauap) then
	      write(*,'(24h WARNING--dtau must be <,g10.4)') dtauap
	      write(*,'(42h Estimate area from 0 to dtau using gampza)')
	      gcdfa=gcdfm*1.01-gamcdf   !add'l area needed to meet gcdfm
	      gampza=2.*gcdfa/dtau   !add to gammpz and try again
	      write(*,*) ' gcdfa =',gcdfa,' gampza =',gampza
	   else
	      dtau=dtau*dtaurf
	   endif
	   goto 90
	else
	   i=i+1   !number of tau steps used to compute Gamma pdf
	   tau=dble(float(i))*dtau
	   gammpd=gampdf(xn,xk,tau)
         if (xn.le.1.d0.and.i.eq.1) then
	      gamcdf=dtau*(gammpo+gammpd)/2.
	   else
	      gamcdf=gamcdf+dtau*gampdf(xn,xk,tau-dtau/2.d0)
	   endif
	   dgpdfe=(gampdf(xn,xk,tau+dtau)-gammpo)/(2.*dtau)	 !est. slope
	   dgpdf=dexp(-tau/xk-xn*dlog(xk)-gammln(xn))	     !actual slope
     &         *((xn-1.d0)*tau**(xn-2.d0)-tau**(xn-1.d0)/xk)
         if (abs(dgpdf).gt.dger.or.gamcdf.lt.gcdfer) then
	      res=0.
	      sadg=0.
	      tauslp=tau   !calc slope error stats for only for tau>tauslp
	      gcdfsl=gamcdf  !gamma cdf at tau=tauslp
	      j=j+1	  !number of tau steps to exclude from error analysis
	   else
	      res=dgpdf-dgpdfe   !slope residual (actual-estimated)
	      sadg=abs(dgpdf)+sadg  !sum of absolute values of slope
	   endif
	   if (abs(res).gt.abs(resmax)) resmax=res   !maximum residual
	   sres=res+sres	  !sum of residuals
	   ssres=res**2.+ssres	!sum-of-squared residuals
	   gammpo=gammpd
	   if (gamcdf.gt.gcdfm.and.taum.lt.0.) then
	      taum=tau
	      rmse=(ssres/(i-j))**0.5
	      aadg=sadg/(i-j)    !average absolute value of slope
	      pra=rmse/aadg*100.
	      avgres=sres/(i-j)	  !average residual
	      if (pra.ge.praer) then
c	         write(*,*) ' RMSE =',rmse
c	         write(*,*) ' Avg abs slope =',aadg
c 	         write(*,*) ' i, j =', i,j
	         write(*,'(30h RMSE as % of avg abs slope = ,g10.4,
     &		       12h for dtau = ,g10.4,10h and tau >,g10.4)') pra,
     &               dtau,tauslp
	         write(*,'(9h P(tau <=,f9.3,4h) = ,f7.5)') taum,gamcdf
	         if (dtau.lt.dtaumi) then
	            write(*,'(25h WARNING2--dtau must be <,g10.4)') dtaumi
	            stop
	         endif
	         gampza=0.
	         dtau=dtau*dtaurf
	         goto 90
	      else
	         write(*,'(74h Statistics concerning actual and estimated 
     &slope of Gamma pdf for dtau = ,g10.4,/,9h and for ,g10.3,
     &9h < tau <=,f9.3)') dtau,tauslp,taum
	         write(*,'(9h P(tau <=,f9.3,4h) = ,f7.5)') tauslp,gcdfsl
	         write(*,'(9h P(tau <=,f9.3,4h) = ,f7.5)') taum,gamcdf
	         write(*,'(26h Root mean square error = ,g10.4)') rmse
	         write(*,'(26h Average absolute slope = ,g10.4)') aadg
	         write(*,'(30h RMSE as % of avg abs slope = ,g10.4)') pra
	         write(*,'(20h Average residual = ,g10.4)') avgres
	         write(*,'(20h Maximum residual = ,g10.4)') resmax
	      endif
	   endif
c	   write(10,'(8(f11.6,3h  ,))')
c	&        xn,xk,tau,gammpd,gamcdf,dgpdfe,dgpdf,res
	endif
	goto 100
c
  110	write(*,*) '+++++++++++++++++++++++++++++++++++++++++++++++++++++'
	write(*,*)
	write(*,*) '+++++++++++++++++++++++++++++++++++++++++++++++++++++'
	write(*,*) 'The following information concerns the memory of the g
     &amma pdf...'
	write(*,'(9h P(tau <=,f9.3,4h) = ,f7.5)') tau,gamcdf
	mem=nint(taum+0.5)
	xmem2=taum+tlagi 
	write(*,'(48h Memory of specified gamma pdf excluding TAUI = ,
     &i5)') mem
      write(*,'(61h Memory is calculated as X, where X defined by P(tau 
     &<= X) = ,f7.5,/,37h and X rounded up to nearest integer.)') gcdfm
	write(*,'(62h Memory (not rounded) of specified gamma pdf includin
     &g TAUI = ,f7.1)') xmem2
c
c  Write gamma pdf to file using final dtau value
c
	if (dtau.gt.1.d-2) then
   	write(10,'(49h xn , xk , tau , gampdf , gamcdf , dgpdfe , dgpdf)')
	gammpo=gammpz
	gamcdf=0.
	do 120 i=1,nint(mem/dtau)
	   tau=dble(float(i))*dtau
	   gammpd=gampdf(xn,xk,tau)
         if (xn.le.1.d0.and.i.eq.1) then
	      gamcdf=dtau*(gammpz+gammpd)/2.
	   else
	      gamcdf=gamcdf+dtau*gampdf(xn,xk,tau-dtau/2.d0)
	   endif
	   dgpdfe=(gampdf(xn,xk,tau+dtau)-gammpo)/(2.*dtau)  !est. slope
	   dgpdf=dexp(-tau/xk-xn*dlog(xk)-gammln(xn))	     !actual slope
     &         *((xn-1.d0)*tau**(xn-2.d0)-tau**(xn-1.d0)/xk)
	   write(10,'(7(f11.6,3h  ,))')
     &        xn,xk,tau,gammpd,gamcdf,dgpdfe,dgpdf
	   gammpo=gammpd
  120	continue	   
	write(*,*) 'Gamma pdf is in file gam.csv'
	write(*,*) '+++++++++++++++++++++++++++++++++++++++++++++++++++++'
	write(*,*)
	endif
c
c  Perform convolution of driving function with transfer function to
c  calculate response function
c
	write(*,*) '+++++++++++++++++++++++++++++++++++++++++++++++++++++'
	write(*,*) 'The following information concerns the convolution of 
     &EI with the gamma pdf and a mass balance check (comparison of area
     &s under EI and RCHIN curves)...'
	ntmem=nint(mem/dtau)   !# tau steps in pdf memory
	write(*,*) 'Number of tau steps in gamma pdf memory =',ntmem
	tmp=tlagi/dtau-dint(tlagi/dtau)
      if (nint(tlagi/dtau).eq.0) then
	   nttli=0					   !# tau steps in initial time lag
	elseif (tmp.gt.1.e-6) then     !tlagi uneven multiple of dtau
c	   write(*,*) 'tlagi/dtau-dint(tlagi/dtau) =',tmp
	   nttli=nint(tlagi/dtau+0.5)  !round up for extra step required
	else					    !tlagi even multiple of dtau	   		
         nttli=nint(tlagi/dtau)  	!don't add extra step
	endif
	write(*,*) 'Number of tau steps in initial time lag =',nttli
	ntmeml=ntmem+nttli
	do 130 m=1,ntsr
	      call convol(ntmeml,nttli,dtau,dtdrf,t(m),drf,ntsd,xn,xk,
     &                  tlagi,gammpz,rsf(m))
c	      write(30,*) ' m, t(m), rsf(m) =',m,t(m),rsf(m)
  130	continue
	write(*,*) ' ntsr =',ntsr,'  ;  m =',m
c
c  Calc and compare area under curves to check for conservation of mass.
c
      rsfa=sngl(dt)*rsf(1)/2.   !area under rsf curve
	do 160 m=2,ntsr
	   rsfa=sngl(dt)*(rsf(m)+rsf(m-1))/2.+rsfa
  160 continue
	write(*,*) 'm =',m,'  ;  Area under rsf for t<=tf:',rsfa
	rsfa1=rsfa
	tmpold=rsf(ntsr)
	mtfml=ntsr+nint(mem/dt)+nint(tlagi/dt)  !rsf time steps
	do 180 m=ntsr+1,mtfml  !calc rsf area for: tf<times<tf+mem+tlagi
	   tgtf=dble(float(m))*dt   !tgtf=t(m) for times > tf
	   call convol(ntmeml,nttli,dtau,dtdrf,tgtf,drf,ntsd,xn,xk,
     &               tlagi,gammpz,tmp)
	   rsfa=sngl(dt)*(tmp+tmpold)/2.+rsfa
	   tmpold=tmp
c	   tgtf=dble(float(m+1))*dt
  180 continue
	tmp=rsfa-rsfa1 
	write(*,*) 'm =',m,'  ;  Area under rsf for t>tf:',tmp
c 
	drfa=(drf(1))/2.   !area under drf curve
	do 210 m=2,ntsd
	   drfa=sngl(dtdrf)*(drf(m)+drf(m-1))/2.+drfa
c	   write(30,'(i4,2x,3(g11.5,2x))') i,drf(i-1),drf(i),drfa
  210	continue
	drfa=(drf(ntsd))/2.+drfa
c
  211	write(*,*) 'Area under the driving function EI curve =',drfa
	write(*,*) 'Area under the response function RCHIN curve =',rsfa
	tmp=rsfa/drfa*100.
	write(*,*) 'Area under the response function curve RCHIN is',tmp,'
     & percent that under the driving function curve EI.'
	if (tmp/100..lt.gcdfm.or.tmp/100..gt.1.) then
	   write(*,*) 'ERROR -- The ratio of the area under the response f
     &unction curve to the area under the driving function curve should 
     &be >= ',gcdfm,' or <= 1',gcdfm
	   stop
	endif
	write(*,*) '+++++++++++++++++++++++++++++++++++++++++++++++++++++'
c
c  Write rsf to file rsffil, rsf has time same time units as unit event
c
c	write(30,'(32hReal-Time, DrF-Time, DrF, RsF:n=,f6.2,3h k=,f6.2,
c	&      7h tlagi=,f6.2,8h memory=,i4,17h, unit-event-time)')
c    &      xn,xk,tlagi,mem
	write(30,'(21hTime, EI, Rch-inst:n=,f6.2,6h TAUi=,f6.2,
     &      3h k=,f6.2,8h TAUmem=,f6.1)')
     &      xn,tlagi,xk,xmem2
	m=1
c	tr1=tri-dt*dble(ntsrd-1)*truc
	tr1=t(m)*truc+tri-dtdrf
	tr=tr1
	do 220 i=1,ntsd
	   tdrf=float(m)/ntsrd*dtdrf
	   do 230 k=1,ntsrd
	     if (rsf(m).gt.0.and.rsf(m-1).gt.0.) then
	        ipowm=ifix(log10(rsf(m)))
	        ipowm1=ifix(log10(rsf(m-1)))
	        dmant=rsf(m)/10.**ipowm-rsf(m-1)/10.**ipowm1  !mantissa difference
	     else
	        dmant=999.
	     endif
	     if (dmant.lt.1.e-6.and.dmant.gt.-1.e-6.and.
     &	     ipowm.eq.ipowm1) then  !mark records if rsf(m)=rsf(m-1)
c	        write(30,'(f9.3,3h , ,f6.0,2(2h, ,g11.5),3h , ,f9.3,3h , ,
c	&              g12.6)') tr,tdrf,drf(i),rsf(m),t(m),dmant
	        write(30,'(f9.3,2(2h, ,g11.5),3h , ,g12.6,
     &             45h  ERROR--two consecutive rsf values identical)')
     &             tr,drf(i),rsf(m),dmant
	     else
c	        write(30,'(f9.3,3h , ,f6.0,2(2h, ,g11.5),3h , ,f9.3)') 
c	&	         tr,tdrf,drf(i),rsf(m),t(m)
	        write(30,'(f9.3,2(2h, ,g11.5))') tr,drf(i),rsf(m)
	     endif
	     tr=t(m)*truc+tr1
	     m=m+1
  230	   continue
  220	continue
      write(*,*)
      write(*,*) 'Instantaneous response function is in file: ',rsffil
c
c  Write output to rffil2 for data averaged over dtravg, averaged rsf
c  has same time units as dtravg
c
	write(40,'(16hTime, Rch-avg:n=,f6.2,6h TAUi=,f6.2,
     &      3h k=,f6.2,8h TAUmem=,f6.1,10h, T-s, T-e)')
     &      xn,tlagi,xk,xmem2
	m=1
	tmp=0.
c	tr1=tri-dt*dble(ntsrd)*truc
	tr1=t(m)*truc-dt+tri-dtdrf
	tra1=tr1   !start time of averaging period
c	write(*,*) 'ntsr/ntsra =',ntsr/ntsra
	do 240 i=1,ntsr/ntsra
	   do 250 k=1,ntsra
            tmp=rsf(m)+tmp
	      m=m+1
  250	   continue
	   tra2=t(m-1)*truc+tr1   !end time of averaging period
	   tra=(tra1+tra2)/2.
	   tmp=tmp/ntsra/dtravg
c	   write(40,'(f9.3,4(3h , ,g13.7))') 
c	&	     tra,tmp,tra1,tra2,t(m-1)
         write(40,'(f9.3,1(3h , ,g13.7),2(3h , ,f9.3))') 
     &	     tra,tmp,tra1,tra2
	   tra1=tra2
	   tmp=0.
  240	continue
c	write(*,*) 'The following time step was used to compute the averag
c	&ed response function: ',dtravg
	write(*,*) 'Averaged response function is in file: ',rffil2
c
c  The END!!
c
   99	close(10)
      close(30)
	close(40)
c
	return
      end
c
c************************************************************
c  Performs convolution.
c************************************************************
      subroutine convol(ntmeml,nttli,dtau,dtdrf,tt,drf,ntsd,xn,xk,
     &tlagi,gammpz,tmp)
 	double precision dtau,dtdrf,tau,tt,xn,xk,taumtl,tlagi,tau2,taui
	dimension drf(ntsd)
c
c      taui=dmod(tlagi,dtau)/2.d0
      taui=(tlagi/dtau-dint(tlagi/dtau))*dtau/2.d0
      tmp=tlagi/dtau-dint(tlagi/dtau)
	if (tmp.le.1.d-6) then
	   taui=dtau/2.d0
	endif
      tmp=0.
	do 150 j=1,ntmeml
	   if (j.eq.1) then
	      tau=taui
	   elseif (j.eq.2) then
	      tau2=2.d0*taui+dtau/2.d0
	      tau=tau2
	   else
            tau=tau2+dble(float(j-2))*dtau
	   endif
       mtmtau=ittau(tt,tau,dtdrf)+1  !drf time step for t-tau
	   if (mtmtau.ge.1.and.mtmtau.le.ntsd) then
	      taumtl=tau-tlagi
 	      if (xn.le.1.d0.and.j.eq.nttli+1) then  !adjust for gammpz
			 g=gampdf(xn,xk,dtau) 
               tmp=drf(mtmtau)*(gammpz+g)/2.*dtau+tmp
c	         ttmp=(gammpz+g)/2.*dtau
c			 write(30,*) j,tau,taumtl,ttmp
	      else
	         tmp=drf(mtmtau)*gampdf(xn,xk,taumtl)*dtau+tmp
	      endif
c	      write(30,'(5(g11.4),2i5,2(g11.4))') tt,tau,taumtl,error,tol,
c     &           ipower,itmtau,drf(itmtau),tmp
	   endif  !skip calc because drf=0 for t<ti & t>tf
  150	continue
      goto 999
c
  999	return
	end
c
c************************************************************
c  Calculates (t-tau), including correction for floating
c  point precision.
c************************************************************
	integer function ittau(tt,tau,dtdrf)
c
	double precision tt,tau,dtdrf,error,tol
c
	error=(tt-dint(tt))-(tau-dint(tau))
	if (tt.gt.tau) then
	   ipower=idint(dlog10(tt))
	else
	   ipower=idint(dlog10(tau))
	endif
	if (ipower.ge.0) then
	   tol=-1.d-15*1.d1**ipower   !prec decreases as value increases
	else
	   tol=-1.d-15   !max precision for double precision numbers
	endif
	if (tau.gt.tt) then
	   ittau=-99   !transfer function extends past start of drf
	elseif (error.lt.0.d0.and.error.gt.tol) then
	   ittau=idint((tt-tau-tol)/dtdrf)   !correct float pt precision
	else
	   ittau=idint((tt-tau)/dtdrf)
	endif
c
	return
	end
c
c************************************************************
c  Calculates the Gamma probability density function -- log
c  transformed to avoid overflow.
c************************************************************
	real function gampdf(xn,xk,ttau)
c
	double precision xn,xk,ttau,gammln,gpdfln
c
      if (ttau.le.1.d-15) then
	   gampdf=0.  !Gamma pdf undefined at ttau<=0
	else
	   gpdfln=(-ttau/xk)+(xn-1.d0)*dlog(ttau/xk)-dlog(xk)-gammln(xn)
	   gampdf=dexp(gpdfln)
	endif
c
	return
	end