ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*sbr detrend2
* Use in conjunction with subroutine rcurve.f to perform quick and dirty regional curve
* standardization 
* Writen by ERC and only slighly modified by pjk.
*
* jdt detrending option
* ndt = alpha freidman smoother, lowess smoother
*
      subroutine detrend2(x,f,n1,n,eq,jdt,mxy,ndt)			!Hakan
c
         integer n1,n,jdt,mxy,ndt
	  real*4  x(mxy),f(mxy),eq(7)
      parameter (i1=1,i2=2,i3=3,i4=4,i5=5,i6=6,i10=10,i21=21,iu11=11)
      parameter (nemax=2500,ndmax=i5,ndmaxp1=i6,r0=0e0,r1=1e0,nmax=5000)
      real*4 pam(ndmaxp1,ndmax),yam(ndmaxp1),v(ndmax)
      real*4 span,alpha,w(nmax),sc(nmax,7)
      real*4 wk2(nmax),wk3(nmax)    !dimention must equal mxy
      real*4 rdt
c      
c
      do i=1,7
        eq(i)=0.0
      enddo
      ls=0
      ndif=0
      if(jdt.eq.-5)then
        do i=1,n
          wk2(i)=float(i)
          wk3(i)=x(i)
          w(i)=1.0
        enddo
        span=0.0
        if(n.lt.40)span=0.3
        iper=1
        alpha=float(ndt)
        call supsmu(n,wk2,wk3,w,iper,span,alpha,f,sc)
      elseif(jdt.eq.-4)then
        do i=1,n
          wk2(i)=x(i)
        enddo
        alpha=float(ndt)/100.
c        rdt=float(ndt)
        call lowess(mxy,wk2,f,n,real(alpha))
      elseif(jdt.eq.-3)then
        do i=1,n
          wk2(i)=x(i)
        enddo
c        do i=2,n
c          wk2(i-1)=x(i)
c        enddo
        call msmoor(wk2,f,wk3,n,ndt,iter,ier,mxy)
      elseif(jdt.eq.-1)then
        do i=2,n
          wk2(i-1)=x(i)
        enddo
        call xmean(x(2),n-1,xm1)
        do i=2,n
          wk3(i-1)=x(i)-x(i-1)						                !first-difference
        enddo
        x(1)=-9.99
        n=n-1
        do i=1,n
          x(i+1)=wk3(i)
        enddo
c        call xmean(x,n,xm)
c        sum=0.0
c        do i=1,n
c          sum=sum+xm
c          f(i)=sum
c        enddo
c        call xmean(f,n,xm2)
c        do i=1,n
c          f(i)=f(i)-xm2+xm1
c        enddo
        ls=n/3
        call spline(n,wk2,f,p,ls)
        do i=1,n
          wk2(i)=wk2(i)/f(i)
        enddo
        call moment(x(2),n,ave1,adev,sdev1,var,skew,curt,sens,ac)
        call moment(wk2,n,ave2,adev,sdev2,var,skew,curt,sens,ac)
        do i=1,n
          x(i+1)=(x(i+1)-ave1)/sdev1
          x(i+1)=x(i+1)*sdev2+ave2
        enddo
      elseif(jdt.eq.1.or.jdt.eq.2)then
        do i=1,7
          eq(i)=0.0
        enddo
c        call curve(n,x,f,eq)				!negative exponential curve
        call curve(n1,n,x,f,eq)				!Hakan
        if(eq(2).le.0.or.eq(5).lt.0)then
          do i=1,7
            eq(i)=0.0
          enddo
c          call trend(n,x,f,eq(3),eq(5))
          call trend(n1,n,x,f,eq(3),eq(5))		!Hakan
          if(jdt.eq.1)then
            jdt=4
          elseif(jdt.eq.2.and.eq(3).le.0)then
            jdt=5
          else
            do i=1,7
              eq(i)=0.0
            enddo
            jdt=5
            call xmean(x,n,eq(5))
c            do i=1,n
            do i=1,n1
              f(i)=eq(5)					!Hakan
            enddo
          endif
        endif
      elseif(jdt.eq.3)then
        iopt=1
        do i=1,n
          wk2(i)=i
        enddo
	 nd=iopt+i2
        pam(1,1)=0.0
        pam(1,2)=0.0
        pam(1,3)=0.0
        pam(1,4)=0.0
        pam(1,5)=0.0
