- 积分
- 4041
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2016-11-9
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
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)
|
-
-
这是水汽通量散度
|