phalf_etkf.f90

Go to the documentation of this file.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Time-stamp: <2016-10-20 10:45:32 pbrowne>
!!!
!!!    Routine to change an ensemble N(0,I) to N(0,P)
!!!    Copyright (C) 2015 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
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


subroutine phalf_etkf(nrhs,x,px)
  use output_empire, only : emp_e
  use pf_control
  use sizes
  use comms
  implicit none
  include 'mpif.h'

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

  integer, intent(in) :: nrhs
  real(kind=rk), dimension(state_dim,nrhs), intent(in) :: x
  real(kind=rk), dimension(state_dim,nrhs), intent(out) :: px

  real(kind=rk), dimension(state_dim,pf%count) :: Xp_p
  real(kind=rk), dimension(obs_dim,pf%count) :: yf_p
  real(kind=rk), dimension(obs_dim,pf%nens) :: yf,Ysf
  real(kind=rk), dimension(obs_dim_g,pf%nens) :: Ysf_g
  real(kind=rk), dimension(obs_dim) :: mean_yf
  integer, dimension(npfs) :: start_var,stop_var
  !real(kind=rk), dimension(stateDimension,N) :: Xp
  real(kind=rk), dimension(:,:), allocatable :: Xp
  real(kind=rk), dimension(:,:), allocatable :: Xa
  integer :: mpi_err
  integer :: i,j,number_gridpoints

  integer :: red_obsdim
  real(kind=rk), allocatable, dimension(:,:) :: Ysf_red

  ! Local variables for the SVD
  integer :: r
  real(kind=rk), dimension(:,:), allocatable :: V
  real(kind=rk), dimension(:), allocatable :: S
  real(kind=rk), dimension(pf%nens,pf%nens) :: UT
  integer :: LWORK,INFO
  real(kind=rk), dimension(:), allocatable :: WORK
  integer :: stateDim
  !!> Forecast ensemble on entry, analysis ensemble on exit
  !real(kind=rk), dimension(stateDimension,N), intent(inout) :: x
  real(kind=rk), dimension(state_dim,pf%count) :: x_p
  
  integer :: ensemble_comm

  !first lets make something N(0,Q):
  call qhalf(nrhs,x,xp_p)

  if(comm_version .eq. 1 .or. comm_version .eq. 2) then
     ensemble_comm = pf_mpi_comm
  elseif(comm_version .eq. 3) then
     ensemble_comm = pf_ens_comm
  else
     write(emp_e,*) 'EMPIRE VERSION ',comm_version,' NOT SUPPORTED IN phalf_etkf'
     write(emp_e,*) 'THIS IS AN ERROR. STOPPING'
     stop '-24'
  end if


  !now go into the LETKF code:

  x_p = xp_p


  ! Calculate forecast observations, split into mean and ensemble
  ! perturbation matrix, scale perturbations by inverse square root of
  ! observation covariance

  ! first apply observation operator only to local state vectors
  ! on this mpi process
  call h(obs_dim,pf%count,x_p,yf_p,pf%timestep)


  ! as yf_p should be much smaller than x_p, send this to mpi processes
  ! need to send round all yf_p and store in yf on all processes
  call mpi_allgatherv(yf_p,pf%count*obs_dim,mpi_double_precision,yf&
       &,gblcount*obs_dim,gbldisp*obs_dim,mpi_double_precision&
       &,ensemble_comm,mpi_err)

  ! compute the mean of yf
  mean_yf = sum(yf,dim=2)/real(pf%nens,rk)
  do i = 1,pf%nens
     yf(:,i) = yf(:,i) - mean_yf
  end do

  ! now make yf the forecast perturbation matrix
  yf = yf/sqrt(real(pf%nens-1,rk))
  ! scale yf to become ysf
  call solve_rhalf(obs_dim,pf%nens,yf,ysf,pf%timestep)


  ! now let us compute which state variables will be analysed on each
  ! MPI process:
  do i = 1,npfs
     start_var(i) = ((i-1)*state_dim)/npfs + 1 
     stop_var(i) = (i*state_dim)/npfs
  end do

  !allocate space for Xp and Xa now that we know how many grid points we consider
  allocate(xp(stop_var(pfrank+1)-start_var(pfrank+1)+1,pf%nens))
  allocate(xa(stop_var(pfrank+1)-start_var(pfrank+1)+1,pf%nens))

  !now we have to get Xp filled with all the values from each process
  do i = 1,npfs
     number_gridpoints = stop_var(i)-start_var(i)+1

     call mpi_gatherv(xp_p(start_var(i):stop_var(i),:),& !sendbuf
          pf%count*number_gridpoints,&                     !sendcount
          mpi_double_precision,&                           !sendtype
          xp,&                                             !recvbuf
          gblcount*number_gridpoints,&                     !recvcounts
          gbldisp*number_gridpoints,&                      !displs
          mpi_double_precision,&                           !recvtype
          i-1,&                                            !root
          ensemble_comm,&                                    !comm
          mpi_err)                                         !ierror
  end do


  !PAB  d = y - mean_yf
  !PAB  call solve_rhalf(obs_dim,1,d,dd,pf%timestep)

  if(comm_version .eq. 1 .or. comm_version .eq. 2) then
     ysf_g = ysf
  elseif(comm_version .eq. 3) then
     call mpi_allgatherv(ysf,obs_dim*pf%nens,mpi_double_precision&
          &,ysf_g,obs_dims*pf%nens,obs_displacements*pf%nens&
          &,mpi_double_precision,pf_member_comm,mpi_err)
  else
     print*,'should not enter here. error. stopping :o('
  end if

  !this is for serial processing
  number_gridpoints = stop_var(pfrank+1)-start_var(pfrank+1)+1
  !$OMP PARALLEL DO &
  !$OMP& PRIVATE(stateDim,i), &
  !$OMP& PRIVATE(Ysf_red,red_obsDim,r,V,S,LWORK,WORK,INFO), &
  !$OMP& PRIVATE(UT)
  do j = 1,number_gridpoints
     statedim = 1

     ! count the total number of observations we shall consider for this
     ! state variable
     red_obsdim = obs_dim_g


     allocate(ysf_red(red_obsdim,pf%nens))
     ! reduce the forecast ensemble to only the observations in range
     ! and scale by the distance function
     do i = 1,pf%nens
        ysf_red(:,i) = ysf_g(:,i)
     end do

     ! Compute the SVD
     r = min(red_obsdim,pf%nens)
     allocate(v(red_obsdim,r))
     allocate(s(r))
     lwork = 2*max( 3*r+max(red_obsdim,pf%nens), 5*r )

     allocate(work(lwork))

     call dgesvd('S','A',red_obsdim,pf%nens,ysf_red,red_obsdim,s,v,red_obsdim,ut,pf%nens,work,lwork,info)
     if(info .ne. 0) then
        write(emp_e,*) 'EMPIRE ERROR in phalf_etkf.f90 with the SVD'
        write(emp_e,*) 'SVD failed with INFO = ',info
        write(emp_e,*) 'FYI WORK(1) = ',work(1)
        write(emp_e,*) 'STOPPING'
        stop
     end if
     deallocate(work)

     ! Compute product of forecast ensemble perturbation matrix and U.  We
     ! store the result in xa, which we here call Xa to distinguish this
     ! secondary use of the variable.
     ! pick out the correct row of Xa and Xp here:
     call dgemm('N','T',statedim,pf%nens,pf%nens,1.0d0,xp(j,:),statedim,ut,pf%nens,0.0d0,xa(j,:),statedim)


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


     !put this back into the full state vector
     !wait it is already
     deallocate(ysf_red)
     deallocate(v)
     deallocate(s)


  end do
  !$OMP END PARALLEL DO



  ! now we must scatter xa back into px
  do i = 1,npfs

     number_gridpoints = stop_var(i)-start_var(i)+1

     call mpi_scatterv(xa,&                                !sendbuf  
          gblcount*number_gridpoints,&                     !sendcounts  
          gbldisp*number_gridpoints,&                      !displs
          mpi_double_precision,&                           !sendtype    
          px(start_var(i):stop_var(i),:),&                 !recvbuf  
          pf%count*number_gridpoints,&                     !recvcount  
          mpi_double_precision,&                           !recvtype    
          i-1,&                                            !root     
          ensemble_comm,&                                    !comm     
          mpi_err)                                         !ierror 
  end do



  deallocate(xp)
  deallocate(xa)



end subroutine phalf_etkf