大神帮帮忙，我用NCEP高度场月平均资料集导出500hpa高度场保存成grd文件时，文件大为0

字俣丞 · 发表于 2018-1-21 11:11:24

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

您需要登录才可以下载或查看，没有帐号？立即注册

x

我用NCEP高度场79-2016年月平均资料集导出500hpa高度场保存成grd文件时，文件大小为0，用计算副高指数时，计算不出来。请问是我的脚本编辑错了吗？还是我的副高计算程序错了？请大神帮忙一下

gs脚本
'reinit'
'sdfopen e:\as\NCEP\h2.nc'
'set fwrite e:\as\NCEP\h2.grd'
'set gxout fwrite'
'set x 1 144'
'set y 1 73'
'set lev 500'
t0=1
while (t0<=467)
'set t ' t0
'd hgt*0.1'
t0=t0+1
endwhile
'disable fwrite'

副高计算程序
program fg_zhishu
implicit none
integer,Parameter :: x=144,y=73,yr=38,mo=12
integer,Parameter :: lon1=37,lon2=73,lat10=41,lat80=69  ! lat10--北纬10度，lat80--北纬80度，LON1--东经90，LON2--东经180
integer i,j,t,m,tem,tt,n,s,kr,kkr
real h(x,y,mo,yr),fgm(x,y,mo,yr),fg(x,y,mo,yr),n588(x,mo,yr),s588(x,mo,yr)
real gm1(mo,yr),gq2(mo,yr),gj3(mo,yr),gx4(mo,yr),gb5(mo,yr)

real a,b,c,d,e,f,nn,ch,du
!读取多个文件
character(len=4) :: year
character(len=1) :: index

open(10,file='e:\as\NCEP\h2.grd',form='binary',status='old')
read(10) ((((h(i,j,m,t),i=1,x),j=1,y),m=1,mo),t=1,yr)
close(10)

!!!!!!!!!!!!!!!!!!!!!!!  判断副高的范围 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do t=1,yr
do m=1,mo
!588线南界
do i=lon1,lon2
s588(i,m,t)=0.0
do j=lat10,lat80
if((h(i,j,m,t)/10.0>=588.0) .and. (h(i,j+1,m,t)/10.0>=588.0))then
s588(i,m,t)=j    !南界位置
exit
end if
end do
end do
!588线北界
do i=lon1,lon2
n588(i,m,t)=0.0
do j=lat80,lat10,-1
if((h(i,j,m,t)/10.0>=588.0) .and. (h(i,j-1,m,t)/10.0>=588.0))then
n588(i,m,t)=j    !北界位置
exit
end if
end do
end do
do i=lon1,lon2
!副高分布
if(s588(i,m,t)/=0)then
do j=s588(i,m,t),lat80
fgm(i,j,m,t)=0.0;fg(i,j,m,t)=0.0
if(h(i,j,m,t)/10.0>=588.0)then
fgm(i,j,m,t)=1
fg(i,j,m,t)=h(i,j,m,t)
end if
end do
end if
end do
end do
end do

!!!!!!!!!!!!!!!!!!!!!!!! 1面积指数 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do t=1,yr
do m=1,mo
gm1(m,t)=0.0
do i=lon1,lon2
if(s588(i,m,t)/=0)then
do j=s588(i,m,t),n588(i,m,t)
if((h(i,j,m,t)/10.0>=588.0))then
gm1(m,t)=gm1(m,t)+1
end if
end do
end if
end do
end do
end do
print*, 'ok gm1'

