Merge pull request #236 from KristenBathmann/release/gfsda.v16.x

GSI Issue #233: Fix indexing issues in correlated error

src/gsi/correlated_obsmod.F90

@@ -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
 
@@ -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 
@@ -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 
@@ -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
@@ -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
@@ -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
@@ -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

util/Correlated_Obs/cov_calc.f90

@@ -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)
 
@@ -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)

util/Correlated_Obs/readsatobs.f90

@@ -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
@@ -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
@@ -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
@@ -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
 
@@ -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)
@@ -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