GitHub Issue #118. EFSOI additions to src/gsi and src/enkf directories.

src/enkf/enkf.f90

@@ -123,14 +123,15 @@ module enkf
                   obtype, oberrvarmean, numobspersat, deltapredx, biaspreds,&
                   oberrvar_orig, probgrosserr, prpgerr,&
                   corrlengthsq,lnsigl,obtimel,obloclat,obloclon,obpress,stattype,&
-                  anal_ob
+                  anal_ob,anal_ob_post
 use constants, only: pi, one, zero
 use params, only: sprd_tol, paoverpb_thresh, datapath, nanals,&
                   iassim_order,sortinc,deterministic,numiter,nlevs,&
                   zhuberleft,zhuberright,varqc,lupd_satbiasc,huber,univaroz,&
                   covl_minfact,covl_efold,nbackgrounds,nhr_anal,fhr_assim,&
-                  iseed_perturbed_obs,lupd_obspace_serial,fso_cycling,&
+                  iseed_perturbed_obs,lupd_obspace_serial,efsoi_cycling,&
                   neigv,vlocal_evecs,denkf
+
 use radinfo, only: npred,nusis,nuchan,jpch_rad,predx
 use radbias, only: apply_biascorr, update_biascorr
 use gridinfo, only: nlevs_pres
@@ -825,24 +826,24 @@ subroutine enkf_update()
 
 ! Gathering analysis perturbations 
 ! in observation space for EFSO
-if(fso_cycling) then  
+if(efsoi_cycling) then  
    if(nproc /= 0) then   
       call mpi_send(anal_obchunk,numobsperproc(nproc+1)*nanals,mpi_real,0, &   
                     1,mpi_comm_world,ierr)   
    else   
-      allocate(anal_ob(1:nanals,nobstot))   
+      allocate(anal_ob_post(1:nanals,nobstot))   
       allocate(buffertmp3(nanals,nobs_max))   
       do np=1,numproc-1   
          call mpi_recv(buffertmp3,numobsperproc(np+1)*nanals,mpi_real,np, &   
                        1,mpi_comm_world,mpi_status,ierr)   
          do nob1=1,numobsperproc(np+1)   
             nob2 = indxproc_obs(np+1,nob1)   
-            anal_ob(:,nob2) = buffertmp3(:,nob1)   
+            anal_ob_post(:,nob2) = buffertmp3(:,nob1)   
          end do   
       end do   
       do nob1=1,numobsperproc(1)   
          nob2 = indxproc_obs(1,nob1)   
-         anal_ob(:,nob2) = anal_obchunk(:,nob1)   
+         anal_ob_post(:,nob2) = anal_obchunk(:,nob1)   
       end do   
       deallocate(buffertmp3)   
    end if   

src/enkf/enkf_main.f90

@@ -76,7 +76,7 @@ program enkf_main
  ! reads namelist parameters.
  use params, only : read_namelist,cleanup_namelist,letkf_flag,readin_localization,lupd_satbiasc,&
                     numiter, nanals, lupd_obspace_serial, write_spread_diag,   &
-                    lobsdiag_forenkf, netcdf_diag, fso_cycling, ntasks_io
+                    lobsdiag_forenkf, netcdf_diag, efsoi_cycling, ntasks_io
  ! mpi functions and variables.
  use mpisetup, only:  mpi_initialize, mpi_initialize_io, mpi_cleanup, nproc, &
                        mpi_wtime
@@ -183,7 +183,7 @@ program enkf_main
 
  ! Initialization for writing
  ! observation sensitivity files
- if(fso_cycling) call init_ob_sens()
+ if(efsoi_cycling) call init_ob_sens()
 
  ! read in vertical profile of horizontal and vertical localization length
  ! scales, set values for each ob.
@@ -216,7 +216,7 @@ program enkf_main
 
  ! Output non-inflated
  ! analyses for FSO
- if(fso_cycling) then
+ if(efsoi_cycling) then
     no_inflate_flag=.true.
     t1 = mpi_wtime()
     call gather_chunks()
@@ -240,7 +240,7 @@ program enkf_main
   endif
 
  ! print EFSO sensitivity i/o on root task.
