- 积分
- 449
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2013-4-19
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
求助各位大神,用GRADS转化nc资料为grd:uwnd.2102.nc和vwnd.2012.nc
程序如下:
'reinit'
'sdfopen g:\lunwen\uwnd.2012.nc'
'sdfopen g:\lunwen\vwnd.2012.nc'
'set fwrite g:\lunwen\u.grd'
'set gxout fwrite'
'set fwrite g:\lunwen\v.grd'
'set gxout fwrite'
'set lon 80 220'
'set lat -10 80'
'set z 3'
it=210
while(it<=217)
'set t 'it''
'd uwnd'
it=it+1
'd vwnd'
endwhile
'disable fwrite'
;
写是这么写的实际是分开来生成grd文件的,然后用FORTRAN编写程序计算涡度:
program main
implicit none
integer i,j,nx,ny,it,n
parameter(nx=57,ny=37,it=8)
real u(nx,ny,it),v(nx,ny,it),lat(ny)
real::a=6371000.0,pi=3.14
real vox1(nx,ny,it),vox2(nx,ny,it),vox(nx,ny,it)
real vox3(nx,ny,it),vox4(nx,ny,it),vox5(nx,ny,it)
open(1,file='g:\lunwen\vox.grd',form='binary')
open(2,file='g:\lunwen\u.grd',form='binary')
open(3,file='g:\lunwen\v.grd',form='binary')
do n=1,it
do j=1,ny
read(2)(u(i,j,n),i=1,nx)
read(3)(v(i,j,n),i=1,nx)
end do
end do
close(2);close(3)
do j=1,ny
lat(j)=10.0+float(j-1)*1.0
enddo
do n=1,it
do j=2,ny-1
do i=2,nx-1
vox1(i,j,n)=(v(i+1,j,n)-v(i-1,j,n))/((pi/180.0)*cosd(lat(j)))
vox2(i,j,n)=(u(i+1,j,n)-u(i-1,j,n))/(pi/180.0)
vox3(i,j,n)=(1.0/(2.0*a))*(vox1(i,j,n)-vox2(i,j,n)+2.0*u(i,j,n)*tand(lat(j)))
vox4(i,j,n)=(-1.0)*(u(i,j,n)*(vox3(i+1,j,n)-vox3(i-1,j,n))/(2.0*a*(pi/180.0)*cosd(lat(j))))
vox5(i,j,n)=(-1.0)*(v(i,j,n)*(vox3(i+1,j,n)-vox3(i-1,j,n))/(2.0*a*(pi/180.0)))
vox(i,j,n)=(vox4(i,j,n)+vox5(i,j,n))*(10**10)
enddo
enddo
enddo
write(1)vox
end
得出了vox.grd然后写ctl画图:
dset g:\lunwen\vox.grd
undef -9.99E+33
title NCEP/NCAR REANALYSIS PROJECT
xdef 57 linear 80.000 2.500
ydef 37 linear -10.000 2.500
zdef 1 levels 850
tdef 8 linear 00:00z28JUL2012 24hr
vars 1
vox 11 99 temprature
endvars
然后画图就只能画出t=1的图。。。。。各位我是哪里出错了呢?帮个忙呀!
|
|