etkf_analysis.f90

Go to the documentation of this file.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Time-stamp: <2014-09-25 11:46:38 pbrowne>
!!!
!!!    Ensemble transform Kalman filter
!!!    Copyright (C) 2014  Philip A. Browne
!!!
!!!    This program is free software: you can redistribute it and/or modify
!!!    it under the terms of the GNU General Public License as published by
!!!    the Free Software Foundation, either version 3 of the License, or
!!!    (at your option) any later version.
!!!
!!!    This program is distributed in the hope that it will be useful,
!!!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!!!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!!!    GNU General Public License for more details.
!!!
!!!    You should have received a copy of the GNU General Public License
!!!    along with this program.  If not, see <http://www.gnu.org/licenses/>.
!!!
!!!    Email: p.browne @ reading.ac.uk
!!!    Mail:  School of Mathematical and Physical Sciences,
!!!           University of Reading,
!!!       Reading, UK
!!!       RG6 6BB
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! $RCSfile: etkf_analysis.f90,v $
! $Revision: 1.6 $
! $Date: 2013/12/12 03:58:54 $

subroutine etkf_analysis(num_hor,num_ver,this_hor,this_ver,boundary,x,N,stateDim,obsDim,rho)

implicit none
integer, parameter :: rk = kind(1.0d0)

integer, intent(in) :: N,stateDim,obsDim
! Forecast ensemble on entry, analysis ensemble on exit
real(kind=rk), dimension(stateDim,N), intent(inout) :: x
! The observation
!real(kind=rk), dimension(obsDim), intent(in) :: y
real(kind=rk), dimension(obsDim) :: y
! Inflation parameter; forecast perturbations will be scaled by 1+rho
! if rho is present
real(kind=rk), intent(in) :: rho
! controls for localisation
integer, intent(in) :: num_hor,num_ver,this_hor,this_ver,boundary

! Local variables for the SVD
integer :: r
real(kind=rk), dimension(:,:), allocatable :: V
real(kind=rk), dimension(:), allocatable :: S
real(kind=rk), dimension(N,N) :: UT
integer :: LWORK,INFO
real(kind=rk), dimension(:), allocatable :: WORK

! Miscellaneous local variables
real(kind=rk), dimension(stateDim) :: mean_x
real(kind=rk), dimension(stateDim,N) :: Xp
real(kind=rk), dimension(obsDim) :: mean_yf,d,dd
real(kind=rk), dimension(obsDim,N) :: yf,Ysf
integer :: i

call get_local_observation_data(num_hor,num_ver,this_hor,this_ver,boundary,obsdim,y)

! Split forecast ensemble into mean and perturbation matrix, inflating
! if necessary
mean_x = sum(x,dim=2)/real(n,rk)
do i = 1,n
   xp(:,i) = x(:,i) - mean_x
end do
!if (present(rho)) then
xp = (1.0_rk + rho) * xp
do i = 1,n
   x(:,i) = mean_x + xp(:,i)
end do
!end if
xp = xp/sqrt(real(n-1,rk))

! Calculate forecast observations, split into mean and ensemble
! perturbation matrix, scale perturbations by inverse square root of
! observation covariance
!call H(N,stateDim,obsDim,x,yf)
call h_local(num_hor,num_ver,this_hor,this_ver,boundary,n,statedim,x&
     &,obsdim,yf)
mean_yf = sum(yf,dim=2)/real(n,rk)
do i = 1,n
   yf(:,i) = yf(:,i) - mean_yf
end do
yf = yf/sqrt(real(n-1,rk))
!call solve_rhalf(N,obsDim,yf,Ysf)
call solve_rhalf_local(num_hor,num_ver,this_hor,this_ver,boundary,n&
     &,obsdim,yf,ysf)

! Compute the SVD
r = min(obsdim,n)
allocate(v(obsdim,r))
allocate(s(r))
lwork = 2*max( 3*r+max(obsdim,n), 5*r )
allocate(work(lwork))
call dgesvd('S','A',obsdim,n,ysf,obsdim,s,v,obsdim,ut,n,work,lwork,info)
if(info .ne. 0) then
   print*,'SVD failed with INFO = ',info
   print*,'FYI WORK(1) = ',work(1)
   stop
end if
deallocate(work)

! Compute product of forecast ensemble perturbation matrix and U.  We
! store the result in x, which we here call X to distinguish this
! secondary use of the variable.
call dgemm('N','T',statedim,n,n,1.0d0,xp,statedim,ut,n,0.0d0,x,statedim)

! Build up analysis ensemble mean (to be stored in mean_x); done now
! because we shall be reusing X shortly
d = y - mean_yf
!call solve_rhalf(1,obsDim,d,dd)
call solve_rhalf_local(num_hor,num_ver,this_hor,this_ver,boundary,1&
     &,obsdim,d,dd)
! Only the first r elements of d are in use from here on
call dgemv('T',obsdim,r,1.0d0,v,obsdim,dd,1,0.0d0,d,1)
d(1:r) = s * d(1:r) / (1.0_rk + s**2)
! Only the first r columns of X are used in the following
call dgemv('N',statedim,r,1.0d0,x,statedim,d,1,1.0d0,mean_x,1)

! Build up analysis ensemble perturbation matrix
do i = 1,r
   x(:,i) = x(:,i) / sqrt(1.0_rk + s(i)**2)
end do
call dgemm('N','N',statedim,n,n,1.0d0,x,statedim,ut,n,0.0d0,xp,statedim)

! Put ensemble mean and perturbation matrix back together
x = sqrt(real(n-1,rk))*xp
do i = 1,n
   x(:,i) = mean_x + x(:,i)
end do

end subroutine etkf_analysis