Github issue #42424242: fix for broken generate_ens=T option

src/gsi/bkerror.f90

@@ -62,6 +62,9 @@ subroutine bkerror(grady)
   use control_vectors, only: mvars,nrf,nrf_var,nrf_3d
   use timermod, only: timer_ini,timer_fnl
   use gsi_bundlemod, only: gsi_bundlegetpointer,gsi_bundlemerge,gsi_bundle,gsi_bundledup,gsi_bundledestroy
+  use general_sub2grid_mod, only: general_sub2grid,general_grid2sub
+  use general_commvars_mod, only: s2g_raf
+  use general_commvars_mod, only: s2g_cv
   use hybrid_ensemble_isotropic, only: sqrt_beta_s_mult
   use hybrid_ensemble_parameters, only: l_hyb_ens
   implicit none
@@ -126,7 +129,7 @@ subroutine bkerror(grady)
      endif
 
 !    Apply variances, as well as vertical & horizontal parts of background error
-     call bkgcov(mbundle)
+     call bkgcov(s2g_raf,mbundle)
 
 !    The following lines test that indeed proper application of cgkcov
 !    reproduces results of bkgcov - left as comments (please do not remove

src/gsi/bkgcov.f90

@@ -1,4 +1,4 @@
-subroutine bkgcov(cstate)
+subroutine bkgcov(grd,cstate)
 !$$$  subprogram documentation block
 !                .      .    .                                       .
 ! subprogram:    bkgcov    perform hor & vert of background error 
@@ -42,18 +42,19 @@ subroutine bkgcov(cstate)
   use gsi_bundlemod, only: gsi_bundle
   use gsi_bundlemod, only: gsi_bundlegetpointer
   use general_sub2grid_mod, only: general_sub2grid,general_grid2sub
-  use general_commvars_mod, only: s2g_raf
+  use general_sub2grid_mod, only: sub2grid_info
   implicit none
+  type(sub2grid_info), intent(in) :: grd
 
 ! Passed Variables
   type(gsi_bundle),intent(inout) :: cstate
 
 ! Local Variables
   integer(i_kind) n,n3d,istatus,nlevs
-  real(r_kind),dimension(nlat*nlon*s2g_raf%nlevs_alloc):: hwork
+  real(r_kind),dimension(nlat*nlon*grd%nlevs_alloc):: hwork
   real(r_kind),pointer,dimension(:,:,:):: ptr3d=>NULL()
 
-  nlevs=s2g_raf%nlevs_loc
+  nlevs=grd%nlevs_loc
   n3d=cstate%n3d
 
 ! Multiply by background error variances, and break up skin temp
@@ -68,13 +69,13 @@ subroutine bkgcov(cstate)
   end do
 
 ! Convert from subdomain to full horizontal field distributed among processors
-  call general_sub2grid(s2g_raf,cstate%values,hwork)
+  call general_sub2grid(grd,cstate%values,hwork)
 
 ! Apply horizontal smoother for number of horizontal scales
   call smoothrf(hwork,nlevs)
 
 ! Put back onto subdomains
-  call general_grid2sub(s2g_raf,hwork,cstate%values)
+  call general_grid2sub(grd,hwork,cstate%values)
 
 ! Apply vertical smoother
 !$omp parallel do  schedule(dynamic,1) private(n,ptr3d,istatus)
@@ -90,7 +91,7 @@ subroutine bkgcov(cstate)
   return
 end subroutine bkgcov
 ! -----------------------------------------------------------------------------
-subroutine ckgcov(z,cstate,nval_lenz)
+subroutine ckgcov(grd,z,cstate,nval_lenz)
 !$$$  subprogram documentation block
 !                .      .    .                                       .
 ! subprogram:    ckgcov   sqrt of bkgcov
@@ -132,29 +133,29 @@ subroutine ckgcov(z,cstate,nval_lenz)
   use gridmod, only: nlat,nlon
   use gsi_bundlemod, only: gsi_bundle
   use gsi_bundlemod, only: gsi_bundlegetpointer
-  use general_sub2grid_mod, only: general_grid2sub
-  use general_commvars_mod, only: s2g_raf
+  use general_sub2grid_mod, only: general_grid2sub, sub2grid_info
   use hybrid_ensemble_parameters, only: l_hyb_ens
   use hybrid_ensemble_isotropic, only: sqrt_beta_s_mult
   implicit none
 
 ! Passed Variables
   integer(i_kind)    ,intent(in   ) :: nval_lenz
+  type(sub2grid_info), intent(in) :: grd
   type(gsi_bundle),intent(inout) :: cstate
   real(r_kind),dimension(nval_lenz),intent(in   ) :: z
 
 ! Local Variables
   integer(i_kind) k,n3d,istatus,nlevs
-  real(r_kind),dimension(nlat*nlon*s2g_raf%nlevs_alloc):: hwork
+  real(r_kind),dimension(nlat*nlon*grd%nlevs_alloc):: hwork
   real(r_kind),dimension(:,:,:),pointer:: ptr3d=>NULL()
 
-  nlevs=s2g_raf%nlevs_loc
+  nlevs=grd%nlevs_loc
 
 ! Apply horizontal smoother for number of horizontal scales
   call sqrt_smoothrf(z,hwork,nlevs)
 
 ! Put back onto subdomains
-  call general_grid2sub(s2g_raf,hwork,cstate%values)
+  call general_grid2sub(grd,hwork,cstate%values)
 
 ! Apply vertical smoother
   n3d=cstate%n3d
@@ -174,7 +175,7 @@ subroutine ckgcov(z,cstate,nval_lenz)
   return
 end subroutine ckgcov
 ! -----------------------------------------------------------------------------
-subroutine ckgcov_ad(z,cstate,nval_lenz)
+subroutine ckgcov_ad(grd,z,cstate,nval_lenz)
 !$$$  subprogram documentation block
 !                .      .    .                                       .
 ! subprogram:    ckgcov_ad  adjoint of ckgcov
@@ -217,11 +218,11 @@ subroutine ckgcov_ad(z,cstate,nval_lenz)
   use gridmod, only: nlat,nlon
   use gsi_bundlemod, only: gsi_bundle
   use gsi_bundlemod, only: gsi_bundlegetpointer
-  use general_sub2grid_mod, only: general_sub2grid
-  use general_commvars_mod, only: s2g_raf
+  use general_sub2grid_mod, only: general_sub2grid,sub2grid_info
   use hybrid_ensemble_parameters, only: l_hyb_ens
   use hybrid_ensemble_isotropic, only: sqrt_beta_s_mult
   implicit none
+  type(sub2grid_info), intent(in) :: grd
 
 ! Passed Variables
   integer(i_kind)    ,intent(in   ) :: nval_lenz
@@ -230,10 +231,10 @@ subroutine ckgcov_ad(z,cstate,nval_lenz)
 
 ! Local Variables
   integer(i_kind) k,n3d,istatus,nlevs
-  real(r_kind),dimension(nlat*nlon*s2g_raf%nlevs_alloc):: hwork
+  real(r_kind),dimension(nlat*nlon*grd%nlevs_alloc):: hwork
   real(r_kind),dimension(:,:,:),pointer:: ptr3d=>NULL()
 
-  nlevs=s2g_raf%nlevs_loc
+  nlevs=grd%nlevs_loc
 
 ! Apply static betas
   if(l_hyb_ens) call sqrt_beta_s_mult(cstate)
@@ -251,7 +252,7 @@ subroutine ckgcov_ad(z,cstate,nval_lenz)
   end do
 
 ! Convert from subdomain to full horizontal field distributed among processors
-  call general_sub2grid(s2g_raf,cstate%values,hwork)
+  call general_sub2grid(grd,cstate%values,hwork)
 
 ! Apply horizontal smoother for number of horizontal scales
   call sqrt_smoothrf_ad(z,hwork,nlevs)

src/gsi/control2model.f90

@@ -173,7 +173,7 @@ subroutine control2model(xhat,sval,bval)
 !  Apply sqrt of variance, as well as vertical & horizontal parts of background
 !  error
 
-   call ckgcov(xhat%step(jj)%values(:),wbundle,size(xhat%step(jj)%values(:)))
+   call ckgcov(grid,xhat%step(jj)%values(:),wbundle,size(xhat%step(jj)%values(:)))
 
 !  Get pointers to required state variables
    call gsi_bundlegetpointer (sval(jj),'u'   ,sv_u,   istatus)

src/gsi/control2model_ad.f90

@@ -236,7 +236,7 @@ subroutine control2model_ad(rval,bval,grad)
    enddo
 
 !  Apply adjoint of sqrt-B
-   call ckgcov_ad(gradz,wbundle,nval_lenz)
+   call ckgcov_ad(grid,gradz,wbundle,nval_lenz)
 
 !  Clean up
    call gsi_bundledestroy(wbundle,istatus)

src/gsi/cplr_gfs_ensmod.f90

@@ -1371,6 +1371,7 @@ subroutine put_gfs_ens(this,grd,member,ntindex,pert,iret)
     use hybrid_ensemble_parameters, only: ensemble_path
     use hybrid_ensemble_parameters, only: sp_ens
     use gridmod, only: use_gfs_nemsio, use_gfs_ncio
+    use netcdfgfs_io, only: write_gfsncatm_pert
 
     implicit none
 
@@ -1402,11 +1403,12 @@ subroutine put_gfs_ens(this,grd,member,ntindex,pert,iret)
        endif
        !call write_nemsatm(grd,...)
     else if ( use_gfs_ncio ) then
-       if ( mype == 0 ) then
-          write(6,*) 'write_gfsncatm is not adapted to write out perturbations yet'
-          iret = 999
-       endif
+       !if ( mype == 0 ) then
+       !   write(6,*) 'write_gfsncatm is not adapted to write out perturbations yet'
+       !   iret = 999
+       !endif
        !call write_gfsncatm(grd,...)
+       call write_gfsncatm_pert(grd,sp_ens,filename,mype_atm,pert,ntindex)
     else
        call general_write_gfsatm(grd,sp_ens,sp_ens,filename,mype_atm, &
             pert,ntindex,inithead,iret)

src/gsi/glbsoi.f90

@@ -150,7 +150,7 @@ subroutine glbsoi
   use zrnmi_mod, only: zrnmi_initialize
   use observermod, only: observer_init,observer_set,observer_finalize,ndata
   use timermod, only: timer_ini, timer_fnl
-  use hybrid_ensemble_parameters, only: l_hyb_ens,destroy_hybens_localization_parameters
+  use hybrid_ensemble_parameters, only: l_hyb_ens,destroy_hybens_localization_parameters,create_hybens_localization_parameters,generate_ens,write_generated_ens
   use hybrid_ensemble_isotropic, only: create_ensemble,load_ensemble,destroy_ensemble, &
        hybens_localization_setup,hybens_grid_setup
   use gfs_stratosphere, only: destroy_nmmb_vcoords,use_gfs_stratosphere
@@ -278,6 +278,7 @@ subroutine glbsoi
 
 ! If l_hyb_ens is true, then read in ensemble perturbations
   if(l_hyb_ens) then
+     call create_hybens_localization_parameters
      call load_ensemble
      call hybens_localization_setup
   end if
@@ -312,6 +313,9 @@ subroutine glbsoi
   if (lsensrecompute) jiterlast=jiterend
   if (l4dvar) jiterlast=jiterstart
   if (ladtest_obs) jiterlast=jiterstart
+  if (l_hyb_ens .and. generate_ens .and. write_generated_ens) then
+     jiterlast=0 ! skip analysis loop
+  endif
 
 ! Main outer analysis loop
   do jiter=jiterstart,jiterlast

src/gsi/gsimod.F90

@@ -144,7 +144,7 @@ module gsimod
   use hybrid_ensemble_parameters,only : l_hyb_ens,uv_hyb_ens,aniso_a_en,generate_ens,&
                          n_ens,nlon_ens,nlat_ens,jcap_ens,jcap_ens_test,oz_univ_static,&
                          regional_ensemble_option,fv3sar_ensemble_opt,merge_two_grid_ensperts, &
-                         full_ensemble,pseudo_hybens,pwgtflg,&
+                         full_ensemble,pseudo_hybens,pwgtflg,write_generated_ens,&
                          beta_s0,s_ens_h,s_ens_v,init_hybrid_ensemble_parameters,&
                          readin_localization,write_ens_sprd,eqspace_ensgrid,grid_ratio_ens,&
                          readin_beta,use_localization_grid,use_gfs_ens,q_hyb_ens,i_en_perts_io, &
@@ -715,8 +715,8 @@ module gsimod
        idmodel,iwrtinc,lwrite4danl,nhr_anal,jiterstart,jiterend,lobserver,lanczosave,llancdone, &
        lferrscale,print_diag_pcg,tsensible,lread_obs_save,lread_obs_skip, &
        use_gfs_ozone,check_gfs_ozone_date,regional_ozone,lwrite_predterms,&
-       lwrite_peakwt,use_gfs_nemsio,use_gfs_ncio,sfcnst_comb,liauon,use_prepb_satwnd,l4densvar,ens_nstarthr,&
-       use_gfs_stratosphere,pblend0,pblend1,step_start,diag_precon,lrun_subdirs,&
+       lwrite_peakwt,use_gfs_nemsio,use_gfs_ncio,sfcnst_comb,liauon,use_prepb_satwnd,l4densvar,&
+       ens_nstarthr,use_gfs_stratosphere,pblend0,pblend1,step_start,diag_precon,lrun_subdirs,&
        use_sp_eqspace,lnested_loops,lsingleradob,thin4d,use_readin_anl_sfcmask,&
        luse_obsdiag,id_drifter,id_ship,verbose,print_obs_para,lsingleradar,singleradar,lnobalance, &
        missing_to_nopcp,minobrangedbz,minobrangedbz,maxobrangedbz,&
@@ -1255,6 +1255,7 @@ module gsimod
 !     oz_univ_static- if true, decouple ozone from other variables and defaults to static B (ozone only)
 !     aniso_a_en - if true, then use anisotropic localization of hybrid ensemble control variable a_en.
 !     generate_ens - if true, then generate internal ensemble based on existing background error
+!     write_generated_ens - if true, then writed generated internal ensemble based on existing background error
 !     n_ens        - number of ensemble members.
 !     nlon_ens     - number of longitudes on ensemble grid (may be different from analysis grid nlon)
 !     nlat_ens     - number of latitudes on ensemble grid (may be different from analysis grid nlat)
@@ -1313,7 +1314,7 @@ module gsimod
   namelist/hybrid_ensemble/l_hyb_ens,uv_hyb_ens,q_hyb_ens,aniso_a_en,generate_ens,n_ens,nlon_ens,nlat_ens,jcap_ens,&
                 pseudo_hybens,merge_two_grid_ensperts,regional_ensemble_option,fv3sar_bg_opt,fv3sar_ensemble_opt,full_ensemble,pwgtflg,&
                 jcap_ens_test,beta_s0,s_ens_h,s_ens_v,readin_localization,eqspace_ensgrid,readin_beta,&
-                grid_ratio_ens, &
+                grid_ratio_ens,write_generated_ens, &
                 oz_univ_static,write_ens_sprd,use_localization_grid,use_gfs_ens, &
                 i_en_perts_io,l_ens_in_diff_time,ensemble_path,ens_fast_read,sst_staticB