求助 等熵位涡时间循环

★翼★

做等熵位涡分析，出现了一些问题，就高手指点一下
第一个问题：'reinit''open e:\toraji\fnl.ctl'
'set lon 0 360'
'set lat -90 90'
'set lev 1000 100'
"define t=T"
"define u=u"
"define v=v"
"define pp=lev"
"define coriol=2*7.29e-5*sin(lat*3.1415/180)"
"define dudy=cdiff(u,y)/(111177*cdiff(lat,y))"
"define dvdx=cdiff(v,x)/(111177*cdiff(lon,x)*cos(lat*3.1415/180))"
"define dt=t(z-1)*pow(1000/PP(z-1),0.286)-t(z+1)*pow(1000/PP(z+1),0.286)"
"define dp=100*(PP(z-1)-PP(z+1))"
"define dtdp=dt/dp"
"define part1="isen(dvdx,t,PP,315)
"define part2="isen(dudy,t,PP,315)
"define part3="isen(dtdp,t,PP,315)
"define pv315=-9.8*(coriol+part1-part2)*part3"
"set z 1"
*'set ccolor 1'
'set lon 80 150'
'set lat 10 60'
'd pv315*1000000'
'gxprint e:\toraji\2.png x1200 y800 white'
*'disable fwrite'

function isen(field,tgrid,pgrid,tlev)



我用这个程序正常出图，用了fwrite之后，就不能出了，ctl如下
dset e:\toraji\is315.dat
title ispv
undef -9.99e33
xdef 71 linear 80 1
ydef 51 linear 10 1
zdef 1 levels 1 1
tdef 1 linear 00z31JuL2001 6hr
vars 1
ipv  0  99 Geopotential vor
endvars

第二个问题，我做了一个时间循环，就是在315-K这个层上的多个时次的PV，但是出来的总是第一个时间的数据，好像while t不起作用，gs如下
*等熵位涡分析
'reinit'
'open e:\toraji\fnl.ctl'
'set lon 0 360'
'set lat -90 90'
*'set lev 1000 100'
*'set t 1 12'
*'set fwrite e:\toraji\is315.dat'
*'set gxout fwrite'
p=1
while(p<13)
'set t 'p''
zz=2
while(zz<3)
*'set z 2 22'
'set lev 1000 100'
"define t=T"
"define u=u"
"define v=v"
"define pp=lev"
"define coriol=2*7.29e-5*sin(lat*3.1415/180)"
"define dudy=cdiff(u,y)/(111177*cdiff(lat,y))"
"define dvdx=cdiff(v,x)/(111177*cdiff(lon,x)*cos(lat*3.1415/180))"
"define dt=t(z-1)*pow(1000/PP(z-1),0.286)-t(z+1)*pow(1000/PP(z+1),0.286)"
"define dp=100*(PP(z-1)-PP(z+1))"
"define dtdp=dt/dp"
"define part1="isen(dvdx,t,PP,315)
"define part2="isen(dudy,t,PP,315)
"define part3="isen(dtdp,t,PP,315)
"define pv330=-9.8*(coriol+part1-part2)*part3"
zz=zz+1
endwhile
"set z 1"
'set lon 80 140'
'set lat 20 55'
'd pv330*1000000'
p=p+1
endwhile
*'gxprint e:\toraji\2.png x1200 y800 white'
*'gxprint e:\toraji\ispv.eps white'
*'disable fwrite'
function isen(field,tgrid,pgrid,tlev)

river

"set z 1"
'set lon 80 140'
'set lat 20 55'
'd pv330*1000000'

你这是怎么个意思，一个时次、一个高度出一张图？直接放在最内层循环出图啊。还有你出图的命令整个在循环外面，那可不就是只出一张图么，也是需要放到最内层循环里面的

★翼★

river 发表于 2015-5-5 07:25
"set z 1"
'set lon 80 140'
'set lat 20 55'

我想做的是一个高度（315-K上的pv），多个时次的，写到一个dat里面。还有，你说的放到内循环，我也试了，出来的都是第一个时次的，如下：
"set z 1"
'set lon 80 140'
'set lat 20 55'
'd pv330*1000000'
zz=zz+1
endwhile
p=p+1
endwhile
里面的'set z 1‘我只是定义到一个维数上了，不然提示错误

river

