用GRADS编写的波活动通量WAF公式，怎么输出结果都是0啊？

jolincai

看了一篇文献，里面有一个波活动通量公式，计算出来画图可以表示波能量的传播方向，所以应该算出来是一个和风场类似的矢量。根据他的公式我自己用grads编写计算过程，不知道哪里有错，是二阶导数的地方吗？求指导！公式如下：

说说我对公式的理解：U V为基本流场，文献里取的2010年4-6月的平均，我做冬季，我取的某年12-2月平均。
                                ψ是扰动流函数，即距平。我算出了某年冬季逐日流函数，用的fish_psi函数，然后用逐日流函数减去冬季平均得到距平值。
                                计算公式之前，我已经准备好了里面 U V 和ψ扰动流函数的数据。
                                以下是我的gs文件，用来计算这个公式，W分为Wx方向和Wy方向两部分；（psi即流函数(这里是其距平值））
'reinit'
'open E:\data\psijp.ctl'
'open E:\data\upj.ctl'
'open E:\data\vpj.ctl'
'set fwrite E:\waf.grd'
'set gxout fwrite'
'set x 1 144'
'set y 1 73'
i=1
while(i<=90)
'set t 'i''
'define dpx=cdiff(psi,x)'
'define dpy=cdiff(psi,y)'
'define dx=cdiff(lon,x)*3.1416/180'
'define dy=cdiff(lat,y)*3.1416/180'
'define psix=dpx/(cos(lat*3.1416/180)*dx)/6.37e6'*扰动流函数对x求导
'define psiy=dpy/dy/6.37e6'*扰动流函数对y求导
'define dpxx=cdiff(psix,x)'
'define psixx=dpxx/(cos(lat*3.1416/180)*dx)/6.37e6'**扰动流函数对x求二阶导
'define dpxy=cdiff(psix,y)'
'define psixy=dpxy/dy/6.37e6'*扰动流函数对x求导后对y求二阶导（不知道求二阶导出错没，应该就是把一阶导的结果再次求导）
'define wx=((u.2)*(psix*psix-psi*psixx)+(v.3)*(psix*psiy-psi*psixy))/(2*(u.2))'
'define wy=((u.2)*(psix*psiy-psi*psixy)+(v.3)*(psiy*psiy-psi*psixy))/(2*(u.2))'
'd wx'
'd wy'
i=i+1
endwhile
'disable fwrite'
求指导啊~~~~~~~~~~~~~~~~~~~~~~~~~~

jolincai

@river 求大神...

jolincai

更正：W公式左边的因子U应该是风速大小。

river

你这个真的好难，我都看不懂这个公式了······
首先你准备好的几个量，是你自己计算的吗？可以正常出图吗？
然后就是公式里的分母那个风速是取绝对值以后的，你计算过程中没有取绝对值
最后，你不要先输出成二进制资料，你直接让它在屏幕上显示，看看出来的图对不对

jolincai

river 发表于 2017-9-28 15:41
你这个真的好难，我都看不懂这个公式了······
首先你准备好的几个量，是你自己计算的吗？可以正常出 ...

你好，，我准备的平均UV场和流函数场是可以画图的，然后计算过程，其实是把流函数的每一个偏导数求出来然后再相乘，二阶导我是连续计算两次一阶导数得到的...文献中公式的意思是，W=左边因子乘上右边括号里上排的因子是W的一个方向矢量，乘上下一排是另一个方向矢量吗？
文献公式，左边的因子U是风速大小，我后来改成了mag(u,v).
我直接画图可以出来，但是写入二进制后就不行。。。不知道为什么
我随意画的某个时次如图，出了局部地区有明显箭头 其他基本都是...啥也看不出来 。。。也不知道出来的东西对不对...

jolincai

......................

jolincai

jolincai 发表于 2017-9-29 10:39
你好，，我准备的平均UV场和流函数场是可以画图的，然后计算过程，其实是把流函数的每一个偏导数求出来然 ...

图在楼上.......

river

jolincai 发表于 2017-9-29 10:39
你好，，我准备的平均UV场和流函数场是可以画图的，然后计算过程，其实是把流函数的每一个偏导数求出来然 ...

我能说我也不知道对不对吗，你这个太难了。但是我觉得你这个求导的时候除的地球半径太多了，应该在最后一个式子里就是wx和wy那里总的除一个就行了吧。你说数据都是0，是怎么看的数据？

jolincai

river 发表于 2017-9-29 16:48
我能说我也不知道对不对吗，你这个太难了。但是我觉得你这个求导的时候除的地球半径太多了，应该在最后一 ...

我写成二进制后用FORTRAN打开，结果数据都是0，但是如果直接画图的话 好像又有图形...想用这个公式画图研究波能量的传播方向...

river

jolincai 发表于 2017-9-29 17:08
我写成二进制后用FORTRAN打开，结果数据都是0，但是如果直接画图的话 好像又有图形...想用这个公式画图研 ...

是不是读取过程有问题？看图感觉还是有数据的