爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 11244|回复: 13

[求助] fortran读netcdf文件求解[已经解决]

[复制链接]

新浪微博达人勋

发表于 2013-3-2 00:47:19 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 infocert 于 2013-3-2 14:41 编辑

有个nc文件,头格式如下:
  1. netcdf VEGET {
  2. dimensions:
  3.         lon = 64 ;
  4.         lat = 32 ;
  5.         time = UNLIMITED ; // (1992 currently)
  6. variables:
  7.         float lon(lon) ;
  8.                 lon:long_name = "longitude coordinate" ;
  9.                 lon:standard_name = "longitude" ;
  10.                 lon:units = "degrees_east" ;
  11.                 lon:axis = "X" ;
  12.         float lat(lat) ;
  13.                 lat:long_name = "latitude coordinate" ;
  14.                 lat:standard_name = "latitude" ;
  15.                 lat:units = "degrees_north" ;
  16.                 lat:axis = "Y" ;
  17.         float time(time) ;
  18.                 time:long_name = "time" ;
  19.                 time:standard_name = "time" ;
  20.                 time:units = "years since 0001-01-01 00:00:00" ;
  21.                 time:time_origin = "01-JAN-0001 00:00:00" ;
  22.         float vegfrac(time, lat, lon) ;
  23.                 vegfrac:long_name = "fraction de culture par maille" ;
  24.                 vegfrac:standard_name = "fraction_culture" ;
  25.                 vegfrac:units = "%" ;
  26.                 vegfrac:missing_value = -99.99f ;
  27. }
复制代码

我读取这个nc所用的fortran如下[已经修改为正确代码]:
  1. PROGRAM crop_grille_ecbilt

  2.   USE NETCDF
  3.   implicit none

  4.   !variable pour l'ouverture ecriture du netcdf
  5.   character(256) :: filenameIn, filenameOut
  6.   integer :: intputID, unlimDimID, varID, ID, nbTime, iyr

  7.   !variable pour le calcul des moyennes, etc...
  8.   real*4, dimension(:), allocatable :: lon, lat, time, ValueVar
  9.   real*8, dimension(:,:,:), allocatable :: vegfrac

  10.   !variables pour le changement de grille
  11.   !real*8 farea(64,32,1)
  12.   integer nlat,nlon

  13.   integer :: outputID
  14.   integer :: outdimid, outvarid

  15.   if(iargc().eq.2) then
  16.     call getarg(1,filenameIn)
  17.     call getarg(2,filenameOut)
  18.   else
  19.     write(*,'(A)') "Usage: ./ConvertVegetYo [input.nc] [output.nc]"
  20.     stop
  21.   endif

  22.   write(*,*) "Open ", trim(filenameIn)
  23.   call check(nf90_open(filenameIn, nf90_nowrite, intputID))
  24.   call check(nf90_inquire(intputID, unlimitedDimId = unlimDimID))

  25.   write(*,*) "Read longitude"
  26.   call check(nf90_inq_varid(intputID, "lon", ID))
  27.   allocate(lon(64))
  28.   call check(nf90_get_var(intputID, ID, lon, (/ 1 /), (/ 64 /)))

  29.   write(*,*) "Read latitude"
  30.   call check(nf90_inq_varid(intputID, "lat", ID))
  31.   allocate(lat(32))
  32.   call check(nf90_get_var(intputID, ID, lat, (/ 1 /), (/ 32 /)))

  33.   call check(nf90_inq_varid(intputID, "time", ID))
  34.   allocate(time(1992))
  35.   call check(nf90_get_var(intputID, ID, time, (/ 1 /), (/ 1992 /)))

  36.   call check(nf90_inq_varid(intputID, "vegfrac", varID))
  37.   allocate(vegfrac(64,32,1992))
  38.   call check(nf90_get_var(intputID, varID, vegfrac, (/ 1, 1, 1 /), (/ 64,32,1992 /)))
  39.   
  40.   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  41.   !!  Create netcdf START
  42.   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  43.   !write target file
  44.   write(*,*) "Create ", trim(filenameOut)
  45.   call check(nf90_create(filenameOut, nf90_clobber, outputID))
  46.   
  47.   !define dimension
  48.   call check(nf90_def_dim(outputID, "lon", 64, outdimid))
  49.   call check(nf90_def_var(outputID, "lon", NF90_FLOAT, (/ 1 /), outvarid))
  50.   call check(nf90_put_att(outputID, outvarid, "long_name", "longitude coordinate"))
  51.   call check(nf90_put_att(outputID, outvarid, "standard_name", "longitude"))
  52.   call check(nf90_put_att(outputID, outvarid, "units", "degrees_east"))
  53.   call check(nf90_put_att(outputID, outvarid, "axis", "X"))

  54.   call check(nf90_def_dim(outputID, "lat", 32, outdimid))
  55.   call check(nf90_def_var(outputID, "lat", NF90_FLOAT, (/ 2 /), outvarid))
  56.   call check(nf90_put_att(outputID, outvarid, "long_name", "latitude coordinate"))
  57.   call check(nf90_put_att(outputID, outvarid, "standard_name", "latitude"))
  58.   call check(nf90_put_att(outputID, outvarid, "units", "degrees_north"))
  59.   call check(nf90_put_att(outputID, outvarid, "axis", "Y"))
  60.   
  61.   !NF90_UNLIMITED
  62.   call check(nf90_def_dim(outputID, "time", NF90_UNLIMITED, outdimid))
  63.   call check(nf90_def_var(outputID, "time", NF90_FLOAT, (/ 3 /), outvarid))
  64.   call check(nf90_put_att(outputID, outvarid, "long_name", "time"))
  65.   call check(nf90_put_att(outputID, outvarid, "standard_name", "time"))
  66.   call check(nf90_put_att(outputID, outvarid, "units", "years since 0001-01-01 00:00:00"))
  67.   call check(nf90_put_att(outputID, outvarid, "time_origin", "01-JAN-0001 00:00:00"))

  68.   !define variable
  69.   call check(nf90_def_var(outputID, "vegfrac", NF90_FLOAT, (/ 1, 2, 3 /), outvarid))
  70.   call check(nf90_put_att(outputID, outvarid, "long_name", "fraction de culture par maille"))
  71.   call check(nf90_put_att(outputID, outvarid, "standard_name", "fraction_culture"))
  72.   call check(nf90_put_att(outputID, outvarid, "units", "%"))
  73.   call check(nf90_put_att(outputID, outvarid, "missing_value", real(-99.99)))

  74.   !finish the configuration of the output file
  75.   call check(nf90_enddef(outputID))
  76.   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  77.   !! Create END
  78.   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  79.   do iyr=1,1992
  80.     do nlat=1,32
  81.       do nlon=1,64
  82.         !farea(nlon,nlat,iyr)=vegfrac(nlon,nlat,1)
  83.         write(*,*) nlon,",",nlat,",",iyr," => ", farea(nlon,nlat,iyr)
  84.       end do
  85.     end do
  86.   end do
  87.   
  88.   call check(nf90_close(intputID))

  89. END PROGRAM crop_grille_ecbilt


  90. SUBROUTINE check(status)
  91.   USE NETCDF
  92.   IMPLICIT NONE

  93.   INTEGER, INTENT (IN) :: status
  94.   if(status /= nf90_noerr) then
  95.     write(*,*)"Error  : ", trim(nf90_strerror(status))
  96.     stop
  97.   end if
  98. END SUBROUTINE check
