fortran计算地转风出错

木易

各位高手好，我用书上的公式计算地转风，资料是北纬10-90度，经度0-360的500hpa的位势高度场计算地转风，不知道为何地转风计算出来的范围在0-200之间，有的风速直接等于0.有的200度，明显不对，哪位高手帮忙看看，这个程序哪里错了？？点击图片无法上传

!计算地转风，公式就是天元书上的公式，ua=-1/f*偏fai/偏y
!数据范围是纬度10-90度，经度0-360 ，格距2.5度
parameter(m=144,n=33,it1=365,unavl=-9999.0,b=10.0) !北半球格点数
      integer m,n,i,j,it
real*8 a,w,r,f1,f2,b1   !地球自转角速度
!地转风u、v分量，h位势高度场，f地转参数
      real ua(m,n,it1),va(m,n,it1),h(m,n,it1)
       real*8 f,fai
real*8 pi

open(1,file='hgt-500hpa.dat',form='binary')
  do it=1,it1
        read(1)((h(i,j,it),i=1,m),j=1,n)
!  write(*,*)((h(i,j,it),i=1,m),j=1,n)
!pause
       enddo
      close(1)

w=0.0000729 !地球自转角速度
pi=3.1415926
r=6371000.0 !地球半径
open(2,file='1ua-500hpa.dat',form='binary')


do it=1,1
!先将所有值都赋值成缺测，然后边界条件不管，直接从i=2,j=2，开始计算，边界赋值为缺测
    do i=1,m
       do j=1,n
           ua(i,j,it)=unavl
           va(i,j,it)=unavl
       enddo !j
    enddo ! i

     !求格点的fai值，便于求f
    do j=2,n-1   
       fai=3.14*(b+2.5*(j-1))/180.0 !每个格点纬度值所对应的弧度值
      f=2*w*sin(fai) !科氏参数
      
   !扇形的弧长=2πr×角度/360  
      b1=r*cos(fai)*2*3.14*2.5/360.0 !纬度上2.5°角度对应的弧长
      b2=r*2*3.14*2.5/360.0 !经度上2.5°角度对应的弧长
        f1=b1*f
        f2=b2*f
   pause     
   do i=2,m-1
pause
    ua(i,j,it)=-180*(h(i,j,it)-h(i,j-1,it))/f2
    va(i,j,it)=180*(h(i,j,it)-h(i-1,j,it))/f1
    write(*,*)va(i,j,it),j,i,f1,f,fai,b1,b,(b+2.5*(j-1))
   enddo ! i

enddo  ! j
!   write(2)((ua(i,j,it),i=1,m),j=1,n)
!   write(*,*)((ua(i,j,it),i=1,m),j=1,n)
   pause
    enddo  !it
!close(2)

end

柿子柿子柿子

1.那个φ你理解是什么？
2.差分公式精度不够
3.Fortran有COSD()、SIND()这一类的函数计算角度的三角函数

木易

柿子柿子柿子 发表于 2015-10-30 18:49
1.那个φ你理解是什么？
2.差分公式精度不够
3.Fortran有COSD()、SIND()这一类的函数计算角度的三角函数

程序中标注的第一行的fai（没打出来符号），我理解的是位势高度，位势米。地转参数中的fai（f=2*w*sin（fai）没打出来符号用拼音代替）我理解的是那个维度的弧度值。

木易

柿子柿子柿子 发表于 2015-10-30 18:49
1.那个φ你理解是什么？
2.差分公式精度不够
3.Fortran有COSD()、SIND()这一类的函数计算角度的三角函数

差分公式的精度不够是？？？这个是格距2.5的数据，我也是按照书上的公式编写的的呢？我把计算弧度的sin（）和计算度数的sind（）都算了一下，结果差不多，都是很大，另外有人说计算地转风要用位势十米去计算，对么？

柿子柿子柿子

木易 发表于 2015-11-4 09:21
程序中标注的第一行的fai（没打出来符号），我理解的是位势高度，位势米。地转参数中的fai（f=2*w*sin（f ...

如果沒記錯的話，位勢和位勢高度是兩個東西

柿子柿子柿子

木易 发表于 2015-11-4 09:50
差分公式的精度不够是？？？这个是格距2.5的数据，我也是按照书上的公式编写的的呢？我把计算弧度的sin（ ...

不，我是說前插後插公式只有一階精度。
不過也有可能不是這個問題。

木易

柿子柿子柿子 发表于 2015-11-5 16:35
如果沒記錯的話，位勢和位勢高度是兩個東西

他们好像是一个东西呢，高度和位势高度不是一个

柿子柿子柿子

木易 发表于 2015-11-5 17:42
他们好像是一个东西呢，高度和位势高度不是一个

我又找了一下書，位勢Φ和位勢高度H還是有差別的
雖然就只差了個9.8

柿子柿子柿子

我還是沒看懂你的那個科氏參數怎麼求的
推導出來的地轉風算式如下：
ug=-(0.67197*10**5/sinφ)*(ΔH/Δy)
vg=(0.67197*10**5/sinφ)*(ΔH/Δx)

φ為緯度，H為位勢高度

qq469015280

想请问 如果是南纬，就是-30N  这个该怎么算，加个绝对值吗？