- 积分
- 2897
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2012-9-1
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看，没有帐号？立即注册
x
本帖最后由 nunu18 于 2016-2-2 11:43 编辑
在论坛里找了很久，但有关水汽净收支的程序少的可怜，而且即便有，但程序有待验证。我是个菜鸟，不怎么会编程，基本上都是靠着大家的无私帮助走过来的。最近一直想计算某个区域的水汽净收支，就是东西南北四个边界的水汽净收支，虽然朋友给了现成的程序最后也能正常运行，但我有几个困惑，想和大家交流以下。首先是从fnl资料中取了地面海平面气压，温度、相对湿度和风场资料，生成了grd资料。
剩下是用fortran读取数据，并进行相关计算，最后输出。
原fortran程序，ix=46，iy=91，我有一直怀疑经纬度写反了，因为按GRADS的存储方式，先经、纬、层次、时间循环，因此x方向应该是91个格点，y方向46个格点才行，因此所有的计算里面经度在最内层循环，所以我把原先程序中的IX IY均改掉了，循环中把所有的IX放最里面循环。我这种改法对不对？忘大家都思考思考
gs文件
'reinit'
'reset'
'open f:/sqsz/fnl_20150831_00_00.ctl '
'set gxout fwrite'
'set fwrite f:/sqsz/20150831_00.grd'
'set lat 15 60 '
'set lon 50 140'
'd PRESsfc/100'
z=1
while(z<=21 )
'set z ' z
'd rhprs'
z=z+1
endwhile
z=1
while(z<=21 )
'set z ' z
'd tmpprs'
z=z+1
endwhile
z=1
while(z<=21 )
'set z ' z
'd ugrdprs'
z=z+1
endwhile
z=1
while(z<=21 )
'd vgrdprs'
z=z+1
endwhile
'disable fwrite'
;
以下是fortran 程序
由grib 数据计算水汽通量、流量、收支
program main
parameter(i_num=91,j_num=46,lev=21)
real ps(i_num,j_num),rh(i_num,j_num,lev)
real u(i_num,j_num,lev),v(i_num,j_num,lev)
real t(i_num,j_num,lev),p(lev)
real td(i_num,j_num,lev),q(i_num,j_num,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(21)
end type str_lev
type(str_lev) levels(4)
data p/1000, 975, 950, 925, 900, 850, 800, 750, 700, 650,600, 550, 500, 450, 400, 350, 300, 250, 200, 150, 100/
data head/'整层','地面-700','700-500','500-300'/
time=6.*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,1/ ! 整层
data levels(2).p/1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0/ ! 地面-700
data levels(3).p/0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0/ ! 700-500
data levels(4).p/0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1/ ! 500-300
! 读入数据
open(806,file='20150829_18.grd',form='unformatted',access='direct',recl=46*91)
irec=1
!rh(i,j,k)读取相对湿度
do k=1,21
read(806,rec=irec) ((rh(i,j,k),i=1,i_num),j=1,j_num)
irec=irec+1
end do
!t(i,j,k)读取温度
do k=1,21
read(806,rec=irec) ((t(i,j,k),i=1,i_num),j=1,j_num)
irec=irec+1
end do
!u(i,j,k)读取u风
do k=1,21
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
do k=1,21
read(806,rec=irec) ((v(i,j,k),i=1,i_num),j=1,j_num)
! print*, v(i,j,k)
irec=irec+1
! print*, i,j,k,v(i,j,k)
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
t=t-273.15
rh=rh/100
u=u*100
v=v*100
! 计算 比湿
call CAL_RH_TD(T,RH,LEV,I_NUM,J_NUM,TD) !调用子程序计算TD(i,j,K)
call CAL_Q(TD,P,LEV,I_NUM,J_NUM,Q) !计算水汽
q=q*1000
print *, 'ps=',ps(2,3)
print *, 'rh=',rh(2,3,6)
print *, 't=',t(2,3,6)
print *, 'u=',u(2,3,6)
print *, 'v=',v(2,3,6)
print *, 'td=',td(2,3,6)
print *, 'q=',q(2,3,6)
open(811,file='20150829_18_1.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(42.),-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(36.),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_RH_TD(T,RH,LEV,I_NUM,J_NUM,TD)
! THIS SUBROUTINE CALCULATE DEW POINT BY USING TEMPERATURE AND RELATIVE HUMIDITY
! INPUT T ========> TEMPERATURE (CELSIUS DEGREE)
! RH =======> RELATIVE HUMIDITY (E/ES)
! LEV ======> LEVEL
! I_NUM ======> 行数
! J_NUM ======> 列数
! OUTPUT TD =======> DEW POINT (CELSIUS DEGREE)
DIMENSION T(I_NUM,J_NUM,LEV),RH(I_NUM,J_NUM,LEV)
DIMENSION TD(I_NUM,J_NUM,LEV)
INTEGER LEV
DO 1101 K=1,LEV
do 1101 j=1,j_NUM
do 1101 i=1,i_NUM
IF(RH(i,j,K)>1.0) RH(i,j,K)=1.0
IF(RH(i,j,K)<0.05) RH(i,j,K)=0.05
RMA=ALOG(RH(i,j,K))*461.*273.15/4.18/597./1000.+T(i,j,K)/(273.15+T(i,j,K))
TD(i,j,K)=273.15*RMA/(1-RMA)
1101 continue
RETURN
END
SUBROUTINE CAL_Q(TD,P,LEV,I_NUM,J_NUM,Q)
! THIS SUBROUTINE CALCULATE SPECIFIC HUMIDITY !!!!!!!!!!!!!!!!!!!!!!!
! NOTICE THIS SUBROUTINE IS SUPPORTED BY SUBROUTINE CAL_E
! INPUT TD ======> DEW POINT TEMPERATURE (CELSIUS DEGREE)
! P =====> PRESSURE (HPA)
! LEV ====> LEVEL
! OUTPUT Q ====> SPECIFIC HUMIDITY
! TEMP ARRAY
! E =====> VAPOR PRESSURE (HPA)
INTEGER LEV
DIMENSION TD(I_NUM,J_NUM,LEV),P(LEV),Q(I_NUM,J_NUM,LEV)
DIMENSION E(I_NUM,J_NUM,LEV)
CALL CAL_E(TD,LEV,I_NUM,J_NUM,E)
DO 2101 K=1,LEV
do 2101 j=1,j_NUM
do 2101 i=1,i_NUM
IF(E(I,J,K).LT.100.)THEN
Q(I,J,K)=0.622*E(I,J,K)/(P(K)-0.378*E(I,J,K))
ELSE
Q(I,J,K)=99999.9
ENDIF
2101 continue
RETURN
END
SUBROUTINE CAL_E(TD,LEV,I_NUM,J_NUM,E) !!!!!!!!!!!!!!!!!!!!!111
! THIS SUBROUTINE CALCULATE VAPOR PRESURE
! NOTICE THE FORMULA ARE DIFFERENT WHEN DEW POINT TEMPERATURE
! IS ABOVE OR BELOW ZERO DEGREE
! INPUT TD =====> DEW POINT (CELSIUS DEGREE)
! LEV ====> LEVEL
! OUTPUT E ====> VAPOR PRESSURE (HPA)
INTEGER LEV
DIMENSION TD(I_NUM,J_NUM,LEV),E(I_NUM,J_NUM,LEV)
DO 3101 K=1,LEV
do 3101 j=1,j_NUM
do 3101 i=1,i_NUM
IF(TD(I,J,K).GT.-150.AND.TD(I,J,K).LT.60.)THEN
IF(TD(I,J,K).GE.0) E(I,J,K)=6.112*10**(7.5*TD(I,J,K)/(237.3+TD(I,J,K)))
IF(TD(I,J,K).LT.0) E(I,J,K)=6.112*10**(9.5*TD(I,J,K)/(265.5+TD(I,J,K)))
ELSE
E(I,J,K)=99999.9
ENDIF
3101 continue
RETURN
END
subroutine cal_vap_flux(u,q,ps,p,ip,lev,I_NUM,J_NUM,uq)
! 计算水汽通量
! INPUT U =====> 风速( cm/s)
! q =====> 比湿( g/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方向
! dis =====> 边界长度
! 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个经度的长度,单位 cm
real num ,pi
pi=3.1415926
dis= Cos(num * pi / 180) * 6371* 2 * pi / 360 * 1000 *100
end function
|
|
窥视卡
雷达卡
发表于 2016-1-4 19:07:07
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
楼主
{:eb303:}高手都去哪儿啦？