- if(fso_cycling) call print_ob_sens()
+ if(efsoi_cycling) call print_ob_sens()
 
  ! print innovation statistics for posterior on root task.
  if (nproc == 0 .and. numiter > 0) then
@@ -268,7 +268,7 @@ program enkf_main
 
  call controlvec_cleanup()
  call loadbal_cleanup()
- if(fso_cycling) call destroy_ob_sens()
+ if(efsoi_cycling) call destroy_ob_sens()
  call cleanup_namelist()
 
  ! write log file (which script can check to verify completion).

src/enkf/enkf_obs_sensitivity.f90

@@ -18,6 +18,7 @@ module enkf_obs_sensitivity
 !   destroy_ob_sens - Deallocate variables
 !
 ! Variable Definitions:
+!   adloc_chunk - Coordinates of observation response
 !   obsense_kin - forecast sensitivity on each observations (kinetic energy)
 !   obsense_dry - forecast sensitivity on each observations (dry total energy)
 !   obsense_moist - forecast sensitivity on each observations (moist total energy)
@@ -32,23 +33,22 @@ module enkf_obs_sensitivity
 use mpisetup, only: mpi_real4,mpi_sum,mpi_comm_io,mpi_in_place,numproc,nproc,&
                 mpi_integer,mpi_wtime,mpi_status,mpi_real8,mpi_max,mpi_realkind
 use kinds, only: r_single,r_kind,r_double,i_kind
-use params, only: fso_calculate,latbound,nlevs,nanals,datestring, &
+use params, only: efsoi_flag,latbound,nlevs,nanals,datestring, &
                   lnsigcutoffsatnh,lnsigcutoffsattr,lnsigcutoffsatsh, &
                   lnsigcutoffpsnh,lnsigcutoffpstr,lnsigcutoffpssh, &
                   lnsigcutoffnh,lnsigcutofftr,lnsigcutoffsh, &
                   corrlengthnh,corrlengthtr,corrlengthsh, &
                   obtimelnh,obtimeltr,obtimelsh,letkf_flag, &
-                  nbackgrounds
+                  nbackgrounds,adrate,eft
 use constants, only: zero,one,half,rearth,pi,deg2rad,rad2deg
 use enkf_obsmod, only: nobstot,nobs_conv,nobs_oz,nobs_sat,obtype,obloclat, &
                        obloclon,obpress,indxsat,oberrvar,stattype,obtime,ob, &
                        ensmean_ob,ensmean_obnobc,obsprd_prior,obfit_prior, &
-                       oberrvar_orig,biaspreds,anal_ob,nobstot,lnsigl, &
+                       oberrvar_orig,biaspreds,anal_ob_post,nobstot,lnsigl, &
                        corrlengthsq,obtimel,oblnp,obloc
 use convinfo, only: convinfo_read,init_convinfo
 use ozinfo, only: ozinfo_read,init_oz
 use radinfo, only: radinfo_read,jpch_rad,nusis,nuchan,npred
-use gridinfo, only: latsgrd,lonsgrd,nlevs_pres,npts
 use loadbal, only: indxproc,grdloc_chunk,numptsperproc,npts_max,kdtree_grid
 use covlocal, only: latval
 use kdtree2_module, only: kdtree2_create
@@ -57,9 +57,10 @@ module enkf_obs_sensitivity
 
 private
 public init_ob_sens,destroy_ob_sens,print_ob_sens,read_ob_sens,&
-       obsense_kin,obsense_dry,obsense_moist
+       obsense_kin,obsense_dry,obsense_moist,adloc_chunk
 
 real(r_kind),allocatable,dimension(:) :: obsense_kin,obsense_dry,obsense_moist
+real(r_single),allocatable,dimension(:,:) :: adloc_chunk
 
 ! Structure for observation sensitivity information output
 type obsense_header
@@ -97,6 +98,7 @@ module enkf_obs_sensitivity
 
 contains
 
+
 subroutine init_ob_sens
 !$$$  subprogram documentation block
 !                .      .    .                                       .
