Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

GSI Issue #233: Fix indexing issues in correlated error #236

Merged
merged 3 commits into from
Oct 18, 2021
Merged

GSI Issue #233: Fix indexing issues in correlated error #236

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
53f781a
updates to fix indexing issue
KristenBathmann Oct 15, 2021
33e46e6
cleaned up and bug fix
KristenBathmann Oct 15, 2021
9a41e68
cleaned up cov_calc.f90
KristenBathmann Oct 15, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
104 changes: 56 additions & 48 deletions src/gsi/correlated_obsmod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -333,7 +333,7 @@ end subroutine ini_
! !INTERFACE:
!
subroutine set_(instrument,fname,mask,method,kreq,kmut,ErrorCov)
use radinfo, only: nusis,iuse_rad,jpch_rad,varch
use radinfo, only: nusis,iuse_rad,jpch_rad,varch,nuchan
use constants, only: zero
implicit none

Expand Down Expand Up @@ -371,8 +371,8 @@ subroutine set_(instrument,fname,mask,method,kreq,kmut,ErrorCov)
character(len=*),parameter :: myname_=myname//'*set'
integer(i_kind) nch_active !number of channels accounted for in the covariance file
integer(i_kind) nctot !the total number of channels (active+passive), according to the covariance file
integer(i_kind) lu,ii,jj,ioflag,iprec,coun,couns,istart,indR
integer(i_kind),dimension(:),allocatable:: indxR !channel indices read in from the covariance file
integer(i_kind) lu,ii,jj,ioflag,iprec,coun,couns,istart,indR,nctotf
integer(i_kind),dimension(:),allocatable:: indxR,indxRf !channel indices read in from the covariance file
real(r_kind),dimension(:,:),allocatable:: Rcov !covariance matrix read in from the covariance file
real(r_single),allocatable, dimension(:,:) :: readR4 ! nch_active x nch_active
real(r_double),allocatable, dimension(:,:) :: readR8 ! nch_active x nch_active
Expand Down Expand Up @@ -400,6 +400,7 @@ subroutine set_(instrument,fname,mask,method,kreq,kmut,ErrorCov)
coun=0
couns=0
istart=0
nctotf=0
do ii=1,jpch_rad
if (nusis(ii)==ErrorCov%instrument) then
if (couns==0) then
Expand All @@ -409,6 +410,7 @@ subroutine set_(instrument,fname,mask,method,kreq,kmut,ErrorCov)
if (iuse_rad(ii)>0) then
coun=coun+1
endif
nctotf=nctotf+1
endif
enddo
! if no data available, turn off Correlated Error
Expand All @@ -420,7 +422,7 @@ subroutine set_(instrument,fname,mask,method,kreq,kmut,ErrorCov)
ErrorCov%nch_active = coun
if (.not.GMAO_ObsErrorCov) ErrorCov%nctot = nctot
call create_(coun,ErrorCov)
allocate(indxR(nch_active),Rcov(nctot,nctot))
allocate(indxRf(nch_active),indxR(nch_active),Rcov(nctot,nctot))

