diff --git a/src/gsi/correlated_obsmod.F90 b/src/gsi/correlated_obsmod.F90 index b477c9ef80..94834bde1d 100644 --- a/src/gsi/correlated_obsmod.F90 +++ b/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 diff --git a/util/Correlated_Obs/cov_calc.f90 b/util/Correlated_Obs/cov_calc.f90 index 87b741a928..38f8ca86f8 100644 --- a/util/Correlated_Obs/cov_calc.f90 +++ b/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) diff --git a/util/Correlated_Obs/readsatobs.f90 b/util/Correlated_Obs/readsatobs.f90 index 55ba064a95..2bf67942ef 100644 --- a/util/Correlated_Obs/readsatobs.f90 +++ b/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,13 +141,15 @@ 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 @@ -153,7 +157,7 @@ subroutine get_chaninfo_nc(filename,chan_choice) 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