希望有大神能帮我看看这个边界水汽收支的程序到底是不是对的

yuejun0026

        program main

        parameter(lon=360,lat=181,lev=19,tim=24,l=1.11133,pi=3.1415926)
        integer::i,j,z,t,bz
        real::fuq(lon,lat,lev,tim),fvq(lon,lat,lev,tim),
     $divf(lon,lat,lev,tim),f(lon,lat,lev,tim),allfuq(lon,lat,tim),
     $allfvq(lon,lat,tim),pres(lon,lat,tim),bj(tim,4)
        real::we(tim),sn(tim),total(tim)
        real::qu,qv,df,sqtl,sum
        integer::p(lev)
        p=(/1000,975,950,925,900,850,800,750,700,650,600,550,500,450,
     $400,350,300,250,200/)
c        读入数据
        open(10,FILE="H:\data\fnl\gs\bjsqsz\sqtl.grd",form="binary")
        do t=1,tim
c        print*,t
        do z=1,lev
                read(10) ((fuq(i,j,z,t),i=1,lon),j=1,lat)
        enddo

        do z=1,lev
                read(10) ((fvq(i,j,z,t),i=1,lon),j=1,lat)
        enddo

        do z=1,lev
                read(10) ((divf(i,j,z,t),i=1,lon),j=1,lat)
        enddo

        do z=1,lev
                read(10) ((f(i,j,z,t),i=1,lon),j=1,lat)
        enddo
       
        read(10) ((pres(i,j,t),i=1,lon),j=1,lat)
       
        enddo
       

c     判断垂直积分的起始层次
        allfuq=0
        allfvq=0
        do i=1,lon
c        print*,'i',i
        do j=1,lat
c        print*,'j',j
        do t=1,tim
c        print*,'t',t
        do z=1,lev
c        print*,'z',z
c        print*,p(z)
        if((pres(i,j,t).gt.1000).and.(pres(i,j,t).lt.200))then
                        print*,'pres_wrong'
        endif
        if(pres(i,j,t).ge.(p(z)))then
                        bz=z
                        print*,'goto'
                        goto 40
        endif
        enddo

       
c        计算水汽通量垂直积分
40        do z=bz,lev-1
c        print*,'zz',z
        allfuq(i,j,t)=allfuq(i,j,t)+0.5*((fuq(i,j,z,t)+fuq(i,j,z+1,
     $t))*(p(z)-p(z+1)))

        allfvq(i,j,t)=allfvq(i,j,t)+0.5*((fvq(i,j,z,t)+fvq(i,j,z+1,
     $t))*(p(z)-p(z+1)))
        enddo
        allfuq(i,j,t)=allfuq(i,j,t)*100*0.001
        allfvq(i,j,t)=allfvq(i,j,t)*100*0.001
        enddo
        enddo
        enddo
c     确定东、西、南、北边界
        alat=24+90
        blat=36+90
        alon=97
        blon=109
c        计算四个边界的水汽收支
        open(20,file="H:\data\fnl\gs\bjsqsz\water_incoming.txt")
        write(20,'(a)')'day,W_boundary,E_boundary,S_boundary,N_boundary,
     $WE_incoming,SN_incoming,total_incoming'
        bj=0
        do t=1,tim
        do j=alat,blat
        bj(t,1)=bj(t,1)+allfuq(alon,j,t)
        bj(t,2)=bj(t,2)+allfuq(blon,j,t)
        enddo
        do i=alon,blon
        bj(t,3)=bj(t,3)+allfvq(i,alat,t)
        bj(t,4)=bj(t,4)+allfvq(i,blat,t)
        enddo
c        考虑边界长度
        bj(t,1)=bj(t,1)*l
        bj(t,2)=bj(t,2)*l
        bj(t,3)=bj(t,3)*l*cos((alat-90)/180*pi)
        bj(t,4)=bj(t,4)*l*cos((blat-90)/180*pi)
        enddo
c        计算东、西、南、北以及总净水气输送量
        do t=1,tim
        we(t)=bj(t,1)-bj(t,2)
        sn(t)=bj(t,3)-bj(t,4)
        total(t)=we(t)+sn(t)
        enddo
c        输出结果
        do t=1,tim
        write(20,'(i8,7f12.2)')t,(bj(t,ib),ib=1,4),we(t),sn(t),total(t)
        enddo
       
        end


这是我根据找到的一个程序自己改的，但是不知道出来的结果是不是对的，希望哪位大神能够帮我看看，或者给我一种检验这个结果是否正确的方法。其中，fuq是纬向水汽通量，fvq是经向水汽通量，divf是水汽通量散度（其实程序中并不会用到），f是水汽通量，pres是海平面气压。

万分感谢了 ！

yuejun0026

哎哟  自己顶一下了

yuejun0026

哎哟  怎么没大神回复一下啊？

虎跃平阳

你好，我最近也在算这种边界的水汽收支，请问你算出来的量级是多少？我看到文献上的单位是kg*s**-1。

lgg00828

我在网上就找到了楼主的这一段计算水汽收支的程序，感谢您的分享，希望大神指导一下，程序没问题吧？

yuejun0026

lgg00828 发表于 2014-5-22 16:32
我在网上就找到了楼主的这一段计算水汽收支的程序，感谢您的分享，希望大神指导一下，程序没问题吧？

不知道哦   现在已经放弃这个程序了

lgg00828

yuejun0026 发表于 2014-5-22 20:21
不知道哦   现在已经放弃这个程序了

请问楼主，计算水汽收支的时候，需要乘以边界长度吗？

yuejun0026

lgg00828 发表于 2014-5-26 10:43
请问楼主，计算水汽收支的时候，需要乘以边界长度吗？

需要哈  是要乘以边界长度的

yuejun0026

lgg00828 发表于 2014-5-26 10:43
请问楼主，计算水汽收支的时候，需要乘以边界长度吗？

需要哈  是要乘以边界长度的

lgg00828

yuejun0026 发表于 2014-5-26 11:17
需要哈  是要乘以边界长度的

不知道楼主有相关的程序不