!!!!!!!!!!!!!!!!!!!!!!!! 2强度指数 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do t=1,yr
do m=1,mo
gq2(m,t)=0.0
do i=lon1,lon2
if(s588(i,m,t)/=0)then
do j=s588(i,m,t),n588(i,m,t)
if((h(i,j,m,t)/10.0>=588.0))then
gq2(m,t)=gq2(m,t)+(h(i,j,m,t)/10.0-587.0)
end if
end do
end if
end do
end do
end do
print*, 'ok gq2'
!!!!!!!!!!!!!!!!!!!!!!!!3 西脊点指数!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do t=1,yr
do m=1,mo
gj3(m,t)=0.0;tt=0
do i=lon1,lon2
do j=lat10,lat80
if(h(i,j,m,t)/10.0>=588 .and. h(i+1,j,m,t)/10.0>=588)then
tt=i
exit
end if
end do
if(tt/=0)then
exit
end if
end do
gj3(m,t)=(tt-1)*2.5
if(gj3(m,t)<=0)then
gj3(m,t)=180
end if
end do
end do
print*, 'ok gj3'
!!!!!!!!!!!!!!!!!!!!!!!! 4 脊线位置指数 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do t=1,yr
do m=1,mo
gx4(m,t)=0.0;tt=0
do i=lon1,lon2
if(n588(i,m,t)/=0 .and. s588(i,m,t)/=0)then
tt=tt+1
gx4(m,t)=gx4(m,t)+(n588(i,m,t)+s588(i,m,t))/2.0
end if
end do
if(gx4(m,t)/=0)then
gx4(m,t)=(gx4(m,t)-1)/float(tt)*2.5-90.0
end if
end do
end do
print*, 'ok gx4'
!!!!!!!!!!!!!!!!!!!!!!!! 5 北界位置指数 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do t=1,yr
do m=1,mo
gb5(m,t)=0.0;tt=0
do i=lon1,lon2
if(n588(i,m,t)/=0.0)then
tt=tt+1
gb5(m,t)=gb5(m,t)+n588(i,m,t)
end if
end do
if(gb5(m,t)/=0)then
gb5(m,t)=(gb5(m,t)-1)/float(tt)*2.5-90.0
end if
end do
end do
print*, 'ok gb5'

!!!!!!!!!!!!!!!!!!!!!!!! 输出结果 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!循环输出文件可根据需要输出数据
do t =1,yr
write( year,'(i4)' ) t+1978             !change 1978为其实年份-1
write(*,*) t+1978,year                   !!change
open (9,file='E:\FGwpsh\'//trim(adjustl(year))//'.1.txt',status='replace')
do m=1,mo
write(9,*) gm1(m,t)
end do
close(9)
open (9,file='E:\FGwpsh\'//trim(adjustl(year))//'.2.txt',status='replace')
do m=1,mo
write(9,*) gq2(m,t)
end do
close(9)
open (9,file='E:\FGwpsh\'//trim(adjustl(year))//'.3.txt',status='replace')
do m=1,mo
write(9,*) gj3(m,t)
end do
close(9)
open (9,file='E:\FGwpsh\'//trim(adjustl(year))//'.4.txt',status='replace')
do m=1,mo
write(9,*) gx4(m,t)
end do
close(9)
open (9,file='E:\FGwpsh\'//trim(adjustl(year))//'.5.txt',status='replace')
do m=1,mo
write(9,*) gb5(m,t)
end do
close(9)

end do

end program

muggle · 发表于 2018-1-21 11:34:58

'set fwrite e:\as\NCEP\h2.grd'
'set gxout fwrite'

调整一下顺序试试？改成
'set gxout fwrite'
'set fwrite e:\as\NCEP\h2.grd'

Masterpiece · 发表于 2018-1-21 11:49:39

个人觉得可能是程序读入数据的时候有问题
而且有一行我有些疑问，因为我用过GFS的数据，里边的位势高度得除98，才能转成位势十米，而楼主是乘以0.1，觉得可能是这个问题，导致后边if判据没生效，几个输出量全是0了

建议先将后边的处理语句注释掉，然后试着输出程序读到的内容，看看人工判断能不能让if执行，之后再查if语句

字俣丞 · 发表于 2018-1-21 12:27:41

Masterpiece 发表于 2018-1-21 11:49
个人觉得可能是程序读入数据的时候有问题
而且有一行我有些疑问，因为我用过GFS的数据，里边的位势高度得 ...

除以10和除以9.8差不多的

字俣丞 · 发表于 2018-1-21 12:30:19

@river @mofangbao 请您帮我看看

tulalang · 发表于 2018-1-21 20:26:41

如果可以的话，可以放点积分，然后上传数据让大家下载后帮忙测试

字俣丞 · 发表于 2018-1-21 22:32:41

tulalang 发表于 2018-1-21 20:26
如果可以的话，可以放点积分，然后上传数据让大家下载后帮忙测试

好的！如果可以的话，帮我看一下

字俣丞 · 发表于 2018-1-21 22:36:05

tulalang 发表于 2018-1-21 20:26
如果可以的话，可以放点积分，然后上传数据让大家下载后帮忙测试

但是再分析资料集太大，没办法上传

		自动登录	找回密码
密码			立即注册

[分享资料] 大神帮帮忙，我用NCEP高度场月平均资料集导出500hpa高度场保存成grd文件时，文件大为0

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