C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C c c c combined with SEIFE input and output routines. c input: a real signal; output: amplitude spectrum linear and in dB c EW March 2012. The FFTPACK routines work well only when the number c of data points is "round" (multiple prime factors 2,3,5, only one c over 5. FFT is inaccurate for large prime factors (up to 19 is OK). c c c PROGRAM seifft use nas_system C implicit real*8 (a-h,o-z) INTEGER I, N, LENSAV, IER, LENWRK, ndim PARAMETER(Ndim=1200000) C PARAMETER(LENSAV= Ndim + INT(LOG(REAL(Ndim))/LOG(2.)) + 16) PARAMETER(LENSAV= 1201000) PARAMETER(LENWRK=Ndim) dimension R(Ndim) dimension WSAVE(LENSAV), WORK(LENWRK) integer fac(15) character infile*48,text*72,pro(48)*60,datform*3,par*56,name*10 logical ncl if(iargc().eq.0) stop 'filename must be supplied as runstring' call getarg(1,infile) lin=0 n=ndim datform='fil' call inpasc(infile,lin,text,pro,npro,n,dt,tmin,tsec,r,ndim, & datform) C C --- INITIALIZE FFT C WRITE(6,*) 'PROGRAM TRFFT1 AND RELATED MESSAGES:' CALL RFFT1I (N,WSAVE,fac,LENSAV,IER) IF (IER.NE.0) THEN WRITE(6,*) 'ERROR ',IER,' IN ROUTINE RFFT1I' STOP END IF C C --- PERFORM FORWARD TRANSFORM C CALL RFFT1F (N,1,R,N, . WSAVE,fac,LENSAV,WORK,LENWRK,IER) IF (IER.NE.0) THEN WRITE(6,*) 'ERROR ',IER,' IN ROUTINE RFFT1F !' STOP END IF C WRITE(6,'(A,/)') ' END PROGRAM TRFFT1 AND RELATED MESSAGES' C --- PERFORM BACKWARD TRANSFORM C c CALL RFFT1B (N,1,R,N,WSAVE,fac,LENSAV,WORK,LENWRK,IER) c c IF (IER.NE.0) THEN c WRITE(6,*) 'ERROR ',IER,' IN ROUTINE RFFT1B !' c STOP c END IF C df=1.d0/(n*dt) par='fre' npro=npro+1 pro(npro)= & ' Multi-Radix Fourier Transform SEIFFT based on FFTPACK' npro=npro+1 pro(npro)= & ' Coefficients are monochromatic time-domain amplitudes' ncl=.false. npro=npro+1 write(pro(npro),*) ' Frequency bin is ',1.d0/(n*dt),' Hz' name='seifft.cpx' write(*,*) &'writing complex spectrum: dc, re1, im1, re2, im2 etc.' npro=npro+1 pro(npro)=' complex amplitudes: dc, re1, im1, re2, im2 etc.' call outform(name,par,text,pro,npro,ncl,n,df,tmin,tsec,r) name='seifft.amp' write(*,*) & 'writing amplitude spectrum: dc, amp1, amp2 etc.' nn=n/2 work(1)=abs(r(1)) do i=2,nn work(i)=dsqrt(r(2*i-2)**2+r(2*i-1)**2) enddo pro(npro)=' absolute amplitudes: dc, amp1, amp2 etc.' call outform(name,par,text,pro,npro,ncl,nn,df,tmin,tsec,work) name='seifft.log' write(*,*) 'writing logamp spectrum in decibels: logdc, logamp1, l &ogamp2 etc.' nn=n/2 work(1)=10.d0*dlog10(max(1d-12,r(1)**2)) do i=2,nn work(i)=10.d0*dlog10(max(1d-12,r(2*i-2)**2+r(2*i-1)**2)) enddo pro(npro)=' log. ampl. in decibels: logdc, logamp1, logamp2 etc. &' call outform(name,par,text,pro,npro,ncl,nn,df,tmin,tsec,work) open(9,file='winplot.par') write(9,7) 0,1,24.,16.,1,nn 7 format(i3,',',i3,",",f5.1,",",f5.1,",",i8,",",i8,", 1.") write(9,'(a)') name close(9) write(6,*) ' Plot-parameter file winplot.par was generated' STOP END c subroutine inpasc(name,lin,text,prolog,npro,n,dt,tmin,tsec,x,ndim, & datform) implicit real*8 (a-h,o-z) dimension x(ndim),y(60),ix(60) character iform*20,name*48,text*72,pro*80,code*4,rest*40 character prolog(48)*60,code1*6,code2*6,datform*3 character zeile*120,b,zarr(120),zeil*40 logical first c b=' ' cpv=1.000d0 do j=1,n x(j)=0.d0 enddo write(6,*) '===' write(6,*) ' Opening file ',name open(7,file=name,status='old') c read ASCII files if(datform.eq.'asc'.or.datform.eq.'sca') then write(6,22) trim(name),0,ndim,dt do i=1,lin read(7,'(a)') text write(6,*) ' header line: ',trim(text) prolog(i)=text enddo npro=lin text=name if(datform.eq.'sca') then c single-column ascii do j=1,n read(7,*,err=232,end=231) x(j) enddo 231 n=j-1 write(6,*) n,' samples read from file ',trim(name) return 232 write(6,*) 'error (maybe a non-numeric character) in data line ' & ,n,' of file ',trim(name) stop else c multi- and variable-column ascii c note: zeile is padded with blanks, therefore there is always a blank after the last digit nz=lin nval=0 do read(7,'(a)',end=131) zeile c read(zeile,'(120a)') zarr nn=0 do i=1,119 c if(zarr(i).ne.b.and.zarr(i+1).eq.b) nn=nn+1 if(zeile(i:i).ne.b.and.zeile(i+1:i+1).eq.b) nn=nn+1 enddo nz=nz+1 read(zeile,*,err=132) (y(j),j=1,nn) do j=1,nn x(nval+j)=y(j) enddo nval=nval+nn enddo 131 write(6,*) nval,' samples read from file ',trim(name) n=min(n,nval) close(7) return 132 write(6,*) 'error (maybe a non-numeric character) in data line ' & ,nz,' of file ',trim(name) stop endif c read data in SEIFE format else if(datform.eq.'fil') then read(7,'(a)') text write(6,*) ' header line: ',trim(text) npro=0 21 read(7,'(a)') pro if(pro(1:1).ne.'%'.and.pro(1:1).ne.'}') goto 20 npro=npro+1 if(npro.gt.48) stop ' cannot handle more than 48 header lines' prolog(npro)=pro(3:62) goto 21 20 read(pro,1) nn,iform,dt,tmin,tsec 1 format(i10,a20,3f10.3) write(6,22) trim(name),nn,ndim,dt 22 format(' reading file ',a,' ',i7,'/',i7,' samples',' dt=', & f10.6) if(n.lt.nn.or.nn.eq.0) then c write(6,'(" number of samples limited to",i7)') n nn=n endif n=nn if(iform(1:3).eq.'fre'.or.index(iform,'i').gt.0) then write(6,*) ' data will be read in free-format' read(7,*,err=25,end=23) (x(j),j=1,n) else read(7,iform,end=23) (x(j),j=1,n) endif c read data in MARS88 format else if(datform.eq.'mar') then read(7,'(a)') text text(61:72)=' ' read(text(15:17),*) idt dt=idt/1000.d0 read(text(53:54),*) thr read(text(56:57),*) tmin read(text(59:60),*) tsec tmin=tmin+60.d0*thr write(6,22) trim(name),n,ndim,dt do j=1,n read(7,*,err=23,end=23) x(j) enddo else if(datform.eq.'asl') then c read data in ASL format (such as written by Quanterra's Cimarron) read(7,'(a)') zeile text=zeile(1:72) write(6,*) 'header: ',trim(text) read(zeile,'(120a)') zarr nn=0 do i=1,119 if(zarr(i).ne.b.and.zarr(i+1).eq.b) then nn=nn+1 ix(nn)=i endif enddo read(zeile(1:ix(1)),'(a)') code1 read(zeile(ix(1)+2:ix(2)),'(a)') code2 read(zeile(ix(2)+2:ix(7)),*) iyear,iday,ithr,itmin,itsec read(zeile(ix(7)+2:ix(11)),*) t1000,srate8,n,cpv srate=srate8 thr=ithr tmin=itmin tsec=itsec dt=nint(1.d6/srate8)/1.d6 tmin=tmin+60.d0*ithr tsec=tsec+t1000/1000.d0 nz=1 nval=0 do read(7,'(a)',end=31) zeile c read(zeile,'(120a)') zarr nn=0 do i=1,119 c if(zarr(i).ne.b.and.zarr(i+1).eq.b) nn=nn+1 if(zeile(i:i).ne.b.and.zeile(i+1:i+1).eq.b) nn=nn+1 enddo read(zeile,*,err=32) (ix(j),j=1,nn) do j=1,nn x(nval+j)=ix(j) enddo nval=nval+nn enddo 31 write(6,*) nval,' samples read from file ',trim(name) n=min(n,nval) close(7) return 32 write(6,*) 'error (maybe a non-numeric character) in data line ' & ,nz,' of file ',trim(name) stop else if(datform.eq.'psn') then c read PSN single-column ascii with multi-line header read(7,'(a)') zeil if(zeil(1:11).ne.'! PSN ASCII') stop 'not a PSN ASCII file' do read(7,'(a)') zeil if(zeil(1:11).eq.'Start Time:') then read(zeil(24:25),*) thr read(zeil(27:28),*) tmin read(zeil(30:35),*) tsec tmin=tmin+60.d0*thr else if(zeil(1:4).eq.'SPS:') then read(zeil(6:9),*) srate8 srate=srate8 dt=nint(1.d6/srate8)/1.d6 else if(zeil(1:18).eq.'Number of Samples:') then read(zeil(20:30),*) n else if(zeil(1:5).eq.'Data:') then exit endif enddo do j=1,n read(7,*,err=233,end=234) x(j) enddo write(6,*) n,' samples read from file ',trim(name) return 234 n=j-1 write(6,*) n,' samples read from file ',trim(name) return 233 write(6,*) 'error (maybe a non-numeric character) in data line ' & ,n,' of file ',trim(name) stop else if(datform.eq.'bdf') then c read data in BDF format first=.true. nnn=0 40 read(7,'(a)',end=23) text code=text(1:4) rest=text(5:44) c write(*,*) code,rest if(code.eq.'DATA') then goto 41 else if(code.eq.'RATE') then read(rest,*) srate8 srate=srate8 dt=nint(1.d6/srate8)/1.d6 else if(first.and.code.eq.'TIME') then read(rest(1:4),*) year read(rest(6:8),*) day read(rest(10:11),*) thr read(rest(13:14),*) tmin read(rest(16:22),*) tsec tmin=tmin+60.d0*thr first=.false. else if(code.eq.'NSAM') then read(rest,*) nsam endif goto 40 41 read(7,*,err=26,end=23) (x(j),j=nnn+1,nnn+nsam) nnn=nnn+nsam goto 40 endif close(7) return 23 n=j-1 write(6,24) n 24 format(' end of file after ',i8,' samples') close(7) return 25 write(6,26) trim(name),j 26 format(' Input error reading file ',a,' at sample # ',i7) stop end subroutine outform(name,par,text,prol,npro,ncl,n,dt,tmin,tsec,x) c formatted output implicit real*8 (a-h,o-z) dimension x(n) character name*10,text*72,par*56,iform*20,fform*22 character*6 code1*6,code2*6,prol(48)*60 logical ncl if(npro.gt.48) stop ' cannot handle more than 48 header lines' iform=' ' xmax=0.d0 do 2 j=1,n 2 xmax=max(xmax,abs(x(j))) if(index(par,'lab').gt.0) then nvor=dlog10(max(xmax/cpv,2.d0))+1.d0 ndec=max(0,10-nvor) else nvor=dlog10(max(xmax,2.d0))+1.d0 ndec=max(0,10-nvor) endif i1=index(par,'(') i2=index(par,')') if(trim(par).eq.'') par(1:3)='fre' if(index(par,'fre').gt.0.) then write(iform,3) ndec 3 format('(5f13.',i1,')') else if(index(par,'lab').gt.0) then write(iform,4) ndec 4 format('(1x,f10.3,1x,f13.',i1,') ') else if(index(par,'nil').gt.0) then return else if(i1.gt.0.and.i2.gt.i1+4) then read(par(i1:i2),'(a)') iform else write(*,*) ' illegal format ',par stop endif write(6,5) trim(name),n,iform 5 format(' writing file ',a/i9,' samples in format ',a) ndf=min(8,8-floor(dlog10(dt))) write(fform,'(a13,i1,a8)') '(i10,a20,f10.',ndf,',2f10.3)' open(8,file=name) write(8,'(a)') text if(.not.ncl) write(8,'("% ",a)') (prol(j),j=1,npro) write(8,fform) n,iform,dt,tmin,tsec 1 format(i10,a20,f10.5,2f10.3) write(8,iform) (x(j),j=1,n) close(8) write(6,'(" finished!")') return 99 stop ' insufficient or wrong parameters for end command' end C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C SUBROUTINE RFFT1B ( N, INC, R, LENR, WSAVE, fac, LENSAV, 1 WORK, LENWRK, IER) implicit real*8 (a-h,o-z) INTEGER N, INC, LENR, LENSAV, LENWRK, IER dimension R(LENR), WSAVE(LENSAV) ,WORK(LENWRK) integer fac(15) C IER = 0 C IF (LENR .LT. INC*(N-1) + 1) THEN IER = 1 CALL XERFFT ('RFFT1B ', 6) ELSEIF (LENSAV .LT. N + INT(LOG(REAL(N))/LOG(2.)) +4) THEN IER = 2 CALL XERFFT ('RFFT1B ', 8) ELSEIF (LENWRK .LT. N) THEN IER = 3 CALL XERFFT ('RFFT1B ', 10) ENDIF C IF (N .EQ. 1) RETURN C CALL RFFTB1 (N,INC,R,WORK,WSAVE,fac) RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C SUBROUTINE RFFT1F ( N, INC, R, LENR, WSAVE, fac, LENSAV, 1 WORK, LENWRK, IER) implicit real*8 (a-h,o-z) INTEGER N, INC, LENR, LENSAV, LENWRK, IER dimension R(LENR), WSAVE(LENSAV), WORK(LENWRK) integer fac(15) C IER = 0 C IF (LENR .LT. INC*(N-1) + 1) THEN IER = 1 CALL XERFFT ('RFFT1F ', 6) ELSEIF (LENSAV .LT. N + INT(LOG(REAL(N))/LOG(2.)) +4) THEN IER = 2 CALL XERFFT ('RFFT1F ', 8) ELSEIF (LENWRK .LT. N) THEN IER = 3 CALL XERFFT ('RFFT1F ', 10) ENDIF C IF (N .EQ. 1) RETURN C CALL RFFTF1 (N,INC,R,WORK,WSAVE,fac) RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C SUBROUTINE RFFT1I ( N, WSAVE,fac, LENSAV, IER ) implicit real*8 (a-h,o-z) INTEGER N, LENSAV, IER dimension WSAVE(LENSAV) integer fac(15) C IER = 0 C IF (LENSAV .LT. N + INT(LOG(REAL(N))/LOG(2.)) +4) THEN IER = 2 CALL XERFFT ('RFFT1I ', 3) ENDIF C IF (N .EQ. 1) RETURN C CALL RFFTI1 (N,WSAVE,fac) RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C SUBROUTINE RFFTB1 (N,IN,C,CH,WA,FAC) implicit real*8 (a-h,o-z) dimension CH(*), C(IN,*), WA(N) integer FAC(15) C NF = FAC(2) NA = 0 DO 10 K1=1,NF IP = FAC(K1+2) NA = 1-NA IF(IP .LE. 5) GO TO 10 IF(K1 .EQ. NF) GO TO 10 NA = 1-NA 10 CONTINUE HALF = .5 HALFM = -.5 MODN = MOD(N,2) NL = N-2 IF(MODN .NE. 0) NL = N-1 IF (NA .EQ. 0) GO TO 120 CH(1) = C(1,1) CH(N) = C(1,N) DO 118 J=2,NL,2 CH(J) = HALF*C(1,J) CH(J+1) = HALFM*C(1,J+1) 118 CONTINUE GO TO 124 120 DO 122 J=2,NL,2 C(1,J) = HALF*C(1,J) C(1,J+1) = HALFM*C(1,J+1) 122 CONTINUE 124 L1 = 1 IW = 1 DO 116 K1=1,NF IP = FAC(K1+2) L2 = IP*L1 IDO = N/L2 IDL1 = IDO*L1 IF (IP .NE. 4) GO TO 103 IX2 = IW+IDO IX3 = IX2+IDO IF (NA .NE. 0) GO TO 101 CALL R1F4KB (IDO,L1,C,IN,CH,1,WA(IW),WA(IX2),WA(IX3)) GO TO 102 101 CALL R1F4KB (IDO,L1,CH,1,C,IN,WA(IW),WA(IX2),WA(IX3)) 102 NA = 1-NA GO TO 115 103 IF (IP .NE. 2) GO TO 106 IF (NA .NE. 0) GO TO 104 CALL R1F2KB (IDO,L1,C,IN,CH,1,WA(IW)) GO TO 105 104 CALL R1F2KB (IDO,L1,CH,1,C,IN,WA(IW)) 105 NA = 1-NA GO TO 115 106 IF (IP .NE. 3) GO TO 109 IX2 = IW+IDO IF (NA .NE. 0) GO TO 107 CALL R1F3KB (IDO,L1,C,IN,CH,1,WA(IW),WA(IX2)) GO TO 108 107 CALL R1F3KB (IDO,L1,CH,1,C,IN,WA(IW),WA(IX2)) 108 NA = 1-NA GO TO 115 109 IF (IP .NE. 5) GO TO 112 IX2 = IW+IDO IX3 = IX2+IDO IX4 = IX3+IDO IF (NA .NE. 0) GO TO 110 CALL R1F5KB (IDO,L1,C,IN,CH,1,WA(IW),WA(IX2), 1 WA(IX3),WA(IX4)) GO TO 111 110 CALL R1F5KB (IDO,L1,CH,1,C,IN,WA(IW),WA(IX2), 1 WA(IX3),WA(IX4)) 111 NA = 1-NA GO TO 115 112 IF (NA .NE. 0) GO TO 113 CALL R1FGKB (IDO,IP,L1,IDL1,C,C,C,IN,CH,CH,1,WA(IW)) GO TO 114 113 CALL R1FGKB (IDO,IP,L1,IDL1,CH,CH,CH,1,C,C,IN,WA(IW)) 114 IF (IDO .EQ. 1) NA = 1-NA 115 L1 = L2 IW = IW+(IP-1)*IDO 116 CONTINUE RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C SUBROUTINE RFFTF1 (N,IN,C,CH,WA,FAC) implicit real*8 (a-h,o-z) dimension CH(*) ,C(IN,*) ,WA(N) integer FAC(15) C NF = FAC(2) NA = 1 L2 = N IW = N DO 111 K1=1,NF KH = NF-K1 IP = FAC(KH+3) L1 = L2/IP IDO = N/L2 IDL1 = IDO*L1 IW = IW-(IP-1)*IDO NA = 1-NA IF (IP .NE. 4) GO TO 102 IX2 = IW+IDO IX3 = IX2+IDO IF (NA .NE. 0) GO TO 101 CALL R1F4KF (IDO,L1,C,IN,CH,1,WA(IW),WA(IX2),WA(IX3)) GO TO 110 101 CALL R1F4KF (IDO,L1,CH,1,C,IN,WA(IW),WA(IX2),WA(IX3)) GO TO 110 102 IF (IP .NE. 2) GO TO 104 IF (NA .NE. 0) GO TO 103 CALL R1F2KF (IDO,L1,C,IN,CH,1,WA(IW)) GO TO 110 103 CALL R1F2KF (IDO,L1,CH,1,C,IN,WA(IW)) GO TO 110 104 IF (IP .NE. 3) GO TO 106 IX2 = IW+IDO IF (NA .NE. 0) GO TO 105 CALL R1F3KF (IDO,L1,C,IN,CH,1,WA(IW),WA(IX2)) GO TO 110 105 CALL R1F3KF (IDO,L1,CH,1,C,IN,WA(IW),WA(IX2)) GO TO 110 106 IF (IP .NE. 5) GO TO 108 IX2 = IW+IDO IX3 = IX2+IDO IX4 = IX3+IDO IF (NA .NE. 0) GO TO 107 CALL R1F5KF (IDO,L1,C,IN,CH,1,WA(IW),WA(IX2), 1 WA(IX3),WA(IX4)) GO TO 110 107 CALL R1F5KF (IDO,L1,CH,1,C,IN,WA(IW),WA(IX2), 1 WA(IX3),WA(IX4)) GO TO 110 108 IF (IDO .EQ. 1) NA = 1-NA IF (NA .NE. 0) GO TO 109 CALL R1FGKF (IDO,IP,L1,IDL1,C,C,C,IN,CH,CH,1,WA(IW)) NA = 1 GO TO 110 109 CALL R1FGKF (IDO,IP,L1,IDL1,CH,CH,CH,1,C,C,IN,WA(IW)) NA = 0 110 L2 = L1 111 CONTINUE SN = 1.d0/N TSN = 2.d0/N TSNM = -TSN MODN = MOD(N,2) NL = N-2 IF(MODN .NE. 0) NL = N-1 IF (NA .NE. 0) GO TO 120 C(1,1) = SN*CH(1) DO 118 J=2,NL,2 C(1,J) = TSN*CH(J) C(1,J+1) = TSNM*CH(J+1) 118 CONTINUE IF(MODN .NE. 0) RETURN C(1,N) = SN*CH(N) RETURN 120 C(1,1) = SN*C(1,1) DO 122 J=2,NL,2 C(1,J) = TSN*C(1,J) C(1,J+1) = TSNM*C(1,J+1) 122 CONTINUE IF(MODN .NE. 0) RETURN C(1,N) = SN*C(1,N) RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C SUBROUTINE RFFTI1 (N,WA,FAC) implicit real*8 (a-h,o-z) dimension WA(N) integer FAC(15) INTEGER NTRYH(4) DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/4,2,3,5/ C NL = N NF = 0 J = 0 101 J = J+1 IF (J-4) 102,102,103 102 NTRY = NTRYH(J) GO TO 104 103 NTRY = NTRY+2 104 NQ = NL/NTRY NR = NL-NTRY*NQ IF (NR) 101,105,101 105 NF = NF+1 if(nf.gt.13) goto 111 FAC(NF+2) = NTRY NL = NQ IF (NTRY .NE. 2) GO TO 107 IF (NF .EQ. 1) GO TO 107 DO 106 I=2,NF IB = NF-I+2 if(ib.gt.13) goto 111 FAC(IB+2) = FAC(IB+1) 106 CONTINUE FAC(3) = 2 107 IF (NL .NE. 1) GO TO 104 FAC(1) = N FAC(2) = NF TPI = 8.D0*DATAN(1.D0) ARGH = TPI/dble(N) IS = 0 NFM1 = NF-1 L1 = 1 IF (NFM1 .EQ. 0) RETURN DO 110 K1=1,NFM1 if(k1.gt.13) goto 111 IP = FAC(K1+2) LD = 0 L2 = L1*IP IDO = N/L2 IPM = IP-1 DO 109 J=1,IPM LD = LD+L1 I = IS ARGLD = dble(LD)*ARGH FI = 0.d0 DO 108 II=3,IDO,2 I = I+2 FI = FI+1.d0 ARG = FI*ARGLD WA(I-1) = DCOS(ARG) WA(I) = DSIN(ARG) 108 CONTINUE IS = IS+IDO 109 CONTINUE L1 = L2 110 CONTINUE write(*,*) n," samples, ",nf," factors (prime or 4)" write(*,*) write(*,'(5i8)') (fac(i),i=3,nf+2) write(*,*) if(nf.gt.1.and.fac(nf+1).gt.5.and.fac(nf+2).gt.5) &write(*,*) " sorry, this program cannot handle two factors larger & than 5" RETURN 111 stop 'too many prime factors - stop' END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C SUBROUTINE XERFFT( SRNAME, INFO) C C .. Scalar Arguments .. CHARACTER*6 SRNAME INTEGER INFO C C .. C C Purpose C ======= C C XERFFT is an error handler for library FFTPACK version 5.1 routines. C It is called by an FFTPACK 5.1 routine if an input parameter has an C invalid value. A message is printed and execution stops. C C Installers may consider modifying the STOP statement in order to C call system-specific exception-handling facilities. C C Arguments C ========= C C SRNAME (input) CHARACTER*6 C The name of the routine which called XERFFT. C C INFO (input) INTEGER C When a single invalid parameter in the parameter list of C the calling routine has been detected, INFO is the position C of that parameter. In the case when an illegal combination C of LOT, JUMP, N, and INC has been detected, the calling C subprogram calls XERFFT with INFO = -1. C C ===================================================================== C C .. Executable Statements .. C IF (INFO .GE. 1) THEN WRITE( *, '(A,A,A,I3,A)') ' ** On entry to ', SRNAME, 1 ' parameter number ', INFO, ' had an illegal value' ELSEIF (INFO .EQ. -1) THEN WRITE( *, '(A,A,A,A)') ' ** On entry to ', SRNAME, 1 ' parameters LOT, JUMP, N and INC are inconsistent' ELSEIF (INFO .EQ. -2) THEN WRITE( *, '(A,A,A,A)') ' ** On entry to ', SRNAME, 1 ' parameter L is greater than LDIM' ELSEIF (INFO .EQ. -3) THEN WRITE( *, '(A,A,A,A)') ' ** On entry to ', SRNAME, 1 ' parameter M is greater than MDIM' ELSEIF (INFO .EQ. -5) THEN WRITE( *, '(A,A,A,A)') ' ** Within ', SRNAME, 1 ' input error returned by lower level routine' ELSEIF (INFO .EQ. -6) THEN WRITE( *, '(A,A,A,A)') ' ** On entry to ', SRNAME, 1 ' parameter LDIM is less than 2*(L/2+1)' ENDIF C STOP C C End of XERFFT C END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C SUBROUTINE R1F2KB (IDO,L1,CC,IN1,CH,IN2,WA1) implicit real*8 (a-h,o-z) dimension CC(IN1,IDO,2,L1), CH(IN2,IDO,L1,2), WA1(IDO) C DO 101 K=1,L1 CH(1,1,K,1) = CC(1,1,1,K)+CC(1,IDO,2,K) CH(1,1,K,2) = CC(1,1,1,K)-CC(1,IDO,2,K) 101 CONTINUE IF (IDO-2) 107,105,102 102 IDP2 = IDO+2 DO 104 K=1,L1 DO 103 I=3,IDO,2 IC = IDP2-I CH(1,I-1,K,1) = CC(1,I-1,1,K)+CC(1,IC-1,2,K) CH(1,I,K,1) = CC(1,I,1,K)-CC(1,IC,2,K) CH(1,I-1,K,2) = WA1(I-2)*(CC(1,I-1,1,K)-CC(1,IC-1,2,K)) 1 -WA1(I-1)*(CC(1,I,1,K)+CC(1,IC,2,K)) CH(1,I,K,2) = WA1(I-2)*(CC(1,I,1,K)+CC(1,IC,2,K))+WA1(I-1) 1 *(CC(1,I-1,1,K)-CC(1,IC-1,2,K)) 103 CONTINUE 104 CONTINUE IF (MOD(IDO,2) .EQ. 1) RETURN 105 DO 106 K=1,L1 CH(1,IDO,K,1) = CC(1,IDO,1,K)+CC(1,IDO,1,K) CH(1,IDO,K,2) = -(CC(1,1,2,K)+CC(1,1,2,K)) 106 CONTINUE 107 RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C SUBROUTINE R1F3KB (IDO,L1,CC,IN1,CH,IN2,WA1,WA2) implicit real*8 (a-h,o-z) dimension CC(IN1,IDO,3,L1) ,CH(IN2,IDO,L1,3), 1 WA1(IDO) ,WA2(IDO) C ARG=8.d0*datan(1.0d0)/3.d0 TAUR=dcos(ARG) TAUI=dsin(ARG) DO 101 K=1,L1 CH(1,1,K,1) = CC(1,1,1,K)+2.d0*CC(1,IDO,2,K) CH(1,1,K,2) = CC(1,1,1,K)+(2.d0*TAUR)*CC(1,IDO,2,K) 1 -(2.d0*TAUI)*CC(1,1,3,K) CH(1,1,K,3) = CC(1,1,1,K)+(2.d0*TAUR)*CC(1,IDO,2,K) 1 +2.d0*TAUI*CC(1,1,3,K) 101 CONTINUE IF (IDO .EQ. 1) RETURN IDP2 = IDO+2 DO 103 K=1,L1 DO 102 I=3,IDO,2 IC = IDP2-I CH(1,I-1,K,1) = CC(1,I-1,1,K)+(CC(1,I-1,3,K)+CC(1,IC-1,2,K)) CH(1,I,K,1) = CC(1,I,1,K)+(CC(1,I,3,K)-CC(1,IC,2,K)) CH(1,I-1,K,2) = WA1(I-2)* 1 ((CC(1,I-1,1,K)+TAUR*(CC(1,I-1,3,K)+CC(1,IC-1,2,K)))- * (TAUI*(CC(1,I,3,K)+CC(1,IC,2,K)))) 2 -WA1(I-1)* 3 ((CC(1,I,1,K)+TAUR*(CC(1,I,3,K)-CC(1,IC,2,K)))+ * (TAUI*(CC(1,I-1,3,K)-CC(1,IC-1,2,K)))) CH(1,I,K,2) = WA1(I-2)* 4 ((CC(1,I,1,K)+TAUR*(CC(1,I,3,K)-CC(1,IC,2,K)))+ 8 (TAUI*(CC(1,I-1,3,K)-CC(1,IC-1,2,K)))) 5 +WA1(I-1)* 6 ((CC(1,I-1,1,K)+TAUR*(CC(1,I-1,3,K)+CC(1,IC-1,2,K)))- 8 (TAUI*(CC(1,I,3,K)+CC(1,IC,2,K)))) CH(1,I-1,K,3) = WA2(I-2)* 7 ((CC(1,I-1,1,K)+TAUR*(CC(1,I-1,3,K)+CC(1,IC-1,2,K)))+ 8 (TAUI*(CC(1,I,3,K)+CC(1,IC,2,K)))) 8 -WA2(I-1)* 9 ((CC(1,I,1,K)+TAUR*(CC(1,I,3,K)-CC(1,IC,2,K)))- 8 (TAUI*(CC(1,I-1,3,K)-CC(1,IC-1,2,K)))) CH(1,I,K,3) = WA2(I-2)* 1 ((CC(1,I,1,K)+TAUR*(CC(1,I,3,K)-CC(1,IC,2,K)))- 8 (TAUI*(CC(1,I-1,3,K)-CC(1,IC-1,2,K)))) 2 +WA2(I-1)* 3 ((CC(1,I-1,1,K)+TAUR*(CC(1,I-1,3,K)+CC(1,IC-1,2,K)))+ 8 (TAUI*(CC(1,I,3,K)+CC(1,IC,2,K)))) 102 CONTINUE 103 CONTINUE RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C SUBROUTINE R1F4KB (IDO,L1,CC,IN1,CH,IN2,WA1,WA2,WA3) implicit real*8 (a-h,o-z) dimension CC(IN1,IDO,4,L1) ,CH(IN2,IDO,L1,4) , 1 WA1(IDO) , WA2(IDO) , WA3(IDO) C SQRT2=dSQRT(2.d0) DO 101 K=1,L1 CH(1,1,K,3) = (CC(1,1,1,K)+CC(1,IDO,4,K)) 1 -(CC(1,IDO,2,K)+CC(1,IDO,2,K)) CH(1,1,K,1) = (CC(1,1,1,K)+CC(1,IDO,4,K)) 1 +(CC(1,IDO,2,K)+CC(1,IDO,2,K)) CH(1,1,K,4) = (CC(1,1,1,K)-CC(1,IDO,4,K)) 1 +(CC(1,1,3,K)+CC(1,1,3,K)) CH(1,1,K,2) = (CC(1,1,1,K)-CC(1,IDO,4,K)) 1 -(CC(1,1,3,K)+CC(1,1,3,K)) 101 CONTINUE IF (IDO-2) 107,105,102 102 IDP2 = IDO+2 DO 104 K=1,L1 DO 103 I=3,IDO,2 IC = IDP2-I CH(1,I-1,K,1) = (CC(1,I-1,1,K)+CC(1,IC-1,4,K)) 1 +(CC(1,I-1,3,K)+CC(1,IC-1,2,K)) CH(1,I,K,1) = (CC(1,I,1,K)-CC(1,IC,4,K)) 1 +(CC(1,I,3,K)-CC(1,IC,2,K)) CH(1,I-1,K,2)=WA1(I-2)*((CC(1,I-1,1,K)-CC(1,IC-1,4,K)) 1 -(CC(1,I,3,K)+CC(1,IC,2,K)))-WA1(I-1) 1 *((CC(1,I,1,K)+CC(1,IC,4,K))+(CC(1,I-1,3,K)-CC(1,IC-1,2,K))) CH(1,I,K,2)=WA1(I-2)*((CC(1,I,1,K)+CC(1,IC,4,K)) 1 +(CC(1,I-1,3,K)-CC(1,IC-1,2,K)))+WA1(I-1) 1 *((CC(1,I-1,1,K)-CC(1,IC-1,4,K))-(CC(1,I,3,K)+CC(1,IC,2,K))) CH(1,I-1,K,3)=WA2(I-2)*((CC(1,I-1,1,K)+CC(1,IC-1,4,K)) 1 -(CC(1,I-1,3,K)+CC(1,IC-1,2,K)))-WA2(I-1) 1 *((CC(1,I,1,K)-CC(1,IC,4,K))-(CC(1,I,3,K)-CC(1,IC,2,K))) CH(1,I,K,3)=WA2(I-2)*((CC(1,I,1,K)-CC(1,IC,4,K)) 1 -(CC(1,I,3,K)-CC(1,IC,2,K)))+WA2(I-1) 1 *((CC(1,I-1,1,K)+CC(1,IC-1,4,K))-(CC(1,I-1,3,K) 1 +CC(1,IC-1,2,K))) CH(1,I-1,K,4)=WA3(I-2)*((CC(1,I-1,1,K)-CC(1,IC-1,4,K)) 1 +(CC(1,I,3,K)+CC(1,IC,2,K)))-WA3(I-1) 1 *((CC(1,I,1,K)+CC(1,IC,4,K))-(CC(1,I-1,3,K)-CC(1,IC-1,2,K))) CH(1,I,K,4)=WA3(I-2)*((CC(1,I,1,K)+CC(1,IC,4,K)) 1 -(CC(1,I-1,3,K)-CC(1,IC-1,2,K)))+WA3(I-1) 1 *((CC(1,I-1,1,K)-CC(1,IC-1,4,K))+(CC(1,I,3,K)+CC(1,IC,2,K))) 103 CONTINUE 104 CONTINUE IF (MOD(IDO,2) .EQ. 1) RETURN 105 CONTINUE DO 106 K=1,L1 CH(1,IDO,K,1) = (CC(1,IDO,1,K)+CC(1,IDO,3,K)) 1 +(CC(1,IDO,1,K)+CC(1,IDO,3,K)) CH(1,IDO,K,2) = SQRT2*((CC(1,IDO,1,K)-CC(1,IDO,3,K)) 1 -(CC(1,1,2,K)+CC(1,1,4,K))) CH(1,IDO,K,3) = (CC(1,1,4,K)-CC(1,1,2,K)) 1 +(CC(1,1,4,K)-CC(1,1,2,K)) CH(1,IDO,K,4) = -SQRT2*((CC(1,IDO,1,K)-CC(1,IDO,3,K)) 1 +(CC(1,1,2,K)+CC(1,1,4,K))) 106 CONTINUE 107 RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C SUBROUTINE R1F5KB (IDO,L1,CC,IN1,CH,IN2, 1 WA1,WA2,WA3,WA4) implicit real*8 (a-h,o-z) dimension CC(IN1,IDO,5,L1) ,CH(IN2,IDO,L1,5), 1 WA1(IDO) ,WA2(IDO) ,WA3(IDO) ,WA4(IDO) C ARG=8.d0*datan(1.0d0)/5.d0 TR11=dcos(ARG) TI11=dsin(ARG) TR12=dcos(2.d0*ARG) TI12=dsin(2.d0*ARG) DO 101 K=1,L1 CH(1,1,K,1) = CC(1,1,1,K)+2.d0*CC(1,IDO,2,K)+2.d0*CC(1,IDO,4,K) CH(1,1,K,2) = (CC(1,1,1,K)+TR11*2.d0*CC(1,IDO,2,K) 1 +TR12*2.d0*CC(1,IDO,4,K))-(TI11*2.d0*CC(1,1,3,K) 1 +TI12*2.d0*CC(1,1,5,K)) CH(1,1,K,3) = (CC(1,1,1,K)+TR12*2.d0*CC(1,IDO,2,K) 1 +TR11*2.d0*CC(1,IDO,4,K))-(TI12*2.d0*CC(1,1,3,K) 1 -TI11*2.d0*CC(1,1,5,K)) CH(1,1,K,4) = (CC(1,1,1,K)+TR12*2.d0*CC(1,IDO,2,K) 1 +TR11*2.d0*CC(1,IDO,4,K))+(TI12*2.d0*CC(1,1,3,K) 1 -TI11*2.d0*CC(1,1,5,K)) CH(1,1,K,5) = (CC(1,1,1,K)+TR11*2.d0*CC(1,IDO,2,K) 1 +TR12*2.d0*CC(1,IDO,4,K))+(TI11*2.d0*CC(1,1,3,K) 1 +TI12*2.d0*CC(1,1,5,K)) 101 CONTINUE IF (IDO .EQ. 1) RETURN IDP2 = IDO+2 DO 103 K=1,L1 DO 102 I=3,IDO,2 IC = IDP2-I CH(1,I-1,K,1) = CC(1,I-1,1,K)+(CC(1,I-1,3,K)+CC(1,IC-1,2,K)) 1 +(CC(1,I-1,5,K)+CC(1,IC-1,4,K)) CH(1,I,K,1) = CC(1,I,1,K)+(CC(1,I,3,K)-CC(1,IC,2,K)) 1 +(CC(1,I,5,K)-CC(1,IC,4,K)) CH(1,I-1,K,2) = WA1(I-2)*((CC(1,I-1,1,K)+TR11* 1 (CC(1,I-1,3,K)+CC(1,IC-1,2,K))+TR12 1 *(CC(1,I-1,5,K)+CC(1,IC-1,4,K)))-(TI11*(CC(1,I,3,K) 1 +CC(1,IC,2,K))+TI12*(CC(1,I,5,K)+CC(1,IC,4,K)))) 1 -WA1(I-1)*((CC(1,I,1,K)+TR11*(CC(1,I,3,K)-CC(1,IC,2,K)) 1 +TR12*(CC(1,I,5,K)-CC(1,IC,4,K)))+(TI11*(CC(1,I-1,3,K) 1 -CC(1,IC-1,2,K))+TI12*(CC(1,I-1,5,K)-CC(1,IC-1,4,K)))) CH(1,I,K,2) = WA1(I-2)*((CC(1,I,1,K)+TR11*(CC(1,I,3,K) 1 -CC(1,IC,2,K))+TR12*(CC(1,I,5,K)-CC(1,IC,4,K))) 1 +(TI11*(CC(1,I-1,3,K)-CC(1,IC-1,2,K))+TI12 1 *(CC(1,I-1,5,K)-CC(1,IC-1,4,K))))+WA1(I-1) 1 *((CC(1,I-1,1,K)+TR11*(CC(1,I-1,3,K) 1 +CC(1,IC-1,2,K))+TR12*(CC(1,I-1,5,K)+CC(1,IC-1,4,K))) 1 -(TI11*(CC(1,I,3,K)+CC(1,IC,2,K))+TI12 1 *(CC(1,I,5,K)+CC(1,IC,4,K)))) CH(1,I-1,K,3) = WA2(I-2) 1 *((CC(1,I-1,1,K)+TR12*(CC(1,I-1,3,K)+CC(1,IC-1,2,K)) 1 +TR11*(CC(1,I-1,5,K)+CC(1,IC-1,4,K)))-(TI12*(CC(1,I,3,K) 1 +CC(1,IC,2,K))-TI11*(CC(1,I,5,K)+CC(1,IC,4,K)))) 1 -WA2(I-1) 1 *((CC(1,I,1,K)+TR12*(CC(1,I,3,K)- 1 CC(1,IC,2,K))+TR11*(CC(1,I,5,K)-CC(1,IC,4,K))) 1 +(TI12*(CC(1,I-1,3,K)-CC(1,IC-1,2,K))-TI11 1 *(CC(1,I-1,5,K)-CC(1,IC-1,4,K)))) CH(1,I,K,3) = WA2(I-2) 1 *((CC(1,I,1,K)+TR12*(CC(1,I,3,K)- 1 CC(1,IC,2,K))+TR11*(CC(1,I,5,K)-CC(1,IC,4,K))) 1 +(TI12*(CC(1,I-1,3,K)-CC(1,IC-1,2,K))-TI11 1 *(CC(1,I-1,5,K)-CC(1,IC-1,4,K)))) 1 +WA2(I-1) 1 *((CC(1,I-1,1,K)+TR12*(CC(1,I-1,3,K)+CC(1,IC-1,2,K)) 1 +TR11*(CC(1,I-1,5,K)+CC(1,IC-1,4,K)))-(TI12*(CC(1,I,3,K) 1 +CC(1,IC,2,K))-TI11*(CC(1,I,5,K)+CC(1,IC,4,K)))) CH(1,I-1,K,4) = WA3(I-2) 1 *((CC(1,I-1,1,K)+TR12*(CC(1,I-1,3,K)+CC(1,IC-1,2,K)) 1 +TR11*(CC(1,I-1,5,K)+CC(1,IC-1,4,K)))+(TI12*(CC(1,I,3,K) 1 +CC(1,IC,2,K))-TI11*(CC(1,I,5,K)+CC(1,IC,4,K)))) 1 -WA3(I-1) 1 *((CC(1,I,1,K)+TR12*(CC(1,I,3,K)- 1 CC(1,IC,2,K))+TR11*(CC(1,I,5,K)-CC(1,IC,4,K))) 1 -(TI12*(CC(1,I-1,3,K)-CC(1,IC-1,2,K))-TI11 1 *(CC(1,I-1,5,K)-CC(1,IC-1,4,K)))) CH(1,I,K,4) = WA3(I-2) 1 *((CC(1,I,1,K)+TR12*(CC(1,I,3,K)- 1 CC(1,IC,2,K))+TR11*(CC(1,I,5,K)-CC(1,IC,4,K))) 1 -(TI12*(CC(1,I-1,3,K)-CC(1,IC-1,2,K))-TI11 1 *(CC(1,I-1,5,K)-CC(1,IC-1,4,K)))) 1 +WA3(I-1) 1 *((CC(1,I-1,1,K)+TR12*(CC(1,I-1,3,K)+CC(1,IC-1,2,K)) 1 +TR11*(CC(1,I-1,5,K)+CC(1,IC-1,4,K)))+(TI12*(CC(1,I,3,K) 1 +CC(1,IC,2,K))-TI11*(CC(1,I,5,K)+CC(1,IC,4,K)))) CH(1,I-1,K,5) = WA4(I-2) 1 *((CC(1,I-1,1,K)+TR11*(CC(1,I-1,3,K)+CC(1,IC-1,2,K)) 1 +TR12*(CC(1,I-1,5,K)+CC(1,IC-1,4,K)))+(TI11*(CC(1,I,3,K) 1 +CC(1,IC,2,K))+TI12*(CC(1,I,5,K)+CC(1,IC,4,K)))) 1 -WA4(I-1) 1 *((CC(1,I,1,K)+TR11*(CC(1,I,3,K)-CC(1,IC,2,K)) 1 +TR12*(CC(1,I,5,K)-CC(1,IC,4,K)))-(TI11*(CC(1,I-1,3,K) 1 -CC(1,IC-1,2,K))+TI12*(CC(1,I-1,5,K)-CC(1,IC-1,4,K)))) CH(1,I,K,5) = WA4(I-2) 1 *((CC(1,I,1,K)+TR11*(CC(1,I,3,K)-CC(1,IC,2,K)) 1 +TR12*(CC(1,I,5,K)-CC(1,IC,4,K)))-(TI11*(CC(1,I-1,3,K) 1 -CC(1,IC-1,2,K))+TI12*(CC(1,I-1,5,K)-CC(1,IC-1,4,K)))) 1 +WA4(I-1) 1 *((CC(1,I-1,1,K)+TR11*(CC(1,I-1,3,K)+CC(1,IC-1,2,K)) 1 +TR12*(CC(1,I-1,5,K)+CC(1,IC-1,4,K)))+(TI11*(CC(1,I,3,K) 1 +CC(1,IC,2,K))+TI12*(CC(1,I,5,K)+CC(1,IC,4,K)))) 102 CONTINUE 103 CONTINUE RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C SUBROUTINE R1FGKB (IDO,IP,L1,IDL1,CC,C1,C2,IN1, 1 CH,CH2,IN2,WA) implicit real*8 (a-h,o-z) dimension CH(IN2,IDO,L1,IP) ,CC(IN1,IDO,IP,L1) , 1 C1(IN1,IDO,L1,IP) ,C2(IN1,IDL1,IP), 2 CH2(IN2,IDL1,IP) ,WA(IDO) C TPI=8.d0*datan(1.0d0) ARG = TPI/dble(IP) DCP = dcos(ARG) DSP = dsin(ARG) IDP2 = IDO+2 NBD = (IDO-1)/2 IPP2 = IP+2 IPPH = (IP+1)/2 IF (IDO .LT. L1) GO TO 103 DO 102 K=1,L1 DO 101 I=1,IDO CH(1,I,K,1) = CC(1,I,1,K) 101 CONTINUE 102 CONTINUE GO TO 106 103 DO 105 I=1,IDO DO 104 K=1,L1 CH(1,I,K,1) = CC(1,I,1,K) 104 CONTINUE 105 CONTINUE 106 DO 108 J=2,IPPH JC = IPP2-J J2 = J+J DO 107 K=1,L1 CH(1,1,K,J) = CC(1,IDO,J2-2,K)+CC(1,IDO,J2-2,K) CH(1,1,K,JC) = CC(1,1,J2-1,K)+CC(1,1,J2-1,K) 1007 CONTINUE 107 CONTINUE 108 CONTINUE IF (IDO .EQ. 1) GO TO 116 IF (NBD .LT. L1) GO TO 112 DO 111 J=2,IPPH JC = IPP2-J DO 110 K=1,L1 DO 109 I=3,IDO,2 IC = IDP2-I CH(1,I-1,K,J) = CC(1,I-1,2*J-1,K)+CC(1,IC-1,2*J-2,K) CH(1,I-1,K,JC) = CC(1,I-1,2*J-1,K)-CC(1,IC-1,2*J-2,K) CH(1,I,K,J) = CC(1,I,2*J-1,K)-CC(1,IC,2*J-2,K) CH(1,I,K,JC) = CC(1,I,2*J-1,K)+CC(1,IC,2*J-2,K) 109 CONTINUE 110 CONTINUE 111 CONTINUE GO TO 116 112 DO 115 J=2,IPPH JC = IPP2-J DO 114 I=3,IDO,2 IC = IDP2-I DO 113 K=1,L1 CH(1,I-1,K,J) = CC(1,I-1,2*J-1,K)+CC(1,IC-1,2*J-2,K) CH(1,I-1,K,JC) = CC(1,I-1,2*J-1,K)-CC(1,IC-1,2*J-2,K) CH(1,I,K,J) = CC(1,I,2*J-1,K)-CC(1,IC,2*J-2,K) CH(1,I,K,JC) = CC(1,I,2*J-1,K)+CC(1,IC,2*J-2,K) 113 CONTINUE 114 CONTINUE 115 CONTINUE 116 AR1 = 1.d0 AI1 = 0.d0 DO 120 L=2,IPPH LC = IPP2-L AR1H = DCP*AR1-DSP*AI1 AI1 = DCP*AI1+DSP*AR1 AR1 = AR1H DO 117 IK=1,IDL1 C2(1,IK,L) = CH2(1,IK,1)+AR1*CH2(1,IK,2) C2(1,IK,LC) = AI1*CH2(1,IK,IP) 117 CONTINUE DC2 = AR1 DS2 = AI1 AR2 = AR1 AI2 = AI1 DO 119 J=3,IPPH JC = IPP2-J AR2H = DC2*AR2-DS2*AI2 AI2 = DC2*AI2+DS2*AR2 AR2 = AR2H DO 118 IK=1,IDL1 C2(1,IK,L) = C2(1,IK,L)+AR2*CH2(1,IK,J) C2(1,IK,LC) = C2(1,IK,LC)+AI2*CH2(1,IK,JC) 118 CONTINUE 119 CONTINUE 120 CONTINUE DO 122 J=2,IPPH DO 121 IK=1,IDL1 CH2(1,IK,1) = CH2(1,IK,1)+CH2(1,IK,J) 121 CONTINUE 122 CONTINUE DO 124 J=2,IPPH JC = IPP2-J DO 123 K=1,L1 CH(1,1,K,J) = C1(1,1,K,J)-C1(1,1,K,JC) CH(1,1,K,JC) = C1(1,1,K,J)+C1(1,1,K,JC) 123 CONTINUE 124 CONTINUE IF (IDO .EQ. 1) GO TO 132 IF (NBD .LT. L1) GO TO 128 DO 127 J=2,IPPH JC = IPP2-J DO 126 K=1,L1 DO 125 I=3,IDO,2 CH(1,I-1,K,J) = C1(1,I-1,K,J)-C1(1,I,K,JC) CH(1,I-1,K,JC) = C1(1,I-1,K,J)+C1(1,I,K,JC) CH(1,I,K,J) = C1(1,I,K,J)+C1(1,I-1,K,JC) CH(1,I,K,JC) = C1(1,I,K,J)-C1(1,I-1,K,JC) 125 CONTINUE 126 CONTINUE 127 CONTINUE GO TO 132 128 DO 131 J=2,IPPH JC = IPP2-J DO 130 I=3,IDO,2 DO 129 K=1,L1 CH(1,I-1,K,J) = C1(1,I-1,K,J)-C1(1,I,K,JC) CH(1,I-1,K,JC) = C1(1,I-1,K,J)+C1(1,I,K,JC) CH(1,I,K,J) = C1(1,I,K,J)+C1(1,I-1,K,JC) CH(1,I,K,JC) = C1(1,I,K,J)-C1(1,I-1,K,JC) 129 CONTINUE 130 CONTINUE 131 CONTINUE 132 CONTINUE IF (IDO .EQ. 1) RETURN DO 133 IK=1,IDL1 C2(1,IK,1) = CH2(1,IK,1) 133 CONTINUE DO 135 J=2,IP DO 134 K=1,L1 C1(1,1,K,J) = CH(1,1,K,J) 134 CONTINUE 135 CONTINUE IF (NBD .GT. L1) GO TO 139 IS = -IDO DO 138 J=2,IP IS = IS+IDO IDIJ = IS DO 137 I=3,IDO,2 IDIJ = IDIJ+2 DO 136 K=1,L1 C1(1,I-1,K,J) = WA(IDIJ-1)*CH(1,I-1,K,J)-WA(IDIJ)* 1 CH(1,I,K,J) C1(1,I,K,J) = WA(IDIJ-1)*CH(1,I,K,J)+WA(IDIJ)* 1 CH(1,I-1,K,J) 136 CONTINUE 137 CONTINUE 138 CONTINUE GO TO 143 139 IS = -IDO DO 142 J=2,IP IS = IS+IDO DO 141 K=1,L1 IDIJ = IS DO 140 I=3,IDO,2 IDIJ = IDIJ+2 C1(1,I-1,K,J) = WA(IDIJ-1)*CH(1,I-1,K,J)-WA(IDIJ)* 1 CH(1,I,K,J) C1(1,I,K,J) = WA(IDIJ-1)*CH(1,I,K,J)+WA(IDIJ)* 1 CH(1,I-1,K,J) 140 CONTINUE 141 CONTINUE 142 CONTINUE 143 RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C SUBROUTINE R1F2KF (IDO,L1,CC,IN1,CH,IN2,WA1) implicit real*8 (a-h,o-z) dimension CH(IN2,IDO,2,L1) ,CC(IN1,IDO,L1,2) , WA1(IDO) C DO 101 K=1,L1 CH(1,1,1,K) = CC(1,1,K,1)+CC(1,1,K,2) CH(1,IDO,2,K) = CC(1,1,K,1)-CC(1,1,K,2) 101 CONTINUE IF (IDO-2) 107,105,102 102 IDP2 = IDO+2 DO 104 K=1,L1 DO 103 I=3,IDO,2 IC = IDP2-I CH(1,I,1,K) = CC(1,I,K,1)+(WA1(I-2)*CC(1,I,K,2)- 1 WA1(I-1)*CC(1,I-1,K,2)) CH(1,IC,2,K) = (WA1(I-2)*CC(1,I,K,2)-WA1(I-1)* 1 CC(1,I-1,K,2))-CC(1,I,K,1) CH(1,I-1,1,K) = CC(1,I-1,K,1)+(WA1(I-2)*CC(1,I-1,K,2)+ 1 WA1(I-1)*CC(1,I,K,2)) CH(1,IC-1,2,K) = CC(1,I-1,K,1)-(WA1(I-2)*CC(1,I-1,K,2)+ 1 WA1(I-1)*CC(1,I,K,2)) 103 CONTINUE 104 CONTINUE IF (MOD(IDO,2) .EQ. 1) RETURN 105 DO 106 K=1,L1 CH(1,1,2,K) = -CC(1,IDO,K,2) CH(1,IDO,1,K) = CC(1,IDO,K,1) 106 CONTINUE 107 RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C SUBROUTINE R1F3KF (IDO,L1,CC,IN1,CH,IN2,WA1,WA2) implicit real*8 (a-h,o-z) dimension CH(IN2,IDO,3,L1) ,CC(IN1,IDO,L1,3) , 1 WA1(IDO) ,WA2(IDO) C ARG=8.d0*datan(1.d0)/3.d0 TAUR=dcos(ARG) TAUI=dsin(ARG) DO 101 K=1,L1 CH(1,1,1,K) = CC(1,1,K,1)+(CC(1,1,K,2)+CC(1,1,K,3)) CH(1,1,3,K) = TAUI*(CC(1,1,K,3)-CC(1,1,K,2)) CH(1,IDO,2,K) = CC(1,1,K,1)+TAUR* 1 (CC(1,1,K,2)+CC(1,1,K,3)) 101 CONTINUE IF (IDO .EQ. 1) RETURN IDP2 = IDO+2 DO 103 K=1,L1 DO 102 I=3,IDO,2 IC = IDP2-I CH(1,I-1,1,K) = CC(1,I-1,K,1)+((WA1(I-2)*CC(1,I-1,K,2)+ 1 WA1(I-1)*CC(1,I,K,2))+(WA2(I-2)*CC(1,I-1,K,3)+WA2(I-1)* 1 CC(1,I,K,3))) CH(1,I,1,K) = CC(1,I,K,1)+((WA1(I-2)*CC(1,I,K,2)- 1 WA1(I-1)*CC(1,I-1,K,2))+(WA2(I-2)*CC(1,I,K,3)-WA2(I-1)* 1 CC(1,I-1,K,3))) CH(1,I-1,3,K) = (CC(1,I-1,K,1)+TAUR*((WA1(I-2)* 1 CC(1,I-1,K,2)+WA1(I-1)*CC(1,I,K,2))+(WA2(I-2)* 1 CC(1,I-1,K,3)+WA2(I-1)*CC(1,I,K,3))))+(TAUI*((WA1(I-2)* 1 CC(1,I,K,2)-WA1(I-1)*CC(1,I-1,K,2))-(WA2(I-2)* 1 CC(1,I,K,3)-WA2(I-1)*CC(1,I-1,K,3)))) CH(1,IC-1,2,K) = (CC(1,I-1,K,1)+TAUR*((WA1(I-2)* 1 CC(1,I-1,K,2)+WA1(I-1)*CC(1,I,K,2))+(WA2(I-2)* 1 CC(1,I-1,K,3)+WA2(I-1)*CC(1,I,K,3))))-(TAUI*((WA1(I-2)* 1 CC(1,I,K,2)-WA1(I-1)*CC(1,I-1,K,2))-(WA2(I-2)* 1 CC(1,I,K,3)-WA2(I-1)*CC(1,I-1,K,3)))) CH(1,I,3,K) = (CC(1,I,K,1)+TAUR*((WA1(I-2)*CC(1,I,K,2)- 1 WA1(I-1)*CC(1,I-1,K,2))+(WA2(I-2)*CC(1,I,K,3)-WA2(I-1)* 1 CC(1,I-1,K,3))))+(TAUI*((WA2(I-2)*CC(1,I-1,K,3)+WA2(I-1)* 1 CC(1,I,K,3))-(WA1(I-2)*CC(1,I-1,K,2)+WA1(I-1)* 1 CC(1,I,K,2)))) CH(1,IC,2,K) = (TAUI*((WA2(I-2)*CC(1,I-1,K,3)+WA2(I-1)* 1 CC(1,I,K,3))-(WA1(I-2)*CC(1,I-1,K,2)+WA1(I-1)* 1 CC(1,I,K,2))))-(CC(1,I,K,1)+TAUR*((WA1(I-2)*CC(1,I,K,2)- 1 WA1(I-1)*CC(1,I-1,K,2))+(WA2(I-2)*CC(1,I,K,3)-WA2(I-1)* 1 CC(1,I-1,K,3)))) 102 CONTINUE 103 CONTINUE RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C SUBROUTINE R1F4KF (IDO,L1,CC,IN1,CH,IN2,WA1,WA2,WA3) implicit real*8 (a-h,o-z) dimension CC(IN1,IDO,L1,4) ,CH(IN2,IDO,4,L1), 1 WA1(IDO) ,WA2(IDO) ,WA3(IDO) C HSQT2=dsqrt(2.d0)/2.d0 DO 101 K=1,L1 CH(1,1,1,K) = (CC(1,1,K,2)+CC(1,1,K,4)) 1 +(CC(1,1,K,1)+CC(1,1,K,3)) CH(1,IDO,4,K) = (CC(1,1,K,1)+CC(1,1,K,3)) 1 -(CC(1,1,K,2)+CC(1,1,K,4)) CH(1,IDO,2,K) = CC(1,1,K,1)-CC(1,1,K,3) CH(1,1,3,K) = CC(1,1,K,4)-CC(1,1,K,2) 101 CONTINUE IF (IDO-2) 107,105,102 102 IDP2 = IDO+2 DO 104 K=1,L1 DO 103 I=3,IDO,2 IC = IDP2-I CH(1,I-1,1,K) = ((WA1(I-2)*CC(1,I-1,K,2)+WA1(I-1)* 1 CC(1,I,K,2))+(WA3(I-2)*CC(1,I-1,K,4)+WA3(I-1)* 1 CC(1,I,K,4)))+(CC(1,I-1,K,1)+(WA2(I-2)*CC(1,I-1,K,3)+ 1 WA2(I-1)*CC(1,I,K,3))) CH(1,IC-1,4,K) = (CC(1,I-1,K,1)+(WA2(I-2)*CC(1,I-1,K,3)+ 1 WA2(I-1)*CC(1,I,K,3)))-((WA1(I-2)*CC(1,I-1,K,2)+ 1 WA1(I-1)*CC(1,I,K,2))+(WA3(I-2)*CC(1,I-1,K,4)+ 1 WA3(I-1)*CC(1,I,K,4))) CH(1,I,1,K) = ((WA1(I-2)*CC(1,I,K,2)-WA1(I-1)* 1 CC(1,I-1,K,2))+(WA3(I-2)*CC(1,I,K,4)-WA3(I-1)* 1 CC(1,I-1,K,4)))+(CC(1,I,K,1)+(WA2(I-2)*CC(1,I,K,3)- 1 WA2(I-1)*CC(1,I-1,K,3))) CH(1,IC,4,K) = ((WA1(I-2)*CC(1,I,K,2)-WA1(I-1)* 1 CC(1,I-1,K,2))+(WA3(I-2)*CC(1,I,K,4)-WA3(I-1)* 1 CC(1,I-1,K,4)))-(CC(1,I,K,1)+(WA2(I-2)*CC(1,I,K,3)- 1 WA2(I-1)*CC(1,I-1,K,3))) CH(1,I-1,3,K) = ((WA1(I-2)*CC(1,I,K,2)-WA1(I-1)* 1 CC(1,I-1,K,2))-(WA3(I-2)*CC(1,I,K,4)-WA3(I-1)* 1 CC(1,I-1,K,4)))+(CC(1,I-1,K,1)-(WA2(I-2)*CC(1,I-1,K,3)+ 1 WA2(I-1)*CC(1,I,K,3))) CH(1,IC-1,2,K) = (CC(1,I-1,K,1)-(WA2(I-2)*CC(1,I-1,K,3)+ 1 WA2(I-1)*CC(1,I,K,3)))-((WA1(I-2)*CC(1,I,K,2)-WA1(I-1)* 1 CC(1,I-1,K,2))-(WA3(I-2)*CC(1,I,K,4)-WA3(I-1)* 1 CC(1,I-1,K,4))) CH(1,I,3,K) = ((WA3(I-2)*CC(1,I-1,K,4)+WA3(I-1)* 1 CC(1,I,K,4))-(WA1(I-2)*CC(1,I-1,K,2)+WA1(I-1)* 1 CC(1,I,K,2)))+(CC(1,I,K,1)-(WA2(I-2)*CC(1,I,K,3)- 1 WA2(I-1)*CC(1,I-1,K,3))) CH(1,IC,2,K) = ((WA3(I-2)*CC(1,I-1,K,4)+WA3(I-1)* 1 CC(1,I,K,4))-(WA1(I-2)*CC(1,I-1,K,2)+WA1(I-1)* 1 CC(1,I,K,2)))-(CC(1,I,K,1)-(WA2(I-2)*CC(1,I,K,3)- 1 WA2(I-1)*CC(1,I-1,K,3))) 103 CONTINUE 104 CONTINUE IF (MOD(IDO,2) .EQ. 1) RETURN 105 CONTINUE DO 106 K=1,L1 CH(1,IDO,1,K) = (HSQT2*(CC(1,IDO,K,2)-CC(1,IDO,K,4)))+ 1 CC(1,IDO,K,1) CH(1,IDO,3,K) = CC(1,IDO,K,1)-(HSQT2*(CC(1,IDO,K,2)- 1 CC(1,IDO,K,4))) CH(1,1,2,K) = (-HSQT2*(CC(1,IDO,K,2)+CC(1,IDO,K,4)))- 1 CC(1,IDO,K,3) CH(1,1,4,K) = (-HSQT2*(CC(1,IDO,K,2)+CC(1,IDO,K,4)))+ 1 CC(1,IDO,K,3) 106 CONTINUE 107 RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C SUBROUTINE R1F5KF (IDO,L1,CC,IN1,CH,IN2,WA1,WA2,WA3,WA4) implicit real*8 (a-h,o-z) dimension CC(IN1,IDO,L1,5) ,CH(IN2,IDO,5,L1) , 1 WA1(IDO) ,WA2(IDO) ,WA3(IDO) ,WA4(IDO) C ARG=8.d0*datan(1.d0)/5.d0 TR11=dcos(ARG) TI11=dsin(ARG) TR12=dcos(2.d0*ARG) TI12=dsin(2.d0*ARG) DO 101 K=1,L1 CH(1,1,1,K) = CC(1,1,K,1)+(CC(1,1,K,5)+CC(1,1,K,2))+ 1 (CC(1,1,K,4)+CC(1,1,K,3)) CH(1,IDO,2,K) = CC(1,1,K,1)+TR11*(CC(1,1,K,5)+CC(1,1,K,2))+ 1 TR12*(CC(1,1,K,4)+CC(1,1,K,3)) CH(1,1,3,K) = TI11*(CC(1,1,K,5)-CC(1,1,K,2))+TI12* 1 (CC(1,1,K,4)-CC(1,1,K,3)) CH(1,IDO,4,K) = CC(1,1,K,1)+TR12*(CC(1,1,K,5)+CC(1,1,K,2))+ 1 TR11*(CC(1,1,K,4)+CC(1,1,K,3)) CH(1,1,5,K) = TI12*(CC(1,1,K,5)-CC(1,1,K,2))-TI11* 1 (CC(1,1,K,4)-CC(1,1,K,3)) 101 CONTINUE IF (IDO .EQ. 1) RETURN IDP2 = IDO+2 DO 103 K=1,L1 DO 102 I=3,IDO,2 IC = IDP2-I CH(1,I-1,1,K) = CC(1,I-1,K,1)+((WA1(I-2)*CC(1,I-1,K,2)+ 1 WA1(I-1)*CC(1,I,K,2))+(WA4(I-2)*CC(1,I-1,K,5)+WA4(I-1)* 1 CC(1,I,K,5)))+((WA2(I-2)*CC(1,I-1,K,3)+WA2(I-1)* 1 CC(1,I,K,3))+(WA3(I-2)*CC(1,I-1,K,4)+ 1 WA3(I-1)*CC(1,I,K,4))) CH(1,I,1,K) = CC(1,I,K,1)+((WA1(I-2)*CC(1,I,K,2)- 1 WA1(I-1)*CC(1,I-1,K,2))+(WA4(I-2)*CC(1,I,K,5)-WA4(I-1)* 1 CC(1,I-1,K,5)))+((WA2(I-2)*CC(1,I,K,3)-WA2(I-1)* 1 CC(1,I-1,K,3))+(WA3(I-2)*CC(1,I,K,4)-WA3(I-1)* 1 CC(1,I-1,K,4))) CH(1,I-1,3,K) = CC(1,I-1,K,1)+TR11* 1 ( WA1(I-2)*CC(1,I-1,K,2)+WA1(I-1)*CC(1,I,K,2) 1 +WA4(I-2)*CC(1,I-1,K,5)+WA4(I-1)*CC(1,I,K,5))+TR12* 1 ( WA2(I-2)*CC(1,I-1,K,3)+WA2(I-1)*CC(1,I,K,3) 1 +WA3(I-2)*CC(1,I-1,K,4)+WA3(I-1)*CC(1,I,K,4))+TI11* 1 ( WA1(I-2)*CC(1,I,K,2)-WA1(I-1)*CC(1,I-1,K,2) 1 -(WA4(I-2)*CC(1,I,K,5)-WA4(I-1)*CC(1,I-1,K,5)))+TI12* 1 ( WA2(I-2)*CC(1,I,K,3)-WA2(I-1)*CC(1,I-1,K,3) 1 -(WA3(I-2)*CC(1,I,K,4)-WA3(I-1)*CC(1,I-1,K,4))) CH(1,IC-1,2,K) = CC(1,I-1,K,1)+TR11* 1 ( WA1(I-2)*CC(1,I-1,K,2)+WA1(I-1)*CC(1,I,K,2) 1 +WA4(I-2)*CC(1,I-1,K,5)+WA4(I-1)*CC(1,I,K,5))+TR12* 1 ( WA2(I-2)*CC(1,I-1,K,3)+WA2(I-1)*CC(1,I,K,3) 1 +WA3(I-2)*CC(1,I-1,K,4)+WA3(I-1)*CC(1,I,K,4))-(TI11* 1 ( WA1(I-2)*CC(1,I,K,2)-WA1(I-1)*CC(1,I-1,K,2) 1 -(WA4(I-2)*CC(1,I,K,5)-WA4(I-1)*CC(1,I-1,K,5)))+TI12* 1 ( WA2(I-2)*CC(1,I,K,3)-WA2(I-1)*CC(1,I-1,K,3) 1 -(WA3(I-2)*CC(1,I,K,4)-WA3(I-1)*CC(1,I-1,K,4)))) CH(1,I,3,K) = (CC(1,I,K,1)+TR11*((WA1(I-2)*CC(1,I,K,2)- 1 WA1(I-1)*CC(1,I-1,K,2))+(WA4(I-2)*CC(1,I,K,5)-WA4(I-1)* 1 CC(1,I-1,K,5)))+TR12*((WA2(I-2)*CC(1,I,K,3)-WA2(I-1)* 1 CC(1,I-1,K,3))+(WA3(I-2)*CC(1,I,K,4)-WA3(I-1)* 1 CC(1,I-1,K,4))))+(TI11*((WA4(I-2)*CC(1,I-1,K,5)+ 1 WA4(I-1)*CC(1,I,K,5))-(WA1(I-2)*CC(1,I-1,K,2)+WA1(I-1)* 1 CC(1,I,K,2)))+TI12*((WA3(I-2)*CC(1,I-1,K,4)+WA3(I-1)* 1 CC(1,I,K,4))-(WA2(I-2)*CC(1,I-1,K,3)+WA2(I-1)* 1 CC(1,I,K,3)))) CH(1,IC,2,K) = (TI11*((WA4(I-2)*CC(1,I-1,K,5)+WA4(I-1)* 1 CC(1,I,K,5))-(WA1(I-2)*CC(1,I-1,K,2)+WA1(I-1)* 1 CC(1,I,K,2)))+TI12*((WA3(I-2)*CC(1,I-1,K,4)+WA3(I-1)* 1 CC(1,I,K,4))-(WA2(I-2)*CC(1,I-1,K,3)+WA2(I-1)* 1 CC(1,I,K,3))))-(CC(1,I,K,1)+TR11*((WA1(I-2)*CC(1,I,K,2)- 1 WA1(I-1)*CC(1,I-1,K,2))+(WA4(I-2)*CC(1,I,K,5)-WA4(I-1)* 1 CC(1,I-1,K,5)))+TR12*((WA2(I-2)*CC(1,I,K,3)-WA2(I-1)* 1 CC(1,I-1,K,3))+(WA3(I-2)*CC(1,I,K,4)-WA3(I-1)* 1 CC(1,I-1,K,4)))) CH(1,I-1,5,K) = (CC(1,I-1,K,1)+TR12*((WA1(I-2)* 1 CC(1,I-1,K,2)+WA1(I-1)*CC(1,I,K,2))+(WA4(I-2)* 1 CC(1,I-1,K,5)+WA4(I-1)*CC(1,I,K,5)))+TR11*((WA2(I-2)* 1 CC(1,I-1,K,3)+WA2(I-1)*CC(1,I,K,3))+(WA3(I-2)* 1 CC(1,I-1,K,4)+WA3(I-1)*CC(1,I,K,4))))+(TI12*((WA1(I-2)* 1 CC(1,I,K,2)-WA1(I-1)*CC(1,I-1,K,2))-(WA4(I-2)* 1 CC(1,I,K,5)-WA4(I-1)*CC(1,I-1,K,5)))-TI11*((WA2(I-2)* 1 CC(1,I,K,3)-WA2(I-1)*CC(1,I-1,K,3))-(WA3(I-2)* 1 CC(1,I,K,4)-WA3(I-1)*CC(1,I-1,K,4)))) CH(1,IC-1,4,K) = (CC(1,I-1,K,1)+TR12*((WA1(I-2)* 1 CC(1,I-1,K,2)+WA1(I-1)*CC(1,I,K,2))+(WA4(I-2)* 1 CC(1,I-1,K,5)+WA4(I-1)*CC(1,I,K,5)))+TR11*((WA2(I-2)* 1 CC(1,I-1,K,3)+WA2(I-1)*CC(1,I,K,3))+(WA3(I-2)* 1 CC(1,I-1,K,4)+WA3(I-1)*CC(1,I,K,4))))-(TI12*((WA1(I-2)* 1 CC(1,I,K,2)-WA1(I-1)*CC(1,I-1,K,2))-(WA4(I-2)* 1 CC(1,I,K,5)-WA4(I-1)*CC(1,I-1,K,5)))-TI11*((WA2(I-2)* 1 CC(1,I,K,3)-WA2(I-1)*CC(1,I-1,K,3))-(WA3(I-2)* 1 CC(1,I,K,4)-WA3(I-1)*CC(1,I-1,K,4)))) CH(1,I,5,K) = (CC(1,I,K,1)+TR12*((WA1(I-2)*CC(1,I,K,2)- 1 WA1(I-1)*CC(1,I-1,K,2))+(WA4(I-2)*CC(1,I,K,5)-WA4(I-1)* 1 CC(1,I-1,K,5)))+TR11*((WA2(I-2)*CC(1,I,K,3)-WA2(I-1)* 1 CC(1,I-1,K,3))+(WA3(I-2)*CC(1,I,K,4)-WA3(I-1)* 1 CC(1,I-1,K,4))))+(TI12*((WA4(I-2)*CC(1,I-1,K,5)+ 1 WA4(I-1)*CC(1,I,K,5))-(WA1(I-2)*CC(1,I-1,K,2)+WA1(I-1)* 1 CC(1,I,K,2)))-TI11*((WA3(I-2)*CC(1,I-1,K,4)+WA3(I-1)* 1 CC(1,I,K,4))-(WA2(I-2)*CC(1,I-1,K,3)+WA2(I-1)* 1 CC(1,I,K,3)))) CH(1,IC,4,K) = (TI12*((WA4(I-2)*CC(1,I-1,K,5)+WA4(I-1)* 1 CC(1,I,K,5))-(WA1(I-2)*CC(1,I-1,K,2)+WA1(I-1)* 1 CC(1,I,K,2)))-TI11*((WA3(I-2)*CC(1,I-1,K,4)+WA3(I-1)* 1 CC(1,I,K,4))-(WA2(I-2)*CC(1,I-1,K,3)+WA2(I-1)* 1 CC(1,I,K,3))))-(CC(1,I,K,1)+TR12*((WA1(I-2)*CC(1,I,K,2)- 1 WA1(I-1)*CC(1,I-1,K,2))+(WA4(I-2)*CC(1,I,K,5)-WA4(I-1)* 1 CC(1,I-1,K,5)))+TR11*((WA2(I-2)*CC(1,I,K,3)-WA2(I-1)* 1 CC(1,I-1,K,3))+(WA3(I-2)*CC(1,I,K,4)-WA3(I-1)* 1 CC(1,I-1,K,4)))) 102 CONTINUE 103 CONTINUE RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C SUBROUTINE R1FGKF (IDO,IP,L1,IDL1,CC,C1,C2,IN1,CH,CH2,IN2,WA) implicit real*8 (a-h,o-z) dimension CH(IN2,IDO,L1,IP) ,CC(IN1,IDO,IP,L1), 1 C1(IN1,IDO,L1,IP) ,C2(IN1,IDL1,IP), 2 CH2(IN2,IDL1,IP) ,WA(IDO) C TPI=8.d0*datan(1.d0) ARG = TPI/dble(IP) DCP = dcos(ARG) DSP = dsin(ARG) IPPH = (IP+1)/2 IPP2 = IP+2 IDP2 = IDO+2 NBD = (IDO-1)/2 IF (IDO .EQ. 1) GO TO 119 DO 101 IK=1,IDL1 CH2(1,IK,1) = C2(1,IK,1) 101 CONTINUE DO 103 J=2,IP DO 102 K=1,L1 CH(1,1,K,J) = C1(1,1,K,J) 102 CONTINUE 103 CONTINUE IF (NBD .GT. L1) GO TO 107 IS = -IDO DO 106 J=2,IP IS = IS+IDO IDIJ = IS DO 105 I=3,IDO,2 IDIJ = IDIJ+2 DO 104 K=1,L1 CH(1,I-1,K,J) = WA(IDIJ-1)*C1(1,I-1,K,J)+WA(IDIJ) 1 *C1(1,I,K,J) CH(1,I,K,J) = WA(IDIJ-1)*C1(1,I,K,J)-WA(IDIJ) 1 *C1(1,I-1,K,J) 104 CONTINUE 105 CONTINUE 106 CONTINUE GO TO 111 107 IS = -IDO DO 110 J=2,IP IS = IS+IDO DO 109 K=1,L1 IDIJ = IS DO 108 I=3,IDO,2 IDIJ = IDIJ+2 CH(1,I-1,K,J) = WA(IDIJ-1)*C1(1,I-1,K,J)+WA(IDIJ) 1 *C1(1,I,K,J) CH(1,I,K,J) = WA(IDIJ-1)*C1(1,I,K,J)-WA(IDIJ) 1 *C1(1,I-1,K,J) 108 CONTINUE 109 CONTINUE 110 CONTINUE 111 IF (NBD .LT. L1) GO TO 115 DO 114 J=2,IPPH JC = IPP2-J DO 113 K=1,L1 DO 112 I=3,IDO,2 C1(1,I-1,K,J) = CH(1,I-1,K,J)+CH(1,I-1,K,JC) C1(1,I-1,K,JC) = CH(1,I,K,J)-CH(1,I,K,JC) C1(1,I,K,J) = CH(1,I,K,J)+CH(1,I,K,JC) C1(1,I,K,JC) = CH(1,I-1,K,JC)-CH(1,I-1,K,J) 112 CONTINUE 113 CONTINUE 114 CONTINUE GO TO 121 115 DO 118 J=2,IPPH JC = IPP2-J DO 117 I=3,IDO,2 DO 116 K=1,L1 C1(1,I-1,K,J) = CH(1,I-1,K,J)+CH(1,I-1,K,JC) C1(1,I-1,K,JC) = CH(1,I,K,J)-CH(1,I,K,JC) C1(1,I,K,J) = CH(1,I,K,J)+CH(1,I,K,JC) C1(1,I,K,JC) = CH(1,I-1,K,JC)-CH(1,I-1,K,J) 116 CONTINUE 117 CONTINUE 118 CONTINUE GO TO 121 119 DO 120 IK=1,IDL1 C2(1,IK,1) = CH2(1,IK,1) 120 CONTINUE 121 DO 123 J=2,IPPH JC = IPP2-J DO 122 K=1,L1 C1(1,1,K,J) = CH(1,1,K,J)+CH(1,1,K,JC) C1(1,1,K,JC) = CH(1,1,K,JC)-CH(1,1,K,J) 122 CONTINUE 123 CONTINUE C AR1 = 1.d0 AI1 = 0.d0 DO 127 L=2,IPPH LC = IPP2-L AR1H = DCP*AR1-DSP*AI1 AI1 = DCP*AI1+DSP*AR1 AR1 = AR1H DO 124 IK=1,IDL1 CH2(1,IK,L) = C2(1,IK,1)+AR1*C2(1,IK,2) CH2(1,IK,LC) = AI1*C2(1,IK,IP) 124 CONTINUE DC2 = AR1 DS2 = AI1 AR2 = AR1 AI2 = AI1 DO 126 J=3,IPPH JC = IPP2-J AR2H = DC2*AR2-DS2*AI2 AI2 = DC2*AI2+DS2*AR2 AR2 = AR2H DO 125 IK=1,IDL1 CH2(1,IK,L) = CH2(1,IK,L)+AR2*C2(1,IK,J) CH2(1,IK,LC) = CH2(1,IK,LC)+AI2*C2(1,IK,JC) 125 CONTINUE 126 CONTINUE 127 CONTINUE DO 129 J=2,IPPH DO 128 IK=1,IDL1 CH2(1,IK,1) = CH2(1,IK,1)+C2(1,IK,J) 128 CONTINUE 129 CONTINUE C IF (IDO .LT. L1) GO TO 132 DO 131 K=1,L1 DO 130 I=1,IDO CC(1,I,1,K) = CH(1,I,K,1) 130 CONTINUE 131 CONTINUE GO TO 135 132 DO 134 I=1,IDO DO 133 K=1,L1 CC(1,I,1,K) = CH(1,I,K,1) 133 CONTINUE 134 CONTINUE 135 DO 137 J=2,IPPH JC = IPP2-J J2 = J+J DO 136 K=1,L1 CC(1,IDO,J2-2,K) = CH(1,1,K,J) CC(1,1,J2-1,K) = CH(1,1,K,JC) 136 CONTINUE 137 CONTINUE IF (IDO .EQ. 1) RETURN IF (NBD .LT. L1) GO TO 141 DO 140 J=2,IPPH JC = IPP2-J J2 = J+J DO 139 K=1,L1 DO 138 I=3,IDO,2 IC = IDP2-I CC(1,I-1,J2-1,K) = CH(1,I-1,K,J)+CH(1,I-1,K,JC) CC(1,IC-1,J2-2,K) = CH(1,I-1,K,J)-CH(1,I-1,K,JC) CC(1,I,J2-1,K) = CH(1,I,K,J)+CH(1,I,K,JC) CC(1,IC,J2-2,K) = CH(1,I,K,JC)-CH(1,I,K,J) 138 CONTINUE 139 CONTINUE 140 CONTINUE RETURN 141 DO 144 J=2,IPPH JC = IPP2-J J2 = J+J DO 143 I=3,IDO,2 IC = IDP2-I DO 142 K=1,L1 CC(1,I-1,J2-1,K) = CH(1,I-1,K,J)+CH(1,I-1,K,JC) CC(1,IC-1,J2-2,K) = CH(1,I-1,K,J)-CH(1,I-1,K,JC) CC(1,I,J2-1,K) = CH(1,I,K,J)+CH(1,I,K,JC) CC(1,IC,J2-2,K) = CH(1,I,K,JC)-CH(1,I,K,J) 142 CONTINUE 143 CONTINUE 144 CONTINUE RETURN END