Merge pull request #260 from AndrewEichmann-NOAA/master

GitHub Issue NOAA-EMC/GSI#259. Add assimilated flag to observations in osense file

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)

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)

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
 
@@ -473,6 +476,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)