水汽收支的相关计算问题

Zhangzh · 发表于 2018-5-8 22:47:28

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

您需要登录才可以下载或查看，没有帐号？立即注册

x

最近刚开始接触大气方面的知识，想计算某地区的水汽净收支，从论坛里找“

” 的程序(输入的是grd文件)修改了一下，我的输入原文件是nc文件，我先使用grads将nc转为dat格式，然后使用fortran读取dat进行计算，目前遇到的问题如下：
1. nc转dat可不可以一次转化多个变量？
2. 我的转化成dat为啥打不开，说是"描述文件错误"

grads转化nc为dat格式代码：

fortran 程序打开dat文件及运行的报错信息：

grads打开转化的dat文件的报错信息：

请大家帮我看看，谢谢大家

Zhangzh · 发表于 2018-5-9 09:23:06

自己顶一下，别沉啊

xin729914245 · 发表于 2018-5-9 10:23:19

我印象中是可以一次性转化多个变量的……你的gs我觉得没有什么问题啊

Zhangzh · 发表于 2018-5-9 10:47:01

xin729914245 发表于 2018-5-9 10:23
我印象中是可以一次性转化多个变量的……你的gs我觉得没有什么问题啊

您好，谢谢您的回答，多个变量我后面再试试，现在是nc转化成dat后，我的在grads中不能打开，程序里调用也不对，感觉应该是“open(806,file='2007_1.dat',form='unformatted',access='direct',recl=41*65)”这个“recl”设置的有问题？附件是我用的Fortran，您能不能帮忙看看，谢谢

