analysis6.F90 5.81 KB
``````subroutine analysis6(A, E, S, d, ndim, nrens, nrobs, verbose)
! Computes the analysed ensemble for A using the square root formulation
! in algorithm A.2 from Evensen 2004.  I.e. the sophisitcated inversion
! algorithm is used, which appears to be more stable with large m.
! This algorith does not use the full R but assumes (nrens-1) R = E E^T and E is
! specified.

use m_multa
implicit none
integer, intent(in) :: ndim             ! dimension of model state
integer, intent(in) :: nrens            ! number of ensemble members
integer, intent(in) :: nrobs            ! number of observations

real, intent(inout) :: A(ndim,nrens)    ! ensemble matrix
real, intent(in)    :: S(nrobs,nrens)   ! matrix holding HA'
real, intent(in)    :: d(nrobs)         ! vector holding d-HA
real, intent(in)    :: E(nrobs,nrens)   ! matrix holding observation perturbations
logical, intent(in) :: verbose

real ave(ndim)       ! ensemble mean

real S0(nrobs,nrens)

real, allocatable :: U0(:,:),sig0(:),VT0(:,:)
real, allocatable :: U1(:,:),sig1(:),VT1(:,:)
real, allocatable :: U2(:,:),sig2(:),VT2(:,:)
real, allocatable :: X0(:,:),X1(:,:),X2(:,:)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real, allocatable, dimension(:)   :: work,isigma

real X3(nrens,nrens)
real, allocatable :: y1(:)
real, allocatable :: y2(:)
real y3(nrobs)
real y4(nrens)

real sigsum,sigsum1
integer ierr,nrsigma,i,j,lwork, nrmin
integer iblkmax

nrmin=min(nrobs,nrens)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Subtract mean from ensemble
ave(:)=A(:,1)
do i=2,nrens
ave(:)=ave(:)+A(:,i)
enddo
ave=(1.0/real(nrens))*ave
do i=1,nrens
A(:,i)=A(:,i)-ave(:)
enddo

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute SVD of S=HA'  ->  U0, sig0
S0=S
allocate (U0(nrobs,nrmin)  )
allocate (sig0(nrmin))
allocate (VT0(1,1))
lwork=2*max(3*nrens+nrobs,5*nrens)
allocate(work(lwork))
sig0=0.0
!\$OMP CRITICAL
call dgesvd('S', 'N', nrobs, nrens, S0, nrobs, sig0, U0, nrobs, VT0, nrens, work, lwork, ierr)
!\$OMP END CRITICAL
deallocate(work,VT0)
if (ierr /= 0) then
print *,'analysis6: ierr from call dgesvd 0= ',ierr; stop
endif

sigsum=sum( sig0(1:nrmin) )
sigsum1=0.0
! Significant eigenvalues.
nrsigma=0
do i=1,nrmin
if (sigsum1/sigsum < 0.9999) then
nrsigma=nrsigma+1
sigsum1=sigsum1+sig0(i)
else
sig0(i:nrmin)=0.0
exit
endif
enddo

if (verbose) then
write(*,'(a,i5,g13.5)') ' dominant sing. values and share ',nrsigma,sigsum1/sigsum
write(*,'(5g11.3)')sig0
endif

do i=1,nrsigma
sig0(i) = 1.0/sig0(i)
enddo

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute X0=sig0^{*T} U0^T E

! X0= U0^T R:q

allocate(X0(nrmin,nrens))
call dgemm('t', 'n', nrmin, nrens, nrobs, 1.0d0, U0, nrobs, E, nrobs, 0.0d0, X0, nrmin)
!X0=matmul(transpose(U0),E)

do j=1,nrens
do i=1,nrmin
X0(i,j)=sig0(i)*X0(i,j)
enddo
enddo

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute singular value decomposition  of X0(nrmin,nrens)
allocate (U1(nrmin,nrmin)  )
allocate (sig1(nrmin))
allocate (VT1(1,1))
allocate (VT0(1,1))
lwork=2*max(3*nrens+nrobs,5*nrens)
allocate(work(lwork))
sig1=0.0
!\$OMP CRITICAL
call dgesvd('S', 'N', nrmin, nrens, X0, nrmin, sig1, U1, nrmin, VT0, nrens, work, lwork, ierr)
!\$OMP END CRITICAL
deallocate(work,VT1)
if (ierr /= 0) then
print *,'analysis6: ierr from call dgesvd 1= ',ierr; stop
endif

do i=1,nrmin
sig1(i)=1.0/(1.0+sig1(i)**2)
enddo

deallocate (X0,VT0)
!print *,'sig1:',nrmin
!print '(5g11.3)',sig1

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! X1 = U0 * sig0^{-1} * U1
do j=1,nrmin
do i=1,nrmin
U1(i,j)=sig0(i)*U1(i,j)
enddo
enddo

allocate(X1(nrobs,nrmin))
call dgemm('n','n',nrobs,nrmin,nrmin, 1.0,U0,nrobs, U1,nrmin, 0.0,X1,nrobs)
! X1=matmul(U0,U1)
deallocate (U0,sig0,U1)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Update mean
allocate(y1(1:nrmin))
call dgemv('t',nrobs,nrmin,1.0,X1,nrobs,d,1,0.0,y1 ,1)
allocate(y2(1:nrmin))
y2=sig1*y1
call dgemv('n',nrobs,nrmin,1.0,X1,nrobs,y2,1,0.0,y3 ,1)
call dgemv('t',nrobs,nrens,1.0,S ,nrobs,y3,1,0.0,y4 ,1)
call dgemv('n',ndim ,nrens,1.0,A ,ndim ,y4,1,1.0,ave,1)
deallocate(y1,y2)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! X2 = (I+sig1^2)^{-0.5} * X1^T * S
allocate(X2(nrmin,nrens))
call dgemm('t','n',nrmin,nrens,nrobs,1.0,X1,nrobs, S,nrobs, 0.0,X2,nrmin)
!X2=matmul(transpose(X1),S)

do j=1,nrens
do i=1,nrmin
X2(i,j)=sqrt(sig1(i))*X2(i,j)
enddo
enddo

deallocate (sig1,X1)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!
! SVD of X2
lwork=2*max(3*nrens+nrens,5*nrens)
allocate (U2(nrmin,nrmin), sig2(nrmin), VT2(nrens,nrens), work(lwork))
sig2=0.0
!\$OMP CRITICAL
call dgesvd('N', 'A', nrmin, nrens, X2, nrmin, sig2, U2, nrmin, VT2, nrens, work, lwork, ierr)
!\$OMP END CRITICAL
deallocate(work)
if (ierr /= 0) then
print *,'ierr from call dgesvd 2 = ',ierr
stop
endif

allocate(isigma(nrmin))
isigma=1.0
do i=1,nrmin
if ( sig2(i) > 1.0 ) print *,'WARNING (analysis 6): sig2 > 1',i,sig2(i)
isigma(i)=sqrt( max(1.0-sig2(i)**2,0.0) )
enddo

do j=1,nrens
X3(:,j)=VT2(j,:)
enddo

do j=1,nrmin
X3(:,j)=X3(:,j)*isigma(j)
enddo

!print '(a)','A5: sig2: '
!print '(5g11.3)',sig2(1:nrmin)

! Final ensemble perturbations
iblkmax=min(ndim,200)
call  multa(A, X3, ndim, nrens, iblkmax )

do i=1,nrens
A(:,i)=A(:,i)+ave(:)
enddo

deallocate (U2,sig2,VT2,isigma)

end subroutine analysis6

``````