      program pcspectrum
      implicit none

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c NOTE :
c Daehyun Kim modified Dr. Matthew C. Wheeler's code for the MJO
c Simulation Metrics calculation system
c Note that any bugs or errors might be caused by modifications.  
c Contact : Daehyun Kim (kim@climate.snu.ac.kr)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc


c Like pcspectrum.f, except input is from a "createdPCs" file instead of
c an EOF-output file.
c
c Spectrum is computed for 6-month segments (parts), padded with zeroes to 256 days.
c The power from each segment is averaged. 
c
c For the red-noise spectrum, which is a function of the lag-1 autocorrelation,
c the lag-1 autocorrelation is an average of the value from all segments (parts).

      integer lrecIn1
      integer MAXTIME,NS,NPC,nyr
      parameter (MAXTIME=10000)   ! 23+ yrs
      parameter (NS=183)          ! Number of days per segment
      parameter (NPC=2)           ! 2 PCs
      parameter (nyr = 27)
      integer MT,num
      integer nyi(MAXTIME),nmi(MAXTIME),ndi(MAXTIME)
      real pc(MAXTIME,NPC)
      integer imd,itseg(MAXTIME), t, is, n_year, iy

      integer year1, leap, linux
      parameter (year1 = 1979, leap = 1)
      parameter (linux = 4)
      integer yy, dd, nod, it, skip, nseg

      integer NT
      parameter (NT=256)  ! enough for 6 month segments
      integer i,mx,ll
      real r1,r1AV(NPC)    ! lag-1 autocorrelations
      real ff(NT/2),vv(NT/2),vredno(NT/2)  ! plot arrays
      real vvAV(NT/2,NPC)
      real ts(NT),a(NT/2),b(NT/2),v(NT/2)
      real mean,wsave(2*NT+15)
      real totvar,totvarAV(NPC)
      real v3080,v3080red
      real sum,pos
      real vpl,vpr,vpb,vpt,ul,ur,ub,ut,cfux,cfuy
      character dd1*1,dd2*2,dd4*4,dd7*7,dd40*40
      character*100 FOUT
      data FOUT
     &/'/home/enver/work/programs/MJOWG/msd/level_2/ceof/sp256.ts'/

      if (NT/2..ne.float(NT/2)) then
       print*,'make NT divisible by 2 - STOP'
       STOP
      endif

c Input
      inquire(IOLENGTH=lrecIn1) ts(1)

      open(1,file=
     &'/home/enver/work/programs/MJOWG/msd/level_2/ceof/'//
     &'ceof.ts.pr',
     *   access='direct',
c     *   form='unformatted',recl=1*linux,status='old')
     *   form='unformatted',recl=lrecIn1,status='old')

