- 积分
- 3876
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2013-4-4
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
本人用fortran读取并处理nc数据,遇到了问题,求助各位大神,程序是在Linux下用ifort编译连接,并生成了可执行文件main.exe,但是执行main.exe后,并没有如愿生成主程序中28行生成的输入结果的文件,而且在读取地理信息数据的子程序的第17行后,也没有输出status,但是前面打开nc文件,查询变量id号都可以输出status为零,说明打开nc数据和查询变量id号都没有问题,但是读取纬度这一变量就不行,导致主程序中没有生成任何输出结果的文件,不知道为什么,我实在不知道自己哪里写的不对,之前这么读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
|
|