fortran

huangyao

我用fortran读取一个grd文件，然后创建一个dat文件，提示的错误是这样的，这是什么原因呢

提示错误

风之精灵

请给出代码

勿忘草

是不是程序代码中没有给出文件类型

huangyao

风之精灵 发表于 2016-7-17 21:52
请给出代码

        这是我的代码，应该看读取文件和写入文件两句就可以了，其他没提示错误

               implicit none
                integer,parameter               :: nx=50,ny=30,nz=1,nt=20,m1=4,m2=8,l=7
        real,parameter                  :: pi=3.1415926
        real,dimension(:,:),allocatable :: yy,ui
        real,dimension(:),allocatable   :: y,x
        integer                         :: mm,it,ix,iy,iz,irec,ll,i
    real                            :: s
        mm=nx*ny*nz
        ll=2*l+1       
    allocate(yy(nt,mm))       
    open(100,file='D:\bb\500hgt\500.grd',form='unformatted',access='direct',recl=nx*ny)
        irec=1
    do it=1,nt
           do iz=1,nz
             read(100,rec=irec)((yy(it,ix+(iy-1)*nx+(iz-1)*nx*ny),ix=1,nx),iy=1,ny)
         irec=irec+1
      end do
        end do
    close(100)  
        do i=1,mm
          s=0.0
          do it=1,nt
            s=s+yy(it,i)
      end do
      s=s/real(nt)
          do it=1,nt
        yy(it,i)=yy(it,i)-s   
      end do
      s=0.0
          do it=1,nt
            s=s+yy(it,i)*yy(it,i)
      end do
      s=sqrt(s/real(nt))      
          do it=1,nt
            yy(it,i)=yy(it,i)/s   
      end do
        end do   
        do i=1,mm
          allocate(y(nt))
      do it=1,nt
        y(it)=yy(it,i)
      end do
          allocate(x(nt))
      call filter(y,nt,m1,m2,x,pi)
          deallocate(y)
      do it=1,nt
        yy(it,i)=x(it)
      end do
          deallocate(x)
    end do  
    allocate(ui(nt,mm))
    call hilbert(yy,ui,mm,nt,l,ll,pi)
     do it=1,nt
          do i=1,mm
        s=yy(it,i)**2+ui(it,i)**2
            yy(it,i)=sqrt(s)
             enddo
    enddo
        deallocate(ui)
    open(400,file='D:\bb\500hgt\bobao.dat',form='unformatted',access='direct',recl=nx*ny)     
    irec=1
        do it=1,nt
          do iz=1,nz
        write(400,rec=irec)((yy(it,ix+(iy-1)*nx+(iz-1)*nx*ny),ix=1,nx),iy=1,ny)
        irec=irec+1
          enddo
        enddo
    close(400)
        deallocate(yy)   
    end
    subroutine filter(x,n,mw1,mw2,bx,pi)
      real,dimension(n)::x,bx,cx,xy
      w1=2.*pi/real(mw1)
      w2=2.*pi/real(mw2)
      w0=sqrt(w1*w2)
      dt=1.0
      ww=2*abs(sin(w1*dt)/(1+cos(w1*dt))-sin(w2*dt)/(1+cos(w2*dt)))
      ww0=4*sin(w1*dt)*sin(w2*dt)/((1+cos(w1*dt))*(1+cos(w2*dt)))
      a0=2*ww/(4+2*ww+ww0)
      b1=2*(ww0-4)/(4+2*ww+ww0)
      b2=(4-2*ww+ww0)/(4+2*ww+ww0)
      do 89 i0=1,n-2
      i=i0+2
      bx(1)=0.0
      bx(2)=0.0
      cx(i)=a0*(x(i)-x(i-2))-b1*bx(i-1)-b2*bx(i-2)
  89  bx(i)=cx(i)
      do 91 i=1,n
  91  x(i)=bx(n+1-i)
      do 92 i0=1,n-2
      i=i0+2
      xy(1)=x(1)
      xy(2)=x(2)
      cx(i)=a0*(x(i)-x(i-2))-b1*xy(i-1)-b2*xy(i-2)
  92  xy(i)=cx(i)
      do 93 i=1,n
  93  bx(i)=xy(n+1-i)
    return
    end
    subroutine hilbert(ur,ui,mm,nt,l,ll,pi)
    implicit none
    real::pi,ss,x
    integer::nt,mm,l,ll
    real::h(ll),ur(nt,mm),ui(nt,mm)
    integer::i,j,j1,j2,j3,jj
    do j=1,nt
      do i=1,mm
        ui(j,i)=0.0
      enddo
    enddo
    h(l+1)=0.0
    do i=1,l
      ss=sin(pi*real(i)/2.0)**2
      h(l+1+i)=ss*2.0/(real(i)*pi)
      h(l+1-i)=-h(l+1+i)
    enddo
    j1=l+1
    j3=nt-l
    do i=1,mm                        
      do j=j1,j3
        x=0.0
        do jj=1,ll
          j2=j-jj+l+1
          x=x+ur(j2,i)*h(jj)
        enddo
        ui(j,i)=x
      enddo
    enddo              
    return
    end

huangyao

勿忘草 发表于 2016-7-18 10:01
是不是程序代码中没有给出文件类型

给出的是无格式文件，直接读取，我的代码已经给出来了，有点长，就看文件读取和写入两行就可以了，请大神指教

huangyao

自己已经解决，是recl应该再乘以一个4

jianglegejiang

{:5_213:}{:5_213:}{:5_213:}

jianglegejiang

{:5_214:}

tonysonic

很好，非常感谢

佳JIA

文件太长了吧，读不进去