求大神帮忙看看怎么没法输出文件

pandadod · 发表于 2019-6-5 13:28:14

F 77 运行时没有出现debug

[Fortran] 纯文本查看 复制代码

            parameter (nn=500000)
      dimension x1(nn),x2(nn),x3(nn)
      dimension xns(nn),xew(nn),xud(nn)
      dimension f(nn),f1(nn),f2(nn),f3(nn),period(nn)
      character title*32,dim*3,idsn*40,odsn*40,idsn_ctl*40
c
      write(*,*) 'input control file name (*.ctl)'
      read(*,*) idsn_ctl
c
      open(10,file=idsn_ctl)
      read(10,*)
      read(10,*) idsn
      read(10,*)
      read(10,*) sf1,nd
      read(10,*)
      read(10,*) band
      read(10,*)
      read(10,*) odsn
c
c

      call wvread_csv(x1,x2,x3,n,nd,title,hz,dim,idsn)
c
c
      nd2=nd/2+1
c     
      dt=1.0/hz
      t=float(nd)*dt
c
      ls=int(sf1*hz)

      do 110 l=1,nd
      ll=l+ls
      x1(l)=x1(ll)
      x2(l)=x2(ll)
      x3(l)=x3(ll)
  110 continue
c
      call fpac(nd,x1,nd,dt,100,f,g,g,nd2,nfold,df)
      do l=1,nfold
      f1(l)=f(l)
      enddo
      call swin(nfold,f1,g,nd2,100,df,band)
      do m=1,nd2
      xns(m)=f1(m)
      enddo

c
      call fpac(nd,x2,nd,dt,100,f,g,g,nd2,nfold,df)
      do l=1,nfold
      f2(l)=f(l)
      enddo
      call swin(nfold,f2,g,nd2,100,df,band)
      do m=1,nd2
      xew(m)=f2(m)
      enddo

c
      call fpac(nd,x3,nd,dt,100,f,g,g,nd2,nfold,df)
      do l=1,nfold
      f3(l)=f(l)
      enddo
      call swin(nfold,f3,g,nd2,100,df,band)
      do m=1,nd2
      xud(m)=f3(m)
      enddo

c   
c

      do k=1,nd2
      period(k)=t/k
      enddo
c
      open(70,file=odsn)
      write(70,200) 
  200 format(2x,'Period(s)',3x,'Amp_NS',6x,'Amp_EW',6x,'Amp_UD')
      do j=1,nd2
      write(70,220) period(j),xns(j),xew(j),xud(j)
  220 format(7e12.4)
      enddo
c
      close(70)
      stop
      end

c---------------------------------------------------------------------
c     subroutine for fourier,power spectra and autocorrelation
c---------------------------------------------------------------------
c
c                                           coded by y.ohasaki
c
      subroutine fpac(n,x,nd1,dt,ind,f,g,r,nd2,nfold,df)
c
      complex c(8192)
      dimension x(nd1),f(nd2),g(nd2),r(nd2)
c
c     initialization

      do 110 m=1,n
      c(m)=cmplx(x(m),0.)
  110 continue
      nt=2
  120 if(nt.ge.n) go to 130
      nt=nt*2
      go to 120
  130 if(nt.eq.n) go to 150
      do 140 m=n+1,nt
      c(m)=(0.,0.)
  140 continue
  150 nfold=nt/2+1
      t=real(nt)*dt
      df=1./t
c
c
c     fourier transform
c
      call fast(nt,c,8192,-1)
c
c     fourier spectrum
c
      if(ind.eq.001) goto 180
      do 151 k=1,nt
      c(k)=c(k)/nt
  151 continue
      do 160 k=1,nfold
      f(k)=cabs(c(k))*dt
  160 continue
      if(ind.eq.100) return
c
c     power spectrum
c
      if(ind.eq.101) go to 180
      g(1)=f(1)**2/t
      do 170 k=2,nfold-1
      g(k)=2.*f(k)**2/t
  170 continue
      g(nfold)=f(nfold)**2/t
      if(mod(ind,10).eq.0) return
c
c     autocorrelation
c
  180 do 190 k=1,nt
      c(k)=c(k)*conjg(c(k))
  190 continue
      call fast(nt,c,8192,+1)
      r0=real(c(1))
      do 200 j=1,nfold
      r(j)=real(c(j))/r0
  200 continue
c
      return
      end
c