c        call gamoeba(pam,yam,v,ndmaxp1,ndmax,nd,iter,iopt,nemax,n,wk2,x,f,sumsq)
        call gamoeba(pam,yam,v,ndmaxp1,ndmax,nd,iter,iopt,nemax,n1,n,wk2,x,f,sumsq)		!Hakan
        slop=pam(i1,i2)
        eq(1)=pam(i1,i1)
        eq(2)=abs(pam(i1,i2))
        eq(3)=r0
        eq(5)=pam(i1,i3)
        eq(6)=r0
        if(slop.ge.0)then
          jdt=4
          do i=1,7
            eq(i)=0.0
          enddo
c          call trend(n,x,f,eq(3),eq(5))
          call trend(n1,n,x,f,eq(3),eq(5))				!Hakan
        endif
      elseif(jdt.eq.4.or.jdt.eq.5)then
c        call trend(n,x,f,eq(3),eq(5))				        !linear regression
        call trend(n1,n,x,f,eq(3),eq(5))				 	!Hakan
        if(jdt.eq.5.and.eq(3).gt.0)then
          jdt=6
          eq(3)=0.0
          eq(5)=0.0
          call xmean(x,n,eq(5))						            !arithmetic mean
c          do i=1,n
          do i=1,n1								!Hakan
            f(i)=eq(5)
          enddo 
        endif       
      elseif(jdt.eq.6)then
        call xmean(x,n,eq(5))						              !arithmetic mean
c        do i=1,n
        do i=1,n1								!Hakan
          f(i)=eq(5)
        enddo
      elseif(jdt.eq.7)then
        do i=1,n
          wk2(i)=i
        enddo
c        call hughdi(wk2,x,n,f,eq,flag)           !hugershoff growth curve
        call hughdi(wk2,x,n1,n,f,eq,flag)         !Hakan
        do i=1,n
          wk2(i)=i
          wk3(i)=f(i)
        enddo
        iopt=3
	       nd=iopt+i2
        pam(1,1)=eq(1)
        pam(1,2)=eq(2)
        pam(1,3)=eq(3)
        pam(1,4)=0.0
        pam(1,5)=0.0
c        call gamoeba(pam,yam,v,ndmaxp1,ndmax,nd,iter,iopt,nemax,n,wk2,x,f,sumsq)
        call gamoeba(pam,yam,v,ndmaxp1,ndmax,nd,iter,iopt,nemax,n1,n,wk2,x,f,sumsq)		!Hakan
        eq(1)=pam(i1,i1)
        eq(2)=pam(i1,i2)
        eq(3)=abs(pam(i1,i3))
        eq(5)=pam(i1,i5)
        eq(6)=pam(i1,i4)
c        if(eq(2).le.0)then
c          jdt=3
c          do i=1,7
c            eq(i)=0.0
c          enddo
c          call trend(n,x,f,eq(3),eq(5))
c        endif
        eq(7)=jdt
      elseif(jdt.eq.8)then
c        call gexp(x,n,f,wk2,wk3,eq,flag)         !general exponential curve
        call gexp(x,n1,n,f,wk2,wk3,eq,flag)       !Hakan
        do i=1,n
          wk2(i)=i
          wk3(i)=f(i)
        enddo
        iopt=2
	       nd=iopt+i2
        pam(1,1)=eq(1)
        pam(1,2)=eq(2)
        pam(1,3)=0.0
        pam(1,4)=0.0
        pam(1,5)=0.0
c        call gamoeba(pam,yam,v,ndmaxp1,ndmax,nd,iter,iopt,nemax,n,wk2,x,f,sumsq)
        call gamoeba(pam,yam,v,ndmaxp1,ndmax,nd,iter,iopt,nemax,n1,n,wk2,x,f,sumsq)		!Hakan
        eq(1)=pam(i1,i1)
        eq(2)=r1
        eq(3)=abs(pam(i1,i2))
        eq(5)=pam(i1,i4)
        eq(6)=pam(i1,i3)
        if(eq(2).le.0)then
          jdt=3
          do i=1,7
            eq(i)=0.0
          enddo
c          call trend(n,x,f,eq(3),eq(5))
          call trend(n1,n,x,f,eq(3),eq(5))			!Hakan
        endif
        eq(7)=jdt
      else
        if(jdt.gt.9)then
          ls=jdt											                      !fixed 50% frequency response
        elseif(jdt.lt.-9)then
          ls=ifix(abs((float(jdt)/100.0)*float(n)))			!%n frequency response
        endif
        do i=1,n
          wk2(i)=1.0
        enddo
ccc        call spl2w(x,n,f,wk2,p,ls,mxy)
ccc        call spl2w(mxy,n,x,wrk2,pp,f,ls)
        call spline(n,x,f,p,ls)						            !cubic spline filter
      endif
      return
      end
