WRF初始场中尺度滤波试验出错求助！

schliezer

晓小 发表于 2015-10-8 21:17
你好，我也正在学滤波，请问可以给我个滤波的程序吗

program main
implicit none
integer,parameter::xx=71,yy=41,zz=1,tt=13,dg=1
!global scale,1*1
real,parameter::C1=4000,C2=16000,G=0.3

real,external::distance
integer::h,i,j,k,l,m,n
real::r
real::hgt(xx,yy,zz,tt),F1(xx,yy,zz,tt),F2(xx,yy,zz,tt),F(xx,yy,zz,tt),D1(xx,yy,zz,tt),D2(xx,yy,zz,tt),E1(xx,yy,zz,tt),E2(xx,yy,zz,tt),F11(xx,yy,zz,tt),F21(xx,yy,zz,tt),F12(xx,yy,zz,tt),F22(xx,yy,zz,tt),t1,t2,w1,w2

!!!!!!!!!!!!!!!!!!!!!!!!!!!    read initial field        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
open(10,file='E:\ScientificResearchDuringMaster\Wrf_out\19-21\wavefiltering\0719PRESwavefilter.dat',access='direct',recl=4,status='old')
open(20,file='E:\ScientificResearchDuringMaster\Wrf_out\19-21\wavefiltering\PRESfiltered.grd',form='binary')
open(30,file='E:\ScientificResearchDuringMaster\Wrf_out\19-21\wavefiltering\PRESditong.grd',form='binary')
  h=0
do m=1,tt
do n=1,zz
  do j=1,yy
   do i=1,xx
    h=h+1
    read(10,rec=h) hgt(i,j,n,m)
enddo
enddo
enddo
enddo

do m=1,tt
do n=1,zz
  do j=1,yy
   do i=1,xx
    w1=0
    t1=0
    w2=0
    t2=0
   do k=1,yy,dg
    do l=1,xx,dg
     r=distance(-90.0+(j-1)*1.0,-90.0+(k-1)*1.0,(i-1)*1.0,(l-1)*1.0)
     if(r<=sqrt(C1))then
      w1=w1+exp(-r**2/(4*C1))
      t1=t1+hgt(l,k,n,m)*exp(-r**2/(4*C1))
     endif
     if(r<=sqrt(C2))then
      w2=w2+exp(-r**2/(4*C2))
      t2=t2+hgt(l,k,n,m)*exp(-r**2/(4*C2))
     endif
enddo
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   F1(i,j,n,m)=t1/w1     !meso scale
   F2(i,j,n,m)=t2/w2     !large scale
enddo
enddo
enddo
enddo

D1=hgt-F1
D2=hgt-F2

do m=1,tt
do n=1,zz
  do j=1,yy
   do i=1,xx
   w1=0
   t1=0
   w2=0
   t2=0
   do k=1,yy,dg
    do l=1,xx,dg
     r=distance(-90.0+(j-1)*1.0,-90.0+(k-1)*1.0,(i-1)*1.0,(l-1)*1.0)
     if(r<=sqrt(C1))then
      w1=w1+exp(-r**2/(4*G*C1))
      t1=t1+D1(l,k,n,m)*exp(-r**2/(4*G*C1))
     endif
     if(r<=sqrt(C2))then
      w2=w2+exp(-r**2/(4*G*C2))
      t2=t2+D2(l,k,n,m)*exp(-r**2/(4*G*C2))
     endif
    enddo
   enddo

   F11(i,j,n,m)=F1(i,j,n,m)+t1/w1
   F21(i,j,n,m)=F2(i,j,n,m)+t2/w2
enddo
enddo
enddo
enddo
   E1=F11-F1
   E2=F21-F2

do m=1,tt
do n=1,zz
  do j=1,yy
   do i=1,xx
   w1=0
   t1=0
   w2=0
   t2=0
   do k=1,yy,dg
    do l=1,xx,dg
     r=distance(-90.0+(j-1)*1.0,-90.0+(k-1)*1.0,(i-1)*1.0,(l-1)*1.0)
     if(r<=sqrt(C1))then
      w1=w1+exp(-r**2/(4*C1))
      t1=t1+E1(l,k,n,m)*exp(-r**2/(4*C1))
     endif
     if(r<=sqrt(C2))then
      w2=w2+exp(-r**2/(4*C2))
      t2=t2+E2(l,k,n,m)*exp(-r**2/(4*C2))
     endif
    enddo
   enddo
   F12(i,j,n,m)=F11(i,j,n,m)+3*(F11(i,j,n,m)-F1(i,j,n,m))/4-t1/w1
   F22(i,j,n,m)=F21(i,j,n,m)+3*(F21(i,j,n,m)-F2(i,j,n,m))/4-t2/w2
   F(i,j,n,m)=F12(i,j,n,m)-F22(i,j,n,m)
   write(20) F(i,j,n,m)
   write(30) F22(i,j,n,m)
!   print *, i,j,n,m
enddo
enddo
enddo
enddo

close(10)
close(20)
close(30)
end program main

! write(*,*) hgt(1:10,1:10)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!   calculate the 'F0' in the Barnes formula!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!  calculate the 'F2' under different 'C' in the Barnes formula (the second revision)  !!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!calculate the spherical distance between (i,j)  and (k,l)!!!!!!!!!!!!!!!!
function distance(lat1,lat2,lon1,lon2)
implicit none
real,parameter::R=6371.004,pi=3.1415926
real::lat1,lat2,lon1,lon2
real::distance
distance=R*acos(cos(lat1*pi/180)*cos(lat2*pi/180)*cos((lon1-lon2)*pi/180)+sin(lat1*pi/180)*sin(lat2*pi/180))
return
end function distance

晓小

谢谢啊~~

寻找费曼的彩虹

不知楼主解决问题没有？能分享些经验吗？
再有，为什么不在wps前滤波？

yeah...

schliezer 发表于 2015-10-9 12:01
program main
implicit none
integer,parameter::xx=71,yy=41,zz=1,tt=13,dg=1

请问这个一定要用dat格式的数据去算么 用grd的数据去算可以吗

schliezer

yeah... 发表于 2017-2-15 16:21
请问这个一定要用dat格式的数据去算么 用grd的数据去算可以吗

应该也可以，这是针对二进制数据。

schliezer

寻找费曼的彩虹 发表于 2016-3-17 23:02
不知楼主解决问题没有？能分享些经验吗？
再有，为什么不在wps前滤波？

wps前是对fnl资料滤波，滤波完成也得是nc格式的才能被wps处理，但当时还不知道怎么用fortran读写nc文件

yeah...

schliezer 发表于 2017-2-16 09:18
应该也可以，这是针对二进制数据。

谢谢回复 我想请教下 您最后的c1 c2 G选择的多少， 我选了几次参数  效果并不理想  还有响应曲线您有计算过吗？

schliezer

yeah... 发表于 2017-2-17 08:37
谢谢回复 我想请教下 您最后的c1 c2 G选择的多少， 我选了几次参数  效果并不理想  还有响应曲线您有计算 ...

这个要多调试，我取了5000，50000和0.3，图放不上来，可惜

yeah...

schliezer 发表于 2017-2-17 09:27
这个要多调试，我取了5000，50000和0.3，图放不上来，可惜

哦哦哦 谢谢您 那请问您有计算响应曲线吗？就是那个r和拉姆打的关系

yeah...

schliezer 发表于 2017-2-17 09:27
这个要多调试，我取了5000，50000和0.3，图放不上来，可惜

还有我想请教下 这个Barnes滤波能实现高通滤波吗？