fortran 读取nc文件后，数据不对，不知道错误出在哪里，求助！

aiai1992

nc数据文件为，hgt.1950.nc
程序：
  program read_netcdf

      include 'netcdf.inc'
        integer i,j,k,l,irec

      integer*4  ncid, status    ! file control
!-------------------------------------------------------------

      real*4      :: lon(144)
      real*4      :: lat(73),level(17)
      integer*4   :: time(365)
      real*4      :: hgt(144,73,17,365)
!-------------------------------------------------------------
      integer*4   :: start(10)
      integer*4   :: count(10)
      integer     :: dimids(5)! allow up to 10 dimensions
      integer     :: dimid, xtype
      character(len = 31) :: dummy

!   Open netCDF file.
      status = nf_open('hgt.1950.nc', nf_nowrite, ncid)
      if ( status /= nf_noerr) call handle_err(status)

!----------------------------------------------------

      status=nf_inq_varid(ncid, 'lon', lonid)
          if (status /= nf_noerr) call handle_err(status)

      status = nf_inq_var(ncid, lonid, dummy, xtype, ndim, dimids, natts)
      if (status /= nf_noerr) call handle_err(status)

     do j = 1, ndim
        status = nf_inq_dim(ncid, dimids(j), dummy, len)
        if (status /= nf_noerr) call handle_err(status)
        start(j) = 1
            count(j) = len
      enddo

      status = nf_get_vara_real(ncid, lonid, start, count, lon)

!----------------------------------------------------
!   Retrieve data for Variable 'lat'
!print*,ncid
      status=nf_inq_varid(ncid, 'lat', latid)
          if (status /= nf_noerr) call handle_err(status)

      status = nf_inq_var(ncid, latid, dummy, xtype, ndim, dimids, natts)
      if (status /= nf_noerr) call handle_err(status)
      do j = 1, ndim
        status = nf_inq_dim(ncid, dimids(j), dummy, len)
        if (status /= nf_noerr) call handle_err(status)
        start(j) = 1
            count(j) = len
      enddo
      status = nf_get_vara_real(ncid, latid, start, count, lat)

  !----------------------------------------------------
  !  Retrieve data for Variable 'level'
!print*,ncid
      status=nf_inq_varid(ncid, 'level', levelid)
          if (status /= nf_noerr) call handle_err(status)

     status = nf_inq_var(ncid, levelid, dummy, xtype, ndim, dimids, natts)
    if (status /= nf_noerr) call handle_err(status)

     do j = 1, ndim
       status = nf_inq_dim(ncid, dimids(j), dummy, len)
       if (status /= nf_noerr) call handle_err(status)
       start(j) = 1
            count(j) = len
     enddo
! print*,ncid

    status = nf_get_vara_real(ncid, levelid, start, count, level)

!----------------------------------------------------
!   Retrieve data for Variable 'time'
    status=nf_inq_varid(ncid,'time',timeid)
          if (status /= nf_noerr) call handle_err(status)

     status = nf_inq_var(ncid, timeid, dummy, xtype, ndim, dimids, natts)
    if (status /= nf_noerr) call handle_err(status)
     do j = 1, ndim
      status = nf_inq_dim(ncid, dimids(j), dummy, len)
     if(status /= nf_noerr) call handle_err(status)
      start(j) = 1
      count(j) = len
     enddo

     status = nf_get_vara_int(ncid, timeid, start, count, time)

!----------------------------------------------------
!   Retrieve data for Variable 'hgt'

      status=nf_inq_varid(ncid, 'hgt', hgtid)
          if (status /= nf_noerr) call handle_err(status)

      status = nf_inq_var(ncid, hgtid, dummy, xtype, ndim, dimids, natts)
      if (status /= nf_noerr) call handle_err(status)
          write(*,*) dummy, xtype , ndim ,dimids, natts
!          pause
      do j = 1,ndim
        status = nf_inq_dim(ncid, dimids(j), dummy, len)
        if (status /= nf_noerr) call handle_err(status)
        start(j) = 1
            count(j) = len
      enddo

      status = nf_get_vara_real(ncid, hgtid, start, count, hgt)

      status = nf_close(ncid)
      if(status /= nf_noerr) call handle_err(status)
!----------------------------------------------------

print*,hgt(1,1,1,1)
end      

!---------------------- End Program  --------------------------
subroutine handle_err(status)

  include 'netcdf.inc'
  integer, intent ( in) :: status

  if(status /= nf_noerr) then
    print *, trim(nf_strerror(status))
  stop "Stopped"
  end if

end subroutine handle_err

结果：：：


第二个图是我用论坛里面那个read_netcdf.f90读出来的，两次结果一样，但是位势高度不应该是这样的数值啊？我想知道我问题出在哪里了！

aiai1992

@mofangbao

longlivehj

应该是数据还没有解包。还需获取hgt的add_offset和scale_factor属性，对hgt作一个简单的运算。具体参考：http://www.esrl.noaa.gov/psd/data/gridded/faq.html#2

aiai1992

longlivehj 发表于 2013-12-19 15:59

                               
登录/注册后可看大图

应该是数据还没有解包。还需获取hgt的add_offset和scale_factor属性，对hgt作一个简单的运算。具体参考：ht ...

十分感谢！

LOSER

您好，add_offset和scale_factor的值是怎么获取的？

aiai1992

LOSER 发表于 2014-6-16 16:43
您好，add_offset和scale_factor的值是怎么获取的？

我是用PanoplyWin这个小软件获取的，应该还有其他办法

addison

netcdf本身就有这个函数，请用我的小程序重新生成读netcdf数据的fortran程序，就ok了。
见：http://bbs.06climate.com/forum.php?mod=viewthread&tid=25183

LOSER

aiai1992 发表于 2014-6-17 07:19
我是用PanoplyWin这个小软件获取的，应该还有其他办法

谢谢啦～用你的方法解决问题啦～

半个瞎子

aiai1992 发表于 2013-12-19 18:38
十分感谢！

我也来补充一下，很多的nc数据为了节省空间会把实型数据转换为整型数据存储，在用GrADS打开这种nc文件时，GrADS会自动转换为实型数据，但用fortran读这种数据时，需要做一个简单的转换： 整型数据*scale_factor+add_offset=实行数据。

爱知天气

想问以下，status = nf_open('hgt.1950.nc', nf_nowrite, ncid)
      if ( status /= nf_noerr)  中的  nf_open，nf_nowrite，if ( status /= nf_noerr) 是什么意思呢，用来干什么的