★翼★ 发表于 2015-5-5 08:49
我想做的是一个高度（315-K上的pv），多个时次的，写到一个dat里面。还有，你说的放到内循环，我也试了， ...

那你的高度不需要循环啊，直接设置Z=2就好了。你说的就出一个时次是什么意思，只出一张图，还是说你出来的dat文件只有一个时次

★翼★

river 发表于 2015-5-5 08:53
那你的高度不需要循环啊，直接设置Z=2就好了。你说的就出一个时次是什么意思，只出一张图，还是说你出来 ...

目的：做一个时间循环（21t），把315k上的位涡写到一个dat里面。做不做高度循环都得到的21个时次，但是画图的时候，不管设定那个时间，出来的图都是第一个时刻的，没有变化。

river

★翼★ 发表于 2015-5-6 17:55
目的：做一个时间循环（21t），把315k上的位涡写到一个dat里面。做不做高度循环都得到的21个时次，但是画 ...

'reinit'
'open e:\toraji\fnl.ctl'
'set lon 0 360'
'set lat -90 90'
*'set lev 1000 100'
*'set t 1 12'
*'set fwrite e:\toraji\is315.dat'
*'set gxout fwrite'
p=1
while(p<13)
'set t 'p''
zz=2
while(zz<3)
*'set z 2 22'
'set lev 1000 100'
"define t=T"
"define u=u"
"define v=v"
"define pp=lev"
"define coriol=2*7.29e-5*sin(lat*3.1415/180)"
"define dudy=cdiff(u,y)/(111177*cdiff(lat,y))"
"define dvdx=cdiff(v,x)/(111177*cdiff(lon,x)*cos(lat*3.1415/180))"
"define dt=t(z-1)*pow(1000/PP(z-1),0.286)-t(z+1)*pow(1000/PP(z+1),0.286)"
"define dp=100*(PP(z-1)-PP(z+1))"
"define dtdp=dt/dp"
"define part1="isen(dvdx,t,PP,315)
"define part2="isen(dudy,t,PP,315)
"define part3="isen(dtdp,t,PP,315)
"define pv330=-9.8*(coriol+part1-part2)*part3"
'set lon 80 140'
'set lat 20 55'
'd pv330*1000000'
zz=zz+1
endwhile
p=p+1
endwhile
function isen(field,tgrid,pgrid,tlev)
;
试一下

★翼★

river 发表于 2015-5-7 19:53
'reinit'
'open e:\toraji\fnl.ctl'
'set lon 0 360'

出现问题：



如果加上’set z 1'就会出现21个时刻同样的图

river

★翼★ 发表于 2015-5-8 13:18
出现问题：

'reinit'
'open e:\toraji\fnl.ctl'
'set fwrite e:\toraji\is315.dat'
'set gxout fwrite'
p=1
while(p<13)
'set t 'p''
'set z 2'
"define t=T"
"define u=u"
"define v=v"
"define pp=lev"
"define coriol=2*7.29e-5*sin(lat*3.1415/180)"
"define dudy=cdiff(u,y)/(111177*cdiff(lat,y))"
"define dvdx=cdiff(v,x)/(111177*cdiff(lon,x)*cos(lat*3.1415/180))"
"define dt=t(z-1)*pow(1000/PP(z-1),0.286)-t(z+1)*pow(1000/PP(z+1),0.286)"
"define dp=100*(PP(z-1)-PP(z+1))"
"define dtdp=dt/dp"
"define part1="isen(dvdx,t,PP,315)
"define part2="isen(dudy,t,PP,315)
"define part3="isen(dtdp,t,PP,315)
"define pv330=-9.8*(coriol+part1-part2)*part3"
'set lon 80 140'
'set lat 20 55'
'd pv330*1000000'
p=p+1
endwhile
function isen(field,tgrid,pgrid,tlev)
;

★翼★

river 发表于 2015-5-8 15:15
'reinit'
'open e:\toraji\fnl.ctl'
'set fwrite e:\toraji\is315.dat'

因为是从p坐标系差值到等熵坐标系，需要不同P层上的值，而先头'set z 2'就限定了只是利用了一层p上的值，可能造成了这这种结果。如果把'set z '放在后面吧，又出现同样的值

Adam-Lee

求等熵位涡的时候是要把t/u/v都放到同一个dat里面去吗？