m_analysis.F90 10.7 KB

module m_analysis
contains
subroutine prepare_and_analysis(mem,obs,nrsamp)
!  This subroutine computes the EnKF analysis. The improved analysed 
!  estimate psia(iens), when psif(iens) is the forcast, is given by
!  the following equation: 
!  psia(iens)=psif(iens)+{[H*(A-Mean)*(A-Mean)^T]^T*scale}*P^(-1)*h
!  H is the observation matrix, A contains the ensemble, Mean contains
!  the ensemble mean, h is the residual between the observations and 
!  the ensemble members at the position of the observations.

   use mod_dimensions
   use mod_states
   use mod_observations
   use m_random
   implicit none

! Interface
   integer,            intent(in) :: nrsamp
   type(observations), intent(in) :: obs(nrobs)
   type(states),       intent(inout) :: mem(nrsamp)

! Locals
   integer iens, i, j, k, m, nrsigma, obstot
   type(observations), allocatable :: d(:)      ! The active observations
   type(observations), allocatable :: tmpobs(:) ! Observations with noise

   type(states) Mean                            ! ensemble mean

   type(states), allocatable :: tmpens(:)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   real,    allocatable  :: h(:,:), E(:,:), d_vec(:)
   real,    allocatable  :: tmp(:,:)
   real,    allocatable  :: P(:,:)
 
  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! For eigen value decomposition
   real,    allocatable  :: ZZ(:,:)  !contains the eigenvectors of P
   real,    allocatable  :: w(:)     !dsyevx workspace
   real,    allocatable  :: ss(:)    !contains the eigenvalues of P
   integer, allocatable  :: iw(:)    !dsyevx workspace
   integer, allocatable  :: ipvt(:)  !dsyevx workspace
   integer                  info     !completion code for dsyevx 
   real                     abstol   !The abs error tolerance for the eigenvalues
   real,   allocatable   :: ap(:)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   real,    allocatable  :: lambda(:)
   integer                  nndim      !parameter in the dgemm subroutine
   type(states), allocatable :: rtd(:) !Damped Transpose of Representer matrix
   integer                  ix,iy,iz,m2,m3,idim !temporary counters
   integer                  toteig   !total number of eigenvectors found from subroutine dsyevx
   real                     VL, VU
   integer                  IL, IU
   integer mp_my_threadnum, mp_numthreads
   real                     summ, summ2
   real,    allocatable  :: L(:,:)
   real                     scale      !1/nrsamp
   character(len=4) tag4
   character(len=3) tag3
   character(len=7) tag7

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute ensemble mean
   scale = 1.0/float(nrsamp)
   Mean=0.0
   do iens=1,nrsamp
      Mean = Mean + mem(iens)
   enddo
   Mean=scale*Mean

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Pick out active observations
   obstot=0
   do m=1,nrobs
      if (obs(m)%action) obstot=obstot+1
   enddo
   print *,'obstot= ',obstot
   allocate(d(obstot))
   obstot=0
   do m=1,nrobs
      if (obs(m)%action) then
         obstot=obstot+1
         d(obstot)=obs(m)
      endif
   enddo

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute the residual between the observations and the value of the
! ensemble members at the position of the observations.
!        "h(:,iens)=d(:)%d-H*mem(iens)"
! and the matrix  "tmp = H*(A(iens) - Mean)" which are used in the analysis.

   allocate(h(obstot,nrsamp))
   allocate(E(obstot,nrsamp))
   allocate(tmp(nrsamp,obstot))  !GE used instead of allocate(tmp(obstot,nrsamp))
   allocate(tmpobs(obstot))
   allocate(d_vec(obstot))

! Compute non-perturbed innovations
   do m=1, obstot
      if (trim(d(m)%ch)=='N') then
         d_vec(m) = d(m)%d - Mean%N(d(m)%k)
      elseif (trim(d(m)%ch)=='P') then
         d_vec(m) = d(m)%d - Mean%P(d(m)%k)
      elseif (trim(d(m)%ch)=='H') then
         d_vec(m) = d(m)%d - Mean%H(d(m)%k)
      else
         stop 'analysis: error computing d_vec'
      endif
   enddo

   do iens=1,nrsamp

