      program spctime
      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--------------------------------------------------------------------
      integer i,j,t,n,pt,pn,NT,NL,NM,itr,m,nlat,hnlat,nskip,nvar
      integer iy, nyr
      parameter (itr=1,NT=96,NL=144,NM=13,nlat=13)
      parameter (nskip=0,nvar=1,nyr=272)
c Make NT and NL even numbers (it's just easier that way!)
c If they are not even, then code will need to be changed.
      real olr(NL,NM,NT)
      complex EE(NL,NT),CEE1(NL),CEE2(NT)
      real WSAVE1(4*NL+15),WSAVE2(4*NT+15)
      real PEE(NL+1,NT/2+1),totvar,globmean,totvar2,globmean2
      real sumPEE(NL+1,NT/2+1)
      real yrPEE(NL+1,NT/2+1)
      real totvarhalf
      real Rd,eif,s,tt,xx,xxm,ll,Pi,g,ps,k,ss(NL+1),ff(NT/2+1)
      real rand
c EE(n,t) = the initial (real) data set that later becomes the 
c         (complex) space-time spectrum
c PEE(n,t) = the (real) power spectrum
c ss() and ff() are arrays of wavenumbers and frequencies (cycles per day)
c          corresponding to PEE(,)
c totvar = the total variance of the dataset (about the global mean)
c globmean = the global time mean (average in space and time)
c NL = Number of longitudes in 'test' dataset
c NT = Number of times in 'test' dataset

c dimensional quantities are in MKS units.
c itr = No. of discrete samples per day
c s = planetary zonal wavenumber
c k = zonal wavenumber = 2 Pi s / L  (m-1)  
c ps = dimensional phase speed w.r.t. the earth 
C ll = distance around latitude circle (m)
C xx = longitude in degrees
c xxm = longitudinal distance from Greenwich in metres.
c id = time in days  (integer)
c tt = time in seconds (real)
c eif = dimensional frequency
c--------------------------------------------------------------------
c Variables for plotting routine
      integer il,iasf(13),li
      real sleft,sright,fftop,ffbot
      integer xmi,xma,mm,ymi,yma,nn
      parameter (sleft=-15.,sright=15.,fftop=.20,ffbot=0.)
c sleft/sright/fftop/ffbot are the zonal wavenumbers and frequencies
c (in cycles per day) to plot between.
      parameter (xmi=nint(sleft)+NL/2+1)
      parameter (xma=nint(sright)+NL/2+1)
      parameter (ymi=nint(ffbot*NT/(itr))+1)
      parameter (yma=nint(fftop*NT/(itr))+1)
      parameter (mm=xma-xmi+1,nn=yma-ymi+1)         
      real vpl,vpr,vpb,vpt,ul,ur,ub,ut,ffl,cfux,cfuy
      character*3 dd3
      character*2 dd2
      character*1 dd1
      character*80 label1
      data iasf/13*1/
c--------------------------------------------------------------------
c------------------------MAIN PROGRAM--------------------------------

c      print*,'xmi',xmi,' xma=',xma,' ymi=',ymi,' yma=',yma

      Pi = 3.1415926
      g = 9.81
      ll = 2.*Pi*6.37E06           ! at the equator
      totvar = 0.
      globmean = 0.

c Input 
      open(1,
     *file='/home/enver/work/programs/MJOWG/msd/level_2/wk99/u200_n1'//
     *'/data/seg96_over60_sym.gdat', access='direct',
     *   form='unformatted',recl=NL*NM*4,status='old')

      open(11,
     *file='/home/enver/work/programs/MJOWG/msd/level_2/wk99/u200_n1'//
     *'/power/power.sym.gdat',access='direct',
     *   form='unformatted',recl=(NL+1)*(NT/2+1)*4,
     *status='unknown')

c
      do t = 1, NT/2+1
      do n = 1, NL+1
        yrPEE(n,t)=0.
      enddo
      enddo
c    
      do 2003 iy=1,nyr

      do t=1,NT
  	 read(1,rec=nskip+(iy-1)*NT+t) ((olr(n,m,t),n=1,NL),m=1,nlat)
      enddo
c
      do t = 1, NT/2+1
      do n = 1, NL+1
        sumPEE(n,t)=0.
      enddo
      enddo
c
c! summed over lat (15S-15N)
c
      do 1000 m=1,nlat

      totvar=0.
      globmean=0.

      do 100 t=1,NT
      do 100 n=1,NL
	 EE(n,t) = olr(n,m,t)     ! symmetric / asymmetric
        globmean = globmean + REAL(EE(n,t))/float(NT*NL)
 100  continue
c      print*,'eif=',eif,' k=',k

      do t=1,NT
      do n=1,NL
       totvar = totvar + ((REAL(EE(n,t))-globmean)**2)/float(NT*NL)
      enddo
      enddo 
c--------------------------------------------------------------------
c------------------COMPUTING SPACE-TIME SPECTRUM---------------------

      call cffti(NL,WSAVE1)            ! Initialize the (complex) FFT
      do 200 t=1,NT
       do 150 n=1,NL
        CEE1(n) = EE(n,t)
c CEE1(n) contains the grid values around a latitude circle
 150   continue
       call cfftf(NL,CEE1,WSAVE1)
       do 170 n=1,NL
        EE(n,t) = CEE1(n)/float(NL)
 170   continue        
 200  continue

c Now the array EE(n,t) contains the Fourier coefficients (in planetary
c wavenumber space) for each time.

      call cffti(NT,WSAVE2)      ! Initialize FFT for a different length. 
      do 300 n=1,NL
       do 250 t=1,NT
        CEE2(t) = EE(n,t)
c CEE2(t) contains a time-series of the coefficients for a single
c planetary zonal wavenumber
 250   continue
       call cfftf(NT,CEE2,WSAVE2)
       do 270 t=1,NT
        EE(n,t) = CEE2(t)/float(NT)
 270   continue
 300  continue

c Now the array EE(n,t) contains the (complex) space-time spectrum.

c Check some quantities such as the total variance.
      globmean2 = EE(1,1)
      totvar2 = 0.
      totvarhalf = 0.
      do t=1,NT
      do n=1,NL
       if(t.eq.1.and.n.eq.1) then
        continue
       else
        totvar2 = totvar2+(ABS(EE(n,t)))**2
        if(t.gt.1.and.t.lt.NT/2+1) then
         totvarhalf = totvarhalf+(ABS(SQRT(2.)*EE(n,t)))**2
        elseif(t.eq.1.or.t.eq.NT/2+1) then
         totvarhalf = totvarhalf+(ABS(EE(n,t)))**2
        endif
c      N.b. summing for totvarhalf only works if NT is even.
       endif
      enddo
      enddo

c     print*,'globmean =',globmean,'  globmean2=',globmean2
c     print*,'totvar =',totvar,'  totvar2=',totvar2,' totvarhalf=',
c    *     totvarhalf

c Create array PEE(NL+1,NT/2+1) which contains the (real) power spectrum.
c In this array, the westward moving waves will be from pn=1 to NL/2;
c The eastward waves will be for pn=NL/2+2 to NL+1.
c Positive frequencies will be from pt=1,NT/2+1.
c Information about zonal mean (wavenumber 0) will be for pn=NL/2+1.
c Information about time mean will be for pt=1.
c Information about the Nyquist Frequency is at pt=NT/2+1

      do pt=1,NT/2+1
       ff(pt)=float(pt-1)/float(NT)*float(itr) 
       do pn=1,NL+1
       ss(pn)=float(pn-1-NL/2)
       if(pn.le.NL/2) then
        n=NL/2+2-pn
        t=pt
       elseif(pn.eq.NL/2+1) then
        n=1
        t=pt
       elseif(pn.ge.NL/2+2) then
        n=pn-NL/2
        if (pt.eq.1) then
         t=pt
        else
         t=NT+2-pt
        endif
       endif
       PEE(pn,pt) = (CABS(EE(n,t)))**2
       enddo
      enddo

      do t = 1, NT/2+1
      do n = 1, NL+1
       if (mod(nlat,2).ne.0) then 
        if(m.eq.nlat/2+1) then
        else
         sumPEE(n,t)=sumPEE(n,t)+PEE(n,t)   ! summed over lat.
        endif
       else
         sumPEE(n,t)=sumPEE(n,t)+PEE(n,t)   ! summed over lat.
       endif
      enddo
      enddo
1000  continue

	do t = 1, NT/2+1
	do n = 1, NL+1
	   yrPEE(n,t)=yrPEE(n,t)+sumPEE(n,t)/float(nyr)
        enddo
	enddo
2003  continue

      write (11,rec=1) yrPEE
                                   
      END
