用fortran计算地转风涡度出错

surelee

用fortran计算地转风涡度，后用grads画图出现了问题，我觉得我写的fortran程序没有问题，附上程序请会的朋友指点一下不胜感激。

1. 
   

复制代码

说明：用中央差分计算的地转风涡度参考公式是
顺便附上从.nc文件读取的位势高度数据的.gs和画地转风涡度的ctl文件

1. 'sdfopen d:\tex2\Apr2015_ght.nc'
   
2. 'set fwrite d:\tex2\Apr2015_ght.grd'
   
3. 'set gxout fwrite'
   
4. 'set x 1 480'
   
5. 'set y 1 241'
   
6. 'set time 12z28Apr2015'
   
7. 'set lev 500'
   
8. 'd z'
   
9. 'disable fwrite'

复制代码

surelee

原谅我以这种方式附上代码：
program ex2
implicit none
   integer,parameter::nx=480,ny=241,nz=1,nt=1
   integer,parameter:: tx=nx-2,ty=ny-2,tz=nz,tt=nt
   real ght(nx,ny,nz,nt)     !位势高度  
   real ksei(2:tx,2:ty,tz,tt)   !地转风涡度   
   integer ix,iy
   integer i,j
!定义常数
   real::r=6371000 !地球半径
   real f !地转参数
   real ::omiga=7.292*10**(-5)
   real ::g=9.8 !重力加速度
   real dx,dy    !差分格距
   real ::dxx=0.75*2*3.1416/360,dyy=0.75*2*3.1416/360 !经纬网格对应的格距即精度
   real lat !纬度
!读取位势高度场数据
   open(1,file='d:/tex2/Apr2015_ght.grd',form='binary')  
   do iy=1,ny
      do ix=1,nx
           read(1)ght(ix,iy,nz,nt)
          enddo
   enddo
   close(1)
        print*,'read done'
!计算地转风涡度并写入文件
   open(2,file='d:/tex2/Apr2015_ksei.grd',form='binary')
   do y=2,ty
          lat=(-90.0+dxx*(j-1))*2*3.1416/360  !计算格点地理纬度
      do i=2,tx
               f=2*omiga*sin(lat)   !计算地转参数
                   ksei(i,j,tz,tt)=g/(f*(r**2))*(  (ght(i+1,j,nz,nt)+ght(i-1,j,nz,nt)-2.0*ght(i,j,nz,nt))/(((dxx)**2)*cos(lat)*cos(lat))&
                                                  +(ght(i,j+1,nz,nt)+ght(i,j-1,nz,nt)-2.0*ght(i,j,nz,nt))/((dxx)**2)&
                                                                                  -(ght(i,j+1,nz,nt)-ght(i,j-1,nz,nt))*tan(lat)/(2.0*dxx)    )     !计算地转涡度
                   write(2) ksei(i,j,tz,tt)
          enddo
   enddo
   close(2)
   print*,'write done'
end

surelee

ctl文件代码

1. dset d:\tex2\Apr2015_ksei.grd
   
2. title ksei
   
3. undef -999.0
   
4. xdef 479 linear 0.75 0.75
   
5. ydef 240 linear -89.25 0.75
   
6. zdef 1   levels 500
   
7. tdef 1   linear 12z28Apr2015 1dy
   
8. vars 1
   
9. ksei 0 99
   
10. endvars

复制代码

lqouc

ksei的ctl空间格点数在算下，我觉得fortran 输出的是477和238。
另外位势高度的缩写一般是HGT

surelee

lqouc 发表于 2016-5-8 22:36
ksei的ctl空间格点数在算下，我觉得fortran 输出的是477和238。
另外位势高度的缩写一般是HGT

好的我改改试试。

querge

楼主这是哪的参考公式，可以给出出处么，我看的咋跟你的不太一样呢

木易

想问下楼主，你的问题解决了吗？我用你的这个f程序算了一下地转风涡度，可是计算的数值都显示-lns，不出图，感觉好像是计算的不太对呢，公式如果对的话，分母那里面的几个三角函数，是弧度还是度数啊，那个拉姆达是弧长么？？感觉程序里是用dxx代替了，这样也可以？

tonysonic

很好，非常感谢

cgccq

学习了，非常好！

王小二

学习以下~~~~