在fortran中正确读取python生成的存有三维数组的二进制文件的方法

Aristuo

      最近用十几年前的fortran代码跑分段位涡反演。由于本人实在不会fortran，于是尝试将输入场在python中整合好，存为文件，直接在fortran中使用。下面是我将三维场（shape=（ks,ns,ms），其中ks为垂直层数，ns为同一经度下的纬度格点数，ms为同一纬度下的经度格点数）存入二进制文件的代码，使用了numpy库中的tofile函数：

    hgt.astype(np.float32).tofile(f'hgt.bin')


      随后，我在fortran中对hgt.bin文件进行读取。读取方法如下：

      real(kind=4), allocatable :: hgt(:,:,:) !声明变量
      allocate(hgt(ks,ns,ms), stat=ierr) !分配内存
      call read_binary_file('hgt.bin', hgt, 'hgt')  !读取hgt.bin文件到hgt变量中。

     划重点！在python保存和在fortran中读取时，必须设置同样的数据精度（标红部分），否则读出来的值会严重错误！我这里统一使用的是单精度，如果想使用双精度，则在python保存时改为astype(np.float64)，fortran读取时改为kind=8。
[size=13.0667px]      统一数据精度后，fortran读取出来的数值是正确的了。然而，数据的维度却出现了问题。这是因为python和fortran的默认存储顺序是不一致的，具体怎么不一样的不用管，我直接给出解决方案。python存储代码不变，照常存储即可。fortran端代码做如下改动：
[size=13.0667px]

      real(kind=4), allocatable :: hgt0(:,:,:), hgt(:,:,:) !声明变量
      allocate(hgt0(ms,ns,ks), stat=ierr) !分配内存
      allocate(hgt(ks,ns,ms), stat=ierr) !分配内存
      call read_binary_file('hgt.bin', hgt0, 'hgt')  !读取hgt.bin文件到hgt0变量中。

      do k=1,ks
        do j=1,ns
          do i=1,ms
            hgt(k,j,i) = hgt0(i,j,k)
          end do
        end do
      end do


[size=13.0667px]      此时的hgt变量，数据结构就跟python中的相同了。hgt0变量只是一个过渡，后续无需再使用。注意分配内存时的维度设置。

      用到的read_binary_file函数如下：

      subroutine read_binary_file(filename, array, varname)
          character(len=*), intent(in) :: filename, varname
          real(kind=4), intent(out) :: array(:,:,:)
          integer :: fileunit, ierr
          integer(kind=4) :: reclen

          inquire(iolength=reclen) array
          open(newunit=fileunit, file=filename, form='unformatted', &
               access='stream', status='old', iostat=ierr) !, recl=reclen

          if (ierr == 0) then
              read(fileunit, iostat=ierr) array !rec=1
              close(fileunit)
              if (ierr /= 0) then
                  print*, 'Error reading 1', trim(varname), ' file'
                  stop
              endif
          else
              print*, 'Error opening 2', trim(varname), ' file: ', trim(filename)
              stop
          endif
      end subroutine

kearon

你可以用numpy的f2py把Fortran代码转成so，直接用python调用逻辑就可以了