复制代码


我想问,这段代码到底对不对?编译成功,运行时,只打印出第一年的数据,第二年的数为空,然后程序就出错
出错信息为:
  1. forrtl: severe (174): SIGSEGV, segmentation fault occurred
  2. Image              PC                Routine            Line        Source
  3. a.out              0000000000401740  Unknown               Unknown  Unknown
  4. a.out              0000000000400E1C  Unknown               Unknown  Unknown
  5. libc.so.6          00002B1417944CDD  Unknown               Unknown  Unknown
  6. a.out              0000000000400D19  Unknown               Unknown  Unknown
复制代码
请论坛中的前辈帮忙解解,多谢了!!!

改变读数据方式后成功,读vegfrac时不在循环里读,直接一次性读出来



密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2013-3-2 01:46:52 | 显示全部楼层
我一直就没安装上netcdf,不知道你有安装的软件可以传我一份吗?多谢了!
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2013-3-2 09:39:44 | 显示全部楼层
第一年的数据处理的正确吗???
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2013-3-2 09:54:55 | 显示全部楼层

我就是直接在大型机上安装的netcdf4.2/intel 12.1.4,你可以试着自己安装看看,不会的地方可以联系我
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2013-3-2 09:56:13 | 显示全部楼层
qxtlyf 发表于 2013-3-2 09:39
第一年的数据处理的正确吗???

没太仔细看,全数据确实是nc里的数据,不像第二年全是0.000000000
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2013-3-2 10:03:47 | 显示全部楼层
本帖最后由 infocert 于 2013-3-2 10:16 编辑

修改第53,54行为:
  1. farea(nlon,nlat,iyr)=vegfrac(nlon,nlat,1)
  2.         write(*,*) nlon,",",nlat,",",iyr," => ", farea(nlon,nlat,1)
复制代码
现在第二年的结果不再是0.000000了,但第二年只打印几十行,就出现1楼原贴的错误信息,求解,请各位前辈支招,小弟初学,非常感谢!!!
还有个问题,nlon,nlat,iyr这三个变量只是lon,lat,time数组的索引,farea(nlon,nlat,iyr)是不是应该变成farea(lon(nlon),lat(nlat),time(iyr)),同样vegfrac(nlon,nlat,1)是不是也应该改变?
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2013-11-26 14:20:42 | 显示全部楼层
你好 请问你这个只是读取netcdf数据的文件吗?  有没有做过修改某些数值后 重新替换掉原来的netcdf文件?
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2013-11-27 12:41:39 | 显示全部楼层
火_旦 发表于 2013-11-26 14:20
你好 请问你这个只是读取netcdf数据的文件吗?  有没有做过修改某些数值后 重新替换掉原来的netcdf文件?

别用fortran做这些,太麻烦,用ncl比较快
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2013-11-27 15:50:40 | 显示全部楼层
infocert 发表于 2013-11-27 12:41
别用fortran做这些,太麻烦,用ncl比较快

好的  谢谢您的建议  试一下NCL
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2015-12-9 16:57:33 | 显示全部楼层
请问你的问题后来是怎么解决的?
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

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