@@ -129,6 +131,7 @@ subroutine init_ob_sens
 
 end subroutine init_ob_sens
 
+
 subroutine read_ob_sens
 !$$$  subprogram documentation block
 !                .      .    .                                       .
@@ -188,8 +191,15 @@ subroutine read_ob_sens
      write(6,*) 'READ_OBSENSE:  number of members is not correct.',nanals,inhead%nanals
      call stop2(26)
   end if
- ! nobstot=nobsgood
+
   if(nproc == 0) write(6,*) 'total number of obs ',nobstot
+  if(nproc == 0) write(6,*) 'total number of conv obs ',nobs_conv
+  if(nproc == 0) write(6,*) 'total number of oz obs',nobs_oz
+  if(nproc == 0) write(6,*) 'total number of sat obs',nobs_sat
+  if(nproc == 0) write(6,*) 'npred=',inhead%npred
+  if(nproc == 0) write(6,*) 'idate=',inhead%idate
+  if(nproc == 0) write(6,*) 'nanals=',inhead%nanals
+
   ! Allocate arrays
   allocate(obfit_prior(nobstot))
   allocate(obsprd_prior(nobstot))
@@ -208,7 +218,7 @@ subroutine read_ob_sens
   allocate(biaspreds(npred+1,nobs_sat))
   allocate(tmpanal_ob(nanals))
   allocate(tmpbiaspreds(npred+1))
-  if(nproc == 0) allocate(anal_ob(nanals,nobstot))
+  if(nproc == 0) allocate(anal_ob_post(nanals,nobstot))
   ! Read loop over conventional observations
   do nob=1,nobs_conv+nobs_oz
      read(iunit) indata,tmpanal_ob
@@ -225,7 +235,7 @@ subroutine read_ob_sens
      oberrvar_orig(nob) = real(indata%oberrvar_orig,r_kind)
      stattype(nob) = indata%stattype
      obtype(nob) = indata%obtype
-     if(nproc == 0) anal_ob(1:nanals,nob) = real(tmpanal_ob(1:nanals),r_kind)
+     if(nproc == 0) anal_ob_post(1:nanals,nob) = real(tmpanal_ob(1:nanals),r_kind)
   end do
   ! Read loop over satellite radiance observations
   nn = 0
@@ -246,7 +256,7 @@ subroutine read_ob_sens
      stattype(nob) = indata%stattype
      obtype(nob) = indata%obtype
      indxsat(nn) = indata%indxsat
-     if(nproc == 0) anal_ob(1:nanals,nob) = real(tmpanal_ob(1:nanals),r_kind)
+     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
   if(nn /= nobs_sat) then
@@ -332,7 +342,7 @@ subroutine print_ob_sens
   type(obsense_info) :: outdata
   iunit = 10
   ! Gather observation sensitivity informations to the root
-  if(fso_calculate) then
+  if(efsoi_flag) then
      allocate(recbuf(nobstot))
      call mpi_reduce(obsense_kin,recbuf,nobstot,mpi_realkind,mpi_sum,0, &
           & mpi_comm_world,ierr)
@@ -420,7 +430,7 @@ subroutine print_ob_sens
         outdata%stattype = stattype(nob)
         outdata%obtype = obtype(nob)
         outdata%indxsat = 0
-        if(fso_calculate) then
+        if(efsoi_flag) then
            outdata%osense_kin = real(obsense_kin(nob),r_single)
            outdata%osense_dry = real(obsense_dry(nob),r_single)
            outdata%osense_moist = real(obsense_moist(nob),r_single)
@@ -429,9 +439,9 @@ subroutine print_ob_sens
            outdata%osense_dry = 9.9e31_r_single
            outdata%osense_moist = 9.9e31_r_single
         end if
-        tmpanal_ob(1:nanals) = real(anal_ob(1:nanals,nob),r_single)
+        tmpanal_ob(1:nanals) = real(anal_ob_post(1:nanals,nob),r_single)
         write(iunit) outdata,tmpanal_ob
