爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 3032|回复: 4

这些直线是啥,求大神们指教!

[复制链接]

新浪微博达人勋

发表于 2016-11-27 19:48:39 | 显示全部楼层 |阅读模式

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

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 GUTIANWEI 于 2016-11-27 20:20 编辑

【fortran代码】
module global
implicit none
integer,parameter::nx=53,ny=28,nz=8,nt=480
real,parameter::a=6371000,pai=3.1416
real,parameter::dfai=2.5*pai/180,dlenta=2.5*pai/180,g=9.8
integer i,j,iz,it
end module
!定义全局变量
program shixi4
use global
implicit none
real u(nx,ny,nz,nt),v(nx,ny,nz,nt),q(nx,ny,nz,nt)
real fai(nx,ny,nz,nt)
real vapor_flux_x(nx,ny,nz,nt),vapor_flux_y(nx,ny,nz,nt)&
,vapor_flux_D(nx-2,ny-2,nz,nt)
!定义变量
open(1,file='D:\tianzhenshixi\shixi4\data\u.dat'&
,form='unformatted',access='direct',recl=nx*ny*4)
open(2,file='D:\tianzhenshixi\shixi4\data\v.dat'&
,form='unformatted',access='direct',recl=nx*ny*4)
open(3,file='D:\tianzhenshixi\shixi4\data\humidity.dat'&
,form='unformatted',access='direct',recl=nx*ny*4)
do it=1,nt
do iz=1,nz
read(1,rec=iz+(it-1)*17)((u(i,j,iz,it),i=1,nx),j=1,ny)
read(2,rec=iz+(it-1)*17)((v(i,j,iz,it),i=1,nx),j=1,ny)
read(3,rec=iz+(it-1)*nz)((q(i,j,iz,it),i=1,nx),j=1,ny)
end do
end do
close(1)
close(2)
close(3)
!读取原始数据
do it=1,nt
do iz=1,nz
do j=1,ny
do i=1,nx
fai(i,j,iz,it)=(10.0+2.5*j)*pai/180
end do
end do
end do
end do
!计算中间参数
do it=1,nt
do iz=1,nz
do j=1,ny
do i=1,nx
call jisuan_vapor_flux_x(u(i,j,iz,it),q(i,j,iz,it)&
,vapor_flux_x(i,j,iz,it))
call jisuan_vapor_flux_y(v(i,j,iz,it),q(i,j,iz,it)&
,vapor_flux_y(i,j,iz,it))
end do
end do
end do
end do
do it=1,nt
do iz=1,nz
do j=2,ny-1
do i=2,nx-1
call jisuan_vapor_flux_D(fai(i,j,iz,it),vapor_flux_x(i+1,j,iz,it)&
,vapor_flux_x(i-1,j,iz,it),vapor_flux_y(i,j+1,iz,it),&
vapor_flux_y(i,j-1,iz,it),vapor_flux_D(i-1,j-1,iz,it))
end do
end do
end do
end do
!调用外部函数
open(4,file='D:\tianzhenshixi\shixi4\data\vapor_flux_x.dat'&
,form='unformatted',access='direct',recl=nx*ny*4)
open(5,file='D:\tianzhenshixi\shixi4\data\vapor_flux_y.dat'&
,form='unformatted',access='direct',recl=nx*ny*4)
open(6,file='D:\tianzhenshixi\shixi4\data\vapor_flux_D.dat'&
,form='unformatted',access='direct',recl=nx*ny*4)
do it=1,nt
do iz=1,nz
write(4,rec=iz+(it-1)*nz)((vapor_flux_x(i,j,iz,it),i=1,nx),j=1,ny)
write(5,rec=iz+(it-1)*nz)((vapor_flux_y(i,j,iz,it),i=1,nx),j=1,ny)
write(6,rec=iz+(it-1)*nz)((vapor_flux_D(i,j,iz,it),i=1,nx-2),j=1,ny-2)
end do
end do
!新建三个存放计算结果的文件
close(4)
close(5)
close(6)
!关闭文件
end program
!关闭主程序
subroutine jisuan_vapor_flux_x(bu,bq,bvapor_flux_x)
use global
implicit none
real bu,bq,bvapor_flux_x
bvapor_flux_x=1.0/g*bu*bq
end subroutine
!计算x方向水汽通量
subroutine jisuan_vapor_flux_y(bv,bq,bvapor_flux_y)
use global
implicit none
real bv,bq,bvapor_flux_y
bvapor_flux_y=1.0/g*bv*bq
end subroutine
!计算y方向水汽通量
subroutine jisuan_vapor_flux_D(bfai,vapor_flux_x1,vapor_flux_x2&
,vapor_flux_y1,vapor_flux_y2,bvapor_flux_D)
use global
implicit none
real bfai,vapor_flux_x1,vapor_flux_x2,vapor_flux_y1&
,vapor_flux_y2,bvapor_flux_D
bvapor_flux_D=1.0/(2*a)*((vapor_flux_x1-vapor_flux_x2)/&
(dlenta*cos(bfai))+(vapor_flux_y1-vapor_flux_y2)/dfai)
end subroutine
!计算水汽通量散度
【ctl】
dset   D:\tianzhenshixi\shixi4\data\vapor_flux_D.dat
undef  1.e+36
xdef   51 linear 32.5   2.5
ydef   26 linear 15 2.5
zdef   8 levels 1000 925 850 700 600 500 400 300
tdef   480 linear 00:00Z1jan2009 6hr
vars   1
vapor_flux_D   8   99   vapor_flux_D
endvars
【gs】
'reinit'
'open D:\tianzhenshixi\shixi4\data\vapor_flux_D.ctl'
'set lev 850'
'set t 73'
'set grads off'
'd vapor_flux_D*1000000000'
'draw title vapor flux divergence of 850hPa for 2009:1:19:0'
'print D:\tianzhenshixi\shixi4\data\vapor_flux_D.eps'
;
【u,v为17个高度层次,q为8个高度层次 humidity.dat (21.74 MB, 下载次数: 0)
line.png

这是水汽通量散度

这是水汽通量散度
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-11-27 20:34:09 | 显示全部楼层
不是网格线吗?set grid off下?
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-11-27 21:11:22 | 显示全部楼层
不是哎,还是直的。。
no grid.png
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-11-27 21:18:43 | 显示全部楼层
50N度左右的数据应该是有问题。
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-11-27 22:04:03 | 显示全部楼层
Mc.Fish 发表于 2016-11-27 21:18
50N度左右的数据应该是有问题。

好吧,谢谢
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表