用Fortran读取nc数据出现问题

xiuyuanyang

本人用fortran读取并处理nc数据，部分程序 如下：
Module constant
    Implicit None
    Save
    Integer,Parameter::nx=152,ny=110,nz=1,nfile=19,nt=72*nfile,nvar=5
    Real(Kind=4),Parameter::pi=3.1415926,r=6371.004
End Module constant
1Program Grapes_Cmaq
2   Use constant
3   Implicit None
4   !Include "/public/home/dengtao/atmosphere/netcdf/include/netcdf.inc"
5  Integer::i,j,k,l,m,n
6    Integer::ii_gz,jj_gz,ii_py,jj_py
7    !Integer,Parameter::nx=152,ny=110,nz=1,nt=72*19,nvar=5
8    !Real(Kind=4),Parameter::pi=3.1415926
9    Real(Kind=4)::lat_obs,lon_obs
10    Real(Kind=4)::lat(nx,ny,nz,1),lon(nx,ny,nz,1) !存储nc数据的经纬度，实际只有两维
11    Real(Kind=4)::p(nx,ny,nz,nt),t(nx,ny,nz,nt),q(nx,ny,nz,nt),ws(nx,ny,nz,nt),wd(nx,ny,nz,nt),rh(nx,ny,nz,nt) !存储nc数据中的气压、温度、混合比、风向、风速值
12    Real(Kind=4)::p_py(1,1,1,nt),t_py(1,1,1,nt),ws_py(1,1,1,nt),wd_py(1,1,1,nt),rh_py(1,1,1,nt)存储观测值
13    Real(Kind=4)::p_gz(1,1,1,nt),t_gz(1,1,1,nt),ws_gz(1,1,1,nt),wd_gz(1,1,1,nt),rh_gz(1,1,1,nt)
14   Real(Kind=4),Dimension(nvar)::meanerror_py,rel_py,rmse_py,meanerror_gz,rel_gz,rmse_gz !存储计算得到的观测地的模拟值与观测值的平均误差、相对湿度、均方根误差（2013年12月至2014年1月共57天）
15    Character(Len=100)::filename,str !存储nc文件名
16    Character,Dimension(nfile)::date
17    Real relative_humidity,average,sd,cc,meanerror,rmse !声明函数类型 此处出现错误：如果使用了implicit none 语句，则在函数和调用程序中都要声明函数的类型
18    date=(/'131204','131207','131210','131213','131216','131219','131222','131225','131228','131231',&
19    &'140103','140106','140109','140112','140115','140118','140121','140124','140127'/)
20    !读取地理信息数据并计算离观测站最近的点
21    filename='GRIDCRO2D_2013120400'
22    !print*,filename
23    Call readdata_from_nc(filename,lat,lon,p,t,q,ws,wd)
24    !open(30,file='lat.txt')
25    !write(30,'(f7.3)') ((lat(i,j,nz,1),i=1,nx),j=1,ny)
26    lat_obs=23.217;lon_obs=113.483
27    Call distance_of_twopoints(lat,lon,lat_obs,lon_obs,ii_gz,jj_gz)
28    Open(10,file='lat_lon.txt',status='new',access='append')
29    Write(*,'(1x,i3,1x,i3,1x,f7.3,1x,f7.3)') '离观测点1最近点的经纬度为：',ii_gz,jj_gz,lat(ii_gz,jj_gz,nz,1),lon(ii_gz,jj_gz,nz,1)
        Close（10）
30    lat_obs=22.933;lon_obs=113.317
31    Call distance_of_twopoints(lat,lon,lat_obs,lon_obs,ii_py,jj_py)
32    Open(10,file='lat_lon.txt',status='old',access='append')
33    Write(*,'(1x,i3,1x,i3,1x,f7.3,1x,f7.3)') '离观测站点2最近点的经纬度为：',ii_py,jj_py,lat(ii_py,jj_py,nz,1),lon(ii_py,jj_py,nz,1)
34   Close(10)
35End
读取nc地理数据的子程序
1Subroutine readgeodata_from_nc(filename,lat,lon)
2   Use constant
3   Implicit None
4   Include "/public/home/dengtao/atmosphere/netcdf/include/netcdf.inc"
5    Character(Len=100),Intent(In)::filename
6    Real(Kind=4),Intent(Out)::lat(:,:,:,:),lon(:,:,:,:)
7    Integer,Dimension(4)::start,count
8    Integer::ncid,latid,lonid,status
9    character(len=4)::str
10    data  start /1,1,1,1/
11    data count1 /nx,ny,nz,1/
12        status=nf_open(filename,nf_nowrite,ncid)
13        !print*,status,ncid
14        status=nf_inq_varid(ncid,'LAT',latid)
15       status=nf_inq_varid(ncid,'LON',lonid)
16        !print*,status
17        status=nf_get_var_real(ncid,latid,start,count,lat)!纬度
18       print*,status
19       status=nf_get_var_real(ncid,lonid,start,count,lon)!经度
20       status=nf_close(ncid)
  21  End Subroutine readgeodata_from_nc
