Fortran Coder

查看: 95|回复: 8

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

[复制链接]

5

帖子

1

主题

0

精华

新人

F 币
30 元
贡献
14 点
发表于 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


回复

使用道具 举报

5

帖子

1

主题

0

精华

新人

F 币
30 元
贡献
14 点
 楼主| 发表于 2019-6-5 13:29:11 | 显示全部楼层
(*.csv)-------------
HRS019_01324wa.csv
秒,个数
0.0 8192
BAND(HZ)
0.1
输出文件名-----------
kadai4.csv

5

帖子

1

主题

0

精华

新人

F 币
30 元
贡献
14 点
 楼主| 发表于 2019-6-5 13:29:55 | 显示全部楼层
2楼是ctl文件

685

帖子

0

主题

0

精华

大师

F 币
638 元
贡献
357 点

规矩勋章元老勋章新人勋章水王勋章

QQ
发表于 2019-6-5 13:42:34 | 显示全部楼层
没有 HRS019_01324wa.csv 文件,无法测试。

5

帖子

1

主题

0

精华

新人

F 币
30 元
贡献
14 点
 楼主| 发表于 2019-6-5 14:30:47 | 显示全部楼层
已上传该文件
麻烦帮忙检查一下

HRS019_01324wa.zip

194.05 KB, 下载次数: 1

685

帖子

0

主题

0

精华

大师

F 币
638 元
贡献
357 点

规矩勋章元老勋章新人勋章水王勋章

QQ
发表于 2019-6-5 14:42:50 | 显示全部楼层
HRS019_01324wa.csv 文件里
NUMBER OF DATA= 14900
而 ctl文件里个数是
0.0 8192
也就是 8192,导致数据不吻合。
导致 wvread_csv 函数里 read(60,*) temp,x1(i),x2(i),x3(i) 数组越界。

685

帖子

0

主题

0

精华

大师

F 币
638 元
贡献
357 点

规矩勋章元老勋章新人勋章水王勋章

QQ
发表于 2019-6-5 14:48:57 | 显示全部楼层
本帖最后由 vvt 于 2019-6-5 14:52 编辑

你的这个代码应该是要求 NUMBER OF DATA= 14900  这里是 2的整幂次方的。
而且必须要求和 ctl 文件里的一致。
此外,代码里也有不少地方写死了 8192,都需要更改成一致的。

5

帖子

1

主题

0

精华

新人

F 币
30 元
贡献
14 点
 楼主| 发表于 2019-6-5 15:05:46 | 显示全部楼层
vvt 发表于 2019-6-5 14:48
你的这个代码应该是要求 NUMBER OF DATA= 14900  这里是 2的整幂次方的。
而且必须要求和 ctl 文件里的一致 ...

试着改了几处还是不行要哭了

685

帖子

0

主题

0

精华

大师

F 币
638 元
贡献
357 点

规矩勋章元老勋章新人勋章水王勋章

QQ
发表于 2019-6-5 16:51:38 | 显示全部楼层
那你就哭吧,如果哭有用的话,你就继续哭。

哭够了,觉得没有用了,就继续说:“你做了什么改动?哪几个地方?现在遇到的问题是什么?”
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

QQ|捐赠本站|Archiver|关于我们 About Us|群聊|Fcode

GMT+8, 2019-6-27 09:59

Powered by Discuz! X3.2

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表