求助关于olr资料的滤波错误

羊羊羊

我用Barnes方法滤波，之前做的高度场和风场滤出来都没有错（老师看过的），但换成olr资料滤波后用grads出来的图就不对，求助各位大神~~
资料：我用的是olr.day.mean.nc，2.5*2.5的
步骤：1、先把nc资料用grads转换成dat格式（经过对生成的dat画图与原始nc出的图对比，验证生成的dat没有错误）
           2、用Barnes滤波程序对dat滤波，之后写ctl出图错误

lqouc

论坛的大神帮不了你，你这么问，只有全知全能的真神能拯救你。

羊羊羊

lqouc 发表于 2015-5-6 13:51
论坛的大神帮不了你，你这么问，只有全知全能的真神能拯救你。

不好意思，我帖子还没发完，要表达的东西太多，显示超出了字数界限，就一直没发后面的东西，稍等一下，我再看看有没有别的办法发出来

羊羊羊

这是原始的图和滤波后错误的图
老师说，错误很可能是由于没有对olr资料的缺省数据进行处理直接滤波，但我把nc转换成的dat用txt打出来看，并没有发现异常数据。

lqouc

羊羊羊 发表于 2015-5-6 14:07
不好意思，我帖子还没发完，要表达的东西太多，显示超出了字数界限，就一直没发后面的东西，稍等 ...

晕，这玩意没啥字数限制吧，你直接编辑这个帖子贴图贴代码就可以了。多出来的在下面自己回复就好了。

羊羊羊

这是滤波程序：
!include 'distance.inc'
program main
implicit none
integer,parameter::xx=144,yy=73,dg=1
real,parameter::C1=10000,C2=150000,G=0.3,pi=3.14

real,external::distance
integer::i,j,k,l
real::r,ii,jj,kk,ll
real::hgt(yy,xx),F1(yy,xx),F2(yy,xx),F(yy,xx),D1(yy,xx),D2(yy,xx),&
    E1(yy,xx),E2(yy,xx),F11(yy,xx),F21(yy,xx),F12(yy,xx),F22(yy,xx),t1,t2,w1,w2
open(unit=1,file='f:\tbbdata\1982_02\olr\olr1.dat',form='unformatted',access='direct',recl=1,status='old')
k=0
do i=1,yy
  do j=1,xx
   k=k+1
   read(unit=1,rec=k) hgt(i,j)
  enddo
enddo
close(1)
!write(*,*) hgt(1:5,1:5)
!method1

do i=1,yy
  do j=1,xx
   w1=0
   t1=0
   w2=0
   t2=0
   do k=1,yy,dg
    do l=1,xx,dg
   r=distance(-90.0+(i-1)*2.5,-90.0+(k-1)*2.5,(j-1)*2.5,(l-1)*2.5)       
     if(r<=sqrt(C1))then
      w1=w1+exp(-r**2/(4*C1))
      t1=t1+hgt(k,l)*exp(-r**2/(4*C1))
     endif
     if(r<=sqrt(C2))then
      w2=w2+exp(-r**2/(4*C2))
      t2=t2+hgt(k,l)*exp(-r**2/(4*C2))
     endif
    enddo
   enddo
   F1(i,j)=t1/w1
   F2(i,j)=t2/w2
  enddo
enddo

print *,'Finish method1'
F=F1-F2

open(unit=20,file='f:\tbbdata\1982_02\olr\olr01.dat',form='binary')
do i=1,yy
  do j=1,xx
   write(10) F(i,j)
  enddo
enddo
!method2
D1=hgt-F1
D2=hgt-F2
do i=1,yy
  do j=1,xx
   w1=0
   t1=0
   w2=0
   t2=0
   do k=1,yy,dg
    do l=1,xx,dg
          r=distance(-90.0+(i-1)*2.5,-90.0+(k-1)*2.5,(j-1)*2.5,(l-1)*2.5)
     if(r<=sqrt(C1))then
      w1=w1+exp(-r**2/(4*G*C1))
      t1=t1+D1(k,l)*exp(-r**2/(4*G*C1))
     endif
     if(r<=sqrt(C2))then
      w2=w2+exp(-r**2/(4*G*C2))
      t2=t2+D2(k,l)*exp(-r**2/(4*G*C2))
     endif
    enddo
   enddo
   F11(i,j)=F1(i,j)+t1/w1
   F21(i,j)=F2(i,j)+t2/w2
  enddo
enddo

print *,'Finish method2'
F=F11-F21
do i=1,yy
  do j=1,xx
   write(10) F(i,j)
  enddo
enddo
!method3
E1=F11-F1
E2=F21-F2
do i=1,yy
  do j=1,xx
   w1=0
   t1=0
   w2=0
   t2=0
   do k=1,yy,dg
    do l=1,xx,dg
          r=distance(-90.0+(i-1)*2.5,-90.0+(k-1)*2.5,(j-1)*2.5,(l-1)*2.5)
     if(r<=sqrt(C1))then
      w1=w1+exp(-r**2/(4*C1))
      t1=t1+E1(k,l)*exp(-r**2/(4*C1))
     endif
     if(r<=sqrt(C2))then
      w2=w2+exp(-r**2/(4*C2))
      t2=t2+E2(k,l)*exp(-r**2/(4*C2))
     endif
    enddo
   enddo
   F12(i,j)=F11(i,j)+3*(F11(i,j)-F1(i,j))/4-t1/w1
   F22(i,j)=F21(i,j)+3*(F21(i,j)-F2(i,j))/4-t2/w2
  enddo
enddo

print *,'Finish method3'
F=F12-F22
! write(*,*) F(1:5,1:5)
do i=1,yy
  do j=1,xx
   write(20) F(i,j)
  enddo
enddo
close(20)
end program main

function distance(lat1,lat2,lon1,lon2)
implicit none
real,parameter::R=6371.004,pi=3.1415926
real::lat1,lat2,lon1,lon2
real::distance,temp
temp=(cos(lat1*pi/180)*cos(lat2*pi/180)*cos((lon1-lon2)*pi/180)+sin(lat1*pi/180)*sin(lat2*pi/180))
distance=R*acos(temp)
return
end function distance

羊羊羊

原始的数据以及滤波后的数据：olr1.dat是原始的，olr01.dat是滤波后的

羊羊羊

这是出图的ctl：
dset f:/tbbdata/1982_02/olr/olr01.dat
undef -9.99e+33
title OLR
xdef 144 linear 0 2.5
ydef 73 linear -90 2.5
zdef 1 linear 0 1
tdef 1  linear 30MAY1982 1dy
vars 1
olr 0 99 olr
endvars

羊羊羊

lqouc 发表于 2015-5-6 14:16
晕，这玩意没啥字数限制吧，你直接编辑这个帖子贴图贴代码就可以了。多出来的在下面自己回复就好了。

不太会发帖，之前加图的时候直接复制粘贴了，导致帖子过大发不出去，鼓捣了半天才觉悟了，但是为什么发附件要扣下载的金钱呢，我请教人不应该扣别人的金钱呀，有没有办法取消掉？

lqouc

羊羊羊 发表于 2015-5-6 14:25
不太会发帖，之前加图的时候直接复制粘贴了，导致帖子过大发不出去，鼓捣了半天才觉悟了，但是 ...

下载需要的金币是取消不了的。
话说你从哪里觉得这个olr的图是错的呢？之前你认为正确的高度场滤波前后是什么样呢？