用Fortran程序编写水汽通量、水汽通量散度程序

小星星~~

最近，遇到了水汽通量和水汽通量散度计算的问题。以前都是用GRADS计算的，这一次由于业务需要，想尝试一下用FORTRAN编程。

找到一个大学天气诊断实习课的程序，计算了各种物理量场！（程序在附件中，有需要的可下载）别的物理量暂且忽略，这里只讨论水汽通量和水汽通量散度。发现，这个程序计算出的水汽通量散度量级都在1e-3左右，网上查了一下，水汽通量散度量级一般在1e-7。
家园里有个讨论水汽通量量级的帖子，但是用的是GRADS编程，地址：http://bbs.06climate.com/forum.php?mod=viewthread&tid=24090 ，大概意思就是，量级不同的原因是单位没统一。于是我统一了单位，但是量级变成1e-5,差距还是很大....什么原因呢？查看了天原书，终于找到了原因。


以下是计算水汽通量散度的程序段：
DATA p    /1000,925,850,700,500,400,300,250,200,150,100/
!--------------------------计算散度-----------------------------------
div=1.e+36
do it=1,nt
do iz=1,nz
   do ii=1+1,nx-1
   do ij=1+1,ny-1
   div(ii,ij,iz,it)=1.0/(2*R)*( (u(ii+1,ij,iz,it)-u(ii-1,ij,iz,it))/(Delta*PI/180*cos(lat(ij)))  +       &
                                  (v(ii,ij+1,iz,it)-v(ii,ij-1,iz,it))/(Delta*PI/180) -                               &
                  2*v(ii,ij,iz,it)*tan(lat(ij))  )
   enddo
   enddo
enddo
enddo
!-------------------------计算水汽通量和水汽通量散度----------------------------------
adq =1.e+36
adqv=1.e+36
qu=1.e+36
qv=1.e+36


do it=1,nt
do iz=1,nz


do ii=1+1,nx-1
do ij=1+1,ny-1
         quv(ii,ij,iz,it)=q(ii,ij,iz,it)*sqrt(u(ii,ij,iz,it)**2+v(ii,ij,iz,it)**2)*(p(iz)-p(iz+1))/9.8
         adq(ii,ij,iz,it)=u(ii,ij,iz,it)*((q(ii+1,ij,iz,it)-q(ii-1,ij,iz,it))/(2*R*cos(lat(ij))*Delta*PI/180))*(p(iz)-p(iz+1))/9.8 +   &
                        v(ii,ij,iz,it)*((q(ii,ij+1,iz,it)-q(ii,ij-1,iz,it))/(2*R*Delta*PI/180))*(p(iz)-p(iz+1))/9.8         !比湿平流
         adv(ii,ij,iz,it)=q(ii,ij,iz,it)*div(ii,ij,iz,it)*(p(iz)-p(iz+1))/9.8                !风散度项
         adqv(ii,ij,iz,it)=(adv(ii,ij,iz,it)+adq(ii,ij,iz,it))*1e5
enddo
enddo


enddo
enddo

！------------------------------------------------------------------------------
这里面quv、adqv分别代表水汽通量、水汽通量散度
quv没什么说的，就是水汽乘以全风速，再除以g=9.8
adqv的计算就比较有趣了，这个程序的计算方法和天原书上的明显不同：adqv=adq(比湿平流)+adv（风散度）
注意，每一项后面都乘以一个(p(iz)-p(iz+1))，这一项实际上就是△P，也是因为这一项，导致该程序计算量级出现错误。看下P的值，不难发现△P的值是很大的，量级1e2左右，1e-5和1e-7刚好差在这里。


那天原书上是怎么讲的呢？
水汽通量（如下图）：可以看到乘以△P和△l，其单位为g/s。当截面积高取1hPa底边长为1cm时，即△P=1，△l=1，这样我们得到的才是每一层的水汽通量，单位为g/(s.hPa.cm)，因此，确定(p(iz)-p(iz+1))这项要去掉。

水汽通量散度（如下图）：前面说了，程序和天原书的计算方法有差别，天原书是先计算qu、qv然后对这两项分别进行差分处理。而程序里的计算方法不知道根据什么公式来的（拜托有了解的小伙伴出来解说下？）



天原书里有据可依，因此我打算采用天原书里的方法，自己重新编写一个程序。

程序段如下：
!-------------------------计算水汽通量和水汽通量散度----------------------------------
    adqv=9.999E+20
    qu=9.999E+20
    qv=9.999E+20
    quv=9.999E+20

   do i=1+1,nx-1
   do j=1+1,ny-1
     if(qq(i,j)/=9.999E+20.and.uwind(i,j)/=9.999E+20.and.vwind(i,j)/=9.999E+20)then
         quv(i,j)=qq(i,j)*sqrt(uwind(i,j)**2+vwind(i,j)**2)/9.8 !!!!!!!水汽通量 单位g/(s.hPa.cm)
     qu(i,j)=(uwind(i,j)*qq(i,j))/9.8 !!!!!单位：g/(hPa.cm.s)
     qv(i,j)=(vwind(i,j)*qq(i,j))/9.8
     endif
   enddo
   enddo
   do i=1+1,nx-1
   do j=1+1,ny-1
   lat=llat(i,j)*PI/180 !!!转换弧度
   if(qu(i+1,j)/=9.999E+20.and.qu(i-1,j)/=9.999E+20.and.qv(i,j+1)/=9.999E+20.and.qv(i,j-1)/=9.999E+20)then
   adqv(i,j)=1.0/(2*R*100)*( ((qu(i+1,j)-qu(i-1,j))/(Delta*PI/180*cos(lat))) + ((qv(i,j+1)-qv(i,j-1))/(Delta*PI/180)) ) !!!!!水汽通量散度 单位：g/(hPa.cm2.s)
   !adqv(i,j)=adqv(i,j)*1e8
   endif
   enddo
   enddo

！-------------------------------------
水汽通量计算比较简单，就是水汽乘以全风速除以9.8
水汽通量散度是：先计算了qu、qv，然后对二者进行差分处理。
计算时注意数据的单位换算：
q单位g/kg
u,v单位m/s
g=9.8单位m/s2
R=6371e+3为地球半径单位m（计算时还要乘以100换算成cm）
lat为纬度（转换成弧度制）
PI=3.1415926
Delta是格点资料的网格距离（计算时也要换算成弧度Delta*PI/180）
根据上述公式和标准单位，不难推导出adqv的单位是g/(hPa.cm2.s)，有兴趣的小伙伴可以推导下，我已经推导过，没问题。
以上就是我对水汽通量散度编程问题的理解，不知道是否正确？
请各位园友批评指正！！
下图为绘制的水汽通量散度图，我用的是EC细网格资料做的：

水汽通量散度图

nunu18

楼主△P是什么？

小星星~~

nunu18 发表于 2018-5-10 10:34
楼主△P是什么？

两层之间的气压差啊，你可以看天原书，里面有一章讲的水汽通量计算，计算单层的话，这项设为1。

nunu18

明白了，谢谢楼主！楼主好细心

葵花笙

{:5_213:}{:5_213:}{:5_213:}{:5_213:}

albert007zh

谢谢楼主分享，赞

GExynos

真的很谢谢楼主，这学期刚好学了天诊。

水月洞天

谢谢楼主分享，赞

yusw

请问一下，计算散度里的这一项是什么意思， -2*v(ii,ij,iz,it)*tan(lat(ij))

良辰

您好楼主，两种算法您做过检验吗？我算出来为什么差异明显