! Add noise to the observations 
      call random(tmpobs(:)%d,obstot) 
      call random(E(:,iens),obstot) 

      do m=1,obstot
         tmpobs(m)%ch=d(m)%ch
         tmpobs(m)%d = d(m)%d + sqrt(d(m)%var)*tmpobs(m)%d
         tmpobs(m)%d=max(0.0,tmpobs(m)%d)  ! ME: Pga neg. obs.
         !E(m,iens) = tmpobs(m)%d - d(m)%d
         E(m,iens) =  E(m,iens) * sqrt(d(m)%var)
      enddo

      do m=1,obstot
! h(:,iens) = obs(:)%d - H*mem(iens)
         if (trim(d(m)%ch)=='N') then
            h(m,iens)=tmpobs(m)%d-mem(iens)%N(d(m)%k)  
         elseif (trim(d(m)%ch)=='P') then
            h(m,iens)=tmpobs(m)%d-mem(iens)%P(d(m)%k)  
         elseif (trim(d(m)%ch)=='H') then
            h(m,iens)=tmpobs(m)%d-mem(iens)%H(d(m)%k)  
         else
            stop 'analysis: errror computing h'
         endif

! tmp = H*(mem(iens) - Mean)
         if (trim(tmpobs(m)%ch)=='N') then
            tmp(iens,m)=mem(iens)%N(d(m)%k)  - Mean%N(d(m)%k)
         elseif (trim(tmpobs(m)%ch)=='P') then
            tmp(iens,m)=mem(iens)%P(d(m)%k)  - Mean%P(d(m)%k)
         elseif (trim(tmpobs(m)%ch)=='H') then
            tmp(iens,m)=mem(iens)%H(d(m)%k)  - Mean%H(d(m)%k)
         else
            stop 'analysis: error computing tmp'
         endif

       enddo
   enddo

   print *,'tmp= done'

   ! CAll standard analysis
    !call analysis(mem,h,E,transpose(tmp),3*kdim,nrsamp,obstot,.false.)
    !return


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! The observations have uncorrelated data errors
   allocate(P(obstot,obstot))
   P=0.0
   do m=1,obstot
      P(m,m)=d(m)%var
!      print '(4f9.4)', P(m,:)
   enddo

   call analysis2(mem,h,P,transpose(tmp),3*kdim,nrsamp,obstot,.false.)
   return

   !call analysis4(mem,P,transpose(tmp),d_vec,3*kdim,nrsamp,obstot,.false.)
   !return
  
   !call analysis4c(mem,P,transpose(tmp),d_vec,3*kdim,nrsamp,obstot,.false.)
   !return

  ! call analysis5(mem,P,transpose(tmp),d_vec,3*kdim,nrsamp,obstot,.false.)
  ! return
 
    !call analysis5c(mem,P,transpose(tmp),d_vec,3*kdim,nrsamp,obstot,.false.)
    !return

   !call analysis6(mem,E,transpose(tmp),d_vec,3*kdim,nrsamp,obstot,.false.)
   !return

   !call analysis6c(mem,E,transpose(tmp),d_vec,3*kdim,nrsamp,obstot,.false.)
   !return
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Construct the P matrix "P + scale*H*P^f*H"
! defined as:  P=scale*matmul(tmp,transpose(tmp))+P
! Here a rewriting is used to save compiler induced memory allocation
! of temporary variables.  Note that tmp has already been transposed
! to access memory with stride one.  The alternative is depressing!

   do j=1,obstot
   do i=j,obstot
      P(i,j)=scale*dot_product(tmp(:,i),tmp(:,j))+P(i,j)
    
      P(j,i)  = P(i,j)
   enddo
   enddo
   print *,'P ok'
   print '(a,i8)','Upper left of P'
 !  print '(10g11.2)',P(1:min(10,obstot),1:min(10,obstot))
   print '(4g11.2)',P(1:4,1:4)
   print *


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Calculate the eigendecomposition of P using eispack. P = ZZ*ss*ZZ

   allocate(ZZ(obstot,obstot))
   allocate(w(8*obstot))
   allocate(ss(obstot))
   allocate(iw(5*obstot))
   allocate(ipvt(obstot))
   allocate(lambda(obstot))

   ZZ=0.0; w=0.0; ss=0.0; iw=0; ipvt=0; lambda=0.0