! Read GSI-like channel numbers used in estimating R for this instrument
read(lu,IOSTAT=ioflag) indxR
Expand All @@ -442,54 +444,60 @@ subroutine set_(instrument,fname,mask,method,kreq,kmut,ErrorCov)
Rcov(1:nch_active,1:nch_active)=readR8
deallocate(readR8)
endif
coun=0
ErrorCov%R=zero
do ii=1,nctot
if (iuse_rad(ii+istart)>0) then
coun=coun+1
ErrorCov%indxR(coun)=ii
endif
enddo
if (GMAO_ObsErrorCov) then
ErrorCov%indxR(1:nch_active)=indxR(1:nch_active)
ErrorCov%nch_active=nch_active
else
coun=0
ErrorCov%R=zero
do ii=1,nctotf
if (iuse_rad(ii+istart)>0) then
coun=coun+1
ErrorCov%indxR(coun)=ii
indxRf(coun)=nuchan(ii+istart)
endif
enddo
!Add rows and columns for active channels in the satinfo that are not in the covariance file
couns=1
do ii=1,ErrorCov%nch_active
indR=0
do jj=couns,nch_active
if (indxR(jj)==ErrorCov%indxR(ii)) then
indR=jj
couns=jj+1
exit
couns=1
do ii=1,ErrorCov%nch_active
indR=0
do jj=couns,nch_active
if (indxR(jj)==indxRf(ii)) then
indR=jj
couns=jj+1
exit
endif
enddo
if (indR==0) then
do jj=nctot-1,ii,-1
Rcov(jj+1,:)=Rcov(jj,:)
Rcov(:,jj+1)=Rcov(:,jj)
enddo
Rcov(ii,:)=zero
Rcov(:,ii)=zero
Rcov(ii,ii)=varch(istart+ErrorCov%indxR(ii))*varch(istart+ErrorCov%indxR(ii))
endif
enddo
if (indR==0) then
do jj=nctot-1,ii,-1
Rcov(jj+1,:)=Rcov(jj,:)
Rcov(:,jj+1)=Rcov(:,jj)
enddo
Rcov(ii,:)=zero
Rcov(:,ii)=zero
Rcov(ii,ii)=varch(istart+ErrorCov%indxR(ii))*varch(istart+ErrorCov%indxR(ii))
endif
enddo
!Remove rows and columns that are in the covariance file, but not in the satinfo
couns=1
do ii=1,nch_active
indR=0
do jj=couns,ErrorCov%nch_active
if (ErrorCov%indxR(jj)==indxR(ii)) then
indR=jj
couns=jj+1
exit
endif
enddo
if (indR==0) then
do jj=ii,nctot-1
Rcov(jj,:)=Rcov(jj+1,:)
Rcov(:,jj)=Rcov(:,jj+1)
couns=1
do ii=1,nch_active
indR=0
do jj=couns,ErrorCov%nch_active
if (indxRf(jj)==indxR(ii)) then
indR=jj
couns=jj+1
exit
endif
enddo
endif
enddo
ErrorCov%R(1:ErrorCov%nch_active,1:ErrorCov%nch_active)=Rcov(1:ErrorCov%nch_active,1:ErrorCov%nch_active)
if (indR==0) then
do jj=ii,nctot-1
Rcov(jj,:)=Rcov(jj+1,:)
Rcov(:,jj)=Rcov(:,jj+1)
enddo
endif
enddo
endif
ErrorCov%R(1:ErrorCov%nch_active,1:ErrorCov%nch_active)=Rcov(1:ErrorCov%nch_active,1:ErrorCov%nch_active)
! Done reading file
close(lu)
else
Expand Down Expand Up @@ -531,7 +539,7 @@ subroutine set_(instrument,fname,mask,method,kreq,kmut,ErrorCov)
endif
deallocate(diag)
endif
deallocate(indxR,Rcov)
deallocate(indxR,Rcov,indxRf)
initialized_=.true.
end subroutine set_
!EOC
Expand Down
6 changes: 2 additions & 4 deletions util/Correlated_Obs/cov_calc.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -332,12 +332,11 @@ program cov_calc
deallocate(Rcovbig,divbig,ges_avebig1,ges_avebig2)
deallocate(n_pair_hl, obs_pairs_hl)
end if

!output
reclen=kind(Rcov(1,1))
open(26,file=trim(cov_file),form='unformatted')
write(26) nch_active, nctot, reclen
write(26) indR
write(26) indRf
write(26) Rcov
close(26)

Expand All @@ -356,9 +355,8 @@ program cov_calc
write(25,rec=1) Rcorr
close(25)
end if

deallocate(Rcov,chaninfo,errout)
deallocate(indR)
deallocate(indR,indRf)
deallocate(divider)
if (out_corr) then
deallocate(Rcorr)
Expand Down
18 changes: 11 additions & 7 deletions util/Correlated_Obs/readsatobs.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,11 +11,11 @@ module readsatobs
implicit none

