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 ccccccccccccc aditional input parameter - number of the level integer ilev ccccccccccc these parameters are needed for binary integer NlAT,kgit 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=1456,dt=0.125D0,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=23) 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,k,ilat,m real tdays(n),hin(2,n),hout(2,n),ainput(n),finput(n),wk2(61,n) real tday0,dolg cccccccc part introduced for reading from binary files cccccccc names are as they were in the "PHI"-program REAl gh_A1(nlat,kgit),gh_F1(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#################################################### do j=1,jtot Aw_ave(j)=0. end do c#################################################### C*************************************************** Wavelet transform C** let the WAVELET subroutine choose the defaults for these: mother = 0 param = 6.D0 open(3,file='input.dat') open(4,file='out.dat') c--------------------------------------------- Read input data dolg=0. open(1,file='wvl_inp.txt') read(1,*)filea1 write(*,*) ccccc in binary files amp's and ph's are in the same file c read(1,*)filef1 cccccccccccccccccccccccccccccccccccccccc read(1,*)m,tday0,ilat,i1,i2,ilev close(1) print *, m, ilat, i1,i2, n open(2,file=filea1,form='unformatted', * access='direct', status='unknown',recl=4*nlat*kgit*2) c that's how it was for input from text-files for amp's c do i = 1, n c read(2, *, end = 100) tdays(i),(wk2(l,i),l=1,61) c ainput(i)=wk2(ilat,i) cc tDays(i) = tDay0+float(i-1)/4. cc tDays(i) = -61.+float(i) c enddo ccccccccccccccccccccccccccccccccccccccccccc CCCCCCC Input from binary files for amp's and ph's do i = 1, n read(2,rec=i) gh_A1,gh_F1 CCCCCCCCCCC filling work arrays form input arrays ainput(i)=gh_A1(ilat,ilev) finput(i)=gh_F1(ilat,ilev) tDays(i) = tday0+float(i-1)/8. end do ccccccccccccccccccccccccccccccccccccccccccccc print *, tdays close(2) 100 continue imax=i-1 do i=1,imax hin(1,i)=ainput(i)*cos(float(m)*(dolg-finput(i))*pi/180.) hin(2,i)=ainput(i)*sin(float(m)*(dolg-finput(i))*pi/180.) c hin(1,i)=10.*cos(float(m)*(dolg-30.)*pi/180.+ c $ 2.*pi*tDays(i)*24./24.) c hin(2,i)=10.*sin(float(m)*(dolg-30.)*pi/180.+ c $ 2.*pi*tDays(i)*24./24.) write(3,*)tDays(i),hin(1,i),hin(2,i) end do imax=i2-i1+1 c----------------------------------------- End of input data DO k=1,2 DO i=1,imax sst(i)=hin(k,i+i1-1) 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,imax DO j=1,jtot wvl(i,j,k)=wave(i,j) END DO END DO END DO C********************************************************** print results open (16, FILE='Dtide_grads.dx' , $ form='unformatted', * access='direct', status='unknown',recl=4*3*jtot) open (17, file = 'Dtide_aver_gh') open (14, file = 'period.dat') write(14,114) period 114 format(10e12.4) PRINT*,' dt=',dt PRINT*,' mother=',mother PRINT*,' param=',param PRINT*,' s0=',s0 PRINT*,' dj=',dj PRINT*,' jtot=',jtot PRINT*,' npad=',npad PRINT'(/,"Let w = wave(n/2,j)",/)' c------------------------------------------------------------------------ do 333 i=1,n hout(1,i)=0. hout(2,i)=0. do 444 j=1,jtot perdays=period(j) c do i=1,n c If(perdays.ge.0..and.perdays.lt.3.)then c If(perdays.gt.20.) go to 333 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) IF(Absw2(k).GT.0.) THEN phase(k)=-ATAN2(DIMAG(wvl(i,j,k)),DBLE(wvl(i,j,k))) ELSE phase(k)=0. END IF END DO c------------- RECONSTRUCTION of time series !!!!!!!!!!!!!!!!!!!!!!!!!!!! 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)) c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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) IF(Awest.GT.0.) THEN phaseW=180.*ATAN2(Aws,Awc)/pi if(phaseW.lt.0.)phaseW=phaseW+360. ELSE phaseW=400. end IF IF(Aeast.GT.0.) THEN phaseE=180.*ATAN2(Aes,Aec)/pi if(phaseE.lt.0.)phaseE=phaseE+360. ELSE phaseE=400. END IF 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 If(tdays(i+i1-1).ge.-31..and.tdays(i+i1-1).le.90.) c $ write(4,213) tdays(i+i1-1),perdays,Awest,Aeast,amp_SPW, c $ ampl_e,ampl_w A_out(1,j)=amp_SPW A_out(2,j)=ampl_e A_out(3,j)=ampl_w c######################################## IF(i.ge.489.and.i.le.736) then Aw_ave(j)=Aw_ave(j)+ampl_w/248. end if 444 continue 213 format(1x,7E11.4) write(16,rec=i) A_out write(4,*)tDays(i),hout(1,i),hout(2,i) 333 continue c end do do j=1,jtot write(17,*)period(j),Aw_ave(j) end do close(4) END