-        if(.not. fso_calculate) cycle
+        if(.not. efsoi_flag) cycle
         ! Sum up
         nob_conv(iobtyp,ireg) = nob_conv(iobtyp,ireg) + 1
         sumsense_conv(iobtyp,ireg,stkin) = sumsense_conv(iobtyp,ireg,stkin) &
@@ -448,7 +458,7 @@ subroutine print_ob_sens
              rate_conv(iobtyp,ireg,stmoist) = rate_conv(iobtyp,ireg,stmoist) + one
      end do
      ! print out
-     if(fso_calculate) then
+     if(efsoi_flag) then
         print *,'observation impact for conventional obs'
         print *,'region, obtype, nobs, dJ, positive rate[%]:'
         do iobtyp=1,8
@@ -489,9 +499,9 @@ subroutine print_ob_sens
         outdata%oberrvar_orig = real(oberrvar_orig(nob),r_single)
         outdata%stattype = stattype(nob)
         outdata%obtype = obtype(nob)
-        outdata%indxsat = nuchan(nchan)
+        outdata%indxsat = nchan
         tmpbiaspreds(1:npred+1) = real(biaspreds(1:npred+1,nn),r_single)
-        if(fso_calculate) then
+        if(efsoi_flag) then
            outdata%osense_kin = real(obsense_kin(nob),r_single)
            outdata%osense_dry = real(obsense_dry(nob),r_single)
            outdata%osense_moist = real(obsense_moist(nob),r_single)
@@ -500,9 +510,9 @@ subroutine print_ob_sens
            outdata%osense_dry = 9.9e31_r_single
            outdata%osense_moist = 9.9e31_r_single
         end if
-        tmpanal_ob(1:nanals) = real(anal_ob(1:nanals,nob),r_single)
+        tmpanal_ob(1:nanals) = real(anal_ob_post(1:nanals,nob),r_single)
         write(iunit) outdata,tmpanal_ob,tmpbiaspreds
-        if(.not. fso_calculate) cycle
+        if(.not. efsoi_flag) cycle
         ! Sum up
         if (oberrvar(nob) < 1.e10_r_kind .and. nchan > 0) then
            nob_sat(nchan) = nob_sat(nchan) + 1
@@ -521,7 +531,7 @@ subroutine print_ob_sens
         end if
      end do
      ! print out
-     if(fso_calculate) then
+     if(efsoi_flag) then
         print *,'observation impact for satellite brightness temp'
         print *,'instrument, channel #, nobs, dJ, positive rate[%]:'
         do nchan=1,jpch_rad
@@ -565,6 +575,7 @@ subroutine destroy_ob_sens
   if(allocated(obsense_kin)) deallocate(obsense_kin)
   if(allocated(obsense_dry)) deallocate(obsense_dry)
   if(allocated(obsense_moist)) deallocate(obsense_moist)
+  if(allocated(adloc_chunk)) deallocate(adloc_chunk)
   return
 end subroutine destroy_ob_sens
 end module enkf_obs_sensitivity

src/enkf/enkf_obsmod.f90

@@ -146,6 +146,9 @@ module enkf_obsmod
 type(c_ptr)                             :: anal_ob_modens_cp    ! C pointer
 integer :: shm_win, shm_win2
 
+! ob-space posterior ensemble, needed for EFSOI
+real(r_single),public,allocatable, dimension(:,:) :: anal_ob_post   ! Fortran pointer
+
 contains
 
 subroutine readobs()
@@ -453,6 +456,7 @@ subroutine obsmod_cleanup()
 if (allocated(probgrosserr)) deallocate(probgrosserr)
 if (allocated(prpgerr)) deallocate(prpgerr)
 if (allocated(diagused)) deallocate(diagused)
+if (allocated(anal_ob_post)) deallocate(anal_ob_post)
 ! free shared memory segement, fortran pointer to that memory.
 nullify(anal_ob)
 call MPI_Barrier(mpi_comm_world,ierr)

src/enkf/loadbal.f90

@@ -444,8 +444,6 @@ subroutine scatter_chunks
 
 end subroutine scatter_chunks
 
-
-
 subroutine gather_chunks
 ! gather chunks into grdin to write out the ensemble members
 use controlvec, only: ncdim, grdin