小弟初学ETKF,求教程序

一个人闯梦

小弟初次学习ETKF方法，想在小模式（Lorenz96）中实验一下，可
更新后的集合扰动预报效果还不如未更新的扰动，特来求教程序该怎么编写，感激不尽！

还未涉及到膨胀因子和局地化等，仅仅是最初的ETKF方法（Bishop，2001）

不知道观测算子和观测误差协方差该怎么给定，程序就胡乱取了个单位矩阵I
初始集合是随机生成的

ensize_inv = 1.0 / real(ensize)
ensize_m1_inv = 1.0 / real(ensize-1)
!========================================================================  
![1] 计算集合平均ensmean，产生集合扰动Zf
!========================================================================
open(30,file='Xf.txt')
ensmean(:)=0.0
      do i=1,dim
          do j=1,ensize
          ensmean(i) = ensmean(i)+xf(i,j)
          end do
          ensmean(i)=ensmean(i)*ensize_inv
      end do
      
      do j = 1, ensize
          xf(:,j) = xf(:,j) - ensmean(:)                  
      enddo
      xf=xf*ensize_m1_inv                                        !Xf
           
     write(30,"(40f12.6)") xf
                                                
!========================================================================  
![2] 计算hzthz(ensize,ensize)
!========================================================================
    call Obs(dim,nobs,ob_H,oberrcovar)       !调用的观测算子H和观测误差协方差R都是单位阵，不知道该怎么取？？
      
    do i=1,nobs
            oberrcovar(i,i)=1.0/sqrt(oberrcovar(i,i))          !R^(-0.5)
    end do
   
   call Matrix_Mul(oberrcovar,ob_H,nobs,nobs,dim,ob_tilda_H)   !R^(-0.5)*H = Tilda(H)


   call Matrix_Mul(ob_tilda_H,xf,nobs,dim,ensize,ens_ob)       !H*Zf=HZf  因为H,R都是单位阵，就把HXf当作了HZf，对吗？
   
   call MatTran_Mul_Mat(ens_ob, hzthz,nobs,ensize)              !(HZf)'(HZf)

!=====================================================================
![3] 计算hzthz的特征值gamma和特征向量C
!=====================================================================
       call eig(hzthz,ensize,EPS,eig_vector1)          ！求实对称矩阵特征值与特征向量的雅克比过关法
      
       do ij=1,ensize
             eig_value1(ij)=hzthz(ij,ij)                        !hzthz主对角线元素为特征值
       end do
      
!------------------------特征值从大到小排列，特征向量对应排列-----------------------------------------   
eig_value(:)=0.0
eig_vector(:,:)=0.0
       do ik=1,ensize
           indx=maxloc(abs(eig_value1))
           eig_value(ik)=eig_value1(indx(1))
           eig_vector(:,ik)=eig_vector1(:,indx(1))         !第I列特征向量对应第I个特征值
  
           eig_value1(indx(1))=0.0   
           eig_vector1(:,indx(1))=0.0
       end do
      
!================================================================================================
!                       转换矩阵T
!================================================================================================
open(35,file='matrix_T.txt')
do ii=1,ensize
    do jj=1,ensize
        T(jj,ii) = eig_vector(jj,ii)*(1.0/sqrt(eig_value(ii)+1.0))
    end do
end do
write(35,"(15f12.6)")T
!==========================================================================================
open(36,file='Xa.txt')
call Matrix_Mul(xf,T,dim,ensize,ensize,xa)     ！xa=xf*T
write(36,"(40f12.6)")xa

最后得到的xa就是更新后扰动，但叠加在集合平均上做出来的预报还不如Xf叠加平均态的预报效果
求高手帮我看看

我最大的疑问在于观测算子H，和观测误差协方差阵R该怎么给予
还有初始的集合随机生成是否可以？
求高手帮我看看
感激不尽！

一个人闯梦

希望能得到你们的帮助

一个人闯梦

没有人会吗？都沉了

一个人闯梦

在那遥远的地方 发表于 2016-1-12 22:42
wrfda3.5以后的手册里面都要现成的，不知道能不能符合你的要求

我还没接触过wrfda，不过还是先谢谢你