c Output
      open(11,file=trim(FOUT)//'01',status='unknown')
      open(12,file=trim(FOUT)//'02',status='unknown')

      open(21,file=trim(FOUT)//'01.var',status='unknown')
      open(22,file=trim(FOUT)//'02.var',status='unknown')

       
c  --------------------
c    Compute spectra of the PCs for each NS-long segment

      DO 88 imd=1,NPC

      nseg  = 0

      call rffti(NT,wsave)   ! initialize

       do i=1,NT/2
        vvAV(i,imd)=0.
       enddo
       r1AV(imd)=0. 

      do 2003 is= 1,2

c winter
        if (is.eq.1) then 
          n_year = nyr-1
c summer
        elseif (is.eq.2) then
          n_year = nyr
        endif

       yy = year1 - 1
       dd = 0

      do 2003 iy=1,n_year
       yy = yy + 1

      if (leap.eq.1) then

       nod = 365
       if (mod(yy,4).eq.0.and.mod(yy,100).ne.0) nod = 366
       if (mod(yy,400).eq.0) nod = 366

       if (is.eq.1) then

         skip = 304
         if (nod.eq.366) skip = 305

       elseif (is.eq.2) then

        skip = 120
        if (nod.eq.366) skip = 121

       endif

      elseif (leap.eq.2) then

       nod = 365

        if (is.eq.1) then
         skip = 304
        elseif (is.eq.2) then
         skip = 120
        endif

      endif

       nseg = nseg + 1

      do 20 t=1,NS

       if (leap.eq.1) then

           it = dd + skip + t

        elseif (leap.eq.2) then

           it = nod*(iy-1)+skip+t

        endif

         read(1,rec=10*(it-1)+imd) ts(t)

  20  continue

        call detrend(ts,NT,1,NS)
c      Apply a tapering to the first and last 10 days by multiplication
c      by a segment of the cosine curve so that the ends of the series
c      taper towards zero (mean has been removed anyway)! Rest is padded
c      with zeroes.
        call tapertozero(ts,NT,1,NS,10)
        call fastftM(ts,NT,wsave,a,b,v,totvar,mean)

        do i=1,NT/2
         vvAV(i,imd)=vvAV(i,imd)+v(i)
        enddo

c      calculate lag-1 autocorrelation coefficient.
        call AUTO(ts,NT,NS,1,r1)
        r1AV(imd)=r1AV(imd)+r1

c for leap year
      if (leap.eq.1) then
       dd = dd + nod
      endif

2003   continue

        totvarAV(imd)=0.
        do i=1,NT/2
         vvAV(i,imd)=vvAV(i,imd)/float(nseg)
         totvarAV(imd)=totvarAV(imd)+vvAV(i,imd)
        enddo
        r1AV(imd)=r1AV(imd)/float(nseg)
        print*,'imd=',imd,'  No. of segments=', nseg,' r1AV=',r1AV(imd)

  88   CONTINUE

      do 109, imd=1,NPC

c      calculate power spectrum of lag-1 autoregressive process.
       call REDNO(vredno,NT/2,r1AV(imd),totvarAV(imd))

c  ***Plot of Spectrum
       v3080=0.
       v3080red=0.          
       do i=1,NT/2
        ff(i)=float(i)/float(NT)
        vv(i)=vvAV(i,imd)
        if (ff(i).ge.(1./80.).and.ff(i).le.(1/30.)) then
         v3080=v3080+vvAV(i,imd)
         v3080red=v3080red+vredno(i)
        endif
       enddo
c      smooth spectrums
c      do i=1,10
c       call smooth121(vv,NT/2,NT/2)
c       call smooth121(vredno,NT/2,NT/2)
c      enddo
c      multiply variance by frequency to have area-conserving....
       do i=1,NT/2
        vv(i)=vv(i)*ff(i)
        vredno(i)=vredno(i)*ff(i)
       enddo
c     calculate fraction variance in the 30 to 80-day band.
       v3080=v3080/totvarAV(imd)
       v3080red=v3080red/totvarAV(imd)
       print*,'variance in 30-80 day band=',v3080
       print*,'variance in 30-80 day band for redno=',v3080red
        do 111 i = 1, NT/2
111     write(10+imd,22) i, ff(i),vv(i),vredno(i)
22      format(i5,1x,3e13.5)
        write(20+imd,32) v3080*100, v3080red*100
        write(20+imd,32) v3080*100, v3080red*100
        write(20+imd,32) v3080*100, v3080red*100
32      format(2f13.2)

109   continue

      END   !main program

c --------------------------------------------------------------------------
      SUBROUTINE fastftM(ts,N,wsave,a,b,v,totvar,mean)
      implicit none

c   Uses the fftpack routines.
c   N.b. If many FFT calls will be made with the same N, then use fastftM.f
c
c   The a and b are the cosine and sine coefficients respectively.
c   v is the variance in each bin. totvar is the total variance.
c   N, the length of the series, can be either even or odd with numco=N/2

      integer N,numco,i
      real ts(N),coef(20000),a(N/2),b(N/2),v(N/2)
      real totvar,mean,wsave(2*N+15)     !>2*N+15

      numco = N/2       ! This works if N is either odd or even.
      totvar = 0.0

      do i=1,N
       coef(i)=ts(i)
      enddo

c Do the FFT.
      call rfftf(N,coef,wsave)

      do i=1,N
         coef(i) = coef(i) / (float(N)/2.)
      enddo

      mean = coef(1)/2.

      do 950  i = 1, numco-1
       a(i) = coef((i)*2)
       b(i) = -coef((i)*2+1)
       v(i) = (a(i)**2 + b(i)**2) / 2.0
c    Note that the variance of a sine or cosine wave is 1/2
       totvar = totvar + v(i)
 950  continue

        if (float(N/2).eq.(float(N)/2.)) then
c get 'a' (cos coef) for the Nyquist frequency, valid if n is even.
         a(numco) = coef(2*numco) / 2.0
         b(numco) = 0.0
         v(numco) = a(numco) ** 2
         totvar = totvar + v(numco)
        else
c get 'a' (cos coef) and 'b' (sin coef) for the Nyquist
c frequency, valid if N is odd.
         a(numco) = coef(2*numco)
         b(numco) = coef(2*numco+1)
         v(numco) = (a(numco)**2 + b(numco)**2) / 2.0
         totvar = totvar + v(numco)
        endif

 970   continue
      return
      END

c-----------------------------------------------------------------
       subroutine detrend(X2,nx,vmi,nv)
       implicit none

c Removes the linear-squares best fit line for the series.
c Only the data from vmi to vmi+nv-1 is deemed useful, and
c the rest is zeroed. (There are nv points of useful data)
c The line of best fit is given by X2' = a + bX1.
c The X1s are assumed to be evenly spaced ( i.e. X1=float(i) )
c
c m1 = mean of X1
c m2 = mean of X2
c m3 = mean of (X1*X2)
c m4 = mean of (X1*X1)
c m5 = mean of (X2*X2)

       integer nv,i,nx,vmi
       real X2(nx),m1,m2,m3,m4,m5,r,a,b

       if (vmi+nv-1.gt.nx) then
        print*,'ERROR as vmi+nv-1 > nx'
        STOP
       endif

       m1=0.
       m2=0.
       m3=0.
       m4=0.
       m5=0.

       do 10 i = vmi,vmi+nv-1
        m1 = m1+float(i)/float(nv)
        m2 = m2+X2(i)/float(nv)
        m3 = m3+(float(i)*X2(i))/float(nv)
        m4 = m4+(float(i)*float(i))/float(nv)
        m5 = m5+(X2(i)*X2(i))/float(nv)
  10   continue

       r=(m3-m1*m2)/(sqrt(m4-m1**2)*sqrt(m5-m2**2))
       b=(m3-m1*m2)/(m4-m1**2)
       a=m2-b*m1

       do 20 i = vmi,vmi+nv-1
        X2(i) = X2(i) - (a + b*float(i))
  20   continue
c      print*,'Trend line has a=',a,' b=',b

       if (vmi.gt.1.) then
        do 15 i = 1,vmi-1
         X2(i) = 0.
  15    continue
       endif
       if (vmi+nv.le.nx) then
        do 25 i = vmi+nv,nx
         X2(i) = 0.
  25    continue
       endif

       return
       end

c-----------------------------------------------------------------
       subroutine tapertozero(ts,N,nmi,nn,tp)
       implicit none

c "taper" the first and last 'tp' members of 'ts' by multiplication by
c a segment of the cosine curve so that the ends of the series
c taper toward zero. This satisfies the
c periodic requirement of the FFT.
c Only the data from nmi to nmi+nn-1 is deemed useful, and
c the rest is set to zero. (There are nn points of useful data)
       integer i,j,N,tp,nmi,nn
       real Pi,ts(N)
       parameter (Pi=3.1415926)   !tp is number to taper on each end.

       if (N.lt.tp*2.or.nn.lt.tp*2) then
        print*,'No use doing the tapering if less than tp*2 values!'
        STOP
       endif
       if (nmi+nn-1.gt.N) then
        print*,'ERROR as nmi+nn-1 > N'
        STOP
       endif

       do i=1,N
        j=i-nmi+1
        if (j.le.0.or.j.gt.nn) then
         ts(i)=0.
        elseif (j.le.tp) then
         ts(i)= ts(i)*.5*(1.-COS((j-1)*Pi/float(tp)))
        elseif (j.gt.(nn-tp).and.j.le.nn) then
         ts(i)= ts(i)*.5*(1.-COS((nn-j)*Pi/float(tp)))
        else
         ts(i)= ts(i)
        endif
       enddo

       return
       end
c-----------------------------------------------------------------

      subroutine smoothrm(vv,vn,nn,pts)
      implicit none
c Smooths vv by passing it through a pts-point running mean.
c Towards the beginning and end of the series the number of points
c is reduced and it is also the case that the total sum is conserved.
c There are 'nn' pieces of useful information, which may be less
c than or equal to 'vn'.
c  Make sure pts is ODD.
      integer nn,pts,i,j,hpts,vn 
      real vv(vn),dum(5000)

      if(float(pts/2).eq.float(pts)/2.) then
       print*,'pts should be odd, STOPPING'
       STOP
      endif
      hpts = (pts-1)/2
      do i=1,hpts
       dum(i)=0.
       do j=1,(hpts-i+1)
        dum(i)=dum(i)+2.*vv(j)/float(pts) 
       enddo
       do j=hpts-i+2,hpts+i
        dum(i)=dum(i)+vv(j)/float(pts) 
       enddo
      enddo
      do i=(hpts+1),(nn-hpts)
       dum(i)=0.
       do j=i-hpts,i+hpts
        dum(i)=dum(i)+vv(j)/float(pts)
       enddo
      enddo
      do i=(nn-hpts+1),nn
       dum(i)=0.
       do j=(i-hpts),(2*nn-i-hpts)
        dum(i)=dum(i)+vv(j)/float(pts)
       enddo
       do j=(2*nn-i-hpts+1),nn
        dum(i)=dum(i)+2.*vv(j)/float(pts)
       enddo
      enddo

      do i=1,nn
       vv(i)=dum(i)
      enddo
      RETURN
      END

c-----------------------------------------------------------------
      subroutine smooth121(vv,vn,nn)
      implicit none
c Smooths vv by passing it through a 1-2-1 filter.
c The first and last points are given 3-1 (1st) or 1-3 (last)
c weightings (Note that this conserves the total sum).
c The routine also skips-over missing data (assigned to be
c a value of 1.E36).
c There are 'nn' pieces of useful information, which may be less
c than or equal to 'vn'.

      integer nn,i,vn
      real spv,vv(vn),dum(5000)

      if (nn.gt.5000) then
       print*,'need to increase 5000 in smooth121.f'
       STOP
      endif

      spv=1.E36
      i=0
 10   continue
      i=i+1
      if (vv(i).eq.spv) then
       dum(i)=spv
      elseif(i.eq.1.or.vv(i-1).eq.spv) then
       dum(i)=(3.*vv(i)+vv(i+1))/4.
      elseif(i.eq.nn.or.vv(i+1).eq.spv) then
       dum(i)=(vv(i-1)+3.*vv(i))/4.
      else
       dum(i)=(1.*vv(i-1)+2.*vv(i)+1.*vv(i+1))/4.
      endif
      if (i.ne.nn) goto 10

      do i=1,nn
       vv(i)=dum(i)
      enddo
      RETURN
      END


c------------------------------------------------------------------
      SUBROUTINE  AUTO  (TS, N, NT, LAG, R)
 
C        THE AUTOCORRELATION COEFFICIENT FOR THE GIVEN LAG IS
C        CALCULATED.  SEE MITCHELL, ET AL., 1966, P. 60.
 
      DIMENSION  TS(N)
 
      REAL*8  SUM1, SUM2, SUM3, SUM4, SUM5
      INTEGER  NTEMP1, NTEMP2, NTEMP3, LAG,NT
 
      NTEMP1 = NT - LAG
 
      SUM1 = 0.0
      SUM2 = 0.0
      SUM3 = 0.0
      SUM4 = 0.0
      SUM5 = 0.0
 
 
      DO 2970  I = 1, NTEMP1
         NTEMP2 = I + LAG
         SUM1 = SUM1 + TS(I) * TS(NTEMP2)
         SUM2 = SUM2 + TS(I)
         SUM4 = SUM4 + TS(I) ** 2
2970  CONTINUE
 
      NTEMP3 = LAG + 1
      DO 2980  I = NTEMP3, NT
         SUM3 = SUM3 + TS(I)
         SUM5 = SUM5 + TS(I) ** 2
2980  CONTINUE
 
      R = ((N-LAG) * SUM1 - SUM2 * SUM3) / (SQRT ((N-LAG) * SUM4 -
     1    SUM2 ** 2) * SQRT ((N-LAG) * SUM5 - SUM3 ** 2))
 
c     write(*,2990) LAG, R
2990  FORMAT ('0AUTOCORRELATION COEF, LAG ',I3,' = ', F10.7)
 
 
      RETURN
      END
 
c------------------------------------------------------------------
 
      SUBROUTINE  REDNO (RN, N, R, totvar)
 
 
C        THE RED NOISE SPECTRUM FOR THE GIVEN LAG 1 AUTOCORRELATION
C        COEFFICIENT IS CALCULATED.  SEE GILMAN, ET AL., JAS, 1963.
 
      DIMENSION  RN (N)
 
      REAL  SBAR
 
      PI = 3.141592654

c     SBAR = 1.0 / FLOAT (N)
      SBAR = totvar / FLOAT (N)
 
      DO 3000  I = 1, N
         RN(I) = SBAR * (1 - R**2) / (1 + R**2 - 2 * R *
     1           COS (PI * FLOAT(I) / FLOAT(N)))
3000  CONTINUE
 
      RETURN
      END
 
c----------------------------end of program---------------------------------

