请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 8427|回复: 19

[经验总结] 有关水汽净收支,程序能正常输出,但有小问题和大家探讨

[复制链接]

新浪微博达人勋

发表于 2016-1-4 19:07:07 | 显示全部楼层 |阅读模式

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

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

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


错误提示.png

20150831_00.grd

1.36 MB, 下载次数: 7, 下载积分: 金钱 -5

密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-1-5 11:13:44 | 显示全部楼层
自己顶自己......这帖子不能沉啊
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-1-5 14:51:43 | 显示全部楼层
帮楼主顶,最近也遇到类似问题
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-1-5 16:53:55 | 显示全部楼层
ziqiangbuxi 发表于 2016-1-5 14:51
帮楼主顶,最近也遇到类似问题

其实我大概知道问题出在读取文件的,就是应该设置个rec=???,但对这个不清楚怎么修改
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-1-5 16:57:00 | 显示全部楼层
{:eb303:}高手都去哪儿啦?
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-1-6 11:54:28 | 显示全部楼层
不要沉啊~~~~~~
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2016-2-2 11:45:01 | 显示全部楼层
上次写代码的时候,粘贴错了,这个程序才是对的,希望大家能回答我的疑惑
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-10-12 00:02:38 | 显示全部楼层
请问楼主这个问题解决了吗?
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-10-12 11:22:54 | 显示全部楼层
liutong 发表于 2016-10-12 00:02
请问楼主这个问题解决了吗?

最后是解决了,但现在又想不起来怎么改的了
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-10-12 13:10:35 | 显示全部楼层
nunu18 发表于 2016-10-12 11:22
最后是解决了,但现在又想不起来怎么改的了

好的吗,非常感谢你能回复,我自己在研究下吧,谢谢哈~
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

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

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

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