请各位大侠帮忙看看我的理查森数程序有木有问题啊！！！多谢了！

血影残尸

因为需要作海雾研究，用到虚位温、理查森数、模拟能见度几个物理量，但是发现和文献做的图有较大差异，请教各位高手呵呵
本人是按照某硕士论文的公式做的，图也以他的作参考：
由于不知道如何判断饱和或非饱和，所以直接用非饱和的虚位温公式
能见度计算公式
空气密度计算，由于湿空气公式太难，就直接用干空气的，不过看了密度表                             ，貌似湿空气干空气的密度差不多
理查森数，这个是Z坐标的，我自己换算成P坐标的如下：

下面是程序：
!!!!!!计算虚位温qe及比湿q
       DO 920 IT=1,NT
       DO 920 K=1,L
       DO 920 I=1,N
       DO 920 J=1,M
       es(I,J,K,IT)=6.112*EXP(17.67*(T(I,J,K,IT)-273.16)/
     &(T(I,J,K,IT)-29.65))!es为水汽压
       q(I,J,K,IT)=rh(I,J,K,IT)*(0.62197*es(I,J,K,IT)/
     &(iP(K)-es(I,J,K,IT)))/100.0 !q为比湿
       TP(i,j,k,it)=T(I,J,K,IT)*((1000.0/iP(K))**0.286)    !tp是位温
         qev(i,j,k,it)=TP(i,j,k,it)*(1+0.61*q(I,J,K,IT)) !qev是虚位温
920    CONTINUE
         print *,qev(122,120,4,nt)-273.16

!!!!!!计算空气密度
       DO 120 IT=1,NT
       DO 120 K=1,L
       DO 120 I=1,N
       DO 120 J=1,M
         ro(i,j,k,it)=ro0*0.27316*IP(k)/1.013/T(I,J,K,IT)
120    CONTINUE
         print *,ro(122,120,4,nt)

!!!!!!计算平均虚位温
       DO 320 IT=1,NT
       DO 320 I=1,N
       DO 320 J=1,M
         qv(i,j,it)=0.0
       DO K=1,L
         qv(i,j,it)=qv(i,j,it)+qev(i,j,k,it)
         enddo
         qv(i,j,it)=qv(i,j,it)/(L*1.0)
320    CONTINUE
         print *,qv(122,120,nt)-273.16
!!!!!!计算理查森数Ri                  
       DO 168 IT=1,NT
       DO 168 K=1,L
ccccc计算dp
       IF(K.EQ.1) THEN
         KK1=K+1
         KK0=K
         DP=IP(KK1)-IP(KK0)       ! 前差
       ELSE IF(K.EQ.L) THEN
         KK1=K
         KK0=K-1
         DP=IP(KK1)-IP(KK0)       ! 后差
       ELSE
         KK1=K+1
         KK0=K-1
         DP=(IP(KK1)-IP(KK0))*2   ! 中央差
       ENDIF
       DO 168 I=1,n
       DO 168 J=1,m
cccccc计算理查森数
         Ri(I,J,K,IT)=(-1.0)*((qev(I,J,KK1,IT)-qev(I,J,KK0,IT))/DP)/
     &   (((V(I,J,KK1,IT)-V(I,J,KK0,IT))/DP)**2+
     &   ((U(I,J,KK1,IT)-U(I,J,KK0,IT))/DP)**2)/
     &   (ro(i,j,k,it)*qv(i,j,it))
168    CONTINUE
         print *,ri(122,120,4,nt)

!!!!!!计算水平能见度Vi                  
       DO 520 IT=1,NT
       DO 520 K=1,L
       DO 520 I=1,N
       DO 520 J=1,M
         VI(I,J,K,IT)=0.02/144.7/((ro(i,j,k,it)*q(I,J,K,IT))**0.88)
         VI(I,J,K,IT)=(-1.0)*log(VI(I,J,K,IT))/(log(e))
520    CONTINUE
         print *,vi(122,120,2,nt)

然后结果如下：

求教各位大侠了，万分感激！！！
PS：怕某些同学无法下载（自动回复：请不要使用迅雷等下载工具，使用IE右键另存为就可以）附件，就把主要的东东贴出来，可以下载附件的同学可以直接看源程序和源文献：
计算公式.doc (403.5 KB, 下载次数: 66)

2013-2-7 11:32 上传

点击文件名下载附件 ri_vi.for (3.26 KB, 下载次数: 25)

2013-2-7 11:32 上传

点击文件名下载附件

jjmesky

你的Ri坐标转换是否正确？！

天然卷の小博

我正好需要看虚位温的算法，希望有用

爱的侍者

jjmesky 发表于 2013-2-18 21:13

                               
登录/注册后可看大图

你的Ri坐标转换是否正确？！

坐标为啥转换 ，直接 计算不行？？
用的什么数据，我用梯度风观测数据直接计算 z坐标

爱的侍者

坐标为啥转换 ，直接 计算不行？？
用的什么数据，我用梯度风观测数据直接计算 z坐标

HZK

能不能用这个公式，直接用z坐标计算，高度就是等压面上的位势高度乘以g，用这个公式可以省去密度和位温的计算

jjmesky

HZK 发表于 2013-3-1 21:50

                               
登录/注册后可看大图

能不能用这个公式，直接用z坐标计算，高度就是等压面上的位势高度乘以g，用这个公式可以省去密度和位温的计 ...

这个公式常用在通量塔的梯度观测的数据计算，上面说的某硕士论文中用的是L波段雷达观测数据。

钱磊037

请问你的程序里密度计算为什么会这么算呢？ ro(i,j,k,it)=ro0*0.27316*IP(k)/1.013/T(I,J,K,IT)
为什么不直接用P=rou*R*T计算呢？