program spctime implicit none c-------------------------------------------------------------------- integer i,j,t,n,pt,pn,NT,NL,NM,itr,m,nlat,hnlat,nskip,nvar integer iy, nyr cccc parameter (itr=2,NT=itr*360,NL=96) !itr=No.per day (1x or 4x) parameter (itr=1,NT=itr*96,NL=128,NM=12,nlat=12) parameter (nskip=0,nvar=1,nyr=49) 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 tmp(NL,NM),tmp2(NL,nlat,NT),asym(NL,nlat,NT), sym(NL,nlat,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 character*8 FOUT data FOUT /'out.pres'/ 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= * '/atmos/milee/dkim/mbmt/work/ST/data/snu_noit/st_in.prcp.gdat', * access='direct', * form='unformatted',recl=NL*NM,status='old') c do t = 1, NT/2+1 do n = 1, NL+1 yrPEE(n,t)=0. enddo enddo c do 2003 iy=1,nyr print*,iy do 20 t=1,NT read(1,rec=nskip+(iy-1)*NT+t) tmp do 30 m=1,nlat do 30 n=1,NL tmp2(n,m,t)=tmp(n,m) 30 continue 20 continue c c call detrend (tmp2,NL,nlat,NT) do 40 t=1,NT do 40 n=1,NL do 50 m=1,nlat/2+1 sym(n,m,t)=tmp2(n,m,t)*0.5+tmp2(n,nlat-m+1,t)*0.5 50 continue do 55 m=1,nlat/2 sym(n,nlat-m+1,t)=sym(n,m,t) 55 continue 40 continue do 60 t=1,NT do 60 n=1,NL do 60 m=1,nlat asym(n,m,t)=tmp2(n,m,t)-sym(n,m,t) 60 continue 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 c print*,m c globmean=0. do 100 t=1,NT do 100 n=1,NL EE(n,t) = sym(n,m,t) ! symmetric / asymmetric c EE(n,t) = tmp2(n,m,t) globmean = globmean + REAL(EE(n,t))/float(NT*NL) 100 continue c print*,'eif=',eif,' k=',k totvar=0. 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) = (ABS(EE(n,t)))**2 enddo enddo c do pt=1,NT/2+1 c do pn=1,NL+1 c write(21,*) 'pt=',pt,' pn=',pn,' PEE=',PEE(pn,pt) c enddo c enddo c do t=1,NT c do n=1,NL c if (CABS(EE(n,t)).gt..05) then c write(21,*) 't=',t,' n=',n,' EE=',EE(n,t),' ABS(EE)', c # CABS(EE(n,t)) c endif c enddo c enddo c do t = 1, NT/2+1 do n = 1, NL+1 c sumPEE(n,t)=sumPEE(n,t)+PEE(n,t)/float(nlat) sumPEE(n,t)=sumPEE(n,t)+PEE(n,t) ! summed over lat. 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 c c-------------------------------------------------------------------- c------------ GRADS OUTPUT ------------------------ c-------------------------------------------------------------------- open(11,file='./gdat/'//FOUT,access='direct', * form='unformatted',recl=(NL+1)*(NT/2+1),status='unknown') write (11,rec=1) yrPEE open(12,file='./gdat/'//FOUT//'.ctl',status='unknown') write (12,1234) FOUT,NL+1,(ss(i),i=1,nl+1) * ,NT/2+1,(ff(i),i=1,nt/2+1) 1234 format ('DSET ^',a8/ & '*'/ & 'UNDEF -999.'/ & 'TITLE SPCTIME OUTPUT'/ & '*'/ & 'XDEF ',i4,1x,'LEVELS ',129(1x,f6.1) / & '*'/ & 'YDEF ',i4,1x,'LEVELS ',49(1x,f6.3) / ) write (12, 1235) '* ' write (12, 1235) 'ZDEF 1 LEVELS 1000. ' write (12, 1235) '* ' write (12, 1235) 'TDEF 1 LINEAR 1mar2001 1dy ' write (12, 1235) '* ' write (12, 1235) 'VARS 1 ' write (12, 1235) 'power 0 99 power ' write (12, 1235) 'ENDVARS ' 1235 format (a30) c c working-out plotting boundaries c print*,'sleft=',ss(xmi),' sright=',ss(xma) c print*,'ffbot=',ff(ymi),' fftop=',ff(yma) c!!!!!!! NO NCAR GRAPHICS : skip these !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! END c----------------------------------------------------------------------- subroutine detrend (tmp2,NL,nlat,NT) implicit none integer NL,nlat,NT,n,m,t, tmean,sum real tmp2(NL,nlat,NT),y(NT) real a, b do 10 n=1,NL do 10 m=1,nlat do 20 t=1,NT y(t)=tmp2(n,m,t) 20 continue tmean=sum(y,NT) call reg (y,NT,a,b) do 30 t=1,NT tmp2(n,m,t)=tmp2(n,m,t)-(a+b*float(t-1)) c tmp2(n,m,t)=tmp2(n,m,t)-tmean/float(NT) 30 continue 10 continue return end c subroutine reg (y,n,a,b) implicit none integer i,n real a,b,xbar,ybar,sum real y(n) real xy(n),xx(n) c xbar=n*(n+1)/2. ybar=sum (y,n)/float(n) do 100 i=1,n xy(i)=(i-xbar)*(y(i)-ybar) xx(i)=(i-xbar)**2 100 continue b=sum(xy,n)/sum(xx,n) a=ybar-b*xbar c return end function sum (x,n) implicit none integer n,i real x(n),sumx, sum sumx=0. do 100 i=1,n sumx=sumx+x(i) 100 continue sum=sumx return end