Skip to content

Commit

Permalink
Jeff's bugix to avoid the mystery crash.
Browse files Browse the repository at this point in the history
  • Loading branch information
ClaraDraper-NOAA committed Aug 21, 2023
1 parent 49e7aec commit 4c3a86e
Showing 1 changed file with 11 additions and 10 deletions.
21 changes: 11 additions & 10 deletions src/enkf/letkf.f90
Original file line number Diff line number Diff line change
Expand Up @@ -470,6 +470,7 @@ subroutine letkf_update()
! is used as workspace and is modified on output), and analysis
! weights for ensemble perturbations represent posterior ens perturbations, not
! analysis increments for ensemble perturbations.

call letkf_core(nobsl2,hxens,obens,dep,&
wts_ensmean,wts_ensperts,pa,&
rdiag,rloc(1:nobsl2),nens,nens/nanals,getkf_inflation,denkf,getkf)
Expand Down Expand Up @@ -729,9 +730,9 @@ subroutine letkf_core(nobsl,hxens,hxens_orig,dep,&
allocate(work3(nanals,nanals),evecs(nanals,nanals))
allocate(rrloc(nobsl),gammapI(nanals),evals(nanals),gamma_inv(nanals))
! for dsyevr
allocate(iwork(10*nanals),work1(70*nanals))
!allocate(iwork(10*nanals),work1(70*nanals))
! for dsyevd
!allocate(iwork(3+5*nanals),work1(1+6*nanals+2*nanals*nanals))
allocate(iwork(3+5*nanals),work1(1+6*nanals+2*nanals*nanals))

! HZ^T = hxens sqrt(Rinv)
rrloc = rdiaginv * rloc
Expand All @@ -754,19 +755,19 @@ subroutine letkf_core(nobsl,hxens,hxens_orig,dep,&
hxens,nanals,0.d0,work3,nanals)
! evecs contains eigenvectors of HZ^T HZ, or left singular vectors of HZ
! evals contains eigenvalues (singular values squared)
call dsyevr('V','A','L',nanals,work3,nanals,vl,vu,1,nanals,-1.d0,nanals,evals,evecs, &
nanals,isuppz,work1,lwork,iwork,liwork,ierr)
!call dsyevr('V','A','L',nanals,work3,nanals,vl,vu,1,nanals,-1.d0,nanals,evals,evecs, &
! nanals,isuppz,work1,lwork,iwork,liwork,ierr)
! use LAPACK dsyevd instead of dsyevr
!evecs = work3
!call dsyevd('V','L',nanals,evecs,nanals,evals,work1,lwork,iwork,liwork,ierr)
evecs = work3
call dsyevd('V','L',nanals,evecs,nanals,evals,work1,lwork,iwork,liwork,ierr)
else ! single precision
call sgemm('n','t',nanals,nanals,nobsl,1.e0,hxens,nanals, &
hxens,nanals,0.e0,work3,nanals)
call ssyevr('V','A','L',nanals,work3,nanals,vl,vu,1,nanals,-1.e0,nanals,evals,evecs, &
nanals,isuppz,work1,lwork,iwork,liwork,ierr)
!call ssyevr('V','A','L',nanals,work3,nanals,vl,vu,1,nanals,-1.e0,nanals,evals,evecs, &
! nanals,isuppz,work1,lwork,iwork,liwork,ierr)
! use LAPACK dsyevd instead of dsyevr
!evecs = work3
!call ssyevd('V','L',nanals,evecs,nanals,evals,work1,lwork,iwork,liwork,ierr)
evecs = work3
call ssyevd('V','L',nanals,evecs,nanals,evals,work1,lwork,iwork,liwork,ierr)
end if
if (ierr .ne. 0) print *,'warning: dsyev* failed, ierr=',ierr
deallocate(work1,iwork,work3) ! no longer needed
Expand Down

0 comments on commit 4c3a86e

Please sign in to comment.