c---------------------------------------------------------------------
c     subroutine for fast fourier transform
c---------------------------------------------------------------------
      subroutine fast(n,x,nd1,ind)
      complex x(nd1),temp,theta
      j=1
      do 140 i=1,n
      if(i.ge.j) go to 110
      temp=x(j)
      x(j)=x(i)
      x(i)=temp
  110 m=n/2
  120 if(j.le.m) go to 130
      j=j-m
      m=m/2
      if(m.ge.2) go to 120
  130 j=j+m
  140 continue
      kmax=1
  150 if(kmax.ge.n) return
      istep=kmax*2
      do 170 k=1,kmax
      theta=cmplx(0.0,3.141593*float(ind*(k-1))/float(kmax))
      do 160 i=k,n,istep
      j=i+kmax
      temp=x(j)*cexp(theta)
      x(j)=x(i)-temp
      x(i)=x(i)+temp
  160 continue
  170 continue
      kmax=istep
      go to 150
      end
c

c
c---------------------------------------------------------------------
c     subroutine for smoothed spectra by parzen's spectral window
c---------------------------------------------------------------------
c                                               coded by y.ohsaki
c
      subroutine swin(nfold,f,g,nd,ind,df,band)
c
      dimension f(nd),g(nd),w(101),g1(4497),g2(4497)
c
c     initialization
c
      t=1./df
      if(band.eq.0.) go to 120
      udf=1.854305/band*df
      if(udf.gt.0.5) go to 230
      lmax=ifix(2./udf)+1
      if(lmax.gt.101) go to 240
c
c     spectral window
c
      w(1)=0.75*udf
      do 110 l=2,lmax
      dif=1.570796*real(l-1)*udf
      w(l)=w(1)*(sin(dif)/dif)**4
  110 continue
c
c     conversion from fourier to power spectrum
c
  120 if(ind.ne.100) go to 140
      g(1)=f(1)**2/t
      do 130 k=2,nfold-1
      g(k)=2.*f(k)**2/t
  130 continue
      g(nfold)=f(nfold)**2/t
c
c     smoothing of power spectrum
c
  140 if(band.eq.0.) go to 210
      ll=lmax*2-1
      ln=ll-1+nfold
      lt=(ll-1)*2+nfold
      le=lt-lmax+1
      do 150 k=1,lt
      g1(k)=0.
  150 continue
      do 160 k=1,nfold
      g1(ll-1+k)=g(k)
  160 continue
      do 180 k=lmax,le
      s=w(1)*g1(k)
      do 170 l=2,lmax
      s=s+w(l)*(g1(k-l+1)+g1(k+l-1))
  170 continue
      g2(k)=s
  180 continue
      do 190 l=2,lmax
      g2(ll+l-1)=g2(ll+l-1)+g2(ll-l+1)
      g2(ln-l+1)=g2(ln-l+1)+g2(ln+l-1)
  190 continue
      do 200 k=1,nfold
      g(k)=g2(ll-1+k)
  200 continue
c
c     smoothed fourier spectrum
c
  210 f(1)=sqrt(g(1)*t)
      do 220 k=2,nfold-1
      f(k)=sqrt(g(k)*t/2.)
  220 continue
      f(nfold)=sqrt(g(nfold)*t)
      return
c
  230 write(30,601)
      stop
  240 write(30,602)
      stop
c
c     format statements
c
  601 format('bandwidth is too narrow')
  602 format('bandwidth is too wide')
      end


c---- 攇宍僨乕僞傪撉傒崬傓僒僽儖乕僠儞 ----------------------------
      subroutine wvread_csv(x1,x2,x3,n,nd1,title,hz,dim,idsn)
      dimension x1(nd1),x2(nd1),x3(nd1)
      character title*32,dim*3,idsn*40,stcode*6,date*22,dim0*19
c
      open(60,file=idsn)
      read(60,601) stcode
  601 format(11x,a6)
      read(60,602) n
  602 format(16x,i5)
      read(60,603) hz
  603 format(15x,f5.1)
      read(60,*) dim0
       if(dim0.eq.'UNIT  = gal(cm/s/s)') dim='ACC'
       if(dim0.eq.'UNIT  = kine(cm/s) ') dim='VEL'
       if(dim0.eq.'UNIT  = cm         ') dim='DIS'
      read(60,604) date
  604 format(15x,a22)
      title=stcode(1:6)//' ,  '//date(1:22)
      read(60,*)
      do i=1,n
      read(60,*) temp,x1(i),x2(i),x3(i)
      enddo
      close(60)
c
      return
      end
c