program main
parameter(i_num=41,j_num=65,lev=20)
real ps(i_num,j_num)
real u(i_num,j_num,lev),v(i_num,j_num,lev)
real q(i_num,j_num,lev),p(lev)
real uq(i_num,j_num),vq(i_num,j_num)
real vapour(10),vap_in,vap_out,time
real vap_inp,vap_outp
integer irec,b(10)
character head(4)*10
type str_lev
integer p(20)
end type str_lev
type (str_lev) levels(4)
data p/1000, 975, 950, 925, 900, 875, 850, 800, 825, 775, 750, 700, 650, 600, 550, 500, 450, 400, 350, 300/
data head/'1000-300','1000-800','800-500','500-300'/
time=30*24*60*60/1000/1000
data levels(1)%p/1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1/ ! 1000-300
data levels(2)%p/1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0/ ! 1000-800
data levels(3)%p/0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0/ ! 800-500
data levels(4)%p/0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1/ ! 500-300
! 读入数据
open(806,file='2007_1.dat',form='unformatted',access='direct',recl=41*65)
irec=1
!u(i,j,k)读取u风
do k=1,20
read(806,rec=irec) ((u(i,j,k),i=1,i_num),j=1,j_num)
irec=irec+1
! print*,u(i,j,k)
end do
!v(i,j,k)读取v风
do k=1,20
read(806,rec=irec) ((v(i,j,k),i=1,i_num),j=1,j_num)
! print*, v(i,j,k)
irec=irec+1
end do
!q(i,j,k)读取比湿(相对湿度)
do k=1,20
read(806,rec=irec) ((q(i,j,k),i=1,i_num),j=1,j_num)
! print*, q(i,j,k)
irec=irec+1
end do
!读取地面气压
do k=1,1
read(806,rec=irec) ((ps(i,j),i=1,i_num),j=1,j_num)
irec=irec+1
end do
open(811,file='result.txt')
write(811,'(10x,7a15)') 'west','north','east','south','all_in','all_out','all'
! 计算水汽通量
do i=1,4
vapour=0
vap_inp=0
vap_outp=0
uq=0
vq=0
call cal_vap_flux(u,q,ps,p,levels(i)%p,lev,I_NUM,J_NUM,uq)
call cal_vap_flux(v,q,ps,p,levels(i)%p,lev,I_NUM,J_NUM,vq)
! 计算水汽收支
call cal_vapour(uq,vq,i_num,j_num,22,28,23,1,dis(0.),1,vapour(1),vap_inp,vap_outp)
print *,'1 ' ,vapour(1),vap_inp,vap_outp
call cal_vapour(uq,vq,i_num,j_num,23,31,28,2,dis(43.),-1,vapour(2),vap_inp,vap_outp)
print *,'2 ' ,vapour(2),vap_inp,vap_outp
call cal_vapour(uq,vq,i_num,j_num,22,28,31,1,dis(0.),-1,vapour(3),vap_inp,vap_outp)
print *,'3 ' ,vapour(3),vap_inp,vap_outp
call cal_vapour(uq,vq,i_num,j_num,23,31,22,2,dis(38.),1,vapour(4),vap_inp,vap_outp)
print *,'4 ' ,vapour(4),vap_inp,vap_outp
write(811,'(a10,7e15.4)') head(i),(vapour(ii),ii=1,4),vap_inp,vap_outp,vap_inp-vap_outp
end do
close(811)
end
subroutine cal_vap_flux(u,q,ps,p,ip,lev,I_NUM,J_NUM,uq)
! 计算水汽通量
! INPUT U =====> 风速( m/s)
! q =====> 比湿( kg/kg)
! ps =====> 地面气压 (hPa)
! p =====> 气压层(hPa)
! Ip =====> 分层计算
! LEV ====> LEVEL
! I_NUM ====> 行数
! J_NUM ====> 列数
! OUTPUT uq ====> 水汽通量
INTEGER LEV,lev_tem,I_NUM,J_NUM
DIMENSION U(I_NUM,J_NUM,LEV),Q(I_NUM,J_NUM,LEV)
dimension PS(I_NUM,J_NUM),P(LEV),Ip(lev)
dimension uq(I_NUM,J_NUM),uq_tem(lev)
uq=0
do 4101 j=1,j_NUM
do 4101 i=1,i_NUM
lev_tem=0
do 4106 k=1,lev
uq_tem(k)=u(i,j,k)*q(i,j,k)
if(ps(i,j).le.p(k)) lev_tem=k
4106 continue
uq(i,j)=uq_tem(lev_tem+1)*(ps(i,j)-p(lev_tem+1))*IP(lev_tem+1)
if(i.eq.2.and.j.eq.3) then
write(*,'(i4,6f10.2)') lev_tem,u(i,j,lev_tem+1),q(i,j,lev_tem+1),uq_tem(lev_tem+1),ps(i,j),p(lev_tem+1),uq(i,j)
endif
do 4111 k= lev_tem+1,lev-1
uq(i,j)=uq(i,j)+(uq_tem(k)+uq_tem(k+1))/2*(p(k)-p(k+1))*IP(k+1)
if(i.eq.2.and.j.eq.3) then
write(*,'(i4,6f10.2)') k,u(i,j,k),q(i,j,k),uq_tem(k),p(k),p(k+1),uq(i,j)
endif
4111 continue
4101 continue
uq=uq/980
print *, uq(2,3)
return
end
subroutine cal_vapour(uq,vq,i_num,j_num,i,j,k,index,dist,bor,vap,vap_in,vap_out)
! 计算水汽收支
! INPUT uq =====> u方向水汽通量
! vq =====> v方向水汽通量
! I_NUM ====> uq,vq的行数
! J_NUM ====> uq,vq的列数
! i =====> 计算时边界的起始点(行或列)
! j =====> 计算时边界的结束点(行或列)
! K =====> 计算时边界行列值中固定不变的点
! index =====> 标示是计算uq方向或vq方向
! dist =====> 边界长度
! bor =====> 边界属性标示
! OUTPUT vap====> 水汽净通过量
! vap_in ====> 流入水汽
! vap_out ====> 流出水汽
integer i_num,j_num,i,j,k,index,bor
real uq(i_num,j_num),vq(i_num,j_num),dist
real vap, vap_in,vap_out
vap=0
time=6.*60*60/1000/1000
do ii=i,j
if(index==2) then
vap=vap+(vq(ii,k)+vq(ii+1,k))*0.5*dist*time
if((vq(ii,k)+vq(ii+1,k))*bor>=0) then
vap_in=vap_in+(vq(ii,k)+vq(ii+1,k))*bor*0.5*dist*time
else
vap_out=vap_out-(vq(ii,k)+vq(ii+1,k))*bor*0.5*dist*time
endif
elseif(index==1) then
vap=vap+(uq(k,ii)+uq(k,ii+1))*0.5*dist*time
if((uq(k,ii)+uq(k,ii+1))*bor>=0) then
vap_in=vap_in+(uq(k,ii)+uq(k,ii+1))*bor*0.5*dist*time
else
vap_out=vap_out-(uq(k,ii)+uq(k,ii+1))*bor*0.5*dist*time
endif
endif
enddo
end
real function dis(num)
!num 度纬圈上，1个经度的长度,单位 m
real num ,pi
pi=3.1415926
dis= Cos(num * pi / 180) * 6371* 2 * pi / 360 * 1000
end function

复制代码

nunu18 · 发表于 2018-5-11 13:16:38

楼主你经度和纬度方向取的格点数并不是41*65，程序不出错才怪

Zhangzh · 发表于 2018-5-11 18:49:49

nunu18 发表于 2018-5-11 13:16
楼主你经度和纬度方向取的格点数并不是41*65，程序不出错才怪

明白了，之前我用的是原数据的格点数，没用转换后dat的格点数，谢谢

曦懿-叶子 · 发表于 2019-12-18 12:44:32

楼主，你有那个程序吗？您方便分享一下吗，我的贡献不够，无法下载（自动回复：请不要使用迅雷等下载工具，点我查看下载帮助），

		自动登录	找回密码
密码			立即注册

[源代码] 水汽收支的相关计算问题

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