读取nc气象数据的子程序（为了查错，只看读取地理数据的部分，因而这一部分注释掉）
!Subroutine readdata_from_nc(filename,p,t,q,ws,wd)
    !Use constant
    !Implicit None
    !Include "/public/home/dengtao/atmosphere/netcdf/include/netcdf.inc"
    !Character(Len=100),Intent(In)::filename
    !Real(Kind=4),Intent(Out)::p(:,:,:,:),t(:,:,:,:),q(:,:,:,:),ws(:,:,:,:),wd(:,:,:,:)
    !Integer,Dimension(4)::start,count
    !Integer::ncid,pid,tid,qid,wsid,wdid,status
    !data start /1,1,1,1/
    !data count /nx,ny,nz,nt/
    !    status=nf_open(filename_d,nf_nowrite,ncid)
    !    status=nf_inq_varid(ncid,'PRSFC',pid)    !气压
    !    status=nf_inq_varid(ncid,'TEMP2',tid)    !温度
    !    status=nf_inq_varid(ncid,'Q2',qid)       !混合比
    !    status=nf_inq_varid(ncid,'WSPD10',wsid)  !风速
    !    status=nf_inq_varid(ncid,'WDIR10',wdid)  !风向
    !    status=nf_get_vara_real(ncid,pid,start,count2,p_d)
    !    status=nf_get_vara_real(ncid,tid,start,count2,t_d)
    !    status=nf_get_vara_real(ncid,qid,start,count2,q_d)
    !    status=nf_get_vara_real(ncid,wsid,start,count2,ws_d)
    !    status=nf_get_vara_real(ncid,wdid,start,count2,wd_d)
    !    status=nf_close(ncid)
    !End Subroutine readdata_from_nc
计算离观测站点最近的点的子程序

        
Subroutine distance_of_twopoints(lat_d1,lon_d1,lat_obs_d,lon_obs_d,ii_d,jj_d)
    Use constant
    Implicit None
    Integer::i,j
    Integer,Intent(Out)::ii_d,jj_d
    Real(Kind=4),Intent(In)::lat_obs_d,lon_obs_d
    Real(Kind=4),Intent(In)::lat_d1(:,:,:,:),lon_d1(:,:,:,:)
    Real(Kind=4)::distance0,c0,distance,c
    c0=sin((90-lat_d1(1,1,1,1))*pi/180)*sin((90-lat_obs_d)*pi/180)*cos((lon_d1(1,1,1,1)-lon_obs_d)*pi/180)+cos((90-lat_d1(1,1,1,1))*pi/180)*cos((90-lat_obs_d)*pi/180)
    distance0=r*acos(c0)*pi/180
    Do j=1,ny
        Do i=1,nx
        c=sin((90-lat_d1(i,j,nz,1))*pi/180)*sin((90-lat_obs_d)*pi/180)*cos((lon_d1(i,j,nz,1)-lon_obs_d)*pi/180)+cos((90-lat_d1(i,j,nz,1))*pi/180)*cos((90-lat_obs_d)*pi/180)
        distance=r*acos(c)*pi/180
        if(distance<distance0)then
            distance0=distance
            ii_d=i
            jj_d=j
        endif
        enddo
    enddo
End Subroutine distance_of_twopoints
子程序是在Linux下用ifort编译连接，并生成了可执行文件main.exe，但是执行main.exe后，并没有如愿生成主程序中28行生成的输入结果的文件，而且在读取地理信息数据的子程序的第17行后，也没有输出status，但是前面打开nc文件，查询变量id号都可以输出status为零，说明打开nc数据和查询变量id号都没有问题，但是读取纬度这一变量就不行，导致主程序中没有生成任何输出结果的文件，不知道为什么，我实在不知道自己哪里写的不对，之前这么读nc数据就没有问题，这回只是把读取过程写成了一个子程序，就一直出现这种情况，哪位前辈可以指点一下，谢谢！！！

zhchlue

感谢楼主的分享

xiuyuanyang

zhchlue 发表于 2015-8-24 16:45
感谢楼主的分享

同学，我这个程序还有问题，你用的时候小心点吧

Richeho

很实用，很有用的东西

tianmeng

data count1 /nx,ny,nz,1/
12        status=nf_open(filename,nf_nowrite,ncid)
13        !print*,status,ncid
14        status=nf_inq_varid(ncid,'LAT',latid)
15       status=nf_inq_varid(ncid,'LON',lonid)
16        !print*,status
17        status=nf_get_var_real(ncid,latid,start,count,lat)!纬度

前面定义是count1，后面用count