From 9de5986aa0482beb9e9491078aa8a0d5f21c2d1b Mon Sep 17 00:00:00 2001 From: "andrew.eichmann" Date: Wed, 1 Dec 2021 18:22:56 +0000 Subject: [PATCH] GitHub Issue NOAA-EMC/GSI#259. Add assimilated flag to observations in osense file. --- src/enkf/enkf.f90 | 21 +++++++++++++++++++-- src/enkf/enkf_obs_sensitivity.f90 | 8 +++++++- src/enkf/enkf_obsmod.f90 | 4 ++++ 3 files changed, 30 insertions(+), 3 deletions(-) diff --git a/src/enkf/enkf.f90 b/src/enkf/enkf.f90 index adf9921e9d..c117e4ba56 100644 --- a/src/enkf/enkf.f90 +++ b/src/enkf/enkf.f90 @@ -123,7 +123,7 @@ module enkf obtype, oberrvarmean, numobspersat, deltapredx, biaspreds,& oberrvar_orig, probgrosserr, prpgerr,& corrlengthsq,lnsigl,obtimel,obloclat,obloclon,obpress,stattype,& - anal_ob,anal_ob_post + anal_ob,anal_ob_post,assimltd_flag use constants, only: pi, one, zero use params, only: sprd_tol, paoverpb_thresh, datapath, nanals,& iassim_order,sortinc,deterministic,numiter,nlevs,& @@ -153,7 +153,8 @@ subroutine enkf_update() ! local variables. integer(i_kind) nob,nob1,nob2,nob3,npob,nf,nf2,ii,nobx,nskip,& - niter,i,nrej,npt,nuse,ncount,ncount_check,nb,np + niter,i,nrej,npt,nuse,ncount,ncount_check,nb,np,& + nuseconvoz,nusesat,nobs_convoz integer(i_kind) indxens1(nanals),indxens2(nanals) integer(i_kind) indxens1_modens(nanals*neigv),indxens2_modens(nanals*neigv) real(r_single) hxpost(nanals),hxprior(nanals),hxinc(nanals),& @@ -757,17 +758,31 @@ subroutine enkf_update() tend = mpi_wtime() if (nproc .eq. 0) then write(6,8003) niter,'timing on proc',nproc,' = ',tend-tbegin,t2,t3,t4,t5,t6,nrej + allocate(assimltd_flag(nobstot)) + assimltd_flag = 99999 if (iassim_order == 2) then ncount_check = ncount else ncount_check = nobstot endif nuse = 0; covl_fact = 0. + nuseconvoz=0; nusesat = 0 + nobs_convoz = nobs_conv + nobs_oz do nob1=1,ncount_check nob = indxassim(nob1) if (iskip(nob) .ne. 1) then covl_fact = covl_fact + sqrt(corrlengthsq(nob)/corrlengthsq_orig(nob)) nuse = nuse + 1 + assimltd_flag(nob) = 1 + if (nob .le. nobs_convoz) then + nuseconvoz = nuseconvoz +1 + else if (nob .gt. nobs_convoz) then + nusesat = nusesat + 1 + else + print *,'nob ', nob ,' falling through' + endif + else + assimltd_flag(nob) = 0 endif enddo nskip = nobstot-nuse @@ -776,6 +791,8 @@ subroutine enkf_update() if (covl_fact < 0.99) print *,'mean covl_fact = ',covl_fact if (nskip > 0) print *,nskip,' out of',nobstot,'obs skipped,',nuse,' used' if (nsame > 0) print *,nsame,' out of', nobstot-nskip,' same lat/long' + if (nuseconvoz > 0 ) print *,nuseconvoz,' out of',nobs_conv + nobs_oz ,'convobs used' + if (nusesat > 0 ) print *,nusesat ,' out of',nobs_sat ,'satobs used' if (nrej > 0) print *,nrej,' obs rejected by varqc' endif 8003 format(i2,1x,a14,1x,i5,1x,a3,6(f7.2,1x),i4) diff --git a/src/enkf/enkf_obs_sensitivity.f90 b/src/enkf/enkf_obs_sensitivity.f90 index afa2e242e8..6c37936f31 100644 --- a/src/enkf/enkf_obs_sensitivity.f90 +++ b/src/enkf/enkf_obs_sensitivity.f90 @@ -45,7 +45,7 @@ module enkf_obs_sensitivity obloclon,obpress,indxsat,oberrvar,stattype,obtime,ob, & ensmean_ob,ensmean_obnobc,obsprd_prior,obfit_prior, & oberrvar_orig,biaspreds,anal_ob_post,nobstot,lnsigl, & - corrlengthsq,obtimel,oblnp,obloc + corrlengthsq,obtimel,oblnp,obloc ,assimltd_flag use convinfo, only: convinfo_read,init_convinfo use ozinfo, only: ozinfo_read,init_oz use radinfo, only: radinfo_read,jpch_rad,nusis,nuchan,npred @@ -91,6 +91,7 @@ module enkf_obs_sensitivity integer(i_kind) :: stattype ! Observation type character(len=20) :: obtype ! Observation element / Satellite name integer(i_kind) :: indxsat ! Satellite index (channel) + integer(i_kind) :: assimltd_flag ! is assimilated? flag real(r_single) :: osense_kin ! Observation sensitivity (kinetic energy) [J/kg] real(r_single) :: osense_dry ! Observation sensitivity (Dry total energy) [J/kg] real(r_single) :: osense_moist ! Observation sensitivity (Moist total energy) [J/kg] @@ -215,6 +216,7 @@ subroutine read_ob_sens allocate(stattype(nobstot)) allocate(obtype(nobstot)) allocate(indxsat(nobs_sat)) + allocate(assimltd_flag(nobstot)) allocate(biaspreds(npred+1,nobs_sat)) allocate(tmpanal_ob(nanals)) allocate(tmpbiaspreds(npred+1)) @@ -235,6 +237,7 @@ subroutine read_ob_sens oberrvar_orig(nob) = real(indata%oberrvar_orig,r_kind) stattype(nob) = indata%stattype obtype(nob) = indata%obtype + assimltd_flag(nob) = indata%assimltd_flag if(nproc == 0) anal_ob_post(1:nanals,nob) = real(tmpanal_ob(1:nanals),r_kind) end do ! Read loop over satellite radiance observations @@ -256,6 +259,7 @@ subroutine read_ob_sens stattype(nob) = indata%stattype obtype(nob) = indata%obtype indxsat(nn) = indata%indxsat + assimltd_flag(nob) = indata%assimltd_flag if(nproc == 0) anal_ob_post(1:nanals,nob) = real(tmpanal_ob(1:nanals),r_kind) biaspreds(1:npred+1,nn) = real(tmpbiaspreds(1:npred+1),r_kind) end do @@ -430,6 +434,7 @@ subroutine print_ob_sens outdata%stattype = stattype(nob) outdata%obtype = obtype(nob) outdata%indxsat = 0 + outdata%assimltd_flag = assimltd_flag(nob) if(efsoi_flag) then outdata%osense_kin = real(obsense_kin(nob),r_single) outdata%osense_dry = real(obsense_dry(nob),r_single) @@ -500,6 +505,7 @@ subroutine print_ob_sens outdata%stattype = stattype(nob) outdata%obtype = obtype(nob) outdata%indxsat = nchan + outdata%assimltd_flag = assimltd_flag(nob) tmpbiaspreds(1:npred+1) = real(biaspreds(1:npred+1,nn),r_single) if(efsoi_flag) then outdata%osense_kin = real(obsense_kin(nob),r_single) diff --git a/src/enkf/enkf_obsmod.f90 b/src/enkf/enkf_obsmod.f90 index 72e196d61e..3a1bbe4363 100644 --- a/src/enkf/enkf_obsmod.f90 +++ b/src/enkf/enkf_obsmod.f90 @@ -152,6 +152,9 @@ module enkf_obsmod ! ob-space posterior ensemble, needed for EFSOI real(r_single),public,allocatable, dimension(:,:) :: anal_ob_post ! Fortran pointer +! is the observation assimilated? logical would be preferable, but that confuses +! Python +integer(i_kind),public,allocatable, dimension(:) :: assimltd_flag ! Fortran pointer contains @@ -466,6 +469,7 @@ subroutine obsmod_cleanup() if (allocated(prpgerr)) deallocate(prpgerr) if (allocated(diagused)) deallocate(diagused) if (allocated(anal_ob_post)) deallocate(anal_ob_post) +if (allocated(assimltd_flag)) deallocate(assimltd_flag) ! free shared memory segement, fortran pointer to that memory. nullify(anal_ob) call MPI_Barrier(mpi_comm_world,ierr)