program spacetimeCross implicit none c c This version set-up for datasets on standard 2.5 degree grid. c c Calculates the Co- and Quadrature spectrum between 2 fields c that have information in both space and time........... c c N.b. this program has been written to use only a small amount c of RAM.....thus the dataset is read-in NLAT times during it. c 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 (separated into eastward and westward) 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 P12(n,t) = co-spectrum | c Q12(n,t) = quadrature spectrum | All real numbers c Coh2(n,t) = coherence-squared statistic | On same grid as PEE(n,t) c Phas(n,t) = Phase (in radians) | c c Like gridspectparts.f this program calculates the space-time c spectra for many parts of a longer (say 10 yr) time series. c Each part is of length NT and is separated by NP time steps. c The longitudinal domain of the S-T spectra may be restricted c using the olr_avs xlonw and xlone. Note that tapering in c space and time is done to make the datasets periodic. c The calculations are also only made for the latitudes c between sout and nout. Also note that all longitudes c are still kept for the calculation so that every planetary c wavenumber will be represented. c c Plotting is done later by the routine frqxwnST.f c c OUTPUT is of the form. c NT,NL,xlonw,xlone,sout,nout !> c (ff(t),t=1,NT+1) !> A HEADER c (ss(n),n=1,NL+1) !> c lat c do lat=sout,nout c ((PEE1sum(n,t),n=1,NL+1),t=1,NT+1) c ((PEE2sum(n,t),n=1,NL+1),t=1,NT+1) c ((P12sum(n,t),n=1,NL+1),t=1,NT+1) c ((Q12sum(n,t),n=1,NL+1),t=1,NT+1) c enddo c WHERE c NT = number of times of input (can be even or odd!) c NT = number of frequency bins of variance c NL = number of spatial bins of variance (probably 144) c xlonw = western longitude that will be forced to zero c xlone = eastern longitude that will be forced to zero c sout to nout are the bounding latitudes for the output. c ff(t) are the frequencies in cycles per day. c ss(n) are the planetary zonal wavenumbers. c P12(n,t,lat) and Q12(n,t,lat) are averaged over the many parts. c c corrections to perform auto-coherence !Enver logical lOpened c VARIABLES ASSOCIATED WITH THE FFT integer NT,NP,PP,itr,NL, NM, ip, m c Make NT and NL EVEN numbers as it is easier that way!! parameter (NT=num_s) ! Length of 'part' to be FFTed parameter (NL=144) ! Remember that NLON repeats the first longitude parameter (NM=13, PP=192) real olr(NL,NM,NT), u850(NL,NM,NT) c VARIABLES ASSOCIATED WITH THE INTERNAL PROCESSING IN THIS PROGRAM integer j,tt,lat,lon,lon2,Pcount,n,t,pn,pt real WSAVE1(4*NL+15),WSAVE2(4*NT+15) real ts1(NT),ts2(NT),ls1(NL),ls2(NL) complex EE1(NL,NT),EEo1(NL,NT),CEEa(NL),CEEb(NT) complex EE2(NL,NT),EEo2(NL,NT) real globmean,globmean2,totvar,totvar2 real frq(NT/2-1) c OUTPUT VARIABLES real ff(NT+1),ss(NL+1) real PEE1sum_sym(NL+1,NT/2+1),PEE2sum_sym(NL+1,NT/2+1) real PEE1sum_asy(NL+1,NT/2+1),PEE2sum_asy(NL+1,NT/2+1) real P12sum_sym(NL+1,NT/2+1),Q12sum_sym(NL+1,NT/2+1) real P12sum_asy(NL+1,NT/2+1),Q12sum_asy(NL+1,NT/2+1) real Coh2_sym(NL+1,NT/2+1),Phas_sym(NL+1,NT/2+1) real Coh2_asy(NL+1,NT/2+1),Phas_asy(NL+1,NT/2+1) real v1_sym(NL+1,NT/2+1),v2_sym(NL+1,NT/2+1) real v1_asy(NL+1,NT/2+1),v2_asy(NL+1,NT/2+1) c INPUT FILES open(11, #file='/home/enver/work/programs/MJOWG/msd/level_2/'// #'coh2/olr_av/data/'// #'seg256_over206_sym_asy.gdat', #recl=NL*NM*recl_linux, access='direct', # form='unformatted',status='old') inquire(file='/home/enver/work/programs/MJOWG/msd/level_2/'// #'coh2/olr_av/data/'// #'seg256_over206_sym_asy.gdat',OPENED=lopened) if ( .not. lopened ) then !corrections for auto-coherence Enver open(12, #file='/home/enver/work/programs/MJOWG/msd/level_2/'// #'coh2/olr_av/data/'// #'seg256_over206_sym_asy.gdat', # recl=NL*NM*recl_linux, access='direct', # form='unformatted',status='old') endif !corrections for auto-coherence Enver c OUTPUT FILE open(21, #file='/home/enver/work/programs/MJOWG/msd/level_2/coh2/'// #'olr_av/power/power.sym.gdat', # recl=(NL+1)*(NT/2+1)*recl_linux, # access='direct', #form='unformatted',status='unknown') open(22, #file='/home/enver/work/programs/MJOWG/msd/level_2/coh2/'// #'olr_av/power/power.asy.gdat', # recl=(NL+1)*(NT/2+1)*recl_linux, #access='direct', #form='unformatted',status='unknown') open(31, #file='/home/enver/work/programs/MJOWG/msd/level_2/coh2/'// #'olr_av/power/coh2.sym.gdat', # recl=(NL+1)*(NT/2+1)*recl_linux, #access='direct', #form='unformatted',status='unknown') open(32, #file='/home/enver/work/programs/MJOWG/msd/level_2/coh2/'// #'olr_av/power/coh2.asy.gdat', # recl=(NL+1)*(NT/2+1)*recl_linux, #access='direct', #form='unformatted',status='unknown') c-------------------------------------------------------------------- c------------------------MAIN PROGRAM-------------------------------- c EE(n,t) = the initial (real) data set that later becomes the c (complex) space-time spectrum c ss() and ff() are arrays of wavenumbers and frequencies (cycles per day) c corresponding to PEE( ) if(float(NT/2).ne.float(NT)/2.) then print*,'NT must be even, it is probably more efficient that way' STOP elseif(float(NL/2).ne.float(NL)/2.) then print*,'NL must be even' STOP endif c WRITE HEADER TO OUTPUT do t=1,NT+1 ff(t) = float(t-1-NT/2)/float(NT) ! in cycles per day. enddo do n=1,NL+1 ss(n) = float(n-1-NL/2) enddo c Initialize SUM arrays c out of bounds subscript: noticed from Dr. Dennis Shea at NCAR c do t=1,NT+1 do t=1,NT/2+1 do n=1,NL+1 PEE1sum_sym(n,t)=0. PEE2sum_sym(n,t)=0. P12sum_sym(n,t)=0. Q12sum_sym(n,t)=0. PEE1sum_asy(n,t)=0. PEE2sum_asy(n,t)=0. P12sum_asy(n,t)=0. Q12sum_asy(n,t)=0. enddo enddo C***************HERE IS THE BIG LOOP********************************** do 2007 ip=1,PP c read do t=1,NT read(11,rec=(ip-1)*NT+t) ((olr(n,m,t),n=1,NL),m=1,NM) if ( .not. lopened ) then !corrections for auto-coherence Enver read(12,rec=(ip-1)*NT+t) ((u850(n,m,t),n=1,NL),m=1,NM) else u850(:,:,t) = olr(:,:,t) endif enddo c Do a loop in latitude!! - average do m = 1, NM do 100 t=1,NT do 100 n=1,NL EE1(n,t) = olr(n,m,t) EE2(n,t) = u850(n,m,t) 100 continue c Initialize FFTs call cffti(NL,WSAVE1) call cffti(NT,WSAVE2) c-------------------------------------------------------------------- c---------------COMPUTING SPACE-TIME SPECTRUM for EE1---------------- do 201 t=1,NT do 151 n=1,NL CEEa(n) = EE1(n,t) c CEEa(n) contains the grid values around a latitude circle 151 continue call cfftf(NL,CEEa,WSAVE1) do 171 n=1,NL EE1(n,t) = CEEa(n)/float(NL) 171 continue 201 continue c Now the array EE(n,t) contains the Fourier coefficients (in planetary c wavenumber space) for each time. do 300 n=1,NL do 251 t=1,NT CEEb(t) = EE1(n,t) c CEEb(t) contains a time-series of the coefficients for a single c planetary zonal wavenumber 251 continue call cfftf(NT,CEEb,WSAVE2) do 270 t=1,NT EE1(n,t) = CEEb(t)/float(NT) 270 continue 300 continue c Now the array EE1(n,t) contains the space-time spectrum. c-------------------------------------------------------------------- c---------------COMPUTING SPACE-TIME SPECTRUM for EE2---------------- do 202 t=1,NT do 152 n=1,NL CEEa(n) = EE2(n,t) c CEEa(n) contains the grid values around a latitude circle 152 continue call cfftf(NL,CEEa,WSAVE1) do 172 n=1,NL EE2(n,t) = CEEa(n)/float(NL) 172 continue 202 continue c Now the array EE(n,t) contains the Fourier coefficients (in planetary c wavenumber space) for each time. do 302 n=1,NL do 252 t=1,NT CEEb(t) = EE2(n,t) c CEEb(t) contains a time-series of the coefficients for a single c planetary zonal wavenumber 252 continue call cfftf(NT,CEEb,WSAVE2) do 272 t=1,NT EE2(n,t) = CEEb(t)/float(NT) 272 continue 302 continue c Now the array EE2(n,t) contains the space-time spectrum. c------------------------------------------------------------------- c Create array PEE(NL+1,NT+1) which contains the (real) power spectrum. c In this array, the negative wavenumbers will be from pn=1 to NL/2; c The positive wavenumbers will be for pn=NL/2+2 to NL+1. c Negative frequencies will be from pt=1 to NT/2. c Positive frequencies will be from pt=NT/2+2 to NT+1. c Information about zonal mean will be for pn=NL/2+1. c Information about time mean will be for pt=NT/2+1. c Information about the Nyquist Frequency is at pt=1 and pt=NT+1 c In PEE, I define the WESTWARD waves to be either +ve frequency c and -ve wavenumber or -ve freq and +ve wavenumber. c EASTWARD waves are either +ve freq and +ve wavenumber OR -ve c freq and -ve wavenumber. c original one c do 191 pt=1,NT+1 c do 189 pn=1,NL+1 c if(pn.le.NL/2) then c n=NL/2+2-pn c if(pt.le.NT/2) then c t=NT/2+pt c else c t=pt-NT/2 c endif c elseif(pn.ge.NL/2+1) then c n=pn-NL/2 c if (pt.le.NT/2+1) then c t=NT/2+2-pt c else c t=NT+NT/2+2-pt c endif c endif do pt=1,NT/2+1 do pn=1,NL+1 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 if (m.ge.1.and.m.le.(NM/2+1)) then PEE1sum_sym(pn,pt)=PEE1sum_sym(pn,pt) #+(CABS(EE1(n,t)))**2 /float(PP*(NM/2+1)) PEE2sum_sym(pn,pt)=PEE2sum_sym(pn,pt) #+(CABS(EE2(n,t)))**2 /float(PP*(NM/2+1)) P12sum_sym(pn,pt)=P12sum_sym(pn,pt) #+REAL(CONJG(EE1(n,t))*EE2(n,t))/float(PP*(NM/2+1)) Q12sum_sym(pn,pt)=Q12sum_sym(pn,pt) #+REAL((0.,1.)*CONJG(EE1(n,t))*EE2(n,t))/float(PP*(NM/2+1)) else PEE1sum_asy(pn,pt)=PEE1sum_asy(pn,pt) #+(CABS(EE1(n,t)))**2 /float(PP*(NM/2)) PEE2sum_asy(pn,pt)=PEE2sum_asy(pn,pt) #+(CABS(EE2(n,t)))**2 /float(PP*(NM/2)) P12sum_asy(pn,pt)=P12sum_asy(pn,pt) #+REAL(CONJG(EE1(n,t))*EE2(n,t))/float(PP*(NM/2)) Q12sum_asy(pn,pt)=Q12sum_asy(pn,pt) #+REAL((0.,1.)*CONJG(EE1(n,t))*EE2(n,t))/float(PP*(NM/2)) endif enddo enddo c do pt=1,NT+1 c do pn=1,NL+1 c Coh2(pn,pt)=(P12sum(pn,pt)**2+Q12sum(pn,pt)**2)/ c # (PEE1sum(pn,pt)*PEE2sum(pn,pt)) c Phas(pn,pt)=ATAN(Q12sum(pn,pt)/P12sum(pn,pt)) c enddo c enddo enddo 2007 continue c***************END BIG LOOP HERE***************************** c write write(21,rec=1) ((PEE1sum_sym(n,t),n=1,NL+1),t=1,NT/2+1) write(21,rec=2) ((PEE2sum_sym(n,t),n=1,NL+1),t=1,NT/2+1) write(21,rec=3) ((P12sum_sym(n,t),n=1,NL+1),t=1,NT/2+1) write(21,rec=4) ((Q12sum_sym(n,t),n=1,NL+1),t=1,NT/2+1) write(22,rec=1) ((PEE1sum_asy(n,t),n=1,NL+1),t=1,NT/2+1) write(22,rec=2) ((PEE2sum_asy(n,t),n=1,NL+1),t=1,NT/2+1) write(22,rec=3) ((P12sum_asy(n,t),n=1,NL+1),t=1,NT/2+1) write(22,rec=4) ((Q12sum_asy(n,t),n=1,NL+1),t=1,NT/2+1) c ********* APLLY SMOOTHING TO THE SPECTRUM ************************ c Apply smoothing to PEE1,PEE2,P12, and Q12 before calculating Coh2 and Phase. c Smoothing in frq only do n=1, NL+1 do t=4,NT/2-1 frq(t-4+1)=PEE1sum_sym(n,t) enddo call smooth121(frq,NT/2-1,NT/2-1) do t=4,NT/2-1 PEE1sum_sym(n,t)=frq(t-4+1) enddo enddo do n=1, NL+1 do t=4,NT/2-1 frq(t-4+1)=PEE2sum_sym(n,t) enddo call smooth121(frq,NT/2-1,NT/2-1) do t=4,NT/2-1 PEE2sum_sym(n,t)=frq(t-4+1) enddo enddo do n=1, NL+1 do t=4,NT/2-1 frq(t-4+1)=P12sum_sym(n,t) enddo call smooth121(frq,NT/2-1,NT/2-1) do t=4,NT/2-1 P12sum_sym(n,t)=frq(t-4+1) enddo enddo do n=1, NL+1 do t=4,NT/2-1 frq(t-4+1)=Q12sum_sym(n,t) enddo call smooth121(frq,NT/2-1,NT/2-1) do t=4,NT/2-1 Q12sum_sym(n,t)=frq(t-4+1) enddo enddo c asymmetric do n=1, NL+1 do t=4,NT/2-1 frq(t-4+1)=PEE1sum_asy(n,t) enddo call smooth121(frq,NT/2-1,NT/2-1) do t=4,NT/2-1 PEE1sum_asy(n,t)=frq(t-4+1) enddo enddo do n=1, NL+1 do t=4,NT/2-1 frq(t-4+1)=PEE2sum_asy(n,t) enddo call smooth121(frq,NT/2-1,NT/2-1) do t=4,NT/2-1 PEE2sum_asy(n,t)=frq(t-4+1) enddo enddo do n=1, NL+1 do t=4,NT/2-1 frq(t-4+1)=P12sum_asy(n,t) enddo call smooth121(frq,NT/2-1,NT/2-1) do t=4,NT/2-1 P12sum_asy(n,t)=frq(t-4+1) enddo enddo do n=1, NL+1 do t=4,NT/2-1 frq(t-4+1)=Q12sum_asy(n,t) enddo call smooth121(frq,NT/2-1,NT/2-1) do t=4,NT/2-1 Q12sum_asy(n,t)=frq(t-4+1) enddo enddo c Calculate the Coherence-squared statistic and the phase. do t=1,NT/2+1 do n=1,NL+1 c if (PEE1sum_sym(n,t).ne.dmiss) then Coh2_sym(n,t)=(P12sum_sym(n,t)**2+Q12sum_sym(n,t)**2)/ # (PEE1sum_sym(n,t)*PEE2sum_sym(n,t)) Phas_sym(n,t)=ATAN(Q12sum_sym(n,t)/P12sum_sym(n,t)) if(n.le.NL/2+1) then v1_sym(n,t)=Q12sum_sym(n,t)/ # sqrt(Q12sum_sym(n,t)**2+P12sum_sym(n,t)**2) else v1_sym(n,t)=-1.*Q12sum_sym(n,t)/ # sqrt(Q12sum_sym(n,t)**2+P12sum_sym(n,t)**2) Q12sum_sym(n,t)=-1.*Q12sum_sym(n,t) endif v2_sym(n,t)=P12sum_sym(n,t)/ # sqrt(Q12sum_sym(n,t)**2+P12sum_sym(n,t)**2) c else c Coh2_sym(n,t)=dmiss c Phas_sym(n,t)=dmiss c endif if (Coh2_sym(n,t).lt..05) then v1_sym(n,t)=0.0 v2_sym(n,t)=0.0 endif enddo enddo do t=1,NT/2+1 do n=1,NL+1 c if (PEE1sum_asy(n,t).ne.dmiss) then Coh2_asy(n,t)=(P12sum_asy(n,t)**2+Q12sum_asy(n,t)**2)/ # (PEE1sum_asy(n,t)*PEE2sum_asy(n,t)) Phas_asy(n,t)=ATAN(Q12sum_asy(n,t)/P12sum_asy(n,t)) if(n.le.NL/2+1) then v1_asy(n,t)=Q12sum_asy(n,t)/ # sqrt(Q12sum_asy(n,t)**2+P12sum_asy(n,t)**2) else v1_asy(n,t)=-1.*Q12sum_asy(n,t)/ # sqrt(Q12sum_asy(n,t)**2+P12sum_asy(n,t)**2) Q12sum_asy(n,t)=-1.*Q12sum_asy(n,t) endif v2_asy(n,t)=P12sum_asy(n,t)/ # sqrt(Q12sum_asy(n,t)**2+P12sum_asy(n,t)**2) c else c Coh2_asy(n,t)=dmiss c Phas_asy(n,t)=dmiss c endif if (Coh2_asy(n,t).lt..05) then v1_asy(n,t)=0.0 v2_asy(n,t)=0.0 endif enddo enddo write(31,rec=1) ((Coh2_sym(n,t),n=1,NL+1),t=1,NT/2+1) write(31,rec=2) ((Phas_sym(n,t),n=1,NL+1),t=1,NT/2+1) write(31,rec=3) ((v1_sym(n,t),n=1,NL+1),t=1,NT/2+1) write(31,rec=4) ((v2_sym(n,t),n=1,NL+1),t=1,NT/2+1) write(32,rec=1) ((Coh2_asy(n,t),n=1,NL+1),t=1,NT/2+1) write(32,rec=2) ((Phas_asy(n,t),n=1,NL+1),t=1,NT/2+1) write(32,rec=3) ((v1_asy(n,t),n=1,NL+1),t=1,NT/2+1) write(32,rec=4) ((v2_asy(n,t),n=1,NL+1),t=1,NT/2+1) STOP END !main program 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=-999. 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