public :: get_satobs_data, get_chaninfo
public :: indR, chaninfo, errout
public :: indR, indRf,chaninfo, errout
public :: nch_active,nctot
public :: RadData

integer(i_kind),dimension(:),allocatable:: indR !indices of the assimlated channels
integer(i_kind),dimension(:),allocatable:: indR,indRf !indices of the assimlated channels
real(r_kind),dimension(:),allocatable:: chaninfo !wavenumbers of assimilated channels
real(r_kind),dimension(:),allocatable:: errout !satinfo obs errors of assimilated channels
integer(i_kind):: nch_active !number of actively assimilated channels
Expand Down Expand Up @@ -86,7 +86,7 @@ subroutine get_chaninfo_bin(filename,chan_choice)
nch_active=nch_active+1
end do
endif
allocate(indR(nch_active),chaninfo(nch_active),errout(nch_active))
allocate(indRf(nch_active),indR(nch_active),chaninfo(nch_active),errout(nch_active))
i=0
do n=1,header_fix0%nchan
if (chan_choice==full_chan) then
Expand All @@ -96,6 +96,7 @@ subroutine get_chaninfo_bin(filename,chan_choice)
indR(i)=n
endif
end do
indRf=indR
do n=1,nch_active
chaninfo(n)=header_chan0(indR(n))%wave
errout(n)=header_chan0(indR(n))%varch
Expand All @@ -112,7 +113,7 @@ subroutine get_chaninfo_nc(filename,chan_choice)
character(len=9), intent(in) :: filename
integer(i_kind), intent(in):: chan_choice
integer(i_kind) iunit, nobs, i,j
integer(i_kind), dimension(:), allocatable ::Use_Flag,chind
integer(i_kind), dimension(:), allocatable ::Use_Flag,chind,schind
real(r_double), dimension(:), allocatable:: Waves,Inv_Errors
real(r_kind), dimension(:), allocatable:: Wave,Inv_Error

Expand All @@ -121,12 +122,13 @@ subroutine get_chaninfo_nc(filename,chan_choice)
nobs = nc_diag_read_get_dim(iunit,'nobs')
if (nobs <= 0) call nc_diag_read_close(filename)
nctot = nc_diag_read_get_dim(iunit,'nchans')
allocate(Use_Flag(nctot),chind(nobs),Wave(nctot),Inv_Error(nctot))
allocate(Use_Flag(nctot),chind(nobs),schind(nctot),Wave(nctot),Inv_Error(nctot))
allocate(Waves(nctot),Inv_Errors(nctot))
call nc_diag_read_get_var(iunit, 'use_flag', Use_Flag)
call nc_diag_read_get_var(iunit, 'Channel_Index', chind)
call nc_diag_read_get_var(iunit, 'wavenumber', Waves)
call nc_diag_read_get_var(iunit, 'error_variance', Inv_Errors)
call nc_diag_read_get_var(iunit, 'sensor_chan', schind)
Wave=real(Waves,r_kind)
Inv_Error=real(Inv_Errors,r_kind)
call nc_diag_read_close(filename)
Expand All @@ -139,21 +141,23 @@ subroutine get_chaninfo_nc(filename,chan_choice)
nch_active=nch_active+1
enddo
endif
allocate(indR(nch_active),chaninfo(nch_active),errout(nch_active))
allocate(indRf(nch_active),indR(nch_active),chaninfo(nch_active),errout(nch_active))
i=0
do j=1,nctot
if (chan_choice==full_chan) then
indRf(j)=schind(j)
indR(j)=j
else if (Use_Flag(chind(j))>0) then
i=i+1
indRf(i)=schind(j)
indR(i)=j
end if
end do
do j=1,nch_active
chaninfo(j)=Wave(chind(indR(j)))
errout(j)=Inv_Error(chind(indR(j)))
end do
deallocate(Use_flag,chind,Wave,Inv_Error,Waves,Inv_Errors)
deallocate(Use_flag,chind,schind,Wave,Inv_Error,Waves,Inv_Errors)
end subroutine get_chaninfo_nc

! read radiance data
Expand Down