C**************************************************************************** C WAVETEST: Example Fortran program for WAVELET, using NINO3 SST dataset C C COMPILE: f90 chisqr.f cfftpack.f wavelet.f wavetest.f -o wavetest C C See "http://paos.colorado.edu/research/wavelets/" C C Copyright (C) 1998, Christopher Torrence and Gilbert P. Compo C This software may be used, copied, or redistributed as long as it is not C sold and this copyright notice is reproduced on each copy made. This C routine is provided as is without any express or implied warranties C whatsoever. C C Modified: November 1999 by Arjan van Dijk to include IMPLICIT NONE and C to convert all routines to DOUBLE precision. C**************************************************************************** PROGRAM wavetest IMPLICIT none INTEGER n,subscale,jtot,nprn,jprn,l,k ccccccccccc these parameters are needed for binary integer NlAT,kgit,lev,lat cccccccccccccccccccccccccccccc DOUBLE PRECISION dt,s0,dj,absw2(2),phase(2), $ perdays,c1,c2,c3,c4,c6,c7 REAL Aws,Awc,Aes,Aec,Awest,Aeast,PhaseW,PhaseE,amp_SPW, $ ampl_w,ampl_e C--------------- these parameters depend on the particular time series PARAMETER (n=1201,dt=0.0833333D0,s0=2.*dt) PARAMETER (subscale=32) PARAMETER (dj=1.D0/subscale,jtot=11*subscale) cccccccc part introduced for reading from binary files cccccccc names are as they were in the "PHI"-program PARAMETER (NlAT=36) parameter (kgit=56) real hout(2,n) cccccccccccccccccccccccccccccccccccccccccccccccccccc C C Note: for accurate reconstruction and wavelet-derived variance C do not pad with zeroes, set s0=dt (for Paul set s0=dt/4), and use C a large "jtot" (even though the extra scales will be within C the cone of influence). C For plotting purposes, it is only necessary to use C s0=2dt (for Morlet) and "jtot" from Eqn(10) Torrence&Compo(1998). INTEGER mother,ibase2,npad integer i1,i2,imax,jmax,lev1,lev2,m real tdays(n),hin(2,n) real tday0,dolg cccccccc part introduced for reading from binary files cccccccc names are as they were in the "PHI"-program REAl gh_A(nlat,kgit),gh_F(nlat,kgit) REAl gh_rec_A(nlat,kgit,n),gh_rec_F(nlat,kgit,n) REAl gh_r1_A(nlat,kgit),gh_r1_F(nlat,kgit) ccccccccccccccccccccccccccccccccccccccccccccccccccccc DOUBLE PRECISION sst(n),recon_sst(n),param,pi,www DOUBLE PRECISION scale(jtot),period(jtot),coi(n) DOUBLE COMPLEX wave(n,jtot),wvl(n,jtot,2) real A_out(3,352),Aw_ave(352) INTEGER i,j,isigtest,javg1,javg2 DOUBLE PRECISION lag1,siglvl,dof(jtot) DOUBLE PRECISION fft_theor(jtot),signif(jtot),ymean,variance DOUBLE PRECISION recon_mean,recon_vari DOUBLE PRECISION Cdelta,psi0 DOUBLE PRECISION global_ws(jtot),global_signif(jtot) DOUBLE PRECISION savg_dof(jtot),savg_signif(jtot),sstENSO(n) c------------------------------------------------------------------ character*30 wk1, filea1, filef1 pi = 4.D0*ATAN(1.D0) ibase2 = NINT(LOG(DBLE(n))/LOG(2.D0))+1 npad = INT(2.D0**ibase2) c npad = n ! this is for no padding with zeroes c#################################################### c do j=1,jtot c Aw_ave(j)=0. c end do dolg=0. c#################################################### C*************************************************** Wavelet transform C** let the WAVELET subroutine choose the defaults for these: mother = 0 param = 6.D0 c--------------------------------------------- Read input data open(1,file='wvl_rec_td_m1.txt') read(1,*)filea1 read(1,*)m,tday0,i1,i2,lev1,lev2 read(1,*)filef1 close(1) print *, m,i1,i2,n,lev1,lev2 open(2,file=filea1,form='unformatted', * access='direct', status='unknown',recl=4*nlat*kgit*2) c------------------------------------- OUTPUT FILE open(3,file=filef1,form='unformatted', * access='direct', status='unknown',recl=4*nlat*kgit*2) do 1000 lev=lev1,lev2 write(*,*) (float(lev)-float(lev1))/(float(lev2)-float(lev1)) do 1001 lat=1,NLAT CCCCCCC Input from binary files for amp's and ph's do i = 1, n read(2,rec=i) gh_A,gh_F hin(1,i)=gh_A(lat,lev)*cos(float(m)*(dolg-gh_F(lat,lev))*pi/180.) hin(2,i)=gh_A(lat,lev)*sin(float(m)*(dolg-gh_F(lat,lev))*pi/180.) tDays(i) = tday0+float(i-1)/12. end do c----------------------------------------- End of input data DO k=1,2 DO i=1,n sst(i)=hin(k,i) end do C********************************************** get the wavelet transform CALL WAVELET(n,sst,dt,mother,param,s0,dj,jtot,npad, & wave,scale,period,coi) DO i=1,n DO j=1,jtot wvl(i,j,k)=wave(i,j) END DO END DO END DO C********************************************************** print results c open (16, FILE='tp_DT_rec_grads.dx' , c $ form='unformatted', c * access='direct', status='unknown',recl=4*3*jtot) open (14, file = 'period.dat') c write(14,114) period 114 format(10e12.4) c PRINT*,' lev=',lev c PRINT*,' dt=',dt c PRINT*,' mother=',mother c PRINT*,' param=',param c PRINT*,' s0=',s0 c PRINT*,' dj=',dj c PRINT*,' jtot=',jtot c PRINT*,' npad=',npad c PRINT'(/,"Let w = wave(n/2,j)",/)' c------------------------------------------------------------------------ do 2000 i=1,n hout(1,i)=0. hout(2,i)=0. do 2001 j=1,jtot c--------------------------------------------------------- power spectrum c Absw2=ABS(wave(i,j))**2 c----------------------------------------------------- amplitude spectrum DO k=1,2 Absw2(k) =ABS(wvl(i,j,k))/sqrt(period(j)/dt) phase(k)=-ATAN2(DIMAG(wvl(i,j,k)),DBLE(wvl(i,j,k))) END DO C-------------------------------------------------reconstruction part c----------------------------------------------- if(period(j).GT.0.7.and.period(j).LT.1.25) then hout(1,i)=hout(1,i)+Absw2(1)*COS(phase(1))* $ dj/0.776*SQRT(SQRT(pi)) hout(2,i)=hout(2,i)+Absw2(2)*COS(phase(2))* $ dj/0.776*SQRT(SQRT(pi)) endif c------------------- CALCULATION OF THE WESTWARD AND EASTWARD c PROPAGATING WAVES Aws= 0.5*(Absw2(1)*cos(pi/2.-phase(1))- $ Absw2(2)*cos( phase(2))) Awc= 0.5*(Absw2(1)*cos( phase(1))+ $ Absw2(2)*cos(pi/2.-phase(2))) Aes=-0.5*(Absw2(1)*cos(pi/2.-phase(1))+ $ Absw2(2)*cos( phase(2))) Aec= 0.5*(Absw2(1)*cos( phase(1))- $ Absw2(2)*cos(pi/2.-phase(2))) Awest=SQRT(Aws*Aws+Awc*Awc) Aeast=SQRT(Aes*Aes+Aec*Aec) phaseW=180.*ATAN2(Aws,Awc)/pi if(phaseW.lt.0.)phaseW=phaseW+360. phaseE=180.*ATAN2(Aes,Aec)/pi if(phaseE.lt.0.)phaseE=phaseE+360. if (Awest.ge.Aeast) then ampl_w=Awest-Aeast ampl_e=0. amp_SPW=2.*Aeast ELSE ampl_e=Aeast-Awest ampl_w=0. amp_SPW=2.*Awest end if c A_out(1,j)=amp_SPW c A_out(2,j)=ampl_e c A_out(3,j)=ampl_w c######################################## c IF(tDays(i).ge.1.and.tDays(i).le.59) then c IF(i.ge.732.and.i.le.1440) then c Aw_ave(j)=Aw_ave(j)+ampl_w/float(nlat)/float(lev2-lev1+1)/708. c end if 2001 continue c------------------reconstruction gh_rec_A(lat,Lev,i)=SQRT(hout(1,i)*hout(1,i)+hout(2,i)*hout(2,i)) c------ the phase - position of ridge gh_rec_F(lat,Lev,i)=-ATAN2(hout(2,i),hout(1,i))*180./PI IF(gh_rec_F(lat,Lev,i).LT.0.)then gh_rec_F(lat,Lev,i)=gh_rec_F(lat,Lev,i)+360. endif gh_rec_F(lat,Lev,i)=gh_rec_F(lat,Lev,i)/FLOAT(m) 2000 continue 1001 continue 1000 continue c open(33,file='td_reconstruction.dx',form='unformatted', c * access='direct', status='unknown',recl=4*nlat*kgit*2) do i = 1, n do lev=1,kgit do lat=1,nlat gh_r1_A(lat,lev)=gh_rec_A(lat,lev,i) gh_r1_F(lat,lev)=gh_rec_F(lat,lev,i) enddo end do write(3,rec=i) gh_r1_A,gh_r1_F end do STOP END