[文件读写] 语法错误读取不了数据，求助

program main 
 implicit none
        integer i,j,k,l,m,n
        Integer,Parameter :: Nxm   = 224, Nxp  = 1057, Nx = Nxm + Nxp
        Integer,Parameter :: Ny    = 161
        Integer,Parameter :: Nz    = 161
        Real*8              U(-2:Nx+2,-2:Ny+2,-2:Nz+2)
        Real*8              V(-2:Nx+2,-2:Ny+2,-2:Nz+2)
        Real*8              W(-2:Nx+2,-2:Ny+2,-2:Nz+2)
        Real*8              P(-2:Nx+2,-2:Ny+2,-2:Nz+2)
        Real*8              T(-2:Nx+2,-2:Ny+2,-2:Nz+2)
        Real*8              dudx(-2:Nx+2,-2:Ny+2,-2:Nz+2)
        Real*8              dudy(-2:Nx+2,-2:Ny+2,-2:Nz+2)
        Real*8              dudz(-2:Nx+2,-2:Ny+2,-2:Nz+2)
        Real*8              dvdx(-2:Nx+2,-2:Ny+2,-2:Nz+2)
        Real*8              dvdy(-2:Nx+2,-2:Ny+2,-2:Nz+2)
        Real*8              dvdz(-2:Nx+2,-2:Ny+2,-2:Nz+2)
        Real*8              dwdx(-2:Nx+2,-2:Ny+2,-2:Nz+2)
        Real*8              dwdy(-2:Nx+2,-2:Ny+2,-2:Nz+2)
        Real*8              dwdz(-2:Nx+2,-2:Ny+2,-2:Nz+2)

        Real*8              dd  (3,3,-2:Nx+2,-2:Ny+2,-2:Nz+2)
        Real*8             Omega(3,  -2:Nx+2,-2:Ny+2,-2:Nz+2)
        integer(kind=4)               :: u_file
        integer(kind=4)               :: v_file
        integer(kind=4)               :: w_file
        Real*8              Levi(1:3,1:3,1:3) 
        Real*8,Parameter :: pi  = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825d0
        Real*8,Parameter :: Ldelta = 2.0d0 * pi
        Real*8,Parameter :: Lx  = 14.4* Ldelta
        Real*8,Parameter :: Ly  = 1.0d0 * Ldelta
        Real*8,Parameter :: Lz  = 1.0d0 * Ldelta
        Real*8,Parameter :: Lxm = 1.2 * Ldelta
        Real*8,Parameter :: Lxp = Lx - Lxm
        Real*8,Parameter :: dy = Ly / ( 161.0d0 - 1.0d0 )
        Real*8,Parameter :: dz = Lz / ( 161.0d0 - 1.0d0 )
        Real*8              dx (-1:Nx)
        Real*8              c1_dx (0:Nx)
        Real*8              c1_dx2(0:Nx)
        Real*8              x(0:Nx)
        Real*8              y(0:Ny)
        Real*8              z(0:Nz)

        Real*8              dxm(0:Nx)
        Real*8              dxhm(0:Nx)
        Real*8              c1_dxm(0:Nx)
        Real*8              c1_dxhm(0:Nx)
        Real*8              dxp(0:Nx)
        Real*8              dxhp(0:Nx)
        Real*8              c1_dxp(0:Nx)
        Real*8              c1_dxhp(0:Nx)

        Real*8              dxm_dxp(Nx-1)
        Real*8              c1_dxhpxm(Nx-1)
        Real*8              dxp_dxm(Nx-1)
        Real*8              c1_dxhmxp(Nx-1)
        Real*8,Parameter :: alp_x1  = 0.9d0
        Real*8,Parameter :: alp_x2  = 0.2d0
        Real*8,Parameter :: alpha_x  = ( Lxp / (1057.0d0 - 1.0d0) ) - ( alp_x1 * Lxm / 224.0d0 )
        Real*8,Parameter :: c1_4       = 1.0d0 /   4.0d0
        Real*8              x_alf, x_Nxm
        Real*8,Parameter :: c1_dz      = 1.0d0 /   dz
        Real*8,Parameter :: c1_dy      = 1.0d0 /   dy
        integer             iter, Count, rl ,Nt
        real*8              time
 
       
        Open(40, File = 'fgg0126Ob_ityp_UVWPT_bu', Form = 'unFormatted', Status = 'unknown')
        read(40) Nt, time
        read(40) U
        read(40) V
        read(40) W
        read(40) P
        read(40) T

        Close(40)
        
        write(*,*) Nt

        x(0) = 0.0d0

         Do i = 1, Nxm, 1
          Select Case(i)
          Case(            : Nxm  /4 )
          x(i) = x(i-1) + ( (alp_x2 / Dble( Nxm/4 )) * Dble(i) + ( 1.0d0 - alp_x2 ) )  &
                        * ( alpha_x * ( dcos(2.0d0 * pi * (Dble(i - 1) )/Dble(Nxm) ) ) **2 + (Lxp/Dble(Nxp -1) - alpha_x) )
          Case( Nxm/4   +1 : Nxm*3/4 )
          x(i) = x(i-1) + ( alpha_x * ( dcos(2.0d0 * pi * (Dble(Nxm/4) )/Dble(Nxm) ) ) **2 + (Lxp/Dble(Nxp -1) - alpha_x) )
          Case( Nxm*3/4 +1 : Nxm     )
          x(i) = x(i-1) + ( alpha_x * ( dcos(2.0d0 * pi * (Dble(i - 1) )/Dble(Nxm) ) ) **2 + (Lxp/Dble(Nxp -1) - alpha_x) )
          End Select   
        End Do

        x_alf = 2.0d0
        x_Nxm = x(Nxm)


        Do i = Nxm, Nx, 1
          x (i) = (i - Dble(Nxm)) * ( Lxp / Dble(Nxp -1) ) + Lxm
        End Do

        dx(0) = x(0) - ( 2.0d0 * x(0) - x(1) )
        Do i = 1, Nx, 1
          dx(i) = x(i) - x(i-1)  
        End Do

        !$OMP PARALLEL DO PRIVATE(i) SHARED(dxm,dxp,c1_dx,c1_dx2)
        Do i = 0, Nx, 1
          dxm(i)    = 0.5d0 * dx(i)
          dxp(i)    = 0.5d0 * dx(i)
          c1_dx (i) = 1.0d0 /  dx(i)
          c1_dx2(i) = 1.0d0 / (dx(i  ) * dx(i))
        End Do
        !$OMP End PARALLEL DO

        !$OMP PARALLEL DO PRIVATE(i) SHARED(dxhm,dxhp)
        Do i = 1, Nx-1, 1
          dxhm(i)     = ( dx(i-1) + dx(i  ) ) * 0.5d0
          dxhp(i)     = ( dx(i  ) + dx(i+1) ) * 0.5d0
        End Do
        !$OMP End PARALLEL DO

        !$OMP PARALLEL DO PRIVATE(i) SHARED(c1_dxm,c1_dxp,c1_dxhm,c1_dxhp)
        Do i = 1, Nx-1, 1
          c1_dxm (i)  = 1.0d0 /   dxm (i)
          c1_dxp (i)  = 1.0d0 /   dxp (i)
          c1_dxhm(i)  = 1.0d0 /   dxhm(i)
          c1_dxhp(i)  = 1.0d0 /   dxhp(i)
        End Do
        !$OMP End PARALLEL DO


