GitHub Issue NOAA-EMC#118. Implement EFSOI process in FV3 GFS workflo…

…w. Changes to src/enkf

src/enkf/enkf.f90

@@ -129,8 +129,9 @@ module enkf
                   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,7 +826,7 @@ 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)   

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,13 +33,13 @@ 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, &
@@ -48,18 +49,19 @@ module enkf_obs_sensitivity
 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 gridinfo_efsoi, only: latsgrd,lonsgrd,nlevs_pres,npts,id_u,id_v
+use loadbal!, only: indxproc,grdloc_chunk,numptsperproc,npts_max,kdtree_grid
 use covlocal, only: latval
 use kdtree2_module, only: kdtree2_create
 
 implicit none
 
 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 +99,7 @@ module enkf_obs_sensitivity
 
 contains
 
+
 subroutine init_ob_sens
 !$$$  subprogram documentation block
 !                .      .    .                                       .
@@ -129,6 +132,7 @@ subroutine init_ob_sens
 
 end subroutine init_ob_sens
 
+
 subroutine read_ob_sens
 !$$$  subprogram documentation block
 !                .      .    .                                       .
@@ -188,8 +192,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))
@@ -332,7 +343,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 +431,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)
@@ -431,7 +442,7 @@ subroutine print_ob_sens
         end if
         tmpanal_ob(1:nanals) = real(anal_ob(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 +459,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 +500,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)
@@ -502,7 +513,7 @@ subroutine print_ob_sens
         end if
         tmpanal_ob(1:nanals) = real(anal_ob(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 +532,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 +576,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/gridio_gfs.f90

@@ -40,7 +40,7 @@ module gridio
 !   language: f95
 !
 !$$$
- use constants, only: zero,one,cp,fv,rd,tiny_r_kind,max_varname_length,t0c,r0_05
+ use constants, only: zero,one,cp,fv,rd,tiny_r_kind,max_varname_length,t0c,r0_05,constants_initialized
  use params, only: nlons,nlats,nlevs,use_gfs_nemsio,pseudo_rh, &
                    cliptracers,datapath,imp_physics,use_gfs_ncio,cnvw_option, &
                    nanals

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