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 integer lrecIn1 integer MAXTIME,NS,NPC,nyr parameter (MAXTIME=10000) ! 23+ yrs parameter (NS=180) ! Number of days per segment parameter (NPC=1) ! 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 = 1) integer yy, dd, nod, it, skip, nseg, nd, imode character*2 number integer NT parameter (NT=180) ! 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, FOUT2 data FOUT &/'/home/enver/work/programs/MJOWG/msd/level_1/u200_n1/eof/'/ data FOUT2 &/'sum/pcps.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_1/u200_n1/eof/'// &'sum/eof.ts.pr', * access='direct', cEnver * form='unformatted',recl=1*linux,status='old') * form='unformatted',recl=lrecIn1,status='old') do imode = 1, 5 write(number,'(i2.2)') imode c Output open(11,file=trim(FOUT)//trim(FOUT2)//number,status='unknown') 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= 2, 2 n_year = nyr yy = year1 - 1 dd = 0 do 2003 iy=1,n_year yy = yy + 1 nseg = nseg + 1 do 20 t=1,NS if (leap.eq.1) then if (is.eq.1) then it = dd + t elseif (is.eq.2) then it = 184*(iy-1)+t endif elseif (leap.eq.0) then if (is.eq.1) then it = 181*(iy-1)+t elseif (is.eq.2) then it = 184*(iy-1)+t endif endif read(1,rec=10*(it-1)+imode) ts(t) 20 continue c for leap year winter season if (leap.eq.1.and.is.eq.1) then nd = 181 if (mod(yy+1,4).eq.0.and.mod(yy+1,100).ne.0) nd = 182 if (mod(yy,400).eq.0) nd = 182 dd = dd + nd endif 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 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 do i=1,NT/2 ff(i)=float(i)/float(NT) vv(i)=vvAV(i,imd) enddo do i=1,NT/2 vv(i)=vv(i)*ff(i) vredno(i)=vredno(i)*ff(i) enddo do 111 i = 1, NT/2 111 write(11,22) i, ff(i),vv(i),vredno(i) 22 format(1x,i5,3e13.5) 109 continue enddo ! imode 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 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---------------------------------