计算视热源的Fortran程序

清宵逝水

在用论坛的程序（http://bbs.06climate.com/forum.php?mod=viewthread&tid=18894）计算高原地区大气视热源时，遇到与前辈们同样的问题，即量级不对或者差三倍。在论坛进一步搜索发现了@LiuJiawei在改编ncl版（http://bbs.06climate.com/forum.php?mod=viewthread&tid=40487&extra=&page=1）时提出了一些问题，根据他的意见我对原始程序进行修改后，数值仍然不对。在单步调试原始程序后，发现最大的问题在平流项的计算，原始代码为：

1. ug(i,j,l)=(u(i,j,l)*(t(i+1,j,l)-t(i-1,j,l))/dx +(v(i,j,l)*(t(i,j+1,l)-t(i,j-1,l))))/dy

复制代码

通过编辑器的括号匹配功能可以发现，x，y方向的差分并不是平级的，正确写法应为

1. (u(i,j,l)*(t(i+1,j,l)-t(i-1,j,l))/dx + (v(i,j,l)*(t(i,j+1,l)-t(i,j-1,l))/dy))

复制代码

或

1. u(i,j,l)*(t(i+1,j,l)-t(i-1,j,l))/dx + v(i,j,l)*(t(i,j+1,l)-t(i,j-1,l))/dy

复制代码

附件中给出了修改后的代码，并增加了一些注释。
另外补充几点
0.！！！如果采用不同时间分辨率和不同来源的资料，需要仔细修改循环变量，同时一定要调整各个变量的单位！！！
1. 我用的是ncep2逐日资料进行计算，最后结果相对文献（ec资料）的结果略小。也许还有一些小问题，希望使用者如果有空仔细读一遍代码
2.编译器问题，源程序应该是在ivf下编译运行的，文件存取格式用的是binary，gfortran不支持这种格式，我统一改为了unformatted格式，但是需要注意recl参数（即记录长度）的问题，ifort的单位是4字节，只需recl=格点数即可，gfortram的recl单位为字节，需要在格点数的基础上*4（近期了解到stream格式可以解决不同编译器的问题，但是需要对程序结构进行一定的改变，故没有采用）

cyc95611

写了一个转grd的gs,这样写可以吗
'reinit'
year=1981
while(year<=2014)
'set gxout fwrite'
'set fwrite D:/ncep/uwnd/uwnd.'year'.grd'
'sdfopen D:/ncep/uwnd/uwnd.'year'.nc'
'set x 1 144'
'set y 1 73'
if ((math_mod(year,4)=0&math_mod(year,100)!=0)|math_mod(year,400)=0)
i=1
while(i<=366)
'set t 'i''
j=1
while(j<=17)
'set z 'j''
'd uwnd'
j=j+1
endwhile
i=i+1
endwhile
else
i=1
while(i<=365)
'set t 'i''
j=1
while(j<=17)
'set z 'j''
'd uwnd'
j=j+1
endwhile
i=i+1
endwhile
endif
'disable fwrite'
'close 1'
year=year+1
endwhile
;

cyc95611

t和z循环的顺序需要改么

cyc95611

因为您说注意nc转化为grd时的存储顺序,每一笔记录为一个时次的一个层，就是先写t的循环再写z的循环的意思对吧

清宵逝水

cyc95611 发表于 2020-9-19 21:16
因为您说注意nc转化为grd时的存储顺序,每一笔记录为一个时次的一个层，就是先写t的循环再写z的循环的意思对 ...

是的，t和z的循环顺序是对的，t在最外层

要自律的朋友

能问下楼主最后算出来的整层视热源的量级大概多大W/m^2吗

清宵逝水

要自律的朋友 发表于 2020-12-1 22:12
能问下楼主最后算出来的整层视热源的量级大概多大W/m^2吗

我算出来的南坡最大不到500，也就10^2的量级，量级上和文献一致，数值有差异

要自律的朋友

清宵逝水 发表于 2020-12-2 19:57
我算出来的南坡最大不到500，也就10^2的量级，量级上和文献一致，数值有差异

好的 谢谢回复

ANNA0217

请问，如果求500-100hPa ,是不是将   do l=ip,nl-6   改成do l=5,nl-6   即可

清宵逝水

ANNA0217 发表于 2022-3-19 20:14
请问，如果求500-100hPa ,是不是将   do l=ip,nl-6   改成do l=5,nl-6   即可

我看了一下程序，感觉从 do l=6,nl-6 就行了吧，如果6代表500的话，你也可以都试试