!       ------ Velocity Gradient------
        
        Do k =1,  Nz-1
         Do j = 1,  Ny-1 
           Do i =   1, Nx-1, 1  
!       ------ (i,j,k) = (1,1)------
          dudx(i,j,k) =             (    - U(i-1,j  ,k) + U(i  ,j  ,k)    )*c1_dx(i)        
!       ------ (i,j,k) = (1,2)------
          dudy(i,j,k) = c1_4 * ( +  (    - U(i-1,j-1,k) + U(i-1,j,k  )    )*c1_dy                 &  
                                 +  (    - U(i-1,j,k  ) + U(i-1,j+1,k)    )*c1_dy                 &     
                                 +  (    - U(i  ,j-1,k) + U(i  ,j,k  )    )*c1_dy                 &
                                 +  (    - U(i  ,j  ,k) + U(i,  j+1,k)    )*c1_dy       )         
!       ------ (i,j,k) = (1,3)------
          dudz(i,j,k) = c1_4 * ( +  (    - U(i-1,j,k-1) + U(i-1,j,k  )    )*c1_dz                 & 
                                 +  (    - U(i-1,j,k  ) + U(i-1,j,k+1)    )*c1_dz                  &           
                                 +  (    - U(i  ,j,k-1) + U(i  ,j,k  )    )*c1_dz           & 
                                 +  (    - U(i  ,j,k  ) + U(i,  j,k+1)    )*c1_dz       )      
