subroutine analysis2(A, D, R, S, ndim, nrens, nrobs, verbose) ! Computes the analysed ensemble for A 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) :: D(nrobs,nrens) ! matrix holding observation innovations real, intent(in) :: S(nrobs,nrens) ! matrix holding HA' real, intent(inout) :: R(nrobs,nrobs) ! Error covariance matrix for observations logical, intent(in) :: verbose real, allocatable, dimension(:,:) :: X1,X2,U,X4,Reps,V real X3(nrobs,nrens) real, allocatable, dimension(:) :: sig,work real sigsum,sigsum1,oneobs(1,1) integer ierr,nrsigma,i,j,lwork, m integer iblkmax character(len=2) tag2 if (nrobs > 1) then ! R=float(nrens)*R+matmul(S,transpose(S)) call dgemm('n','t', nrobs, nrobs, nrens, & 1.0, S, nrobs, & S, nrobs, & float(nrens), R, nrobs) allocate (U(nrobs,nrobs) ) allocate (V(nrobs,nrobs) ) allocate (sig(nrobs) ) lwork=2*max(3*nrobs+nrobs,5*nrobs) allocate(work(lwork)) sig=0.0 call dgesvd('A', 'A', nrobs, nrobs, R, nrobs, sig, U, nrobs, V, nrobs, work, lwork, ierr) deallocate(work) if (ierr /= 0) then print *,'ierr from call dgesvd= ',ierr stop endif open(10,file='sigma.dat') do i=1,nrobs write(10,'(i5,g12.3)')i,sig(i) enddo close(10) sigsum=sum( sig(1:nrobs) ) sigsum1=0.0 ! Significant eigenvalues. nrsigma=0 do i=1,nrobs ! singular values are in descending order if (sigsum1/sigsum < 0.999) then nrsigma=nrsigma+1 sigsum1=sigsum1+sig(i) sig(i) = 1.0/sig(i) else sig(i:nrobs)=0.0 exit endif enddo if (verbose) then write(*,'(a,i5,g13.5)') ' dominant sing. values and share ',nrsigma,sigsum1/sigsum write(*,'(8g12.3)')1./sig(1:nrsigma), sig(nrsigma+1:nrobs) endif allocate (X1(nrobs,nrobs)) do i=1,nrobs do j=1,nrobs X1(i,j)=sig(i)*U(j,i) enddo enddo deallocate(sig) allocate (X2(nrobs,nrens)) ! X2=matmul(X1,D) call dgemm('n','n',nrobs,nrens,nrobs,1.0,X1,nrobs,D ,nrobs,0.0,X2,nrobs) deallocate(X1) ! X3=matmul(V,X2) call dgemm('t','n',nrobs,nrens,nrobs,1.0,V ,nrobs,X2,nrobs,0.0,X3,nrobs) deallocate(V) deallocate(X2) else oneobs=matmul(S,transpose(S))+R print *,'oneobs: ',oneobs(1,1) X3=D/oneobs(1,1) endif if (2_8*ndim*nrobs < 1_8*nrens*(nrobs+ndim)) then ! Code for few observations ( m