! Eigen value factorization using rsm:
!  call rsm(obstot,obstot,P,ss,obstot,ZZ,w,ipvt,info)


   allocate (ap(obstot*(obstot+1)/2) )
   k=0
   do j=1,obstot
   do i=1,j
      k=k+1
    ap(k)=P(i,j) 
    enddo
    enddo



   call dspev(21,ap,ss,ZZ,obstot,obstot,w,2*obstot)
   deallocate(ap)



   print *,'Eigenvalues of P:'
   print '(20g10.2)',ss(obstot:1:-1)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Write the eigenvalues to a file eig.dat
   summ=sum(ss)
   print *,'sum(eigenvalues)=',summ

   open(10,file='eig.dat')
      do i=obstot,1,-1
!         write(10,'(i4,g13.5)')i,ss(i)/summ
         write(10,'(i4,g13.5)')obstot-i+1,ss(i)/summ
      enddo
   close(10)
   print *,'eig.dat done'


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Find the number of significant eigenvalues 
!   To avoid a singular invers of P=ZZ*ss*ZZ and a noisy solution, we
!   truncate the spectrum of P. lambda contains the truncated inverse of
!   the matrix diag(ss)

   print *,'computing the number of significant eigen values'
   summ2=0.0
   do i=obstot,1,-1
      summ2=summ2+ss(i)
      if (summ2/summ>0.99) then
         nrsigma=obstot-i+1
         exit
      endif
   enddo
   lambda=0.0
   do k=obstot,nrsigma,-1
      lambda(k)=1.0/ss(k)
   enddo
   print *,'Number of significant eigenvalues are: ',nrsigma


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! The analysis increment is given as 
!     scale*(A-Mean)*[H*(A-Mean)]^T*{ZZ*lambda*ZZ^T*h}
! which also can be written as
!     scale*(A-Mean)*tmp^T*{ZZ*lambda*ZZ^T*h}
! or
!     scale*(A-Mean)*L

! ZZ^T*h  (rotated measurement residuals, or residuals represented in the ZZ space)
   h=matmul(transpose(ZZ),h)
   print *,'made ZZ^T*h'

! lambda*(ZZ^T*h)
   do j=1,nrsamp
      h(:,j)=lambda(:)*h(:,j)
   enddo
   print *,'made lambda*(ZZ^T*h)='
   print '(15E10.2)',h(:,2)

! ZZ*(lambda*(ZZ^T*h))
   h=matmul(ZZ,h)
   print *,'made ZZ*(lambda*(ZZ^T*h))'

! tmp^T*(ZZ*lambda*ZZ^T*h)            ! L: Not in use now!
   allocate(L(nrsamp,nrsamp))
   L=matmul(tmp,h)

   print *,'made coefficients L=tmp^T*(ZZ*lambda*ZZ^T*h)'


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Final analysis (A-Mean)*L. 
! ndim = nx*ny*nz*(number of 3D fields)+nx*ny*(number of 2D fields)
   nndim=kdim*antvar

   allocate(tmpens(nrsamp));  print *,'allocated tmpens'

!$OMP PARALLEL DO PRIVATE(iens)
   do iens=1,nrsamp
      tmpens(iens)=mem(iens)-Mean
   enddo
   

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Analysis using representer approach

! - compute represeters, ``rtd = scale*tmpens*tmp''.
   allocate (rtd(obstot))

   do i = 1, obstot
      rtd(i)=0.0 ! Initialize representers
   enddo
   print *,'Computing the representer matrix rtd ...'
   call dgemm('n','n',nndim,obstot,nrsamp,scale,tmpens,nndim,tmp,nrsamp,1.0,rtd,nndim)
   print *,'... done representers'

   deallocate (tmpens)
   deallocate (tmp)

! - compute analysis ``A = rtd*h + A''.

   print *,'Computing the new analysis ...'
   call dgemm('n','n',nndim,nrsamp,obstot,1.0,rtd,nndim,h,obstot,1.0,mem,nndim)
   print *,'... done'
   
   deallocate (rtd)

   do iens=1,nrsamp
     mem(iens)%N=max(0.0,mem(iens)%N) 
     mem(iens)%P=max(0.0,mem(iens)%P)
     mem(iens)%H=max(0.0,mem(iens)%H)
   enddo

end subroutine prepare_and_analysis

end module m_analysis