!       ------ (i,j,k) = (2,1)------
          dvdx(i,j,k) = c1_4 * ( +  (    - V(i-1,j-1,k) + V(i  ,j-1,k)    )*c1_dxhm(i)            &
                                 +  (    - V(i  ,j-1,k) + V(i+1,j-1,k)    )*c1_dxhp(i)            &
                                 +  (    - V(i-1,j  ,k) + V(i  ,j  ,k)    )*c1_dxhm(i)            &
                                 +  (    - V(i  ,j  ,k) + V(i+1,j  ,k)    )*c1_dxhp(i)  )         
!       ------ (i,j,k) = (2,2)------
          dvdy(i,j,k) =             (    - V(i  ,j-1,k) + V(i  ,j  ,k)    )*c1_dy  
!       ------ (i,j,k) = (2,3)------
          dvdz(i,j,k) = c1_4 * ( +  (    - V(i,j-1,k-1) + V(i,j-1,k  )    )*c1_dz          &
                                 +  (    - V(i,j-1,k  ) + V(i,j-1,k+1)    )*c1_dz               &        
                                 +  (    - V(i,j  ,k-1) + V(i,j  ,k  )    )*c1_dz              &
                                 +  (    - V(i,j  ,k  ) + V(i,j  ,k+1)    )*c1_dz       )          

!       ------ (i,j,k) = (3,1)------
          dwdx(i,j,k) =    c1_4 * ( +  (    - W(i-1,j,k-1) + W(i  ,j,k-1)    )*c1_dxhm(i)            &
                                 +  (    - W(i  ,j,k-1) + W(i+1,j,k-1)    )*c1_dxhp(i)            &
                                 +  (    - W(i-1,j,k  ) + W(i  ,j,k  )    )*c1_dxhm(i)            &
                                 +  (    - W(i  ,j,k  ) + W(i+1,j,k  )    )*c1_dxhp(i)  )  

!       ------ (i,j,k) = (3,2)------
          dwdy(i,j,k) =    c1_4 * ( +  (    - W(i,j-1,k-1) + W(i,j  ,k-1)    )*c1_dy          &
                                 +  (    - W(i,j  ,k-1) + W(i,j+1,k-1)    )*c1_dy              &
                                 +  (    - W(i,j-1,k  ) + W(i,j  ,k  )    )*c1_dy              &
                                 +  (    - W(i,j  ,k  ) + W(i,j+1,k  )    )*c1_dy       )           
!       ------ (i,j,k) = (3,3)------                 
          dwdz(i,j,k) =                (    - W(i  ,j,k-1) + W(i  ,j  ,k)    )*c1_dz
        
       End Do
         End Do 
        End Do

        
        Do k = 145,  176
          Do j = 145,  176 
            Do i =   1, Nx-1, 1       
          dd(1,1,i,j,k) = dudx(i,j,k)
          dd(1,2,i,j,k) = dudy(i,j,k)
          dd(1,3,i,j,k) = dudz(i,j,k)
          dd(2,1,i,j,k) = dvdx(i,j,k)
          dd(2,2,i,j,k) = dvdy(i,j,k)
          dd(2,3,i,j,k) = dvdz(i,j,k)
          dd(3,1,i,j,k) = dwdx(i,j,k)
          dd(3,2,i,j,k) = dwdy(i,j,k)
          dd(3,3,i,j,k) = dwdz(i,j,k)
            End Do
          End Do
        End Do

!       ------ Levi-Civita ------
    Do k = 1, 3, 1
    Do j = 1, 3, 1
        Do i = 1, 3, 1
          levi(i,j,k) = 0.0d0  
    End Do; End Do; End Do  

        levi(1,2,3) = 1.0d0 
        levi(3,1,2) = 1.0d0 
        levi(2,3,1) = 1.0d0 

        levi(1,3,2) = -1.0d0 
        levi(2,1,3) = -1.0d0 
        levi(3,2,1) = -1.0d0      



        Omega(:,:,:,:) = 0.0d0
        
        Do k = 1, Nz-1
          Do j = 1, Ny-1 
            Do i = 1, Nx-1, 1     
          Do l = 1, 3       
          Do m = 1, 3 
          Do n = 1, 3  
            Omega(l,i,j,k) = Omega(l,i,j,k) + levi(l,m,n) *  dd(m,n,i,j,k) 
          End Do; End Do; End Do        
        End Do
          End Do
        End Do
end program