diff --git a/src/gsi/class_get_fv3_regional_ensperts.f90 b/src/gsi/class_get_fv3_regional_ensperts.f90
index 6c51a4fd88..bda4a47cd3 100644
--- a/src/gsi/class_get_fv3_regional_ensperts.f90
+++ b/src/gsi/class_get_fv3_regional_ensperts.f90
@@ -31,7 +31,7 @@ subroutine get_fv3_regional_ensperts(this,en_perts,nelen,ps_bar)
     import abstract_get_fv3_regional_ensperts_class
     implicit none
     class(abstract_get_fv3_regional_ensperts_class),intent(inout) :: this
-    type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:)
+    type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:,:)
     integer(i_kind), intent(in   ):: nelen
     real(r_single),dimension(:,:,:),allocatable, intent(inout):: ps_bar
 
diff --git a/src/gsi/class_get_pseudo_ensperts.f90 b/src/gsi/class_get_pseudo_ensperts.f90
index 703a62df3f..beee90f5e9 100644
--- a/src/gsi/class_get_pseudo_ensperts.f90
+++ b/src/gsi/class_get_pseudo_ensperts.f90
@@ -12,7 +12,7 @@ subroutine get_pseudo_ensperts(this,en_perts,nelen)
     import abstract_get_pseudo_ensperts_class
     implicit none
     class(abstract_get_pseudo_ensperts_class), intent(inout) :: this
-    type(gsi_bundle),allocatable, intent(in   ) :: en_perts(:,:)
+    type(gsi_bundle),allocatable, intent(in   ) :: en_perts(:,:,:)
     integer(i_kind),                   intent(in   ) :: nelen
   end subroutine get_pseudo_ensperts
   end interface
diff --git a/src/gsi/class_get_wrf_mass_ensperts.f90 b/src/gsi/class_get_wrf_mass_ensperts.f90
index 902b14ad21..b392bcd174 100644
--- a/src/gsi/class_get_wrf_mass_ensperts.f90
+++ b/src/gsi/class_get_wrf_mass_ensperts.f90
@@ -12,7 +12,7 @@ subroutine get_wrf_mass_ensperts(this,en_perts,nelen,ps_bar)
     import abstract_get_wrf_mass_ensperts_class
     implicit none
     class(abstract_get_wrf_mass_ensperts_class), intent(inout) :: this
-    type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:)
+    type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:,:)
     integer(i_kind), intent(in   ):: nelen
     real(r_single),dimension(:,:,:),allocatable:: ps_bar
   end subroutine get_wrf_mass_ensperts
@@ -25,7 +25,7 @@ subroutine ens_spread_dualres_regional(this,mype,en_perts,nelen,en_bar)
     implicit none
     class(abstract_get_wrf_mass_ensperts_class), intent(inout) :: this
     integer(i_kind),intent(in):: mype
-    type(gsi_bundle),allocatable, intent(in   ) :: en_perts(:,:)
+    type(gsi_bundle),allocatable, intent(in   ) :: en_perts(:,:,:)
     integer(i_kind), intent(in   ):: nelen
     type(gsi_bundle),OPTIONAL,intent(in):: en_bar
   end subroutine ens_spread_dualres_regional
diff --git a/src/gsi/class_get_wrf_nmm_ensperts.f90 b/src/gsi/class_get_wrf_nmm_ensperts.f90
index 19170033e2..9f496c7b03 100644
--- a/src/gsi/class_get_wrf_nmm_ensperts.f90
+++ b/src/gsi/class_get_wrf_nmm_ensperts.f90
@@ -13,7 +13,7 @@ subroutine get_wrf_nmm_ensperts(this,en_perts,nelen,region_lat_ens,region_lon_en
     import abstract_get_wrf_nmm_ensperts_class
     implicit none
     class(abstract_get_wrf_nmm_ensperts_class),intent(inout) :: this
-    type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:)
+    type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:,:)
     integer(i_kind), intent(in   ):: nelen
     real(r_kind),allocatable, intent(inout):: region_lat_ens(:,:),region_lon_ens(:,:)
     real(r_single),dimension(:,:,:),allocatable, intent(inout):: ps_bar
diff --git a/src/gsi/control_vectors.f90 b/src/gsi/control_vectors.f90
index a1c0a61ed9..0847257777 100644
--- a/src/gsi/control_vectors.f90
+++ b/src/gsi/control_vectors.f90
@@ -144,7 +144,7 @@ module control_vectors
    type(GSI_Bundle), pointer :: step(:)
    type(GSI_Bundle), pointer :: motley(:)
    type(GSI_Grid)  :: grid_aens
-   type(GSI_Bundle), pointer :: aens(:,:)
+   type(GSI_Bundle), pointer :: aens(:,:,:)
    real(r_kind), pointer :: predr(:) => NULL()
    real(r_kind), pointer :: predp(:) => NULL()
    real(r_kind), pointer :: predt(:) => NULL()
@@ -296,6 +296,7 @@ subroutine init_anacv
 !   language: f90
 !   machine:  ibm rs/6000 sp
 !
+use hybrid_ensemble_parameters,only:idaen3d,idaen2d
 implicit none
 !character(len=*),parameter:: rcname='anavinfo.txt'
 character(len=*),parameter:: rcname='anavinfo'  ! filename should have extension
@@ -345,6 +346,7 @@ subroutine init_anacv
 allocate(cvarsmd(mvars))
 allocate(atsfc_sdv(mvars))
 allocate(an_amp0(nvars))
+allocate(idaen3d(nc3d),idaen2d(nc2d))
 allocate(incvars_to_zero(nvars))
 allocate(incvars_zero_strat(nvars))
 incvars_to_zero(:) = 'NONE'
@@ -371,12 +373,22 @@ subroutine init_anacv
          cvars2d(nc2d)=trim(adjustl(var))
          nrf2_loc(nc2d)=ii  ! rid of soon
          as2d(nc2d)=aas
+         if(itracer>10) then
+            idaen2d(nc2d)=2
+         else
+            idaen2d(nc2d)=1
+         endif
       else
          nc3d=nc3d+1
          cvars3d(nc3d)=trim(adjustl(var))
          nrf3_loc(nc3d)=ii  ! rid of soon
          nrf_3d(ii)=.true.
          as3d(nc3d)=aas
+         if(itracer>10) then
+            idaen3d(nc3d)=2
+         else
+            idaen3d(nc3d)=1
+         endif
       endif
    endif
    nrf_var(ii)=trim(adjustl(var))
@@ -446,9 +458,11 @@ subroutine allocate_cv(ycv)
 !$$$ end documentation block
 
   use hybrid_ensemble_parameters, only: grd_ens
+  use hybrid_ensemble_parameters, only: naensgrp
   implicit none
   type(control_vector), intent(  out) :: ycv
   integer(i_kind) :: ii,jj,nn,ndim,ierror,n_step,n_aens
+  integer(i_kind) :: ig
   character(len=256)::bname
   character(len=max_varname_length)::ltmp(1) 
   type(gsi_grid) :: grid_motley
@@ -481,7 +495,7 @@ subroutine allocate_cv(ycv)
 ! If so, define grid of ensemble control vector
   n_aens=0
   if (l_hyb_ens) then
-      ALLOCATE(ycv%aens(nsubwin,n_ens))
+      ALLOCATE(ycv%aens(nsubwin,naensgrp,n_ens))
       call GSI_GridCreate(ycv%grid_aens,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig)
          if (lsqrtb) then
             n_aens=nval_lenz_en
@@ -524,17 +538,19 @@ subroutine allocate_cv(ycv)
      if (l_hyb_ens) then
 
          ltmp(1)='a_en'
-         do nn=1,n_ens
-            ycv%aens(jj,nn)%values => ycv%values(ii+1:ii+n_aens)
-            write(bname,'(a,i3.3,a,i4.4)') 'Ensemble Control Bundle subwin-',jj,' and member-',nn
-            call GSI_BundleSet(ycv%aens(jj,nn),ycv%grid_aens,bname,ierror,names3d=ltmp,bundle_kind=r_kind)
-            if (ierror/=0) then
-                write(6,*)'allocate_cv: error alloc(ensemble bundle)'
-                call stop2(109)
-            endif
-            ndim=ndim+ycv%aens(jj,nn)%ndim
-
-            ii=ii+n_aens
+         do ig=1,naensgrp
+            do nn=1,n_ens
+               ycv%aens(jj,ig,nn)%values => ycv%values(ii+1:ii+n_aens)
+               write(bname,'(a,i3.3,a,i4.4)') 'Ensemble Control Bundle subwin-',jj,' and member-',nn
+               call GSI_BundleSet(ycv%aens(jj,ig,nn),ycv%grid_aens,bname,ierror,names3d=ltmp,bundle_kind=r_kind)
+               if (ierror/=0) then
+                  write(6,*)'allocate_cv: error alloc(ensemble bundle)'
+                  call stop2(109)
+               endif
+               ndim=ndim+ycv%aens(jj,ig,nn)%ndim
+
+               ii=ii+n_aens
+            enddo
          enddo
 
      endif
@@ -620,15 +636,19 @@ subroutine deallocate_cv(ycv)
 !
 !$$$ end documentation block
 
+  use hybrid_ensemble_parameters, only: naensgrp
   implicit none
   type(control_vector), intent(inout) :: ycv
   integer(i_kind) :: ii,nn,ierror
+  integer(i_kind) :: ig
 
   if (ycv%lallocated) then
      do ii=1,nsubwin
         if (l_hyb_ens) then
-           do nn=n_ens,1,-1
-              call GSI_BundleUnset(ycv%aens(ii,nn),ierror)
+           do ig=1,naensgrp
+              do nn=n_ens,1,-1
+                 call GSI_BundleUnset(ycv%aens(ii,ig,nn),ierror)
+              enddo
            enddo
         endif
 !       if (beta_s0>tiny_r_kind) then
@@ -845,9 +865,11 @@ real(r_quad) function qdot_prod_sub(xcv,ycv)
 !
 !$$$ end documentation block
 
+  use hybrid_ensemble_parameters, only: naensgrp
   implicit none
   type(control_vector), intent(in   ) :: xcv, ycv
   integer(i_kind) :: ii,nn,m3d,m2d,i,j,itot
+  integer(i_kind) :: ig,nigtmp
   real(r_quad),allocatable,dimension(:) :: partsum
 
   qdot_prod_sub=zero_quad
@@ -858,9 +880,11 @@ real(r_quad) function qdot_prod_sub(xcv,ycv)
         qdot_prod_sub=qdot_prod_sub+qdot_product( xcv%step(ii)%values(:) ,ycv%step(ii)%values(:) )
      end do
      if(l_hyb_ens) then
-        do nn=1,n_ens
-           do ii=1,nsubwin
-              qdot_prod_sub=qdot_prod_sub+qdot_product( xcv%aens(ii,nn)%values(:) ,ycv%aens(ii,nn)%values(:) )
+        do ig=1,naensgrp
+           do nn=1,n_ens
+              do ii=1,nsubwin
+                 qdot_prod_sub=qdot_prod_sub+qdot_product( xcv%aens(ii,ig,nn)%values(:) ,ycv%aens(ii,ig,nn)%values(:) )
+              end do
            end do
         end do
      endif
@@ -869,7 +893,7 @@ real(r_quad) function qdot_prod_sub(xcv,ycv)
         m3d=xcv%step(ii)%n3d
         m2d=xcv%step(ii)%n2d
         itot=max(m3d,0)+max(m2d,0)
-        if(l_hyb_ens)itot=itot+n_ens
+        if(l_hyb_ens)itot=itot+naensgrp*n_ens
         allocate(partsum(itot))
 !$omp parallel do  schedule(dynamic,1) private(i)
         do i = 1,m3d
@@ -880,9 +904,12 @@ real(r_quad) function qdot_prod_sub(xcv,ycv)
            partsum(m3d+i) = dplevs(xcv%step(ii)%r2(i)%q,ycv%step(ii)%r2(i)%q,ihalo=1)
         enddo
         if(l_hyb_ens) then
+           do ig=1,naensgrp
+              nigtmp=n_ens*(ig-1)
 !$omp parallel do  schedule(dynamic,1) private(i)
-           do i = 1,n_ens
-              partsum(m3d+m2d+i) = dplevs(xcv%aens(ii,i)%r3(1)%q,ycv%aens(ii,i)%r3(1)%q,ihalo=1)
+              do i = 1,n_ens
+                 partsum(m3d+m2d+nigtmp+i) = dplevs(xcv%aens(ii,ig,i)%r3(1)%q,ycv%aens(ii,ig,i)%r3(1)%q,ihalo=1)
+              end do
            end do
         end if
         do i=1,itot
@@ -933,6 +960,7 @@ subroutine qdot_prod_vars_eb(xcv,ycv,prods,eb)
 !
 !$$$ end documentation block
 
+  use hybrid_ensemble_parameters, only: naensgrp
   implicit none
   type(control_vector), intent(in   ) :: xcv, ycv
   character(len=*)    , intent(in   ) :: eb
@@ -941,6 +969,8 @@ subroutine qdot_prod_vars_eb(xcv,ycv,prods,eb)
   real(r_quad) :: zz(nsubwin)
   integer(i_kind) :: ii,i,nn,m3d,m2d
   real(r_quad),allocatable,dimension(:) :: partsum
+  integer(i_kind) :: ig
+  integer(i_kind) ::ngtmp,nn0
 
   prods(:)=zero_quad
   zz(:)=zero_quad
@@ -953,9 +983,11 @@ subroutine qdot_prod_vars_eb(xcv,ycv,prods,eb)
         end do
      endif
      if(trim(eb) == 'cost_e') then
-        do nn=1,n_ens
-           do ii=1,nsubwin
-              zz(ii)=zz(ii)+qdot_product( xcv%aens(ii,nn)%values(:) ,ycv%aens(ii,nn)%values(:) )
+        do ig=1,naensgrp
+           do nn=1,n_ens
+              do ii=1,nsubwin
+                 zz(ii)=zz(ii)+qdot_product( xcv%aens(ii,ig,nn)%values(:) ,ycv%aens(ii,ig,nn)%values(:) )
+              end do
            end do
         end do
      endif
@@ -981,21 +1013,25 @@ subroutine qdot_prod_vars_eb(xcv,ycv,prods,eb)
      end if
      if(trim(eb) == 'cost_e') then
         do ii=1,nsubwin ! RTod: somebody could work in opt/zing this ...
-           allocate(partsum(n_ens))
+           allocate(partsum(n_ens*naensgrp))
+           do ig=1,naensgrp
+              ngtmp=(ig-1)*n_ens
 !$omp parallel do  schedule(dynamic,1) private(nn,m3d,m2d)
-           do nn=1,n_ens
-              partsum(nn) = zero_quad
-              m3d=xcv%aens(ii,nn)%n3d
-              do i = 1,m3d
-                 partsum(nn)= partsum(nn) + dplevs(xcv%aens(ii,nn)%r3(i)%q,ycv%aens(ii,nn)%r3(i)%q,ihalo=1)
+              do nn=1,n_ens
+                 nn0=nn+ngtmp
+                 partsum(nn0) = zero_quad
+                 m3d=xcv%aens(ii,ig,nn)%n3d
+                 do i = 1,m3d
+                    partsum(nn0)= partsum(nn0) + dplevs(xcv%aens(ii,ig,nn)%r3(i)%q,ycv%aens(ii,ig,nn)%r3(i)%q,ihalo=1)
+                 enddo
+                 m2d=xcv%aens(ii,ig,nn)%n2d
+                 do i = 1,m2d
+                    partsum(nn0)= partsum(nn0) + dplevs(xcv%aens(ii,ig,nn)%r2(i)%q,ycv%aens(ii,ig,nn)%r2(i)%q,ihalo=1)
+                 enddo
               enddo
-              m2d=xcv%aens(ii,nn)%n2d
-              do i = 1,m2d
-                 partsum(nn)= partsum(nn) + dplevs(xcv%aens(ii,nn)%r2(i)%q,ycv%aens(ii,nn)%r2(i)%q,ihalo=1)
-              enddo
-           enddo
-           do nn=1,n_ens
-             zz(ii)=zz(ii)+partsum(nn)
+           end do
+           do nn=1,n_ens*naensgrp
+              zz(ii)=zz(ii)+partsum(nn)
            end do
            deallocate(partsum)
         end do
@@ -1348,6 +1384,7 @@ subroutine random_cv(ycv,kseed)
 !
 !$$$ end documentation block
 
+use hybrid_ensemble_parameters, only: naensgrp
 implicit none
 type(control_vector)     , intent(inout) :: ycv
 integer(i_kind), optional, intent(in   ) :: kseed
@@ -1355,6 +1392,7 @@ subroutine random_cv(ycv,kseed)
 integer(i_kind):: ii,jj,nn,iseed
 integer, allocatable :: nseed(:) ! Intentionaly default integer
 real(r_kind), allocatable :: zz(:)
+integer(i_kind) :: ig
 
 iseed=iadatebgn
 if (present(kseed)) iseed=iseed+kseed
@@ -1381,10 +1419,12 @@ subroutine random_cv(ycv,kseed)
 if (nval_lenz_en>0) then
    allocate(zz(nval_lenz_en))
    do nn=1,n_ens
-      do jj=1,nsubwin
-         call random_number(zz)
-         do ii=1,nval_lenz_en
-            ycv%aens(jj,nn)%values(ii) = two*zz(ii)-one
+      do ig=1,naensgrp
+         do jj=1,nsubwin
+            call random_number(zz)
+            do ii=1,nval_lenz_en
+               ycv%aens(jj,ig,nn)%values(ii) = two*zz(ii)-one
+            enddo
          enddo
       enddo
    enddo
diff --git a/src/gsi/cplr_get_fv3_regional_ensperts.f90 b/src/gsi/cplr_get_fv3_regional_ensperts.f90
index d64f266026..f99ca39790 100644
--- a/src/gsi/cplr_get_fv3_regional_ensperts.f90
+++ b/src/gsi/cplr_get_fv3_regional_ensperts.f90
@@ -73,7 +73,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar)
 
      implicit none
      class(get_fv3_regional_ensperts_class), intent(inout) :: this
-     type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:)
+     type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:,:)
      integer(i_kind), intent(in   ):: nelen
      real(r_single),dimension(:,:,:),allocatable,intent(inout):: ps_bar
  
@@ -248,7 +248,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar)
         en_bar(m)%values=zero
  
         do n=imem_start,n_ens
-           en_perts(n,m)%valuesr4 = zero
+           en_perts(n,1,m)%valuesr4 = zero
         enddo
  
         mm1=mype+1
@@ -311,7 +311,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar)
  ! SAVE ENSEMBLE MEMBER DATA IN COLUMN VECTOR
            do ic3=1,nc3d
  
-              call gsi_bundlegetpointer(en_perts(n,m),trim(cvars3d(ic3)),w3,istatus)
+              call gsi_bundlegetpointer(en_perts(n,1,m),trim(cvars3d(ic3)),w3,istatus)
               if(istatus/=0) then
                  write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for ensemble member ',n
                  call stop2(9992)
@@ -471,7 +471,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar)
  
            do ic2=1,nc2d
     
-              call gsi_bundlegetpointer(en_perts(n,m),trim(cvars2d(ic2)),w2,istatus)
+              call gsi_bundlegetpointer(en_perts(n,1,m),trim(cvars2d(ic2)),w2,istatus)
               if(istatus/=0) then
                  write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for ensemble member ',n
                  call stop2(9994)
@@ -540,7 +540,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar)
  
         do n=imem_start,n_ens
            do i=1,nelen
-              en_perts(n,m)%valuesr4(i)=(en_perts(n,m)%valuesr4(i)-en_bar(m)%values(i))*sig_norm
+              en_perts(n,1,m)%valuesr4(i)=(en_perts(n,1,m)%valuesr4(i)-en_bar(m)%values(i))*sig_norm
            end do
         end do
 
@@ -880,7 +880,7 @@ subroutine ens_spread_dualres_regional_fv3_regional(this,mype,en_perts,nelen)
 
     class(get_fv3_regional_ensperts_class), intent(inout) :: this
     integer(i_kind),intent(in):: mype
-    type(gsi_bundle),allocatable, intent(in   ) :: en_perts(:,:)
+    type(gsi_bundle),allocatable, intent(in   ) :: en_perts(:,:,:)
     integer(i_kind), intent(in   ):: nelen
   
     type(gsi_bundle):: sube,suba
@@ -922,7 +922,7 @@ subroutine ens_spread_dualres_regional_fv3_regional(this,mype,en_perts,nelen)
          do n=1,n_ens
             do i=1,nelen
                sube%values(i)=sube%values(i) &
-                 +(en_perts(n,1)%valuesr4(i))*(en_perts(n,1)%valuesr4(i))
+                 +(en_perts(n,1,1)%valuesr4(i))*(en_perts(n,1,1)%valuesr4(i))
             end do
          end do
      
diff --git a/src/gsi/cplr_get_pseudo_ensperts.f90 b/src/gsi/cplr_get_pseudo_ensperts.f90
index 8c59cd0be8..9d12165409 100644
--- a/src/gsi/cplr_get_pseudo_ensperts.f90
+++ b/src/gsi/cplr_get_pseudo_ensperts.f90
@@ -57,7 +57,7 @@ subroutine get_pseudo_ensperts_wrf(this,en_perts,nelen)
     implicit none
 !   Declare passed variables
     class(get_pseudo_ensperts_class),  intent(inout) :: this
-    type(gsi_bundle),allocatable, intent(in   ) :: en_perts(:,:)
+    type(gsi_bundle),allocatable, intent(in   ) :: en_perts(:,:,:)
     integer(i_kind),                   intent(in   ) :: nelen
     
 
@@ -600,12 +600,12 @@ subroutine get_pseudo_ensperts_wrf(this,en_perts,nelen)
   ! NOTE:  because library perturbation bundle structure is same as ensemble perturbation, 
   !        use same ipc3d and ipc2d indices for lib_perts and en_perts bundles.
   
-    call gsi_bundlegetpointer (en_perts(1,1),cvars3d,ipc3d,istatus)
+    call gsi_bundlegetpointer (en_perts(1,1,1),cvars3d,ipc3d,istatus)
     if(istatus/=0) then
        write(6,*) 'get_pseudo_enperts: cannot find 3d pointers'
        call stop2(999)
     endif
-    call gsi_bundlegetpointer (en_perts(1,1),cvars2d,ipc2d,istatus)
+    call gsi_bundlegetpointer (en_perts(1,1,1),cvars2d,ipc2d,istatus)
     if(istatus/=0) then
        write(6,*) 'get_pseudo_enperts: cannot find 2d pointers'
        call stop2(999)
@@ -618,9 +618,9 @@ subroutine get_pseudo_ensperts_wrf(this,en_perts,nelen)
           do k=1,grd_ens%nsig
              do j=1,grd_ens%lon2
                 do i=1,grd_ens%lat2
-                   en_perts(n,1)%r3(ipc3d(ic3))%qr4(i,j,k) = wgt(i,j) * &
+                   en_perts(n,1,1)%r3(ipc3d(ic3))%qr4(i,j,k) = wgt(i,j) * &
                                 lib_perts(n)%r3(ipc3d(ic3))%qr4(i,j,k)  + &
-                                (one-wgt(i,j))*en_perts(n,1)%r3(ipc3d(ic3))%qr4(i,j,k)
+                                (one-wgt(i,j))*en_perts(n,1,1)%r3(ipc3d(ic3))%qr4(i,j,k)
                 end do
              end do
           end do
@@ -631,9 +631,9 @@ subroutine get_pseudo_ensperts_wrf(this,en_perts,nelen)
   
           do j=1,grd_ens%lon2
              do i=1,grd_ens%lat2
-                en_perts(n,1)%r2(ipc2d(ic2))%qr4(i,j) = wgt(i,j)  * &
+                en_perts(n,1,1)%r2(ipc2d(ic2))%qr4(i,j) = wgt(i,j)  * &
                                lib_perts(n)%r2(ipc2d(ic2))%qr4(i,j) + &
-                               (one-wgt(i,j))*en_perts(n,1)%r2(ipc2d(ic2))%qr4(i,j)
+                               (one-wgt(i,j))*en_perts(n,1,1)%r2(ipc2d(ic2))%qr4(i,j)
              end do
           end do
   
@@ -667,12 +667,12 @@ subroutine get_pseudo_ensperts_wrf(this,en_perts,nelen)
   
           write(filename,"('enpert',i3.3)") n
   
-          call gsi_bundlegetpointer(en_perts(n,1),cvars3d,ipc3d,istatus)
+          call gsi_bundlegetpointer(en_perts(n,1,1),cvars3d,ipc3d,istatus)
           if(istatus/=0) then
             write(6,*) 'mtong: cannot find 3d pointers'
             call stop2(999)
           endif
-          call gsi_bundlegetpointer(en_perts(n,1),cvars2d,ipc2d,istatus)
+          call gsi_bundlegetpointer(en_perts(n,1,1),cvars2d,ipc2d,istatus)
           if(istatus/=0) then
             write(6,*) 'mtong: cannot find 2d pointers'
             call stop2(999)
@@ -686,7 +686,7 @@ subroutine get_pseudo_ensperts_wrf(this,en_perts,nelen)
                    do k=1,grd_ens%nsig
                       do j=1,grd_ens%lon2
                          do i=1,grd_ens%lat2
-                            u(i,j,k) = en_perts(n,1)%r3(ipc3d(ic3))%qr4(i,j,k)
+                            u(i,j,k) = en_perts(n,1,1)%r3(ipc3d(ic3))%qr4(i,j,k)
                          end do
                       end do
                    end do
@@ -696,7 +696,7 @@ subroutine get_pseudo_ensperts_wrf(this,en_perts,nelen)
                    do k=1,grd_ens%nsig
                       do j=1,grd_ens%lon2
                          do i=1,grd_ens%lat2
-                            v(i,j,k) = en_perts(n,1)%r3(ipc3d(ic3))%qr4(i,j,k)
+                            v(i,j,k) = en_perts(n,1,1)%r3(ipc3d(ic3))%qr4(i,j,k)
                          end do
                       end do
                    end do
@@ -706,7 +706,7 @@ subroutine get_pseudo_ensperts_wrf(this,en_perts,nelen)
                    do k=1,grd_ens%nsig
                       do j=1,grd_ens%lon2
                          do i=1,grd_ens%lat2
-                            tv(i,j,k) = en_perts(n,1)%r3(ipc3d(ic3))%qr4(i,j,k)
+                            tv(i,j,k) = en_perts(n,1,1)%r3(ipc3d(ic3))%qr4(i,j,k)
                          end do
                       end do
                    end do
@@ -716,7 +716,7 @@ subroutine get_pseudo_ensperts_wrf(this,en_perts,nelen)
                    do k=1,grd_ens%nsig
                       do j=1,grd_ens%lon2
                          do i=1,grd_ens%lat2
-                            rh(i,j,k) = en_perts(n,1)%r3(ipc3d(ic3))%qr4(i,j,k)
+                            rh(i,j,k) = en_perts(n,1,1)%r3(ipc3d(ic3))%qr4(i,j,k)
                          end do
                       end do
                    end do
@@ -731,7 +731,7 @@ subroutine get_pseudo_ensperts_wrf(this,en_perts,nelen)
   
                    do j=1,grd_ens%lon2
                       do i=1,grd_ens%lat2
-                         ps(i,j) = en_perts(n,1)%r2(ipc2d(ic2))%qr4(i,j)
+                         ps(i,j) = en_perts(n,1,1)%r2(ipc2d(ic2))%qr4(i,j)
                       end do
                    end do
   
diff --git a/src/gsi/cplr_get_wrf_mass_ensperts.f90 b/src/gsi/cplr_get_wrf_mass_ensperts.f90
index eb94c6c0f3..775c9c4195 100644
--- a/src/gsi/cplr_get_wrf_mass_ensperts.f90
+++ b/src/gsi/cplr_get_wrf_mass_ensperts.f90
@@ -63,7 +63,7 @@ subroutine get_wrf_mass_ensperts_wrf(this,en_perts,nelen,ps_bar)
 
       implicit none
       class(get_wrf_mass_ensperts_class), intent(inout) :: this
-      type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:)
+      type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:,:)
       integer(i_kind), intent(in   ):: nelen
       real(r_single),dimension(:,:,:),allocatable:: ps_bar
   
@@ -148,7 +148,7 @@ subroutine get_wrf_mass_ensperts_wrf(this,en_perts,nelen,ps_bar)
          en_bar%values=zero
   
          do n=1,n_ens
-            en_perts(n,it)%valuesr4 = zero
+            en_perts(n,1,it)%valuesr4 = zero
          enddo
   
          if ( .not. l_use_dbz_directDA ) then !direct reflectivity DA option does not employ this 
@@ -253,7 +253,7 @@ subroutine get_wrf_mass_ensperts_wrf(this,en_perts,nelen,ps_bar)
   ! SAVE ENSEMBLE MEMBER DATA IN COLUMN VECTOR
             member_data_loop: do ic3=1,nc3d
   
-               call gsi_bundlegetpointer(en_perts(n,it),trim(cvars3d(ic3)),w3,istatus)
+               call gsi_bundlegetpointer(en_perts(n,1,it),trim(cvars3d(ic3)),w3,istatus)
                if(istatus/=0) then
                   write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for ensemble member ',n
                   call stop2(999)
@@ -442,7 +442,7 @@ subroutine get_wrf_mass_ensperts_wrf(this,en_perts,nelen,ps_bar)
   
             member_mass_loop: do ic2=1,nc2d
      
-               call gsi_bundlegetpointer(en_perts(n,it),trim(cvars2d(ic2)),w2,istatus)
+               call gsi_bundlegetpointer(en_perts(n,1,it),trim(cvars2d(ic2)),w2,istatus)
                if(istatus/=0) then
                   write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for ensemble member ',n
                   call stop2(999)
@@ -515,7 +515,7 @@ subroutine get_wrf_mass_ensperts_wrf(this,en_perts,nelen,ps_bar)
   
          do n=1,n_ens
             do i=1,nelen
-               en_perts(n,it)%valuesr4(i)=(en_perts(n,it)%valuesr4(i)-en_bar%values(i))*sig_norm
+               en_perts(n,1,it)%valuesr4(i)=(en_perts(n,1,it)%valuesr4(i)-en_bar%values(i))*sig_norm
             end do
          end do
 
@@ -2147,7 +2147,7 @@ subroutine ens_spread_dualres_regional_wrf(this,mype,en_perts,nelen,en_bar)
     class(get_wrf_mass_ensperts_class), intent(inout) :: this
     type(gsi_bundle),OPTIONAL,intent(in):: en_bar
     integer(i_kind),intent(in):: mype
-    type(gsi_bundle),allocatable, intent(in   ) :: en_perts(:,:)
+    type(gsi_bundle),allocatable, intent(in   ) :: en_perts(:,:,:)
     integer(i_kind), intent(in   ):: nelen
   
     type(gsi_bundle):: sube,suba
@@ -2201,7 +2201,7 @@ subroutine ens_spread_dualres_regional_wrf(this,mype,en_perts,nelen,en_bar)
        do n=1,n_ens
           do i=1,nelen
              sube%values(i)=sube%values(i) &
-               +en_perts(n,1)%valuesr4(i)*en_perts(n,1)%valuesr4(i)
+               +en_perts(n,1,1)%valuesr4(i)*en_perts(n,1,1)%valuesr4(i)
           end do
        end do
   
@@ -2212,7 +2212,7 @@ subroutine ens_spread_dualres_regional_wrf(this,mype,en_perts,nelen,en_bar)
        do n=1,n_ens
           do i=1,nelen
              sube%values(i)=sube%values(i) &
-               +(en_perts(n,1)%valuesr4(i)-en_bar%values(i))*(en_perts(n,1)%valuesr4(i)-en_bar%values(i))
+               +(en_perts(n,1,1)%valuesr4(i)-en_bar%values(i))*(en_perts(n,1,1)%valuesr4(i)-en_bar%values(i))
           end do
        end do
    
diff --git a/src/gsi/cplr_get_wrf_nmm_ensperts.f90 b/src/gsi/cplr_get_wrf_nmm_ensperts.f90
index 9de32ea384..3e30822639 100644
--- a/src/gsi/cplr_get_wrf_nmm_ensperts.f90
+++ b/src/gsi/cplr_get_wrf_nmm_ensperts.f90
@@ -72,7 +72,7 @@ subroutine get_wrf_nmm_ensperts_wrf(this,en_perts,nelen,region_lat_ens,region_lo
  
       implicit none
       class(get_wrf_nmm_ensperts_class), intent(inout) :: this
-      type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:)
+      type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:,:)
       integer(i_kind), intent(in   ):: nelen
       real(r_kind),allocatable, intent(inout):: region_lat_ens(:,:),region_lon_ens(:,:)
       real(r_single),dimension(:,:,:),allocatable, intent(inout):: ps_bar
@@ -689,7 +689,7 @@ subroutine get_wrf_nmm_ensperts_wrf(this,en_perts,nelen,region_lat_ens,region_lo
   ! SAVE ENSEMBLE MEMBER DATA IN COLUMN VECTOR
          do ic3=1,nc3d
   
-            call gsi_bundlegetpointer(en_perts(n,1),trim(cvars3d(ic3)),w3,istatus)
+            call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars3d(ic3)),w3,istatus)
             if(istatus/=0) then
                write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for ensemble member ',n
                call stop2(999)
@@ -791,7 +791,7 @@ subroutine get_wrf_nmm_ensperts_wrf(this,en_perts,nelen,region_lat_ens,region_lo
   
          do ic2=1,nc2d
   
-            call gsi_bundlegetpointer(en_perts(n,1),trim(cvars2d(ic2)),w2,istatus)
+            call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars2d(ic2)),w2,istatus)
             if(istatus/=0) then
                write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for ensemble member ',n
                call stop2(999)
@@ -828,7 +828,7 @@ subroutine get_wrf_nmm_ensperts_wrf(this,en_perts,nelen,region_lat_ens,region_lo
       if (use_gfs_stratosphere)then
           do n=1,n_ens
              do ic3=1,nc3d
-                call gsi_bundlegetpointer(en_perts(n,1),trim(cvars3d(ic3)),w3,istatus)
+                call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars3d(ic3)),w3,istatus)
                 if(istatus/=0) then
                    write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' &
                    for ensemble member ',n
@@ -893,7 +893,7 @@ subroutine get_wrf_nmm_ensperts_wrf(this,en_perts,nelen,region_lat_ens,region_lo
   
       do n=1,n_ens
          do i=1,nelen
-            en_perts(n,1)%valuesr4(i)=(en_perts(n,1)%valuesr4(i)-en_bar%values(i))*sig_norm
+            en_perts(n,1,1)%valuesr4(i)=(en_perts(n,1,1)%valuesr4(i)-en_bar%values(i))*sig_norm
          end do
       end do
   
@@ -2917,7 +2917,7 @@ subroutine ens_member_mean_dualres_regional(this,en_bar,mype,en_perts,nelen)
     class(get_wrf_nmm_ensperts_class), intent(inout) :: this
     type(gsi_bundle),intent(in):: en_bar
     integer(i_kind),intent(in):: mype
-    type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:)
+    type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:,:)
     integer(i_kind), intent(in   ):: nelen
   
     type(gsi_bundle):: sube,suba
@@ -2961,7 +2961,7 @@ subroutine ens_member_mean_dualres_regional(this,en_bar,mype,en_perts,nelen)
   
        do i=1,nelen
           if(n <= n_ens)then
-             sube%values(i)=en_perts(n,1)%valuesr4(i)*sig_norm_inv+en_bar%values(i)
+             sube%values(i)=en_perts(n,1,1)%valuesr4(i)*sig_norm_inv+en_bar%values(i)
           else
              sube%values(i)=en_bar%values(i)
           end if
diff --git a/src/gsi/cplr_gfs_ensmod.f90 b/src/gsi/cplr_gfs_ensmod.f90
index 582574b52f..5c9934919b 100644
--- a/src/gsi/cplr_gfs_ensmod.f90
+++ b/src/gsi/cplr_gfs_ensmod.f90
@@ -440,7 +440,7 @@ subroutine move2bundle_(grd3d,en_loc3,atm_bundle,m_cvars2d,m_cvars3d,iret)
 !      if(trim(cvars2d(m))=='sst') sst=en_loc3(:,:,m_cvars2d(m)) !no sst for now
     enddo
 
-    km = en_perts(1,1)%grid%km
+    km = en_perts(1,1,1)%grid%km
 !$omp parallel do  schedule(dynamic,1) private(m) 
     do m=1,nc3d
        if(trim(cvars3d(m))=='sf')then
diff --git a/src/gsi/en_perts_io.f90 b/src/gsi/en_perts_io.f90
index 1a9dc58a50..f873a8520f 100644
--- a/src/gsi/en_perts_io.f90
+++ b/src/gsi/en_perts_io.f90
@@ -124,7 +124,7 @@ subroutine en_perts_get_from_save_fulldomain
 
      do ic3=1,nc3d
 
-        call gsi_bundlegetpointer(en_perts(n,1),trim(cvars3d(ic3)),w3,istatus)
+        call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars3d(ic3)),w3,istatus)
         if(istatus/=0) then
            write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for ensemble member ',n
            call stop2(999)
@@ -142,7 +142,7 @@ subroutine en_perts_get_from_save_fulldomain
 
      do ic2=1,nc2d
 
-        call gsi_bundlegetpointer(en_perts(n,1),trim(cvars2d(ic2)),w2,istatus)
+        call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars2d(ic2)),w2,istatus)
         if(istatus/=0) then
            write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for ensemble member ',n
            call stop2(999)
@@ -221,7 +221,7 @@ subroutine en_perts_get_from_save
 !
      do ic3=1,nc3d
 
-        call gsi_bundlegetpointer(en_perts(n,1),trim(cvars3d(ic3)),w3,istatus)
+        call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars3d(ic3)),w3,istatus)
         if(istatus/=0) then
            write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for ensemble member ',n
            call stop2(999)
@@ -238,7 +238,7 @@ subroutine en_perts_get_from_save
 
      do ic2=1,nc2d
 
-        call gsi_bundlegetpointer(en_perts(n,1),trim(cvars2d(ic2)),w2,istatus)
+        call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars2d(ic2)),w2,istatus)
         if(istatus/=0) then
            write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for ensemble member ',n
            call stop2(999)
@@ -289,6 +289,7 @@ subroutine en_perts_save
   use kinds, only: r_kind,i_kind,r_single
   use gsi_bundlemod, only: GSI_BundleGetPointer
   use mpeu_util, only: die
+  use hybrid_ensemble_parameters, only: nsclgrp
   implicit none
 
   real(r_single),pointer,dimension(:,:,:):: w3
@@ -300,6 +301,10 @@ subroutine en_perts_save
 
   integer(i_kind) ic3,ic2
   integer(i_kind) iunit
+  if(nsclgrp>1) then
+    write(6,*)"nsclgrp >1 is not considerred in this part, stop"
+    stop
+  endif
 
   iunit=20
   write(filename,'(a,I4.4)') 'saved_en_perts.pe',mype
@@ -312,7 +317,7 @@ subroutine en_perts_save
 !
      do ic3=1,nc3d
 
-        call gsi_bundlegetpointer(en_perts(n,1),trim(cvars3d(ic3)),w3,istatus)
+        call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars3d(ic3)),w3,istatus)
         if(istatus/=0) then
            write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for ensemble member ',n
            call stop2(999)
@@ -324,7 +329,7 @@ subroutine en_perts_save
      end do
      do ic2=1,nc2d
 
-        call gsi_bundlegetpointer(en_perts(n,1),trim(cvars2d(ic2)),w2,istatus)
+        call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars2d(ic2)),w2,istatus)
         if(istatus/=0) then
            write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for ensemble member ',n
            call stop2(999)
diff --git a/src/gsi/ens_spread_mod.f90 b/src/gsi/ens_spread_mod.f90
index afab9737ac..da70d31ab0 100644
--- a/src/gsi/ens_spread_mod.f90
+++ b/src/gsi/ens_spread_mod.f90
@@ -90,7 +90,7 @@ subroutine ens_spread_dualres_regional(en_bar)
      do n=1,n_ens
         do i=1,nelen
            sube%values(i)=sube%values(i) &
-             +en_perts(n,1)%valuesr4(i)*en_perts(n,1)%valuesr4(i)
+             +en_perts(n,1,1)%valuesr4(i)*en_perts(n,1,1)%valuesr4(i)
         end do
      end do
 
@@ -101,7 +101,7 @@ subroutine ens_spread_dualres_regional(en_bar)
      do n=1,n_ens
         do i=1,nelen
            sube%values(i)=sube%values(i) &
-             +(en_perts(n,1)%valuesr4(i)-en_bar%values(i))*(en_perts(n,1)%valuesr4(i)-en_bar%values(i))
+             +(en_perts(n,1,1)%valuesr4(i)-en_bar%values(i))*(en_perts(n,1,1)%valuesr4(i)-en_bar%values(i))
         end do
      end do
  
diff --git a/src/gsi/ensctl2model.f90 b/src/gsi/ensctl2model.f90
index 525bf33af6..12e1fe374e 100644
--- a/src/gsi/ensctl2model.f90
+++ b/src/gsi/ensctl2model.f90
@@ -43,6 +43,7 @@ subroutine ensctl2model(xhat,mval,eval)
 use gsi_metguess_mod, only: gsi_metguess_get
 use mod_strong, only: tlnmc_option
 use timermod, only: timer_ini,timer_fnl
+use hybrid_ensemble_parameters,only: naensgrp
 implicit none
 
 ! Declare passed variables
@@ -58,7 +59,7 @@ subroutine ensctl2model(xhat,mval,eval)
 integer(i_kind), parameter :: ncvars = 5
 integer(i_kind) :: icps(ncvars)
 type(gsi_bundle):: wbundle_c ! work bundle
-type(gsi_bundle),allocatable :: ebundle(:)
+type(gsi_bundle),allocatable :: ebundle(:,:)
 character(len=3), parameter :: mycvars(ncvars) = (/  &  ! vars from CV needed here
                                'sf ', 'vp ', 'ps ', 't  ',    &
                                'q  '/)
@@ -77,6 +78,7 @@ subroutine ensctl2model(xhat,mval,eval)
 
 logical :: do_getprs_tl,do_normal_rh_to_q,do_tv_to_tsen,do_getuv,lstrong_bk_vars
 logical :: do_tlnmc,do_q_copy
+integer(i_kind) :: ig
 ! ****************************************************************************
 
 ! Initialize timer
@@ -117,17 +119,21 @@ subroutine ensctl2model(xhat,mval,eval)
    do_tlnmc = lstrong_bk_vars .and. ( (tlnmc_option==3) .or. &
          (jj==ibin_anl .and. tlnmc_option==2) )
 
-   allocate(ebundle(n_ens))
-   do nn=1,n_ens
-      call gsi_bundlecreate (ebundle(nn),xhat%aens(1,1),'c2m ensemble work',istatus)
-      if(istatus/=0) then
-         write(6,*) trim(myname), ': trouble creating work ens-bundle'
-         call stop2(999)
-      endif
+   allocate(ebundle(naensgrp,n_ens))
+   do ig=1,naensgrp
+      do nn=1,n_ens
+         call gsi_bundlecreate (ebundle(ig,nn),xhat%aens(1,1,1),'c2m ensemble work',istatus)
+         if(istatus/=0) then
+            write(6,*) trim(myname), ': trouble creating work ens-bundle'
+            call stop2(999)
+         endif
+      enddo
    enddo
 
 !  Apply square-root of ensemble error covariance
-   call ckgcov_a_en_new_factorization(xhat%aens(jj,1)%values(:),ebundle)
+   do ig=1,naensgrp
+      call ckgcov_a_en_new_factorization(ig,xhat%aens(jj,ig,1)%values(:),ebundle(ig,:))
+   enddo
    call sqrt_beta_e_mult(ebundle)
 
 ! Initialize ensemble contribution to zero
@@ -230,12 +236,14 @@ subroutine ensctl2model(xhat,mval,eval)
       call stop2(999)
    endif
 
-   do nn=n_ens,1,-1 ! first in; last out
-      call gsi_bundledestroy(ebundle(nn),istatus)
-      if(istatus/=0) then
-         write(6,*) trim(myname), ': trouble destroying work ens bundle, ', istatus
-         call stop2(999)
-      endif
+   do ig=1,naensgrp
+      do nn=n_ens,1,-1 ! first in; last out
+         call gsi_bundledestroy(ebundle(ig,nn),istatus)
+         if(istatus/=0) then
+            write(6,*) trim(myname), ': trouble destroying work ens bundle, ', istatus
+            call stop2(999)
+         endif
+      enddo
    enddo
    deallocate(ebundle)
 
diff --git a/src/gsi/ensctl2model_ad.f90 b/src/gsi/ensctl2model_ad.f90
index 8c491b6de8..706dafc59c 100644
--- a/src/gsi/ensctl2model_ad.f90
+++ b/src/gsi/ensctl2model_ad.f90
@@ -42,6 +42,7 @@ subroutine ensctl2model_ad(eval,mval,grad)
 use gsi_metguess_mod, only: gsi_metguess_get
 use mod_strong, only: tlnmc_option
 use timermod, only: timer_ini,timer_fnl
+use hybrid_ensemble_parameters,only: naensgrp
 implicit none
 
 ! Declare passed variables
@@ -57,7 +58,7 @@ subroutine ensctl2model_ad(eval,mval,grad)
 integer(i_kind), parameter :: ncvars = 5
 integer(i_kind) :: icps(ncvars)
 type(gsi_bundle):: wbundle_c ! work bundle
-type(gsi_bundle),allocatable :: ebundle(:)
+type(gsi_bundle),allocatable :: ebundle(:,:)
 real(r_kind) :: grade(nval_lenz_en)
 character(len=3), parameter :: mycvars(ncvars) = (/  &  ! vars from CV needed here
                                'sf ', 'vp ', 'ps ', 't  ',    &
@@ -77,7 +78,7 @@ subroutine ensctl2model_ad(eval,mval,grad)
 
 logical :: do_getuv,do_tv_to_tsen_ad,do_normal_rh_to_q_ad,do_getprs_ad
 logical :: do_tlnmc,lstrong_bk_vars,do_q_copy
-
+integer(i_kind) :: ig
 !****************************************************************************
 
 ! Initialize timer
@@ -122,13 +123,15 @@ subroutine ensctl2model_ad(eval,mval,grad)
          (jj==ibin_anl .and. tlnmc_option==2) )
 
    !allocate(grade(nval_lenz_en))
-   allocate(ebundle(n_ens))
-   do nn=1,n_ens
-      call gsi_bundlecreate (ebundle(nn),grad%aens(1,1),'m2c ensemble work',istatus)
-      if(istatus/=0) then
-         write(6,*) trim(myname), ': trouble creating work ens-bundle'
-         call stop2(999)
-      endif
+   allocate(ebundle(naensgrp,n_ens))
+   do ig=1,naensgrp
+      do nn=1,n_ens
+         call gsi_bundlecreate (ebundle(ig,nn),grad%aens(1,1,1),'m2c ensemble work',istatus)
+         if(istatus/=0) then
+            write(6,*) trim(myname), ': trouble creating work ens-bundle'
+            call stop2(999)
+         endif
+      enddo
    enddo
 
 !  Create a temporary bundle similar to grad, and copy contents of grad into it
@@ -210,8 +213,10 @@ subroutine ensctl2model_ad(eval,mval,grad)
       if(do_getprs_ad) call getprs_ad(cv_ps,cv_tv,rv_prse)
    end if
 
-   do nn=1,n_ens
-      ebundle(nn)%values=grad%aens(jj,nn)%values
+   do ig=1,naensgrp
+      do nn=1,n_ens
+         ebundle(ig,nn)%values=grad%aens(jj,ig,nn)%values
+      enddo
    enddo
    if(dual_res) then
       call ensemble_forward_model_ad_dual_res(wbundle_c,ebundle,jj)
@@ -222,7 +227,23 @@ subroutine ensctl2model_ad(eval,mval,grad)
 
 !  Apply square-root of ensemble error covariance
    call sqrt_beta_e_mult(ebundle)
-   call ckgcov_a_en_new_factorization_ad(grade,ebundle)
+   do ig=1,naensgrp
+      call ckgcov_a_en_new_factorization_ad(ig,grade,ebundle(ig,:))
+
+      do ii=1,n_ens
+         grad%aens(jj,ig,ii)%values=grad%aens(jj,ig,ii)%values + grade
+      enddo
+   enddo
+
+   do ig=1,naensgrp
+      do nn=n_ens,1,-1 ! first in; last out
+         call gsi_bundledestroy(ebundle(ig,nn),istatus)
+         if(istatus/=0) then
+            write(6,*) trim(myname), ': trouble destroying work ens bundle', ig, nn, istatus
+            call stop2(999)
+         endif
+      enddo
+   enddo
 
    call gsi_bundledestroy(wbundle_c,istatus)
    if (istatus/=0) then
@@ -230,20 +251,6 @@ subroutine ensctl2model_ad(eval,mval,grad)
       call stop2(999)
    endif
 
-   do ii=1,n_ens
-      grad%aens(jj,ii)%values=grad%aens(jj,ii)%values + grade
-   enddo
-!  do ii=1,nval_lenz_en
-!     grad%aens(jj,1)%values(ii)=grad%aens(jj,1)%values(ii)+grade(ii)
-!  enddo
-
-   do nn=n_ens,1,-1 ! first in; last out
-      call gsi_bundledestroy(ebundle(nn),istatus)
-      if(istatus/=0) then
-         write(6,*) trim(myname), ': trouble destroying work ens bundle', nn, istatus
-         call stop2(999)
-      endif
-   enddo
    deallocate(ebundle)
 !  deallocate(grade)
 
diff --git a/src/gsi/ensctl2state.f90 b/src/gsi/ensctl2state.f90
index 4214b01d85..0d6d3042c5 100644
--- a/src/gsi/ensctl2state.f90
+++ b/src/gsi/ensctl2state.f90
@@ -165,9 +165,9 @@ subroutine ensctl2state(xhat,mval,eval)
 !  For 4densvar, this is the "3D/Time-invariant contribution from static B"
 
    if(dual_res) then
-      call ensemble_forward_model_dual_res(wbundle_c,xhat%aens(1,:),jj)
+      call ensemble_forward_model_dual_res(wbundle_c,xhat%aens(1,:,:),jj)
    else
-      call ensemble_forward_model(wbundle_c,xhat%aens(1,:),jj)
+      call ensemble_forward_model(wbundle_c,xhat%aens(1,:,:),jj)
    end if
 
 !  Get pointers to required state variables
diff --git a/src/gsi/ensctl2state_ad.f90 b/src/gsi/ensctl2state_ad.f90
index 212d0312c9..4c038c8c6e 100644
--- a/src/gsi/ensctl2state_ad.f90
+++ b/src/gsi/ensctl2state_ad.f90
@@ -260,9 +260,9 @@ subroutine ensctl2state_ad(eval,mval,grad)
 !$omp end parallel sections
 
    if(dual_res) then
-      call ensemble_forward_model_ad_dual_res(wbundle_c,grad%aens(1,:),jj)
+      call ensemble_forward_model_ad_dual_res(wbundle_c,grad%aens(1,:,:),jj)
    else
-      call ensemble_forward_model_ad(wbundle_c,grad%aens(1,:),jj)
+      call ensemble_forward_model_ad(wbundle_c,grad%aens(1,:,:),jj)
    end if
 
 end do
diff --git a/src/gsi/get_gefs_ensperts_dualres.f90 b/src/gsi/get_gefs_ensperts_dualres.f90
index c547a9f062..fa3d0ecbdd 100644
--- a/src/gsi/get_gefs_ensperts_dualres.f90
+++ b/src/gsi/get_gefs_ensperts_dualres.f90
@@ -93,32 +93,32 @@ subroutine get_gefs_ensperts_dualres
   type(sub2grid_info) :: grd_tmp
 
 ! Create perturbations grid and get variable names from perturbations
-  if(en_perts(1,1)%grid%im/=grd_ens%lat2.or. &
-     en_perts(1,1)%grid%jm/=grd_ens%lon2.or. &
-     en_perts(1,1)%grid%km/=grd_ens%nsig ) then
+  if(en_perts(1,1,1)%grid%im/=grd_ens%lat2.or. &
+     en_perts(1,1,1)%grid%jm/=grd_ens%lon2.or. &
+     en_perts(1,1,1)%grid%km/=grd_ens%nsig ) then
      if (mype==0) then
         write(6,*) 'get_gefs_ensperts_dualres: grd_ens ', grd_ens%lat2,grd_ens%lon2,grd_ens%nsig
-        write(6,*) 'get_gefs_ensperts_dualres: pertgrid', en_perts(1,1)%grid%im, en_perts(1,1)%grid%jm, en_perts(1,1)%grid%km
+        write(6,*) 'get_gefs_ensperts_dualres: pertgrid', en_perts(1,1,1)%grid%im, en_perts(1,1,1)%grid%jm, en_perts(1,1,1)%grid%km
         write(6,*) 'get_gefs_ensperts_dualres: inconsistent dims, aborting ...'
      endif
      call stop2(999)
  endif
 
-  call gsi_bundlegetpointer (en_perts(1,1),cvars3d,ipc3d,istatus)
+  call gsi_bundlegetpointer (en_perts(1,1,1),cvars3d,ipc3d,istatus)
   if(istatus/=0) then
     write(6,*) ' get_gefs_ensperts_dualres',': cannot find 3d pointers'
     call stop2(999)
   endif
-  call gsi_bundlegetpointer (en_perts(1,1),cvars2d,ipc2d,istatus)
+  call gsi_bundlegetpointer (en_perts(1,1,1),cvars2d,ipc2d,istatus)
   if(istatus/=0) then
     write(6,*) ' get_gefs_ensperts_dualres',': cannot find 2d pointers'
     call stop2(999)
   endif
 
 
-  im=en_perts(1,1)%grid%im
-  jm=en_perts(1,1)%grid%jm
-  km=en_perts(1,1)%grid%km
+  im=en_perts(1,1,1)%grid%im
+  jm=en_perts(1,1,1)%grid%jm
+  km=en_perts(1,1,1)%grid%km
   bar_norm = one/float(n_ens)
   sig_norm=sqrt(one/max(one,n_ens-one))
 
@@ -126,14 +126,14 @@ subroutine get_gefs_ensperts_dualres
   call gsi_enscoupler_create_sub2grid_info(grd_tmp,km,npe,grd_ens)
 
   ! Allocate bundle to hold mean of ensemble members
-  call gsi_bundlecreate(en_bar,en_perts(1,1)%grid,'ensemble',istatus,names2d=cvars2d,names3d=cvars3d)
+  call gsi_bundlecreate(en_bar,en_perts(1,1,1)%grid,'ensemble',istatus,names2d=cvars2d,names3d=cvars3d)
   if ( istatus /= 0 ) &
      call die('get_gefs_ensperts_dualres',': trouble creating en_bar bundle, istatus =',istatus)
 
   ! Allocate bundle used for reading members
   allocate(en_read(n_ens))
   do n=1,n_ens
-     call gsi_bundlecreate(en_read(n),en_perts(1,1)%grid,'ensemble member',istatus,names2d=cvars2d,names3d=cvars3d)
+     call gsi_bundlecreate(en_read(n),en_perts(1,1,1)%grid,'ensemble member',istatus,names2d=cvars2d,names3d=cvars3d)
      if ( istatus /= 0 ) &
         call die('get_gefs_ensperts_dualres',': trouble creating en_read bundle, istatus =',istatus)
   end do
@@ -147,7 +147,7 @@ subroutine get_gefs_ensperts_dualres
 !$omp parallel do schedule(dynamic,1) private(m,n)
   do m=1,ntlevs_ens
      do n=1,n_ens
-       en_perts(n,m)%valuesr4=zero_single
+       en_perts(n,1,m)%valuesr4=zero_single
      end do
   end do
 
@@ -254,7 +254,7 @@ subroutine get_gefs_ensperts_dualres
 
        end do !c3d
        do i=1,nelen
-          en_perts(n,m)%valuesr4(i)=en_read(n)%values(i)
+          en_perts(n,1,m)%valuesr4(i)=en_read(n)%values(i)
           en_bar%values(i)=en_bar%values(i)+en_read(n)%values(i)
        end do
 
@@ -295,7 +295,7 @@ subroutine get_gefs_ensperts_dualres
 !$omp parallel do schedule(dynamic,1) private(n,i,ic3,ipic,k,j)
      do n=1,n_ens
         do i=1,nelen
-           en_perts(n,m)%valuesr4(i)=en_perts(n,m)%valuesr4(i)-en_bar%values(i)
+           en_perts(n,1,m)%valuesr4(i)=en_perts(n,1,m)%valuesr4(i)-en_bar%values(i)
         end do
         if(.not. q_hyb_ens) then
           do ic3=1,nc3d
@@ -304,8 +304,8 @@ subroutine get_gefs_ensperts_dualres
                 do k=1,km
                    do j=1,jm
                       do i=1,im
-                         en_perts(n,m)%r3(ipic)%qr4(i,j,k) = min(en_perts(n,m)%r3(ipic)%qr4(i,j,k),limqens)
-                         en_perts(n,m)%r3(ipic)%qr4(i,j,k) = max(en_perts(n,m)%r3(ipic)%qr4(i,j,k),-limqens)
+                         en_perts(n,1,m)%r3(ipic)%qr4(i,j,k) = min(en_perts(n,1,m)%r3(ipic)%qr4(i,j,k),limqens)
+                         en_perts(n,1,m)%r3(ipic)%qr4(i,j,k) = max(en_perts(n,1,m)%r3(ipic)%qr4(i,j,k),-limqens)
                       end do
                    end do
                 end do
@@ -313,7 +313,7 @@ subroutine get_gefs_ensperts_dualres
           end do
         end if
         do i=1,nelen
-           en_perts(n,m)%valuesr4(i)=en_perts(n,m)%valuesr4(i)*sig_norm
+           en_perts(n,1,m)%valuesr4(i)=en_perts(n,1,m)%valuesr4(i)*sig_norm
         end do
      end do
   end do  ntlevs_ens_loop !end do over bins
@@ -455,7 +455,7 @@ subroutine ens_spread_dualres(en_bar,ibin)
   do n=1,n_ens
      do i=1,nelen
         sube%values(i)=sube%values(i) &
-           +(en_perts(n,ibin)%valuesr4(i)-en_bar%values(i))*(en_perts(n,ibin)%valuesr4(i)-en_bar%values(i))
+           +(en_perts(n,1,ibin)%valuesr4(i)-en_bar%values(i))*(en_perts(n,1,ibin)%valuesr4(i)-en_bar%values(i))
      end do
   end do
 
diff --git a/src/gsi/get_gefs_for_regional.f90 b/src/gsi/get_gefs_for_regional.f90
index 958849a177..a076f0ccfd 100644
--- a/src/gsi/get_gefs_for_regional.f90
+++ b/src/gsi/get_gefs_for_regional.f90
@@ -1333,9 +1333,9 @@ subroutine get_gefs_for_regional
      do ic3=1,nc3d
 
         if(ntlevs_ens > 1) then
-           call gsi_bundlegetpointer(en_perts(n,it),trim(cvars3d(ic3)),w3,istatus)
+           call gsi_bundlegetpointer(en_perts(n,1,it),trim(cvars3d(ic3)),w3,istatus)
         else
-           call gsi_bundlegetpointer(en_perts(n,1),trim(cvars3d(ic3)),w3,istatus)
+           call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars3d(ic3)),w3,istatus)
         endif
         if(istatus/=0) then
            write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for ensemble member ',n
@@ -1422,9 +1422,9 @@ subroutine get_gefs_for_regional
      do ic2=1,nc2d
 
         if(ntlevs_ens > 1) then
-           call gsi_bundlegetpointer(en_perts(n,it),trim(cvars2d(ic2)),w2,istatus)
+           call gsi_bundlegetpointer(en_perts(n,1,it),trim(cvars2d(ic2)),w2,istatus)
         else
-           call gsi_bundlegetpointer(en_perts(n,1),trim(cvars2d(ic2)),w2,istatus)
+           call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars2d(ic2)),w2,istatus)
         endif
         if(istatus/=0) then
            write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for ensemble member ',n
diff --git a/src/gsi/get_nmmb_ensperts.f90 b/src/gsi/get_nmmb_ensperts.f90
index 93d23c837d..ece1780c03 100644
--- a/src/gsi/get_nmmb_ensperts.f90
+++ b/src/gsi/get_nmmb_ensperts.f90
@@ -59,7 +59,7 @@ subroutine get_nmmb_ensperts
    endif
    
    do n=1,n_ens
-      en_perts(n,1)%valuesr4=zero
+      en_perts(n,1,1)%valuesr4=zero
    end do
 
    en_bar%values=zero
@@ -116,7 +116,7 @@ subroutine get_nmmb_ensperts
 
       do ic3=1,nc3d
 
-         call gsi_bundlegetpointer(en_perts(n,1),trim(cvars3d(ic3)),w3,istatus)
+         call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars3d(ic3)),w3,istatus)
          if(istatus/=0) then
             write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for ensemble member ',n,' in get_nmmb_ensperts'
             call stop2(999)
@@ -277,7 +277,7 @@ subroutine get_nmmb_ensperts
 
       do ic2=1,nc2d
 
-         call gsi_bundlegetpointer(en_perts(n,1),trim(cvars2d(ic2)),w2,istatus)
+         call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars2d(ic2)),w2,istatus)
          if(istatus/=0) then
             write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for ensemble member ',n, ' in get_nmmb_ensperts'
             call stop2(999)
@@ -344,7 +344,7 @@ subroutine get_nmmb_ensperts
     
    do n=1,n_ens
       do i=1,nelen
-         en_perts(n,1)%valuesr4(i)=(en_perts(n,1)%valuesr4(i)-en_bar%values(i))*sig_norm      
+         en_perts(n,1,1)%valuesr4(i)=(en_perts(n,1,1)%valuesr4(i)-en_bar%values(i))*sig_norm      
       end do
    end do
         
diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90
index ff636bd70f..859d4cedb1 100644
--- a/src/gsi/gsimod.F90
+++ b/src/gsi/gsimod.F90
@@ -149,7 +149,9 @@ module gsimod
                          beta_s0,beta_e0,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, &
-                         l_ens_in_diff_time,ensemble_path,ens_fast_read,sst_staticB,limqens
+                         l_ens_in_diff_time,ensemble_path,ens_fast_read,sst_staticB,limqens, &
+                         ntotensgrp,nsclgrp,naensgrp,ngvarloc,ntlevs_ens,naensloc, &
+                         i_ensloccov4tim,i_ensloccov4var,i_ensloccov4scl,l_timloc_opt
   use hybrid_ensemble_parameters,only : l_both_fv3sar_gfs_ens,n_ens_gfs,n_ens_fv3sar
   use rapidrefresh_cldsurf_mod, only: init_rapidrefresh_cldsurf, &
                             dfi_radar_latent_heat_time_period,metar_impact_radius,&
@@ -497,6 +499,7 @@ module gsimod
 !                           2. fv3_regional =      .true.  
 !                           3. fv3_cmaq_regional = .true. 
 !                           4. berror_fv3_cmaq_regional = .true. 
+!  09-15-2022 yokota  - add scale/variable/time-dependent localization
 !
 !EOP
 !-------------------------------------------------------------------------
@@ -1309,9 +1312,19 @@ module gsimod
 !     beta_e0 - default weight given to ensemble background error covariance
 !               (if .not. readin_beta). if beta_e0<0, then it is set to
 !               1.-beta_s0 (this is the default)
-!     s_ens_h             - homogeneous isotropic horizontal ensemble localization scale (km)
-!     s_ens_v             - vertical localization scale (grid units for now)
-!                              s_ens_h, s_ens_v, and beta_s0 are tunable parameters.
+!     s_ens_h - horizontal localization correlation length of Gaussian exp(-0.5*(r/L)**2)
+!               (units of km), default = 2828.0
+!     s_ens_v - vertical localization correlation length of Gaussian exp(-0.5*(r/L)**2)
+!               (grid units if s_ens_v>=0, or units of ln(p) if s_ens_v<0), default = 30.0
+!                  in scale/variable/time-dependent localization (SDL/VDL/TDL),
+!                  localization length for i-th scale, j-th variable, and k-th time is
+!                     s_ens_[hv]( i + nsclgrp*(j-1) + nsclgrp*ngvarloc*(k-1) )
+!                        in SDL(nsclgrp>1),         i = 1(largest scale)  .. nsclgrp(smallest scale)
+!                        in VDL(ngvarloc=2),        j = 1(itracer<=10)    .. 2(itracer>=11)
+!                        in TDL(l_timloc_opt=true), k = 1(first time bin) .. ntlevs_ens(last time bin)
+!                  in SDL, scale separation length for i-th scale is also set here as
+!                     s_ens_[hv]( naensgrp+i ) - naensgrp is the total number of localization lengths for SDL/VDL/TDL
+!                        in applying SDL only horizontally, set s_ens_v(naensgrp+i)=0.0
 !     use_gfs_ens  - controls use of global ensemble: .t. use GFS (default); .f. uses user-defined ens
 !     readin_localization - flag to read (.true.)external localization information file
 !     readin_beta         - flag to read (.true.) the vertically varying beta parameters beta_s and beta_e
@@ -1349,14 +1362,27 @@ module gsimod
 !     ensemble_path - path to ensemble members; default './'
 !     ens_fast_read - read ensemble in parallel; default '.false.'
 !     sst_staticB - use only static background error covariance for SST statistic
-!              
-!                         
+!     nsclgrp - number of scale-dependent localization lengths
+!     l_timloc_opt - if true, then turn on time-dependent localization
+!     ngvarloc - number of variable-dependent localization lengths
+!     naensloc - total number of spatial localization lengths and scale separation lengths (should be naensgrp+nsclgrp-1)
+!     i_ensloccov4tim - flag of cross-temporal localization
+!                         =0: cross-temporal covariance is retained
+!                         =1: cross-temporal covariance is zero
+!     i_ensloccov4var - flag of cross-variable localization
+!                         =0: cross-variable covariance is retained
+!                         =1: cross-variable covariance is zero
+!     i_ensloccov4scl - flag of cross-scale localization
+!                         =0: cross-scale covariance is retained
+!                         =1: cross-scale covariance is zero
+!
   namelist/hybrid_ensemble/l_hyb_ens,uv_hyb_ens,q_hyb_ens,aniso_a_en,generate_ens,n_ens,l_both_fv3sar_gfs_ens,n_ens_gfs,n_ens_fv3sar,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,beta_e0,s_ens_h,s_ens_v,readin_localization,eqspace_ensgrid,readin_beta,&
                 grid_ratio_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,limqens
+                i_en_perts_io,l_ens_in_diff_time,ensemble_path,ens_fast_read,sst_staticB,limqens, &
+                nsclgrp,l_timloc_opt,ngvarloc,naensloc,i_ensloccov4tim,i_ensloccov4var,i_ensloccov4scl
 
 ! rapidrefresh_cldsurf (options for cloud analysis and surface 
 !                             enhancement for RR appilcation  ):
@@ -1807,6 +1833,16 @@ subroutine gsimain_initialize
      call stop2(137)
   endif
 
+  ntotensgrp=nsclgrp*ngvarloc
+  if(l_timloc_opt) then
+     naensgrp=ntotensgrp*ntlevs_ens
+  else
+     naensgrp=ntotensgrp
+  endif
+  if(naensloc<naensgrp+nsclgrp-1) naensloc=naensgrp+nsclgrp-1
+  if(mype==0) write(6,*) 'in gsimod: naensgrp,ntotensgrp,nsclgrp,ngvarloc,ntlevs_ens= ', &
+                          naensgrp,ntotensgrp,nsclgrp,ngvarloc,ntlevs_ens
+
   call gsi_4dcoupler_setservices(rc=ier)
   if(ier/=0) call die(myname_,'gsi_4dcoupler_setServices(), rc =',ier)
 
diff --git a/src/gsi/hybrid_ensemble_isotropic.F90 b/src/gsi/hybrid_ensemble_isotropic.F90
index 717815751b..3fa827f1f6 100644
--- a/src/gsi/hybrid_ensemble_isotropic.F90
+++ b/src/gsi/hybrid_ensemble_isotropic.F90
@@ -48,6 +48,7 @@ module hybrid_ensemble_isotropic
 !   2015-04-07  carley  - bug fix to allow grd_loc%nlat=grd_loc%nlon
 !   2016-05-13  parrish - remove beta12mult
 !   2018-02-15  wu      - add code for fv3_regional option
+!   2022-09-15  yokota  - add scale/variable/time-dependent localization
 !
 ! subroutines included:
 !   sub init_rf_z                         - initialize localization recursive filter (z direction)
@@ -158,19 +159,19 @@ module hybrid_ensemble_isotropic
 !      following variables are filter parameters for isotropic
 !       homogeneous localization of hybrid control variable a_en
 
-  real(r_kind),allocatable:: fmatz(:,:,:,:)
-  real(r_kind),allocatable:: fmat0z(:,:,:)
-  real(r_kind),allocatable:: fmatx(:,:,:,:,:)
-  real(r_kind),allocatable:: fmat0x(:,:,:,:)
-  real(r_kind),allocatable:: fmaty(:,:,:,:)
-  real(r_kind),allocatable:: fmat0y(:,:,:)
-  real(r_kind),allocatable:: znorm_new(:,:)
-  real(r_kind),allocatable:: xnorm_new(:,:,:)
-  real(r_kind),allocatable:: ynorm_new(:,:)
+  real(r_kind),allocatable:: fmatz(:,:,:,:,:)
+  real(r_kind),allocatable:: fmat0z(:,:,:,:)
+  real(r_kind),allocatable:: fmatx(:,:,:,:,:,:)
+  real(r_kind),allocatable:: fmat0x(:,:,:,:,:)
+  real(r_kind),allocatable:: fmaty(:,:,:,:,:)
+  real(r_kind),allocatable:: fmat0y(:,:,:,:)
+  real(r_kind),allocatable:: znorm_new(:,:,:)
+  real(r_kind),allocatable:: xnorm_new(:,:,:,:)
+  real(r_kind),allocatable:: ynorm_new(:,:,:)
   real(r_kind),allocatable:: psbar(:)
 
 ! Other local variables for horizontal/spectral localization
-  real(r_kind),allocatable,dimension(:,:)  :: spectral_filter,sqrt_spectral_filter
+  real(r_kind),allocatable,dimension(:,:,:)  :: spectral_filter,sqrt_spectral_filter
   integer(i_kind),allocatable,dimension(:) :: k_index
 
 !    following is for special subdomain to slab variables used when internally generating ensemble members
@@ -231,12 +232,13 @@ subroutine init_rf_z(z_len)
   use constants, only: half,one,rd_over_cp,zero,one_tenth,ten,two
   use hybrid_ensemble_parameters, only: grd_ens
   use hybrid_ensemble_parameters, only: ps_bar
+  use hybrid_ensemble_parameters, only: naensloc
 
   implicit none
 
-  real(r_kind)   ,intent(in) :: z_len(grd_ens%nsig)
+  real(r_kind)   ,intent(in) :: z_len(grd_ens%nsig,naensloc)
 
-  integer(i_kind) k,nxy,i,ii,jj,j,l
+  integer(i_kind) k,nxy,i,ii,jj,j,l,ig
   real(r_kind) aspect(nsig),p_interface(nsig+1),ln_p_int(nsig+1)
   real(r_kind) dlnp,kap1,kapr,d1,rnsig
   real(r_kind),dimension(:,:,:),allocatable:: fmatz_tmp
@@ -249,15 +251,17 @@ subroutine init_rf_z(z_len)
 
 !    use new factorization:
 
-  if(.not.allocated(fmatz))  allocate(fmatz(nxy,2,nsig,2))
-  if(.not.allocated(fmat0z)) allocate(fmat0z(nxy,nsig,2))
+  if(.not.allocated(fmatz))  allocate(fmatz(nxy,2,nsig,2,naensloc))
+  if(.not.allocated(fmat0z)) allocate(fmat0z(nxy,nsig,2,naensloc))
   allocate(fmatz_tmp(2,nsig,2),fmat0z_tmp(nsig,2))
+
+  do ig=1,naensloc
 !   for z_len < zero, use abs val z_len and assume localization scale is in units of ln(p)
-  if(maxval(z_len) > zero) then
+  if(maxval(z_len(:,ig)) > zero) then
 
 !  z_len is in grid units
      do k=1,nsig
-        aspect(k)=z_len(k)**2
+        aspect(k)=z_len(k,ig)**2
      enddo
 
      do i=1,nxy
@@ -265,13 +269,13 @@ subroutine init_rf_z(z_len)
        do l=1,2
          do k=1,nsig
            do j=1,2
-             fmatz(i,j,k,l)=fmatz_tmp(j,k,l)
+             fmatz(i,j,k,l,ig)=fmatz_tmp(j,k,l)
            enddo
          enddo
        enddo
        do l=1,2
          do k=1,nsig
-           fmat0z(i,k,l)=fmat0z_tmp(k,l)
+           fmat0z(i,k,l,ig)=fmat0z_tmp(k,l)
          enddo
        enddo
      enddo
@@ -320,7 +324,7 @@ subroutine init_rf_z(z_len)
 
            do k=1,nsig
               dlnp=abs(ln_p_int(k)-ln_p_int(k+1))
-              d1=abs(z_len(k))/dlnp
+              d1=abs(z_len(k,ig))/dlnp
               d1=min(rnsig,d1)
               aspect(k)=d1**2
 !!            if(mype == 0) write(400,'(" k, vertical localization in grid units for ln(p) scaling =",i4,f10.2,f10.2,f10.2)') &
@@ -331,18 +335,19 @@ subroutine init_rf_z(z_len)
            do l=1,2
              do k=1,nsig
                do j=1,2
-                 fmatz(i,j,k,l)=fmatz_tmp(j,k,l)
+                 fmatz(i,j,k,l,ig)=fmatz_tmp(j,k,l)
                enddo
              enddo
            enddo
            do l=1,2
              do k=1,nsig
-               fmat0z(i,k,l)=fmat0z_tmp(k,l)
+               fmat0z(i,k,l,ig)=fmat0z_tmp(k,l)
              enddo
            enddo
         enddo
      enddo
   end if
+  enddo !ig loop
   deallocate(fmatz_tmp,fmat0z_tmp)
   return
 
@@ -376,37 +381,41 @@ subroutine init_rf_x(x_len,kl)
 
   use hybrid_ensemble_parameters, only: grd_loc
   use hybrid_ensemble_parameters, only: region_dx_ens,region_dy_ens
+  use hybrid_ensemble_parameters, only: naensloc
   use constants, only: half
 
   implicit none
 
-  real(r_kind),intent(in   ) :: x_len(kl)
+  real(r_kind),intent(in   ) :: x_len(kl,naensloc)
   integer(i_kind),intent(in   ) :: kl
 
-  integer(i_kind) i,j,k,l,kk
+  integer(i_kind) i,j,k,l,kk,ig
   real(r_kind) aspect(grd_loc%nlon)
   real(r_kind) fmatc(2,grd_loc%nlon,2),fmat0c(grd_loc%nlon,2)
 
 !    use new factorization:
   if(allocated(fmatx)) deallocate(fmatx)
   if(allocated(fmat0x)) deallocate(fmat0x)
-  allocate(fmatx(grd_loc%nlat,2,grd_loc%nlon,2,kl),fmat0x(grd_loc%nlat,grd_loc%nlon,2,kl))
-  do k=1,kl
-     do i=1,grd_loc%nlat
-        do j=1,grd_loc%nlon
-           aspect(j)=(x_len(k)*region_dy_ens(grd_loc%nlat/2,grd_loc%nlon/2)/region_dx_ens(i,j))**2 ! only works for rotated lat-lon grids
-        enddo
-        call get_new_alpha_beta(aspect,grd_loc%nlon,fmatc,fmat0c)
-        do kk=1,2
+  allocate(fmatx(grd_loc%nlat,2,grd_loc%nlon,2,kl,naensloc),fmat0x(grd_loc%nlat,grd_loc%nlon,2,kl,naensloc))
+  do ig=1,naensloc
+     do k=1,kl
+        do i=1,grd_loc%nlat
            do j=1,grd_loc%nlon
-              do l=1,2
-                 fmatx(i,l,j,kk,k)=fmatc(l,j,kk)
+              ! only works for rotated lat-lon grids
+              aspect(j)=(x_len(k,ig)*region_dy_ens(grd_loc%nlat/2,grd_loc%nlon/2)/region_dx_ens(i,j))**2
+           enddo
+           call get_new_alpha_beta(aspect,grd_loc%nlon,fmatc,fmat0c)
+           do kk=1,2
+              do j=1,grd_loc%nlon
+                 do l=1,2
+                    fmatx(i,l,j,kk,k,ig)=fmatc(l,j,kk)
+                 enddo
+                 fmat0x(i,j,kk,k,ig)=fmat0c(j,kk)
               enddo
-              fmat0x(i,j,kk,k)=fmat0c(j,kk)
            enddo
         enddo
      enddo
-  enddo
+  enddo !ig loop
   return
 
 end subroutine init_rf_x
@@ -438,31 +447,34 @@ subroutine init_rf_y(y_len,kl)
 !$$$
 
   use hybrid_ensemble_parameters, only: grd_loc
+  use hybrid_ensemble_parameters, only: naensloc
   use constants, only: half
 
   implicit none
 
-  real(r_kind),intent(in   ) :: y_len(kl)
+  real(r_kind),intent(in   ) :: y_len(kl,naensloc)
   integer(i_kind),intent(in   ) :: kl
 
   real(r_kind) aspect(grd_loc%nlat)
-  integer(i_kind) i,k
+  integer(i_kind) i,k,ig
 
 !    use new factorization:
   if(allocated(fmaty)) deallocate(fmaty)
   if(allocated(fmat0y)) deallocate(fmat0y)
-  allocate(fmaty(2,grd_loc%nlat,2,kl),fmat0y(grd_loc%nlat,2,kl))
-  do k=1,kl
-     do i=1,grd_loc%nlat
-        aspect(i)=y_len(k)**2
+  allocate(fmaty(2,grd_loc%nlat,2,kl,naensloc),fmat0y(grd_loc%nlat,2,kl,naensloc))
+  do ig=1,naensloc
+     do k=1,kl
+        do i=1,grd_loc%nlat
+           aspect(i)=y_len(k,ig)**2
+        enddo
+        call get_new_alpha_beta(aspect,grd_loc%nlat,fmaty(1,1,1,k,ig),fmat0y(1,1,k,ig))
      enddo
-     call get_new_alpha_beta(aspect,grd_loc%nlat,fmaty(1,1,1,k),fmat0y(1,1,k))
-  enddo
+  enddo !ig loop
   return
 
 end subroutine init_rf_y
 
-subroutine new_factorization_rf_z(f,iadvance,iback)
+subroutine new_factorization_rf_z(f,iadvance,iback,ig)
 !$$$  subprogram documentation block
 !                .      .    .
 ! subprogram:    new_factorization_rf_z
@@ -480,6 +492,7 @@ subroutine new_factorization_rf_z(f,iadvance,iback)
 !     f        - input field to be filtered
 !     iadvance - =1  for forward operator, =2 for adjoint operator
 !     iback    - =2  for forward operator, =1 for adjoint operator
+!     ig       - number for smoothing scales
 !
 !   output argument list:
 !     f        - filtered output
@@ -493,7 +506,7 @@ subroutine new_factorization_rf_z(f,iadvance,iback)
   use hybrid_ensemble_parameters, only: grd_ens
   implicit none
 
-  integer(i_kind),intent(in   ) :: iadvance,iback
+  integer(i_kind),intent(in   ) :: iadvance,iback,ig
   real(r_kind)   ,intent(inout) :: f(grd_ens%latlon11,grd_ens%nsig)
 
   integer(i_kind) i,k,l,nxy,nz
@@ -502,34 +515,34 @@ subroutine new_factorization_rf_z(f,iadvance,iback)
   if(iadvance == 1) then
      do k=1,nz
         do i=1,nxy
-           f(i,k)=znorm_new(i,k)*f(i,k)
+           f(i,k)=znorm_new(i,k,ig)*f(i,k)
         enddo
      enddo
   end if
   do k=1,nz
      do l=1,min(2,k-1)
         do i=1,nxy
-           f(i,k)=f(i,k)-fmatz(i,l,k,iadvance)*f(i,k-l)
+           f(i,k)=f(i,k)-fmatz(i,l,k,iadvance,ig)*f(i,k-l)
         enddo
      enddo
      do i=1,nxy
-        f(i,k)=fmat0z(i,k,iadvance)*f(i,k)
+        f(i,k)=fmat0z(i,k,iadvance,ig)*f(i,k)
      enddo
   enddo
   do k=nz,1,-1
      do l=1,min(2,nz-k)
         do i=1,nxy
-           f(i,k)=f(i,k)-fmatz(i,l,k+l,iback)*f(i,k+l)
+           f(i,k)=f(i,k)-fmatz(i,l,k+l,iback,ig)*f(i,k+l)
         enddo
      enddo
      do i=1,nxy
-        f(i,k)=fmat0z(i,k,iback)*f(i,k)
+        f(i,k)=fmat0z(i,k,iback,ig)*f(i,k)
      enddo
   enddo
   if(iadvance == 2) then
      do k=1,nz
         do i=1,nxy
-           f(i,k)=znorm_new(i,k)*f(i,k)
+           f(i,k)=znorm_new(i,k,ig)*f(i,k)
         enddo
      enddo
   end if
@@ -537,7 +550,7 @@ subroutine new_factorization_rf_z(f,iadvance,iback)
 
 end subroutine new_factorization_rf_z
 
-subroutine new_factorization_rf_x(f,iadvance,iback,nlevs)
+subroutine new_factorization_rf_x(f,iadvance,iback,nlevs,ig)
 !$$$  subprogram documentation block
 !                .      .    .
 ! subprogram:    new_factorization_rf_x
@@ -556,6 +569,8 @@ subroutine new_factorization_rf_x(f,iadvance,iback,nlevs)
 !     f        - input field to be filtered
 !     iadvance - =1  for forward operator, =2 for adjoint operator
 !     iback    - =2  for forward operator, =1 for adjoint operator
+!     nlevs    - number of vertical levels for smoothing
+!     ig       - number for smoothing scales
 !
 !   output argument list:
 !     f        - filtered output
@@ -568,7 +583,7 @@ subroutine new_factorization_rf_x(f,iadvance,iback,nlevs)
   use hybrid_ensemble_parameters, only: grd_loc,vvlocal
   implicit none
 
-  integer(i_kind),intent(in   ) :: iadvance,iback,nlevs
+  integer(i_kind),intent(in   ) :: iadvance,iback,nlevs,ig
   real(r_kind)   ,intent(inout) :: f(grd_loc%nlat,grd_loc%nlon,max(nlevs,1))
 
   integer(i_kind) i,j,k,l,ny,nx,nz
@@ -582,7 +597,7 @@ subroutine new_factorization_rf_x(f,iadvance,iback,nlevs)
         if(iadvance == 1) then
            do j=1,nx
               do i=1,ny
-                 f(i,j,k)=xnorm_new(i,j,k)*f(i,j,k)
+                 f(i,j,k)=xnorm_new(i,j,k,ig)*f(i,j,k)
               enddo
            enddo
         end if
@@ -590,29 +605,29 @@ subroutine new_factorization_rf_x(f,iadvance,iback,nlevs)
         do j=1,nx
            do l=1,min(2,j-1)
               do i=1,ny
-                 f(i,j,k)=f(i,j,k)-fmatx(i,l,j,iadvance,k)*f(i,j-l,k)
+                 f(i,j,k)=f(i,j,k)-fmatx(i,l,j,iadvance,k,ig)*f(i,j-l,k)
               enddo
            enddo
            do i=1,ny
-              f(i,j,k)=fmat0x(i,j,iadvance,k)*f(i,j,k)
+              f(i,j,k)=fmat0x(i,j,iadvance,k,ig)*f(i,j,k)
            enddo
         enddo
 
         do j=nx,1,-1
            do l=1,min(2,nx-j)
               do i=1,ny
-                 f(i,j,k)=f(i,j,k)-fmatx(i,l,j+l,iback,k)*f(i,j+l,k)
+                 f(i,j,k)=f(i,j,k)-fmatx(i,l,j+l,iback,k,ig)*f(i,j+l,k)
               enddo
            enddo
            do i=1,ny
-              f(i,j,k)=fmat0x(i,j,iback,k)*f(i,j,k)
+              f(i,j,k)=fmat0x(i,j,iback,k,ig)*f(i,j,k)
            enddo
         enddo
 
         if(iadvance == 2) then
            do j=1,nx
               do i=1,ny
-                 f(i,j,k)=xnorm_new(i,j,k)*f(i,j,k)
+                 f(i,j,k)=xnorm_new(i,j,k,ig)*f(i,j,k)
               enddo
            enddo
         end if
@@ -625,7 +640,7 @@ subroutine new_factorization_rf_x(f,iadvance,iback,nlevs)
         if(iadvance == 1) then
            do j=1,nx
               do i=1,ny
-                 f(i,j,k)=xnorm_new(i,j,1)*f(i,j,k)
+                 f(i,j,k)=xnorm_new(i,j,1,ig)*f(i,j,k)
               enddo
            enddo
         end if
@@ -633,29 +648,29 @@ subroutine new_factorization_rf_x(f,iadvance,iback,nlevs)
         do j=1,nx
            do l=1,min(2,j-1)
               do i=1,ny
-                 f(i,j,k)=f(i,j,k)-fmatx(i,l,j,iadvance,1)*f(i,j-l,k)
+                 f(i,j,k)=f(i,j,k)-fmatx(i,l,j,iadvance,1,ig)*f(i,j-l,k)
               enddo
            enddo
            do i=1,ny
-              f(i,j,k)=fmat0x(i,j,iadvance,1)*f(i,j,k)
+              f(i,j,k)=fmat0x(i,j,iadvance,1,ig)*f(i,j,k)
            enddo
         enddo
 
         do j=nx,1,-1
            do l=1,min(2,nx-j)
               do i=1,ny
-                 f(i,j,k)=f(i,j,k)-fmatx(i,l,j+l,iback,1)*f(i,j+l,k)
+                 f(i,j,k)=f(i,j,k)-fmatx(i,l,j+l,iback,1,ig)*f(i,j+l,k)
               enddo
            enddo
            do i=1,ny
-              f(i,j,k)=fmat0x(i,j,iback,1)*f(i,j,k)
+              f(i,j,k)=fmat0x(i,j,iback,1,ig)*f(i,j,k)
            enddo
         enddo
 
         if(iadvance == 2) then
            do j=1,nx
               do i=1,ny
-                 f(i,j,k)=xnorm_new(i,j,1)*f(i,j,k)
+                 f(i,j,k)=xnorm_new(i,j,1,ig)*f(i,j,k)
               enddo
            enddo
         end if
@@ -665,7 +680,7 @@ subroutine new_factorization_rf_x(f,iadvance,iback,nlevs)
   return
 end subroutine new_factorization_rf_x
 
-subroutine new_factorization_rf_y(f,iadvance,iback,nlevs)
+subroutine new_factorization_rf_y(f,iadvance,iback,nlevs,ig)
 !$$$  subprogram documentation block
 !                .      .    .
 ! subprogram:    new_factorization_rf_y
@@ -684,6 +699,8 @@ subroutine new_factorization_rf_y(f,iadvance,iback,nlevs)
 !     f        - input field to be filtered
 !     iadvance - =1  for forward operator, =2 for adjoint operator
 !     iback    - =2  for forward operator, =1 for adjoint operator
+!     nlevs    - number of vertical levels for smoothing
+!     ig       - number for smoothing scales
 !
 !   output argument list:
 !     f        - filtered output
@@ -697,7 +714,7 @@ subroutine new_factorization_rf_y(f,iadvance,iback,nlevs)
                       !                  use mpimod, only: mype
   implicit none
 
-  integer(i_kind),intent(in   ) :: iadvance,iback,nlevs
+  integer(i_kind),intent(in   ) :: iadvance,iback,nlevs,ig
   real(r_kind)   ,intent(inout) :: f(grd_loc%nlat,grd_loc%nlon,max(nlevs,1))
 
   integer(i_kind) i,j,k,l,nx,ny,nz
@@ -710,27 +727,27 @@ subroutine new_factorization_rf_y(f,iadvance,iback,nlevs)
 
            if(iadvance == 1) then
               do i=1,ny
-                 f(i,j,k)=ynorm_new(i,k)*f(i,j,k)
+                 f(i,j,k)=ynorm_new(i,k,ig)*f(i,j,k)
               enddo
            end if
 
            do i=1,ny
               do l=1,min(2,i-1)
-                 f(i,j,k)=f(i,j,k)-fmaty(l,i,iadvance,k)*f(i-l,j,k)
+                 f(i,j,k)=f(i,j,k)-fmaty(l,i,iadvance,k,ig)*f(i-l,j,k)
               enddo
-              f(i,j,k)=fmat0y(i,iadvance,k)*f(i,j,k)
+              f(i,j,k)=fmat0y(i,iadvance,k,ig)*f(i,j,k)
            enddo
 
            do i=ny,1,-1
               do l=1,min(2,ny-i)
-                 f(i,j,k)=f(i,j,k)-fmaty(l,i+l,iback,k)*f(i+l,j,k)
+                 f(i,j,k)=f(i,j,k)-fmaty(l,i+l,iback,k,ig)*f(i+l,j,k)
               enddo
-              f(i,j,k)=fmat0y(i,iback,k)*f(i,j,k)
+              f(i,j,k)=fmat0y(i,iback,k,ig)*f(i,j,k)
            enddo
 
            if(iadvance == 2) then
               do i=1,ny
-                 f(i,j,k)=ynorm_new(i,k)*f(i,j,k)
+                 f(i,j,k)=ynorm_new(i,k,ig)*f(i,j,k)
               enddo
            end if
 
@@ -742,27 +759,27 @@ subroutine new_factorization_rf_y(f,iadvance,iback,nlevs)
 
            if(iadvance == 1) then
               do i=1,ny
-                 f(i,j,k)=ynorm_new(i,1)*f(i,j,k)
+                 f(i,j,k)=ynorm_new(i,1,ig)*f(i,j,k)
               enddo
            end if
 
            do i=1,ny
               do l=1,min(2,i-1)
-                 f(i,j,k)=f(i,j,k)-fmaty(l,i,iadvance,1)*f(i-l,j,k)
+                 f(i,j,k)=f(i,j,k)-fmaty(l,i,iadvance,1,ig)*f(i-l,j,k)
               enddo
-              f(i,j,k)=fmat0y(i,iadvance,1)*f(i,j,k)
+              f(i,j,k)=fmat0y(i,iadvance,1,ig)*f(i,j,k)
            enddo
 
            do i=ny,1,-1
               do l=1,min(2,ny-i)
-                 f(i,j,k)=f(i,j,k)-fmaty(l,i+l,iback,1)*f(i+l,j,k)
+                 f(i,j,k)=f(i,j,k)-fmaty(l,i+l,iback,1,ig)*f(i+l,j,k)
               enddo
-              f(i,j,k)=fmat0y(i,iback,1)*f(i,j,k)
+              f(i,j,k)=fmat0y(i,iback,1,ig)*f(i,j,k)
            enddo
 
            if(iadvance == 2) then
               do i=1,ny
-                 f(i,j,k)=ynorm_new(i,1)*f(i,j,k)
+                 f(i,j,k)=ynorm_new(i,1,ig)*f(i,j,k)
               enddo
            end if
 
@@ -798,34 +815,37 @@ subroutine normal_new_factorization_rf_z
 
   use kinds, only: r_kind,i_kind
   use hybrid_ensemble_parameters, only: grd_ens
+  use hybrid_ensemble_parameters, only: naensgrp,naensloc
   use constants, only: zero,one
   implicit none
 
-  integer(i_kind) k,iadvance,iback,nxy
+  integer(i_kind) k,iadvance,iback,nxy,ig
   real(r_kind) f(grd_ens%latlon11,grd_ens%nsig),diag(grd_ens%latlon11,grd_ens%nsig)
 
   if(allocated(znorm_new)) deallocate(znorm_new)
-  allocate(znorm_new(grd_ens%latlon11,grd_ens%nsig))
+  allocate(znorm_new(grd_ens%latlon11,grd_ens%nsig,naensloc))
 
   nxy=grd_ens%latlon11
 
   znorm_new=one
 
-  do k=1,grd_ens%nsig
-      f=zero
-      f(:,k)=one
-
-      iadvance=1 ; iback=2
-      call new_factorization_rf_z(f,iadvance,iback)
-      iadvance=2 ; iback=1
-      call new_factorization_rf_z(f,iadvance,iback)
+  do ig=1,naensgrp
+     do k=1,grd_ens%nsig
+        f=zero
+        f(:,k)=one
 
-      diag(:,k)=sqrt(one/f(:,k))
-  enddo
+        iadvance=1 ; iback=2
+        call new_factorization_rf_z(f,iadvance,iback,ig)
+        iadvance=2 ; iback=1
+        call new_factorization_rf_z(f,iadvance,iback,ig)
+        
+        diag(:,k)=sqrt(one/f(:,k))
+     enddo
 
-  do k=1,grd_ens%nsig
-      znorm_new(:,k)=diag(:,k)
-  enddo
+     do k=1,grd_ens%nsig
+        znorm_new(:,k,ig)=diag(:,k)
+     enddo
+  enddo !ig loop
 
 ! Check result:
   if(debug)then
@@ -834,9 +854,9 @@ subroutine normal_new_factorization_rf_z
        f(:,k)=one
 
        iadvance=1 ; iback=2
-       call new_factorization_rf_z(f,iadvance,iback)
+       call new_factorization_rf_z(f,iadvance,iback,1)
        iadvance=2 ; iback=1
-       call new_factorization_rf_z(f,iadvance,iback)
+       call new_factorization_rf_z(f,iadvance,iback,1)
 
        diag(:,k)=sqrt(one/f(:,k))
     enddo
@@ -876,11 +896,12 @@ subroutine normal_new_factorization_rf_x
 
   use kinds, only: r_kind,i_kind
   use hybrid_ensemble_parameters, only: grd_loc,vvlocal
+  use hybrid_ensemble_parameters, only: naensgrp,naensloc
   use constants, only: zero,one
 
   implicit none
 
-  integer(i_kind) i,j,k,iadvance,iback,kl
+  integer(i_kind) i,j,k,iadvance,iback,kl,ig
   real(r_kind) f(grd_loc%nlat,grd_loc%nlon,grd_loc%kend_alloc+1-grd_loc%kbegin_loc)
   real(r_kind),allocatable:: diag(:,:,:)
 !  real(r_kind) diag(grd_loc%nlat,grd_loc%nlon)
@@ -896,35 +917,37 @@ subroutine normal_new_factorization_rf_x
      kl=1
   endif
   if(allocated(xnorm_new)) deallocate(xnorm_new)
-  allocate(xnorm_new(grd_loc%nlat,grd_loc%nlon,kl))
+  allocate(xnorm_new(grd_loc%nlat,grd_loc%nlon,kl,naensloc))
   if(allocated(diag)) deallocate(diag)
   allocate(diag(grd_loc%nlat,grd_loc%nlon,kl))
   xnorm_new=one
 
-  do j=1,grd_loc%nlon
-     f=zero
-     do k=1,kl
-        do i=1,grd_loc%nlat
-           f(i,j,k)=one
+  do ig=1,naensgrp
+     do j=1,grd_loc%nlon
+        f=zero
+        do k=1,kl
+           do i=1,grd_loc%nlat
+              f(i,j,k)=one
+           enddo
         enddo
-     enddo
-     iadvance=1 ; iback=2
-     call new_factorization_rf_x(f,iadvance,iback,kl)
-     iadvance=2 ; iback=1
-     call new_factorization_rf_x(f,iadvance,iback,kl)
-     do k=1,kl
-        do i=1,grd_loc%nlat
-           diag(i,j,k)=sqrt(one/f(i,j,k))
+        iadvance=1 ; iback=2
+        call new_factorization_rf_x(f,iadvance,iback,kl,ig)
+        iadvance=2 ; iback=1
+        call new_factorization_rf_x(f,iadvance,iback,kl,ig)
+        do k=1,kl
+           do i=1,grd_loc%nlat
+              diag(i,j,k)=sqrt(one/f(i,j,k))
+           enddo
         enddo
      enddo
-  enddo
-  do k=1,kl
-     do j=1,grd_loc%nlon
-        do i=1,grd_loc%nlat
-           xnorm_new(i,j,k)=diag(i,j,k)
+     do k=1,kl
+        do j=1,grd_loc%nlon
+           do i=1,grd_loc%nlat
+              xnorm_new(i,j,k,ig)=diag(i,j,k)
+           enddo
         enddo
      enddo
-  enddo
+  enddo !ig loop
 !           check accuracy of xnorm
   if(debug) then
      do j=1,grd_loc%nlon
@@ -935,9 +958,9 @@ subroutine normal_new_factorization_rf_x
            enddo
         enddo
         iadvance=1 ; iback=2
-        call new_factorization_rf_x(f,iadvance,iback,kl)
+        call new_factorization_rf_x(f,iadvance,iback,kl,1)
         iadvance=2 ; iback=1
-        call new_factorization_rf_x(f,iadvance,iback,kl)
+        call new_factorization_rf_x(f,iadvance,iback,kl,1)
         do k=1,kl
            do i=1,grd_loc%nlat
               diag(i,j,k)=f(i,j,k)
@@ -977,10 +1000,11 @@ subroutine normal_new_factorization_rf_y
 
   use kinds, only: r_kind,i_kind
   use hybrid_ensemble_parameters, only: grd_loc,vvlocal
+  use hybrid_ensemble_parameters, only: naensgrp,naensloc
   use constants, only: zero,one
   implicit none
 
-  integer(i_kind) i,k,lend,lcount,iadvance,iback,kl,loop,ll,iend
+  integer(i_kind) i,k,lend,lcount,iadvance,iback,kl,loop,ll,iend,ig
 !  real(r_kind) f(grd_loc%nlat,grd_loc%nlon*(grd_loc%kend_alloc+1-grd_loc%kbegin_loc)),diag(grd_loc%nlat)
   real(r_kind) f(grd_loc%nlat,grd_loc%nlon,grd_loc%kend_alloc+1-grd_loc%kbegin_loc)
   real(r_kind),allocatable:: diag(:,:)
@@ -997,7 +1021,7 @@ subroutine normal_new_factorization_rf_y
   endif
 
   if(allocated(ynorm_new)) deallocate(ynorm_new)
-  allocate(ynorm_new(grd_loc%nlat,kl))
+  allocate(ynorm_new(grd_loc%nlat,kl,naensloc))
 
   if(allocated(diag)) deallocate(diag)
   allocate(diag(grd_loc%nlat,kl))
@@ -1012,31 +1036,33 @@ subroutine normal_new_factorization_rf_y
     iend=grd_loc%nlon 
   endif
               
-  do loop=1,lend
-     ll=(loop-1)*iend
-     f=zero
-     do k=1,kl
-        do i=1,iend
-        lcount=ll+i
-        f(lcount,i,k)=one
-        if(lcount == grd_loc%nlat) exit
+  do ig=1,naensgrp
+     do loop=1,lend
+        ll=(loop-1)*iend
+        f=zero
+        do k=1,kl
+           do i=1,iend
+              lcount=ll+i
+              f(lcount,i,k)=one
+              if(lcount == grd_loc%nlat) exit
+           enddo
         enddo
-     enddo
 
-     iadvance=1 ; iback=2
-     call new_factorization_rf_y(f,iadvance,iback,kl)
-     iadvance=2 ; iback=1
-     call new_factorization_rf_y(f,iadvance,iback,kl)
+        iadvance=1 ; iback=2
+        call new_factorization_rf_y(f,iadvance,iback,kl,ig)
+        iadvance=2 ; iback=1
+        call new_factorization_rf_y(f,iadvance,iback,kl,ig)
 
-     do k=1,kl
-        do i=1,iend
-           lcount=ll+i
-           diag(lcount,k)=sqrt(one/f(lcount,i,k))
-           ynorm_new(lcount,k)=diag(lcount,k)
-           if(lcount == grd_loc%nlat) exit
+        do k=1,kl
+           do i=1,iend
+              lcount=ll+i
+              diag(lcount,k)=sqrt(one/f(lcount,i,k))
+              ynorm_new(lcount,k,ig)=diag(lcount,k)
+              if(lcount == grd_loc%nlat) exit
+           enddo
         enddo
      enddo
-  enddo
+  enddo !ig loop
 !               check that ynorm is corect
   if(debug) then
      do loop=1,lend
@@ -1051,9 +1077,9 @@ subroutine normal_new_factorization_rf_y
         enddo
 
         iadvance=1 ; iback=2
-        call new_factorization_rf_y(f,iadvance,iback,kl)
+        call new_factorization_rf_y(f,iadvance,iback,kl,1)
         iadvance=2 ; iback=1
-        call new_factorization_rf_y(f,iadvance,iback,kl) 
+        call new_factorization_rf_y(f,iadvance,iback,kl,1) 
 
         do k=1,kl
            do i=1,iend
@@ -1094,37 +1120,40 @@ subroutine create_ensemble
 !$$$
     use hybrid_ensemble_parameters, only: n_ens,grd_ens,ntlevs_ens
     use hybrid_ensemble_parameters, only: nelen,en_perts,ps_bar
+    use hybrid_ensemble_parameters, only: ntotensgrp
 
     implicit none
 
     type(gsi_grid)  :: grid_ens
 
-    integer(i_kind) n,istatus,m
+    integer(i_kind) n,istatus,m,ig
     character(len=*),parameter::myname_=trim(myname)//'*create_ensemble'
 
     nelen=grd_ens%latlon11*(max(0,nc3d)*grd_ens%nsig+max(0,nc2d))
 !   create ensemble perturbations bundles (using newly added r_single capability
 
-    allocate(en_perts(n_ens,ntlevs_ens))
+    allocate(en_perts(n_ens,ntotensgrp,ntlevs_ens))
     call gsi_gridcreate(grid_ens,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig)
  
     do m=1,ntlevs_ens
-       do n=1,n_ens
-          call gsi_bundlecreate(en_perts(n,m),grid_ens,'ensemble perts',istatus, &
-                                names2d=cvars2d,names3d=cvars3d,bundle_kind=r_single)
-          if(istatus/=0) then
-             write(6,*)trim(myname_),': trouble creating en_perts bundle'
-             call stop2(999)
-          endif
+       do ig=1,ntotensgrp
+          do n=1,n_ens
+             call gsi_bundlecreate(en_perts(n,ig,m),grid_ens,'ensemble perts',istatus, &
+                                   names2d=cvars2d,names3d=cvars3d,bundle_kind=r_single)
+             if(istatus/=0) then
+                write(6,*)trim(myname_),': trouble creating en_perts bundle'
+                call stop2(999)
+             endif
+          enddo
        enddo
     enddo
 
 
     allocate(ps_bar(grd_ens%lat2,grd_ens%lon2,ntlevs_ens) )
     if(debug) then
-       write(6,*)' in create_ensemble, grd_ens%latlon11,grd_ens%latlon1n,n_ens,ntlevs_ens=', &
-                                 grd_ens%latlon11,grd_ens%latlon1n,n_ens,ntlevs_ens
-       write(6,*)' in create_ensemble, total bytes allocated=',4*nelen*n_ens*ntlevs_ens
+       write(6,*)' in create_ensemble, grd_ens%latlon11,grd_ens%latlon1n,n_ens,ntotensgrp,ntlevs_ens=', &
+                                 grd_ens%latlon11,grd_ens%latlon1n,n_ens,ntotensgrp,ntlevs_ens
+       write(6,*)' in create_ensemble, total bytes allocated=',4*nelen*n_ens*ntotensgrp*ntlevs_ens
     end if
     return
 
@@ -1245,7 +1274,7 @@ subroutine load_ensemble
           do n=1,n_ens
              call generate_one_ensemble_perturbation(bundle_anl,bundle_ens,seed)
              do ii=1,nelen
-                en_perts(n,m)%valuesr4(ii)=bundle_ens%values(ii)
+                en_perts(n,1,m)%valuesr4(ii)=bundle_ens%values(ii)
                 en_bar(m)%values(ii)=en_bar(m)%values(ii)+bundle_ens%values(ii)
              enddo
           enddo
@@ -1276,9 +1305,9 @@ subroutine load_ensemble
        do m=1,ntlevs_ens
           do n=1,n_ens
              do ii=1,nelen
-                en_perts(n,m)%valuesr4(ii)=(en_perts(n,m)%valuesr4(ii)-en_bar(m)%values(ii)*bar_norm)*sig_norm
+                en_perts(n,1,m)%valuesr4(ii)=(en_perts(n,1,m)%valuesr4(ii)-en_bar(m)%values(ii)*bar_norm)*sig_norm
              enddo
-             call gsi_enscoupler_put_gsi_ens(grd_ens,n,m,en_perts(n,m),istatus)
+             call gsi_enscoupler_put_gsi_ens(grd_ens,n,m,en_perts(n,1,m),istatus)
              if(istatus/=0) then
                  write(6,*)trim(myname_),': trouble writing perts'
                  call stop2(999)
@@ -1634,15 +1663,16 @@ subroutine rescale_ensemble_rh_perturbations
     use gridmod, only: regional
     use hybrid_ensemble_parameters, only: n_ens,grd_ens,grd_anl,grd_a1,grd_e1,p_e2a,ntlevs_ens
     use hybrid_ensemble_parameters, only: en_perts
+    use hybrid_ensemble_parameters, only: ntotensgrp
     use general_sub2grid_mod, only: general_suba2sube
     use berror, only: qvar3d
     implicit none
 
-    integer(i_kind) i,j,k,n,istatus,m
+    integer(i_kind) i,j,k,n,istatus,m,ig
     real(r_kind) qvar3d_ens(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig,1)
     real(r_single),pointer,dimension(:,:,:):: w3=>NULL()
 
-    call gsi_bundlegetpointer(en_perts(1,1),'q',w3,istatus)
+    call gsi_bundlegetpointer(en_perts(1,1,1),'q',w3,istatus)
     if(istatus/=0) then
        write(6,*)' rh variable not available, skip subroutine rescale_ensemble_rh_perturbations'
        return
@@ -1655,17 +1685,19 @@ subroutine rescale_ensemble_rh_perturbations
             reshape(qvar3d,(/size(qvar3d,1),size(qvar3d,2),size(qvar3d,3),1/)),qvar3d_ens,regional)
     end if
     do m=1,ntlevs_ens
+       do ig=1,ntotensgrp
 !$omp parallel do schedule(dynamic,1) private(n,i,j,k,w3,istatus)
-       do n=1,n_ens
-          call gsi_bundlegetpointer(en_perts(n,m),'q',w3,istatus)
-          if(istatus/=0) then
-             write(6,*)' error retrieving pointer to rh variable for ensemble number ',n
-             call stop2(999)
-          end if
-          do k=1,grd_ens%nsig
-             do j=1,grd_ens%lon2
-                do i=1,grd_ens%lat2
-                   w3(i,j,k)=qvar3d_ens(i,j,k,1)*w3(i,j,k)
+          do n=1,n_ens
+             call gsi_bundlegetpointer(en_perts(n,ig,m),'q',w3,istatus)
+             if(istatus/=0) then
+                write(6,*)' error retrieving pointer to rh variable for ensemble number ',n
+                call stop2(999)
+             end if
+             do k=1,grd_ens%nsig
+                do j=1,grd_ens%lon2
+                   do i=1,grd_ens%lat2
+                      w3(i,j,k)=qvar3d_ens(i,j,k,1)*w3(i,j,k)
+                   enddo
                 enddo
              enddo
           enddo
@@ -1699,18 +1731,21 @@ subroutine destroy_ensemble
 !$$$
     use hybrid_ensemble_parameters, only: l_hyb_ens,n_ens,ntlevs_ens
     use hybrid_ensemble_parameters, only: en_perts,ps_bar
+    use hybrid_ensemble_parameters, only: ntotensgrp
     implicit none
 
-    integer(i_kind) istatus,n,m
+    integer(i_kind) istatus,n,m,ig
 
     if(l_hyb_ens) then
        do m=1,ntlevs_ens
-          do n=1,n_ens
-             call gsi_bundleunset(en_perts(n,m),istatus)
-             if(istatus/=0) then
-                write(6,*)'in destroy_ensemble: trouble destroying en_perts bundle'
-                call stop2(999)
-             endif
+          do ig=1,ntotensgrp
+             do n=1,n_ens
+                call gsi_bundleunset(en_perts(n,ig,m),istatus)
+                if(istatus/=0) then
+                   write(6,*)'in destroy_ensemble: trouble destroying en_perts bundle'
+                   call stop2(999)
+                endif
+             enddo
           enddo
        enddo
        deallocate(ps_bar)
@@ -1759,24 +1794,26 @@ subroutine ensemble_forward_model(cvec,a_en,ibin)
 !$$$
     use hybrid_ensemble_parameters, only: n_ens,pwgtflg,pwgt
     use hybrid_ensemble_parameters, only: en_perts
+    use hybrid_ensemble_parameters, only: ntotensgrp,naensgrp
+    use hybrid_ensemble_parameters, only: ensgrp2aensgrp
     use constants, only: zero
 
     implicit none
     type(gsi_bundle),intent(inout) :: cvec
-    type(gsi_bundle),intent(in)    :: a_en(n_ens)
+    type(gsi_bundle),intent(in)    :: a_en(naensgrp,n_ens)
     integer,intent(in)             :: ibin
 
     character(len=*),parameter :: myname_=trim(myname)//'*ensemble_forward_model'
     logical :: nogood
     integer(i_kind) :: i,j,k,n,im,jm,km,ic2,ic3,ipic,ipx,km_tmp
-    integer(i_kind) :: ipc3d(nc3d),ipc2d(nc2d),istatus
+    integer(i_kind) :: ipc3d(nc3d),ipc2d(nc2d),istatus,ig,iaens
 
     im=cvec%grid%im
     jm=cvec%grid%jm
     km=cvec%grid%km
 
 !   Check resolution consistency between static and ensemble components
-    nogood=im/=a_en(1)%grid%im.or.jm/=a_en(1)%grid%jm.or.km/=a_en(1)%grid%km
+    nogood=im/=a_en(1,1)%grid%im.or.jm/=a_en(1,1)%grid%jm.or.km/=a_en(1,1)%grid%km
     if (nogood) then
        write(6,*) myname_,': static&ensemble vectors have inconsistent dims'
        call stop2(999)
@@ -1798,7 +1835,7 @@ subroutine ensemble_forward_model(cvec,a_en,ibin)
  
     ipx=1
 
-!$omp parallel do schedule(dynamic,1) private(j,n,ic3,k,i,ipic)
+!$omp parallel do schedule(dynamic,1) private(j,n,ic3,k,i,ipic,ig,iaens)
     do k=1,km
        do ic3=1,nc3d
           ipic=ipc3d(ic3)
@@ -1807,18 +1844,23 @@ subroutine ensemble_forward_model(cvec,a_en,ibin)
                 cvec%r3(ipic)%q(i,j,k)=zero
              enddo
           enddo
-          do n=1,n_ens
-             do j=1,jm
-                do i=1,im
-                   cvec%r3(ipic)%q(i,j,k)=cvec%r3(ipic)%q(i,j,k) &
-                         +a_en(n)%r3(ipx)%q(i,j,k)*en_perts(n,ibin)%r3(ipic)%qr4(i,j,k)
+          do ig=1,ntotensgrp
+             iaens=ensgrp2aensgrp(ig,ic3,ibin)
+             if(iaens>0) then
+                do n=1,n_ens
+                   do j=1,jm
+                      do i=1,im
+                         cvec%r3(ipic)%q(i,j,k)=cvec%r3(ipic)%q(i,j,k) &
+                              +a_en(iaens,n)%r3(ipx)%q(i,j,k)*en_perts(n,ig,ibin)%r3(ipic)%qr4(i,j,k)
+                      enddo
+                   enddo
                 enddo
-             enddo
+             endif ! iaens>0
           enddo
        enddo
     enddo
 
-!$omp parallel do schedule(dynamic,1) private(j,n,k,i,ic2,ipic)
+!$omp parallel do schedule(dynamic,1) private(j,n,k,i,ic2,ipic,ig,iaens)
     do ic2=1,nc2d
        ipic=ipc2d(ic2)
        do j=1,jm
@@ -1837,27 +1879,37 @@ subroutine ensemble_forward_model(cvec,a_en,ibin)
                 km_tmp = 1
              endif
 
-             do n=1,n_ens
-                do j=1,jm
-                   do k=1,km_tmp
-                      do i=1,im
-                         cvec%r2(ipic)%q(i,j)=cvec%r2(ipic)%q(i,j) &
-                            +a_en(n)%r3(ipx)%q(i,j,k)*en_perts(n,ibin)%r2(ipic)%qr4(i,j)*pwgt(i,j,k)
+             do ig=1,ntotensgrp
+                iaens=ensgrp2aensgrp(ig,ic2+nc3d,ibin)
+                if(iaens>0) then
+                   do n=1,n_ens
+                      do j=1,jm
+                         do k=1,km_tmp
+                            do i=1,im
+                               cvec%r2(ipic)%q(i,j)=cvec%r2(ipic)%q(i,j) &
+                                    +a_en(iaens,n)%r3(ipx)%q(i,j,k)*en_perts(n,ig,ibin)%r2(ipic)%qr4(i,j)*pwgt(i,j,k)
+                            enddo
+                         enddo
                       enddo
-                   enddo
-                enddo
-             enddo ! enddo n_ens
+                   enddo ! enddo n_ens
+                endif ! iaens>0
+             enddo
 
           case('SST')
  
-             do n=1,n_ens
-                do j=1,jm
-                   do i=1,im
-                      cvec%r2(ipic)%q(i,j)=cvec%r2(ipic)%q(i,j) &
-                         +a_en(n)%r3(ipx)%q(i,j,1)*en_perts(n,ibin)%r2(ipic)%qr4(i,j)
-                   enddo
-                enddo
-             enddo ! enddo n_ens
+             do ig=1,ntotensgrp
+                iaens=ensgrp2aensgrp(ig,ic2+nc3d,ibin)
+                if(iaens>0) then
+                   do n=1,n_ens
+                      do j=1,jm
+                         do i=1,im
+                            cvec%r2(ipic)%q(i,j)=cvec%r2(ipic)%q(i,j) &
+                                 +a_en(iaens,n)%r3(ipx)%q(i,j,1)*en_perts(n,ig,ibin)%r2(ipic)%qr4(i,j)
+                         enddo
+                      enddo
+                   enddo ! enddo n_ens
+                endif ! iaens>0
+             enddo
  
        end select
 
@@ -1909,20 +1961,22 @@ subroutine ensemble_forward_model_dual_res(cvec,a_en,ibin)
     use hybrid_ensemble_parameters, only: n_ens,pwgtflg,pwgt
     use hybrid_ensemble_parameters, only: grd_ens,grd_anl,p_e2a
     use hybrid_ensemble_parameters, only: en_perts
+    use hybrid_ensemble_parameters, only: ntotensgrp,naensgrp
+    use hybrid_ensemble_parameters, only: ensgrp2aensgrp
     use general_sub2grid_mod, only: general_sube2suba
     use gridmod,only: regional
     use constants, only: zero
     implicit none
 
     type(gsi_bundle),intent(inout) :: cvec
-    type(gsi_bundle),intent(in)    :: a_en(n_ens)
+    type(gsi_bundle),intent(in)    :: a_en(naensgrp,n_ens)
     integer,intent(in)             :: ibin
 
     character(len=*),parameter::myname_=trim(myname)//'*ensemble_forward_model_dual_res'
     type(gsi_grid)   :: grid_ens,grid_anl
     type(gsi_bundle) :: work_ens,work_anl
     integer(i_kind) :: i,j,k,n,im,jm,km,ic2,ic3,ipic,ipx,km_tmp
-    integer(i_kind) :: ipc2d(nc2d),ipc3d(nc3d),istatus
+    integer(i_kind) :: ipc2d(nc2d),ipc3d(nc3d),istatus,ig,iaens
 
 !   Request ensemble-corresponding fields from control vector
 !    NOTE:  because ensemble perturbation bundle structure is same as control vector, use same ipc3d and
@@ -1958,7 +2012,7 @@ subroutine ensemble_forward_model_dual_res(cvec,a_en,ibin)
     im=work_ens%grid%im
     jm=work_ens%grid%jm
     km=work_ens%grid%km
-!$omp parallel do schedule(dynamic,1) private(j,n,ic3,k,i,ipic)
+!$omp parallel do schedule(dynamic,1) private(j,n,ic3,k,i,ipic,ig,iaens)
     do k=1,km
        do ic3=1,nc3d
           ipic=ipc3d(ic3)
@@ -1967,17 +2021,22 @@ subroutine ensemble_forward_model_dual_res(cvec,a_en,ibin)
                 work_ens%r3(ipic)%q(i,j,k)=zero
              enddo
           enddo
-          do n=1,n_ens
-             do j=1,jm
-                do i=1,im
-                   work_ens%r3(ipic)%q(i,j,k)=work_ens%r3(ipic)%q(i,j,k) &
-                      +a_en(n)%r3(ipx)%q(i,j,k)*en_perts(n,ibin)%r3(ipic)%qr4(i,j,k)
+          do ig=1,ntotensgrp
+             iaens=ensgrp2aensgrp(ig,ic3,ibin)
+             if(iaens>0) then
+                do n=1,n_ens
+                   do j=1,jm
+                      do i=1,im
+                         work_ens%r3(ipic)%q(i,j,k)=work_ens%r3(ipic)%q(i,j,k) &
+                              +a_en(iaens,n)%r3(ipx)%q(i,j,k)*en_perts(n,ig,ibin)%r3(ipic)%qr4(i,j,k)
+                      enddo
+                   enddo
                 enddo
-             enddo
+             endif ! iaens>0
           enddo
        enddo
     enddo
-!$omp parallel do schedule(dynamic,1) private(j,n,k,i,ic2,ipic)
+!$omp parallel do schedule(dynamic,1) private(j,n,k,i,ic2,ipic,ig,iaens)
     do ic2=1,nc2d
        ipic=ipc2d(ic2)
        do j=1,jm
@@ -1996,27 +2055,37 @@ subroutine ensemble_forward_model_dual_res(cvec,a_en,ibin)
                 km_tmp = 1
              endif
 
-             do n=1,n_ens
-                do k=1,km_tmp
-                   do j=1,jm
-                      do i=1,im
-                         work_ens%r2(ipic)%q(i,j)=work_ens%r2(ipic)%q(i,j) &
-                            +a_en(n)%r3(ipx)%q(i,j,k)*en_perts(n,ibin)%r2(ipic)%qr4(i,j)*pwgt(i,j,k)
+             do ig=1,ntotensgrp
+                iaens=ensgrp2aensgrp(ig,ic2+nc3d,ibin)
+                if(iaens>0) then
+                   do n=1,n_ens
+                      do k=1,km_tmp
+                         do j=1,jm
+                            do i=1,im
+                               work_ens%r2(ipic)%q(i,j)=work_ens%r2(ipic)%q(i,j) &
+                                    +a_en(iaens,n)%r3(ipx)%q(i,j,k)*en_perts(n,ig,ibin)%r2(ipic)%qr4(i,j)*pwgt(i,j,k)
+                            enddo
+                         enddo
                       enddo
-                   enddo
-                enddo
-             enddo ! enddo n_ens
+                   enddo ! enddo n_ens
+                endif ! iaens>0
+             enddo
 
           case('SST')
 
-             do n=1,n_ens
-                do j=1,jm
-                   do i=1,im
-                      work_ens%r2(ipic)%q(i,j)=work_ens%r2(ipic)%q(i,j) &
-                         +a_en(n)%r3(ipx)%q(i,j,1)*en_perts(n,ibin)%r2(ipic)%qr4(i,j)
-                   enddo
-                enddo
-             enddo ! enddo n_ens
+             do ig=1,ntotensgrp
+                iaens=ensgrp2aensgrp(ig,ic2+nc3d,ibin)
+                if(iaens>0) then
+                   do n=1,n_ens
+                      do j=1,jm
+                         do i=1,im
+                            work_ens%r2(ipic)%q(i,j)=work_ens%r2(ipic)%q(i,j) &
+                                 +a_en(iaens,n)%r3(ipx)%q(i,j,1)*en_perts(n,ig,ibin)%r2(ipic)%qr4(i,j)
+                         enddo
+                      enddo
+                   enddo ! enddo n_ens
+                endif ! iaens>0
+             enddo
 
        end select
 
@@ -2082,22 +2151,24 @@ subroutine ensemble_forward_model_ad(cvec,a_en,ibin)
 
     use hybrid_ensemble_parameters, only: n_ens,pwgtflg,pwgt
     use hybrid_ensemble_parameters, only: en_perts
+    use hybrid_ensemble_parameters, only: ntotensgrp,naensgrp
+    use hybrid_ensemble_parameters, only: ensgrp2aensgrp
     implicit none
 
     type(gsi_bundle),intent(inout) :: cvec
-    type(gsi_bundle),intent(inout) :: a_en(n_ens)
+    type(gsi_bundle),intent(inout) :: a_en(naensgrp,n_ens)
     integer,intent(in)             :: ibin
 
     character(len=*),parameter :: myname_=trim(myname)//'*ensemble_forward_model_ad'
     logical :: nogood
     integer(i_kind) :: i,j,k,n,im,jm,km,ic2,ic3,ipx,ipic,km_tmp
-    integer(i_kind) :: ipc3d(nc3d),ipc2d(nc2d),istatus
+    integer(i_kind) :: ipc3d(nc3d),ipc2d(nc2d),istatus,ig,iaens
 
     im=cvec%grid%im
     jm=cvec%grid%jm
     km=cvec%grid%km
 !   Check resolution consistency between static and ensemble components
-    nogood=im/=a_en(1)%grid%im.or.jm/=a_en(1)%grid%jm.or.km/=a_en(1)%grid%km
+    nogood=im/=a_en(1,1)%grid%im.or.jm/=a_en(1,1)%grid%jm.or.km/=a_en(1,1)%grid%km
     if (nogood) then
        write(6,*) myname_,': static/ensemble vectors have inconsistent dims'
        call stop2(999)
@@ -2118,51 +2189,58 @@ subroutine ensemble_forward_model_ad(cvec,a_en,ibin)
     endif
 
     ipx=1
-!$omp parallel do schedule(dynamic,1) private(j,n,ic3,k,i,ic2,ipic)
+!$omp parallel do schedule(dynamic,1) private(j,n,ic3,k,i,ic2,ipic,ig,iaens)
     do n=1,n_ens
-       do ic3=1,nc3d
-          ipic=ipc3d(ic3)
-          do k=1,km
-             do j=1,jm
-                do i=1,im
-                      a_en(n)%r3(ipx)%q(i,j,k)=a_en(n)%r3(ipx)%q(i,j,k) &
-                            +cvec%r3(ipic)%q(i,j,k)*en_perts(n,ibin)%r3(ipic)%qr4(i,j,k)
-                enddo
-             enddo
-          enddo
-       enddo
-       do ic2=1,nc2d
-
-          ipic=ipc2d(ic2)
-          select case ( trim(StrUpCase(cvars2d(ic2))) )
- 
-             case('PS')
- 
-                if ( pwgtflg ) then
-                   km_tmp = km
-                else
-                   km_tmp = 1
-                endif
-
-                do k=1,km_tmp
+       do ig=1,ntotensgrp
+          do ic3=1,nc3d
+             iaens=ensgrp2aensgrp(ig,ic3,ibin)
+             if(iaens>0) then
+                ipic=ipc3d(ic3)
+                do k=1,km
                    do j=1,jm
                       do i=1,im
-                         a_en(n)%r3(ipx)%q(i,j,k)=a_en(n)%r3(ipx)%q(i,j,k) &
-                            +cvec%r2(ipic)%q(i,j)*en_perts(n,ibin)%r2(ipic)%qr4(i,j)*pwgt(i,j,k)
+                         a_en(iaens,n)%r3(ipx)%q(i,j,k)=a_en(iaens,n)%r3(ipx)%q(i,j,k) &
+                              +cvec%r3(ipic)%q(i,j,k)*en_perts(n,ig,ibin)%r3(ipic)%qr4(i,j,k)
                       enddo
                    enddo
                 enddo
+             endif ! iaens>0
+          enddo
+          do ic2=1,nc2d
+             iaens=ensgrp2aensgrp(ig,ic2+nc3d,ibin)
+             if(iaens>0) then
+                ipic=ipc2d(ic2)
+                select case ( trim(StrUpCase(cvars2d(ic2))) )
+ 
+                   case('PS')
+ 
+                      if ( pwgtflg ) then
+                         km_tmp = km
+                      else
+                         km_tmp = 1
+                      endif
+
+                      do k=1,km_tmp
+                         do j=1,jm
+                            do i=1,im
+                               a_en(iaens,n)%r3(ipx)%q(i,j,k)=a_en(iaens,n)%r3(ipx)%q(i,j,k) &
+                                    +cvec%r2(ipic)%q(i,j)*en_perts(n,ig,ibin)%r2(ipic)%qr4(i,j)*pwgt(i,j,k)
+                            enddo
+                         enddo
+                      enddo
   
-             case('SST')
+                   case('SST')
   
-                do j=1,jm
-                   do i=1,im
-                      a_en(n)%r3(ipx)%q(i,j,1)=a_en(n)%r3(ipx)%q(i,j,1) &
-                         +cvec%r2(ipic)%q(i,j)*en_perts(n,ibin)%r2(ipic)%qr4(i,j)
-                   enddo
-                enddo
+                      do j=1,jm
+                         do i=1,im
+                            a_en(iaens,n)%r3(ipx)%q(i,j,1)=a_en(iaens,n)%r3(ipx)%q(i,j,1) &
+                                 +cvec%r2(ipic)%q(i,j)*en_perts(n,ig,ibin)%r2(ipic)%qr4(i,j)
+                         enddo
+                      enddo
  
-          end select
+                end select
+             endif ! iaens>0
+          enddo
        enddo
     enddo ! enddo n_ens
     return
@@ -2212,20 +2290,22 @@ subroutine ensemble_forward_model_ad_dual_res(cvec,a_en,ibin)
     use hybrid_ensemble_parameters, only: n_ens,pwgtflg,pwgt
     use hybrid_ensemble_parameters, only: n_ens,grd_ens,grd_anl,p_e2a
     use hybrid_ensemble_parameters, only: en_perts
+    use hybrid_ensemble_parameters, only: ntotensgrp,naensgrp
+    use hybrid_ensemble_parameters, only: ensgrp2aensgrp
     use general_sub2grid_mod, only: general_sube2suba_ad
     use gridmod,only: regional
     use constants, only: zero
     implicit none
 
     type(gsi_bundle),intent(inout) :: cvec
-    type(gsi_bundle),intent(inout) :: a_en(n_ens)
+    type(gsi_bundle),intent(inout) :: a_en(naensgrp,n_ens)
     integer,intent(in)             :: ibin
 
     character(len=*),parameter::myname_=trim(myname)//'*ensemble_forward_model_ad_dual_res'
     type(gsi_grid)   :: grid_ens,grid_anl
     type(gsi_bundle) :: work_ens,work_anl
     integer(i_kind) :: i,j,k,n,im,jm,km,ic2,ic3,ipx,ipic,km_tmp
-    integer(i_kind) :: ipc2d(nc2d),ipc3d(nc3d),istatus
+    integer(i_kind) :: ipc2d(nc2d),ipc3d(nc3d),istatus,ig,iaens
 
 !   Request ensemble-corresponding fields from control vector
 !    NOTE:  because ensemble perturbation bundle structure is same as control vector, use same ipc3d and
@@ -2271,54 +2351,61 @@ subroutine ensemble_forward_model_ad_dual_res(cvec,a_en,ibin)
     endif
 
     ipx=1
-    im=a_en(1)%grid%im
-    jm=a_en(1)%grid%jm
-    km=a_en(1)%grid%km
-!$omp parallel do schedule(dynamic,1) private(j,n,ic3,k,i,ic2,ipic)
+    im=a_en(1,1)%grid%im
+    jm=a_en(1,1)%grid%jm
+    km=a_en(1,1)%grid%km
+!$omp parallel do schedule(dynamic,1) private(j,n,ic3,k,i,ic2,ipic,ig,iaens)
     do n=1,n_ens
-       do ic3=1,nc3d
-          ipic=ipc3d(ic3)
-          do k=1,km
-             do j=1,jm
-                do i=1,im
-                   a_en(n)%r3(ipx)%q(i,j,k)=a_en(n)%r3(ipx)%q(i,j,k) &
-                            +work_ens%r3(ipic)%q(i,j,k)*en_perts(n,ibin)%r3(ipic)%qr4(i,j,k)
-                enddo
-             enddo
-          enddo
-       enddo
-       do ic2=1,nc2d
-
-          ipic=ipc2d(ic2)
-          select case ( trim(StrUpCase(cvars2d(ic2))) )
-
-             case('PS')
-
-                if ( pwgtflg ) then
-                   km_tmp = km
-                else
-                   km_tmp = 1
-                endif
-
-                do k=1,km_tmp
+       do ig=1,ntotensgrp
+          do ic3=1,nc3d
+             iaens=ensgrp2aensgrp(ig,ic3,ibin)
+             if(iaens>0) then
+                ipic=ipc3d(ic3)
+                do k=1,km
                    do j=1,jm
                       do i=1,im
-                         a_en(n)%r3(ipx)%q(i,j,k)=a_en(n)%r3(ipx)%q(i,j,k) &
-                            +work_ens%r2(ipic)%q(i,j)*en_perts(n,ibin)%r2(ipic)%qr4(i,j)*pwgt(i,j,k)
+                         a_en(iaens,n)%r3(ipx)%q(i,j,k)=a_en(iaens,n)%r3(ipx)%q(i,j,k) &
+                              +work_ens%r3(ipic)%q(i,j,k)*en_perts(n,ig,ibin)%r3(ipic)%qr4(i,j,k)
                       enddo
                    enddo
                 enddo
+             endif ! iaens>0
+          enddo
+          do ic2=1,nc2d
+             iaens=ensgrp2aensgrp(ig,ic2+nc3d,ibin)
+             if(iaens>0) then
+                ipic=ipc2d(ic2)
+                select case ( trim(StrUpCase(cvars2d(ic2))) )
+
+                   case('PS')
+
+                      if ( pwgtflg ) then
+                         km_tmp = km
+                      else
+                         km_tmp = 1
+                      endif
+
+                      do k=1,km_tmp
+                         do j=1,jm
+                            do i=1,im
+                               a_en(iaens,n)%r3(ipx)%q(i,j,k)=a_en(iaens,n)%r3(ipx)%q(i,j,k) &
+                                    +work_ens%r2(ipic)%q(i,j)*en_perts(n,ig,ibin)%r2(ipic)%qr4(i,j)*pwgt(i,j,k)
+                            enddo
+                         enddo
+                      enddo
 
-             case('SST')
+                   case('SST')
 
-                do j=1,jm
-                   do i=1,im
-                      a_en(n)%r3(ipx)%q(i,j,1)=a_en(n)%r3(ipx)%q(i,j,1) &
-                         +work_ens%r2(ipic)%q(i,j)*en_perts(n,ibin)%r2(ipic)%qr4(i,j)
-                   enddo
-                enddo
+                      do j=1,jm
+                         do i=1,im
+                            a_en(iaens,n)%r3(ipx)%q(i,j,1)=a_en(iaens,n)%r3(ipx)%q(i,j,1) &
+                                 +work_ens%r2(ipic)%q(i,j)*en_perts(n,ig,ibin)%r2(ipic)%qr4(i,j)
+                         enddo
+                      enddo
 
-          end select
+                end select
+             endif ! iaens>0
+          enddo
        enddo
     enddo ! enddo n_ens
     call gsi_bundledestroy(work_ens,istatus)
@@ -2758,6 +2845,7 @@ subroutine sqrt_beta_e_mult_cvec(grady)
   use kinds, only: r_kind,i_kind
   use gsi_4dvar, only: nsubwin
   use hybrid_ensemble_parameters, only: n_ens,sqrt_beta_e,grd_ens
+  use hybrid_ensemble_parameters, only: naensgrp
   use control_vectors,only: control_vector
   use timermod, only: timer_ini,timer_fnl
 
@@ -2770,19 +2858,21 @@ subroutine sqrt_beta_e_mult_cvec(grady)
 
 ! Declare local variables
   character(len=*),parameter::myname_=myname//'*sqrt_beta_e_mult'
-  integer(i_kind) :: i,j,k,ii,nn
+  integer(i_kind) :: i,j,k,ii,nn,ig
 
   ! Initialize timer
   call timer_ini('sqrt_beta_e_mult')
 
   ! multiply by sqrt_beta_e
-!$omp parallel do schedule(dynamic,1) private(nn,k,j,i,ii)
+!$omp parallel do schedule(dynamic,1) private(nn,k,j,i,ii,ig)
   do j=1,grd_ens%lon2
      do ii=1,nsubwin
-        do nn=1,n_ens
-           do k=1,nsig
-              do i=1,grd_ens%lat2
-                 grady%aens(ii,nn)%r3(1)%q(i,j,k) = sqrt_beta_e(k)*grady%aens(ii,nn)%r3(1)%q(i,j,k)
+        do ig=1,naensgrp
+           do nn=1,n_ens
+              do k=1,nsig
+                 do i=1,grd_ens%lat2
+                    grady%aens(ii,ig,nn)%r3(1)%q(i,j,k) = sqrt_beta_e(k)*grady%aens(ii,ig,nn)%r3(1)%q(i,j,k)
+                 enddo
               enddo
            enddo
         enddo
@@ -2823,6 +2913,7 @@ subroutine sqrt_beta_e_mult_bundle(aens)
 !$$$ end documentation block
   use kinds, only: r_kind,i_kind
   use hybrid_ensemble_parameters, only: n_ens,sqrt_beta_e,grd_ens
+  use hybrid_ensemble_parameters, only: naensgrp
   use gsi_bundlemod, only: gsi_bundle
   use timermod, only: timer_ini,timer_fnl
   use gridmod, only: nsig
@@ -2830,22 +2921,24 @@ subroutine sqrt_beta_e_mult_bundle(aens)
   implicit none
 
 ! Declare passed variables
-  type(gsi_bundle),intent(inout) :: aens(n_ens)
+  type(gsi_bundle),intent(inout) :: aens(naensgrp,n_ens)
 
 ! Declare local variables
   character(len=*),parameter::myname_=myname//'*sqrt_beta_e_mult'
-  integer(i_kind) :: i,j,k,nn
+  integer(i_kind) :: i,j,k,nn,ig
 
   ! Initialize timer
   call timer_ini('sqrt_beta_e_mult')
 
   ! multiply by sqrt_beta_e
-!$omp parallel do schedule(dynamic,1) private(nn,k,j,i)
+!$omp parallel do schedule(dynamic,1) private(nn,k,j,i,ig)
   do j=1,grd_ens%lon2
-     do nn=1,n_ens
-        do k=1,nsig
-           do i=1,grd_ens%lat2
-              aens(nn)%r3(1)%q(i,j,k) = sqrt_beta_e(k)*aens(nn)%r3(1)%q(i,j,k)
+     do ig=1,naensgrp
+        do nn=1,n_ens
+           do k=1,nsig
+              do i=1,grd_ens%lat2
+                 aens(ig,nn)%r3(1)%q(i,j,k) = sqrt_beta_e(k)*aens(ig,nn)%r3(1)%q(i,j,k)
+              enddo
            enddo
         enddo
      enddo
@@ -2885,6 +2978,7 @@ subroutine init_sf_xy(jcap_in)
   use kinds, only: r_kind,i_kind,r_single
   use hybrid_ensemble_parameters,only: s_ens_hv,sp_loc,grd_ens,grd_loc,sp_ens,n_ens,p_sploc2ens,grd_sploc
   use hybrid_ensemble_parameters,only: use_localization_grid
+  use hybrid_ensemble_parameters,only: naensgrp,naensloc
   use gridmod,only: use_sp_eqspace
   use general_specmod, only: general_init_spec_vars
   use constants, only: zero,half,one,two,three,rearth,pi
@@ -2899,7 +2993,7 @@ subroutine init_sf_xy(jcap_in)
 
   integer(i_kind),intent(in   ) :: jcap_in
 
-  integer(i_kind) i,ii,j,k,l,n,jcap,kk,nsigend
+  integer(i_kind) i,ii,j,k,l,n,jcap,kk,nsigend,ig
   real(r_kind),allocatable::g(:),gsave(:)
   real(r_kind) factor
   real(r_kind),allocatable::rkm(:),f(:,:),f0(:,:)
@@ -2925,16 +3019,18 @@ subroutine init_sf_xy(jcap_in)
 
   s_ens_h_min=pi*rearth*.001_r_kind/jcap_in
   do k=1,grd_ens%nsig
-     if(s_ens_hv(k) <  s_ens_h_min) then
-        if(mype == 0) write(6,*)' s_ens_hv(',k,') = ',s_ens_hv(k),' km--too small, min value = ', &
-                                        s_ens_h_min,' km.'
-        if(mype == 0) write(6,*)' s_ens_hv(',k,') reset to min value'
-        s_ens_hv(k)=s_ens_h_min
-     else if(s_ens_hv(k) >  5500._r_kind) then
-        if(mype == 0) write(6,*)' s_ens_hv(',k,') = ',s_ens_hv(k),' km--too large, max value = 5500 km.'
-        if(mype == 0) write(6,*)' s_ens_hv(',k,') reset to max value'
-        s_ens_hv(k)=5500._r_kind
-     end if
+     do ig=1,naensloc
+        if(s_ens_hv(k,ig) <  s_ens_h_min) then
+           if(mype == 0) write(6,*)' s_ens_hv(',k,') = ',s_ens_hv(k,ig),' km--too small, min value = ', &
+                                           s_ens_h_min,' km.'
+           if(mype == 0) write(6,*)' s_ens_hv(',k,') reset to min value'
+           s_ens_hv(k,ig)=s_ens_h_min
+        else if(s_ens_hv(k,ig) >  5500._r_kind) then
+           if(mype == 0) write(6,*)' s_ens_hv(',k,') = ',s_ens_hv(k,ig),' km--too large, max value = 5500 km.'
+           if(mype == 0) write(6,*)' s_ens_hv(',k,') reset to max value'
+           s_ens_hv(k,ig)=5500._r_kind
+        end if
+     enddo
   enddo
 
 
@@ -3059,99 +3155,107 @@ subroutine init_sf_xy(jcap_in)
            rkm(1+(grd_sploc%nlat-2)/2), &
           -rkm(grd_sploc%nlat-(grd_sploc%nlat-2)/2)+rkm(1+(grd_sploc%nlat-2)/2),' km'
 
-  if(.not.allocated(spectral_filter)) allocate(spectral_filter(sp_loc%nc,grd_sploc%nsig))
-  if(.not.allocated(sqrt_spectral_filter)) allocate(sqrt_spectral_filter(sp_loc%nc,grd_sploc%nsig))
+  if(.not.allocated(spectral_filter)) allocate(spectral_filter(naensloc,sp_loc%nc,grd_sploc%nsig))
+  if(.not.allocated(sqrt_spectral_filter)) allocate(sqrt_spectral_filter(naensloc,sp_loc%nc,grd_sploc%nsig))
   allocate(g(sp_loc%nc),gsave(sp_loc%nc))
   allocate(pn0_npole(0:sp_loc%jcap))
   allocate(ksame(grd_sploc%nsig))
-  ksame=.false.
-  do k=2,grd_sploc%nsig
-     if(s_ens_hv(k) == s_ens_hv(k-1))ksame(k)=.true.
-  enddo
-  spectral_filter=zero
-  do k=1,grd_sploc%nsig
-     if(ksame(k))then
-        spectral_filter(:,k)=spectral_filter(:,k-1)
-     else
-        do i=1,grd_sploc%nlat
-           f0(i,1)=exp(-half*(rkm(i)/s_ens_hv(k))**2)
-        enddo
-
-        do j=2,grd_sploc%nlon
+  do ig=1,naensloc
+     ksame=.false.
+     do k=2,grd_sploc%nsig
+        if(s_ens_hv(k,ig) == s_ens_hv(k-1,ig))ksame(k)=.true.
+     enddo
+     spectral_filter(ig,:,:)=zero
+     do k=1,grd_sploc%nsig
+        if(ksame(k))then
+           spectral_filter(ig,:,k)=spectral_filter(ig,:,k-1)
+        else
            do i=1,grd_sploc%nlat
-              f0(i,j)=f0(i,1)
+              f0(i,1)=exp(-half*(rkm(i)/s_ens_hv(k,ig))**2)
            enddo
-        enddo
-
-        call general_g2s0(grd_sploc,sp_loc,g,f0)
-
-        call general_s2g0(grd_sploc,sp_loc,g,f)
 
-!       adjust so value at np = 1
-        f=f/f(grd_sploc%nlat,1)
-        f0=f
-        call general_g2s0(grd_sploc,sp_loc,g,f)
-        call general_s2g0(grd_sploc,sp_loc,g,f)
-        if(mype == 0)then
-           nsigend=k
-           do kk=k+1,grd_sploc%nsig
-              if(s_ens_hv(kk) /= s_ens_hv(k))exit
-              nsigend=nsigend+1
+           do j=2,grd_sploc%nlon
+              do i=1,grd_sploc%nlat
+                 f0(i,j)=f0(i,1)
+              enddo
            enddo
-           write(6,900)k,nsigend,sp_loc%jcap,s_ens_hv(k),maxval(abs(f0-f))
-  900      format(' in init_sf_xy, jcap,s_ens_hv(',i5,1x,'-',i5,'), max diff(f0-f)=', &
-                                        i10,f10.2,e20.10)
-        end if
 
-!            correct spectrum by dividing by pn0_npole
-        gsave=g
+           call general_g2s0(grd_sploc,sp_loc,g,f0)
 
-!       obtain pn0_npole
-        do n=0,sp_loc%jcap
-           g=zero
-           g(2*n+1)=one
            call general_s2g0(grd_sploc,sp_loc,g,f)
-           pn0_npole(n)=f(grd_sploc%nlat,1)
-        enddo
+
+!          adjust so value at np = 1
+           f=f/f(grd_sploc%nlat,1)
+           f0=f
+           call general_g2s0(grd_sploc,sp_loc,g,f)
+           call general_s2g0(grd_sploc,sp_loc,g,f)
+           if(mype == 0)then
+              nsigend=k
+              do kk=k+1,grd_sploc%nsig
+                 if(s_ens_hv(kk,ig) /= s_ens_hv(k,ig))exit
+                 nsigend=nsigend+1
+              enddo
+              write(6,900)k,nsigend,sp_loc%jcap,s_ens_hv(k,ig),maxval(abs(f0-f))
+900           format(' in init_sf_xy, jcap,s_ens_hv(',i5,1x,'-',i5,'), max diff(f0-f)=', &
+                                           i10,f10.2,e20.10)
+           end if
+
+!          correct spectrum by dividing by pn0_npole
+           gsave=g
+
+!          obtain pn0_npole
+           do n=0,sp_loc%jcap
+              g=zero
+              g(2*n+1)=one
+              call general_s2g0(grd_sploc,sp_loc,g,f)
+              pn0_npole(n)=f(grd_sploc%nlat,1)
+           enddo
    
-        g=zero
-        do n=0,sp_loc%jcap
-           g(2*n+1)=gsave(2*n+1)/pn0_npole(n)
-        enddo
+           g=zero
+           do n=0,sp_loc%jcap
+              g(2*n+1)=gsave(2*n+1)/pn0_npole(n)
+           enddo
 
-!       obtain spectral_filter
+!          obtain spectral_filter
 
-        ii=0
-        do l=0,sp_loc%jcap
-           factor=one
-           if(l >  0) factor=half
-           do n=l,sp_loc%jcap
-              ii=ii+1
-              if(sp_loc%factsml(ii)) then
-                 spectral_filter(ii,k)=zero
-              else
-                 spectral_filter(ii,k)=factor*g(2*n+1)
-              end if
-              ii=ii+1
-              if(l == 0 .or. sp_loc%factsml(ii)) then
-                 spectral_filter(ii,k)=zero
+           ii=0
+           do l=0,sp_loc%jcap
+              if(ig>naensgrp) then
+                 factor=one/g(1)
               else
-                 spectral_filter(ii,k)=factor*g(2*n+1)
+                 factor=one
+                 if(l>0) factor=half
               end if
+              do n=l,sp_loc%jcap
+                 ii=ii+1
+                 if(sp_loc%factsml(ii)) then
+                    spectral_filter(ig,ii,k)=zero
+                 else
+                    spectral_filter(ig,ii,k)=factor*g(2*n+1)
+                 end if
+                 ii=ii+1
+                 if(l == 0 .or. sp_loc%factsml(ii)) then
+                    spectral_filter(ig,ii,k)=zero
+                 else
+                    spectral_filter(ig,ii,k)=factor*g(2*n+1)
+                 end if
+              enddo
            enddo
-        enddo
-     end if
-  enddo
+        end if
+     enddo
+  enddo !ig loop
   deallocate(g,gsave,pn0_npole,ksame)
 
 ! Compute sqrt(spectral_filter).  Ensure spectral_filter >=0 zero
 !$omp parallel do schedule(dynamic,1) private(k,i)
-  do k=1,grd_sploc%nsig
-     do i=1,sp_loc%nc
-        if (spectral_filter(i,k) < zero) spectral_filter(i,k)=zero
-        sqrt_spectral_filter(i,k) = sqrt(spectral_filter(i,k))
+  do ig=1,naensloc
+     do k=1,grd_sploc%nsig
+        do i=1,sp_loc%nc
+           if (spectral_filter(ig,i,k) < zero) spectral_filter(ig,i,k)=zero
+           sqrt_spectral_filter(ig,i,k) = sqrt(spectral_filter(ig,i,k))
+        end do
      end do
-  end do
+  enddo !ig loop
 
 !  assign array k_index for each processor, based on grd_loc%kbegin_loc,grd_loc%kend_loc
 
@@ -3167,7 +3271,7 @@ subroutine init_sf_xy(jcap_in)
    do k=grd_loc%kbegin_loc,grd_loc%kend_loc
       ftest(grd_ens%nlat/2,grd_ens%nlon/2,k)=one
    enddo
-   call sf_xy(ftest,grd_loc%kbegin_loc,grd_loc%kend_loc)
+   call sf_xy(1,ftest,grd_loc%kbegin_loc,grd_loc%kend_loc)
    if(mype==0) then
       do j=1,grd_ens%nlon
         do i=1,grd_ens%nlat
@@ -3183,7 +3287,7 @@ subroutine init_sf_xy(jcap_in)
 
 end subroutine init_sf_xy
 
-subroutine sf_xy(f,k_start,k_end)
+subroutine sf_xy(ig,f,k_start,k_end)
 !$$$  subprogram documentation block
 !                .      .    .
 ! subprogram:    sf_xy       spectral isotropic localization for global domain
@@ -3198,6 +3302,7 @@ subroutine sf_xy(f,k_start,k_end)
 !   2010-03-11  parrish - adjust dimensions for f to allow for nlevs=0
 !
 !   input argument list:
+!     ig       - number for smoothing scales
 !     f        - input field to be filtered
 !     k_start  - starting horizontal slab index
 !     k_end    - ending horizontal slab index    (k_end can be less than k_start, meaning there is
@@ -3216,39 +3321,60 @@ subroutine sf_xy(f,k_start,k_end)
   use kinds, only: r_kind,i_kind
   use hybrid_ensemble_parameters, only: grd_ens,sp_loc,p_sploc2ens,grd_sploc
   use hybrid_ensemble_parameters,only: use_localization_grid
-  use egrid2agrid_mod,only: g_egrid2agrid,g_egrid2agrid_ad  
+  use hybrid_ensemble_parameters,only: naensgrp
+  use egrid2agrid_mod,only: g_egrid2agrid,g_egrid2agrid_ad,g_agrid2egrid
   implicit none
 
+  integer(i_kind),intent(in   ) :: ig
   integer(i_kind),intent(in   ) :: k_start,k_end
   real(r_kind)   ,intent(inout) :: f(grd_ens%nlat,grd_ens%nlon,k_start:max(k_start,k_end))
 
+  real(r_kind) g(sp_loc%nc)
   real(r_kind) work(grd_sploc%nlat,grd_sploc%nlon,1)
   integer(i_kind) k
   logical vector(k_start:max(k_start,k_end))
 
   if(.not.use_localization_grid) then
 
+    if(ig>naensgrp) then
+       do k=k_start,k_end
+          call general_g2s0(grd_ens,sp_loc,g,f(:,:,k))
+          g(:)=g(:)*spectral_filter(ig,:,k_index(k))
+          call general_s2g0(grd_ens,sp_loc,g,f(:,:,k))
+       enddo
+    else
 !$omp parallel do schedule(dynamic,1) private(k)
-    do k=k_start,k_end
-       call sfilter(grd_ens,sp_loc,spectral_filter(1,k_index(k)),f(1,1,k))
-    enddo
+       do k=k_start,k_end
+          call sfilter(grd_ens,sp_loc,spectral_filter(ig,:,k_index(k)),f(1,1,k))
+       enddo
+    endif
 
   else
 
     vector=.false.
+    if(ig>naensgrp) then
+       do k=k_start,k_end
+          call g_agrid2egrid(p_sploc2ens,work,f(:,:,k:k),k,k,vector(k:k))
+          call general_g2s0(grd_ens,sp_loc,g,f(:,:,k))
+          g(:)=g(:)*spectral_filter(ig,:,k_index(k))
+          call general_s2g0(grd_ens,sp_loc,g,f(:,:,k))
+          call g_egrid2agrid(p_sploc2ens,work,f(:,:,k:k),k,k,vector(k:k))
+       enddo
+    else
 !$omp parallel do schedule(dynamic,1) private(k,work)
-    do k=k_start,k_end
-       call g_egrid2agrid_ad(p_sploc2ens,work,f(:,:,k:k),k,k,vector(k:k))
-       call sfilter(grd_ens,sp_loc,spectral_filter(:,k_index(k)),f(1,1,k))
-       call g_egrid2agrid(p_sploc2ens,work,f(:,:,k:k),k,k,vector(k:k))
-    enddo
+       do k=k_start,k_end
+          call g_egrid2agrid_ad(p_sploc2ens,work,f(:,:,k:k),k,k,vector(k:k))
+          call sfilter(grd_ens,sp_loc,spectral_filter(ig,:,k_index(k)),f(1,1,k))
+          call g_egrid2agrid(p_sploc2ens,work,f(:,:,k:k),k,k,vector(k:k))
+       enddo
+    end if
 
   endif
   return
 
 end subroutine sf_xy
 
-subroutine sqrt_sf_xy(z,f,k_start,k_end)
+subroutine sqrt_sf_xy(ig,z,f,k_start,k_end)
 !$$$  subprogram documentation block
 !                .      .    .
 ! subprogram:    sqrt_sf_xy  sqrt(sf_xy)
@@ -3261,6 +3387,7 @@ subroutine sqrt_sf_xy(z,f,k_start,k_end)
 !   2011-06-28  parrish  initial documentation
 !
 !   input argument list:
+!     ig       - number for smoothing scales
 !     z        - input spectral space variable
 !     k_start  - starting horizontal slab index
 !     k_end    - ending horizontal slab index    (k_end can be less than k_start, meaning there is
@@ -3282,6 +3409,7 @@ subroutine sqrt_sf_xy(z,f,k_start,k_end)
   use egrid2agrid_mod,only: g_egrid2agrid
   implicit none
 
+  integer(i_kind),intent(in   ) :: ig
   integer(i_kind),intent(in   ) :: k_start,k_end
   real(r_kind)   ,intent(in   ) :: z(sp_loc%nc,k_start:max(k_start,k_end))
   real(r_kind)   ,intent(  out) :: f(grd_ens%nlat,grd_ens%nlon,k_start:max(k_start,k_end))
@@ -3294,7 +3422,7 @@ subroutine sqrt_sf_xy(z,f,k_start,k_end)
   if(.not.use_localization_grid) then
 
     do k=k_start,k_end
-       g(:)=z(:,k)*sqrt_spectral_filter(:,k_index(k))
+       g(:)=z(:,k)*sqrt_spectral_filter(ig,:,k_index(k))
        call general_s2g0(grd_ens,sp_loc,g,f(:,:,k))
     enddo
 
@@ -3302,7 +3430,7 @@ subroutine sqrt_sf_xy(z,f,k_start,k_end)
 
      vector=.false.
      do k=k_start,k_end
-        g(:)=z(:,k)*sqrt_spectral_filter(:,k_index(k))
+        g(:)=z(:,k)*sqrt_spectral_filter(ig,:,k_index(k))
         call general_s2g0(grd_sploc,sp_loc,g,work)
         call g_egrid2agrid(p_sploc2ens,work,f(:,:,k:k),k,k,vector(k:k))
      enddo
@@ -3312,7 +3440,7 @@ subroutine sqrt_sf_xy(z,f,k_start,k_end)
 
 end subroutine sqrt_sf_xy
 
-subroutine sqrt_sf_xy_ad(z,f,k_start,k_end)
+subroutine sqrt_sf_xy_ad(ig,z,f,k_start,k_end)
 !$$$  subprogram documentation block
 !                .      .    .
 ! subprogram:    sqrt_sf_xy_ad  adjoint of sqrt_sf_xy
@@ -3325,6 +3453,7 @@ subroutine sqrt_sf_xy_ad(z,f,k_start,k_end)
 !   2011-06-28  parrish  initial documentation
 !
 !   input argument list:
+!     ig       - number for smoothing scales
 !     f        - grid space variable
 !     k_start  - starting horizontal slab index
 !     k_end    - ending horizontal slab index    (k_end can be less than k_start, meaning there is
@@ -3347,6 +3476,7 @@ subroutine sqrt_sf_xy_ad(z,f,k_start,k_end)
   use egrid2agrid_mod,only: g_egrid2agrid_ad
   implicit none
 
+  integer(i_kind),intent(in   ) :: ig
   integer(i_kind),intent(in   ) :: k_start,k_end
   real(r_kind)   ,intent(inout) :: z(sp_loc%nc,k_start:max(k_start,k_end))
   real(r_kind)   ,intent(inout) :: f(grd_ens%nlat,grd_ens%nlon,k_start:max(k_start,k_end))
@@ -3360,7 +3490,7 @@ subroutine sqrt_sf_xy_ad(z,f,k_start,k_end)
 
     do k=k_start,k_end
        call general_s2g0_ad(grd_ens,sp_loc,g,f(:,:,k))
-       z(:,k)=g(:)*sqrt_spectral_filter(:,k_index(k))
+       z(:,k)=g(:)*sqrt_spectral_filter(ig,:,k_index(k))
     enddo
 
   else
@@ -3369,7 +3499,7 @@ subroutine sqrt_sf_xy_ad(z,f,k_start,k_end)
      do k=k_start,k_end
         call g_egrid2agrid_ad(p_sploc2ens,work,f(:,:,k:k),k,k,vector(k:k))
         call general_s2g0_ad(grd_sploc,sp_loc,g,work)
-        z(:,k)=g(:)*sqrt_spectral_filter(:,k_index(k))
+        z(:,k)=g(:)*sqrt_spectral_filter(ig,:,k_index(k))
      enddo
 
   end if
@@ -3447,6 +3577,7 @@ subroutine bkerror_a_en(grady)
 ! program history log:
 !   2009-09-17  parrish  initial creation of code from a copy of bkerror
 !   2010-05-20  todling  update to use bundle
+!   2022-09-15  yokota   add scale/variable/time-dependent localization
 !
 !   input argument list:
 !     grady    - input field  
@@ -3460,10 +3591,14 @@ subroutine bkerror_a_en(grady)
 !
 !$$$ end documentation block
   use kinds, only: r_kind,i_kind
+  use constants, only: zero
   use gsi_4dvar, only: nsubwin, lsqrtb
   use control_vectors, only: control_vector
   use timermod, only: timer_ini,timer_fnl
   use hybrid_ensemble_parameters, only: n_ens
+  use hybrid_ensemble_parameters, only: naensgrp
+  use hybrid_ensemble_parameters, only: alphacvarsclgrpmat
+  use hybrid_ensemble_parameters, only: nval_lenz_en
   use gsi_bundlemod,only: gsi_bundlegetpointer
   implicit none
 
@@ -3471,7 +3606,9 @@ subroutine bkerror_a_en(grady)
   type(control_vector),intent(inout) :: grady
 
 ! Declare local variables
-  integer(i_kind) ii,ip,istatus
+  integer(i_kind) ii,ip,istatus,k,ig,ig2
+  real(r_kind),allocatable,dimension(:,:) :: z
+  real(r_kind),allocatable,dimension(:) :: ztmp
 
   if (lsqrtb) then
      write(6,*)'bkerror_a_en: not for use with lsqrtb'
@@ -3482,7 +3619,7 @@ subroutine bkerror_a_en(grady)
   call timer_ini('bkerror_a_en')
 
 ! Put things in grady first since operations change input variables
-  call gsi_bundlegetpointer ( grady%aens(1,1),'a_en',ip,istatus)
+  call gsi_bundlegetpointer ( grady%aens(1,1,1),'a_en',ip,istatus)
   if(istatus/=0) then
      write(6,*)'bkerror_a_en: trouble getting pointer to ensemble CV'
      call stop2(317)
@@ -3493,16 +3630,31 @@ subroutine bkerror_a_en(grady)
 
 ! Apply variances, as well as vertical & horizontal parts of background error
   do ii=1,nsubwin
-    !if(test_sqrt_localization) then
-    !        write(6,*)' using ckgcov_a_en_new_factorization'
-    !   allocate(z(nval_lenz_en))
-    !   call ckgcov_a_en_new_factorization_ad(z,grady%aens(ii,1:n_ens))
-    !   call ckgcov_a_en_new_factorization   (z,grady%aens(ii,1:n_ens))
-    !   deallocate(z)
-    !else
-    !        write(6,*)' using bkgcov_a_en_new_factorization'
-        call bkgcov_a_en_new_factorization(grady%aens(ii,1:n_ens))
-    !end if
+     if (naensgrp==1) then
+        call bkgcov_a_en_new_factorization(1,grady%aens(ii,1,1:n_ens))
+     else
+        allocate(z(naensgrp,nval_lenz_en))
+        do ig=1,naensgrp
+           call ckgcov_a_en_new_factorization_ad(ig,z(ig,:),grady%aens(ii,ig,1:n_ens))
+        enddo
+        allocate(ztmp(naensgrp))
+        do k=1,nval_lenz_en
+           ztmp=zero
+           do ig=1,naensgrp
+              do ig2=1,naensgrp
+                 ztmp(ig) = ztmp(ig) + z(ig2,k) * alphacvarsclgrpmat(ig,ig2)  
+              enddo
+           enddo
+           do ig=1,naensgrp
+              z(ig,k) = ztmp(ig)
+           enddo
+        enddo
+        deallocate(ztmp)
+        do ig=1,naensgrp
+           call ckgcov_a_en_new_factorization(ig,z(ig,:),grady%aens(ii,ig,1:n_ens))
+        enddo
+        deallocate(z)
+     endif
   enddo
 
 !  multiply by sqrt_beta_e_mult
@@ -3514,7 +3666,7 @@ subroutine bkerror_a_en(grady)
   return
 end subroutine bkerror_a_en
 
-subroutine bkgcov_a_en_new_factorization(a_en)
+subroutine bkgcov_a_en_new_factorization(ig,a_en)
 !$$$  subprogram documentation block
 !                .      .    .                                       .
 ! subprogram:    bkgcov_a_en copy of bkgcov for hybrid ens var a_en 
@@ -3531,6 +3683,7 @@ subroutine bkgcov_a_en_new_factorization(a_en)
 !                 of horizontal localization length
 !
 !   input argument list:
+!     ig       - number for smoothing scales
 !     a_en     - control variable for ensemble contribution to background error
 !     nlevs    - number of vertical levels for smoothing
 !
@@ -3553,6 +3706,7 @@ subroutine bkgcov_a_en_new_factorization(a_en)
 
 ! Passed Variables
 ! real(r_kind),dimension(grd_loc%latlon1n,n_ens),intent(inout) :: a_en
+  integer(i_kind),intent(in   ) :: ig
   type(gsi_bundle),intent(inout) :: a_en(n_ens)
 
 ! Local Variables
@@ -3579,7 +3733,7 @@ subroutine bkgcov_a_en_new_factorization(a_en)
   iadvance=1 ; iback=2
 !$omp parallel do schedule(dynamic,1) private(k,ii,is,ie)
   do k=1,n_ens
-     call new_factorization_rf_z(a_en(k)%r3(ipnt)%q,iadvance,iback)
+     call new_factorization_rf_z(a_en(k)%r3(ipnt)%q,iadvance,iback,ig)
      ii=(k-1)*a_en(1)%ndim
      is=ii+1
      ie=ii+a_en(1)%ndim
@@ -3592,13 +3746,13 @@ subroutine bkgcov_a_en_new_factorization(a_en)
 ! Apply horizontal smoother for number of horizontal scales
   if(regional) then
      iadvance=1 ; iback=2
-     call new_factorization_rf_x(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
-     call new_factorization_rf_y(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
+     call new_factorization_rf_x(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc,ig)
+     call new_factorization_rf_y(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc,ig)
      iadvance=2 ; iback=1
-     call new_factorization_rf_y(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
-     call new_factorization_rf_x(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
+     call new_factorization_rf_y(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc,ig)
+     call new_factorization_rf_x(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc,ig)
   else
-     call sf_xy(hwork,grd_loc%kbegin_loc,grd_loc%kend_loc)
+     call sf_xy(ig,hwork,grd_loc%kbegin_loc,grd_loc%kend_loc)
   end if
 
 ! Put back onto subdomains
@@ -3613,14 +3767,14 @@ subroutine bkgcov_a_en_new_factorization(a_en)
      is=ii+1
      ie=ii+a_en(1)%ndim
      a_en(k)%values(1:a_en(k)%ndim)=a_en_work(is:ie)
-     call new_factorization_rf_z(a_en(k)%r3(ipnt)%q,iadvance,iback)
+     call new_factorization_rf_z(a_en(k)%r3(ipnt)%q,iadvance,iback,ig)
   enddo
   deallocate(a_en_work)
 
   return
 end subroutine bkgcov_a_en_new_factorization
 
-subroutine ckgcov_a_en_new_factorization(z,a_en)
+subroutine ckgcov_a_en_new_factorization(ig,z,a_en)
 !$$$  subprogram documentation block
 !                .      .    .                                       .
 ! subprogram:    ckgcov_a_en_new_factorization sqrt(bkgcov_a_en_new_factorization)
@@ -3632,6 +3786,7 @@ subroutine ckgcov_a_en_new_factorization(z,a_en)
 !   2011-06-27  parrish, initial documentation
 !
 !   input argument list:
+!     ig       - number for smoothing scales
 !     z        - long vector containing sqrt control vector for ensemble extended control variable
 !
 !   output argument list:
@@ -3654,6 +3809,7 @@ subroutine ckgcov_a_en_new_factorization(z,a_en)
   implicit none
 
 ! Passed Variables
+  integer(i_kind),intent(in   ) :: ig
   type(gsi_bundle),intent(inout) :: a_en(n_ens)
   real(r_kind),dimension(nval_lenz_en),intent(in   ) :: z
 
@@ -3685,11 +3841,11 @@ subroutine ckgcov_a_en_new_factorization(z,a_en)
 ! Make a copy of input variable z to hwork
         hwork=z
         iadvance=2 ; iback=1
-        call new_factorization_rf_y(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
-        call new_factorization_rf_x(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
+        call new_factorization_rf_y(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc,ig)
+        call new_factorization_rf_x(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc,ig)
      else
 #ifdef LATER
-        call sqrt_sf_xy(z,hwork,grd_loc%kbegin_loc,grd_loc%kend_loc)
+        call sqrt_sf_xy(ig,z,hwork,grd_loc%kbegin_loc,grd_loc%kend_loc)
 #else
         write(6,*) ' problem with ibm compiler with "use hybrid_ensemble_isotropic, only: sqrt_sf_xy"'
 #endif /*LATER*/
@@ -3718,14 +3874,14 @@ subroutine ckgcov_a_en_new_factorization(z,a_en)
   do k=1,n_ens
 
      iadvance=2 ; iback=1
-     call new_factorization_rf_z(a_en(k)%r3(ipnt)%q,iadvance,iback)
+     call new_factorization_rf_z(a_en(k)%r3(ipnt)%q,iadvance,iback,ig)
 
   enddo
 
   return
 end subroutine ckgcov_a_en_new_factorization
 
-subroutine ckgcov_a_en_new_factorization_ad(z,a_en)
+subroutine ckgcov_a_en_new_factorization_ad(ig,z,a_en)
 !$$$  subprogram documentation block
 !                .      .    .                                       .
 ! subprogram:    ckgcov_a_en_new_factorization_ad adjoint of ckgcov_a_en_new_factorization
@@ -3739,6 +3895,7 @@ subroutine ckgcov_a_en_new_factorization_ad(z,a_en)
 !   2011-06-27  parrish, initial documentation
 !
 !   input argument list:
+!     ig       - number for smoothing scales
 !     z        - long vector containing sqrt control vector for ensemble extended control variable
 !     a_en     - bundle containing intermediate control variable after multiplication by sqrt(S), the
 !                    ensemble localization correlation.
@@ -3764,6 +3921,7 @@ subroutine ckgcov_a_en_new_factorization_ad(z,a_en)
   implicit none
 
 ! Passed Variables
+  integer(i_kind),intent(in   ) :: ig
   type(gsi_bundle),intent(inout) :: a_en(n_ens)
   real(r_kind),dimension(nval_lenz_en),intent(inout) :: z
 
@@ -3788,7 +3946,7 @@ subroutine ckgcov_a_en_new_factorization_ad(z,a_en)
   do k=1,n_ens
 
      iadvance=1 ; iback=2
-     call new_factorization_rf_z(a_en(k)%r3(ipnt)%q,iadvance,iback)
+     call new_factorization_rf_z(a_en(k)%r3(ipnt)%q,iadvance,iback,ig)
  
   enddo
 
@@ -3819,11 +3977,11 @@ subroutine ckgcov_a_en_new_factorization_ad(z,a_en)
 ! Apply horizontal smoother for number of horizontal scales
      if(regional) then
         iadvance=1 ; iback=2
-        call new_factorization_rf_x(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
-        call new_factorization_rf_y(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc)
+        call new_factorization_rf_x(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc,ig)
+        call new_factorization_rf_y(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc,ig)
         z=hwork
      else
-        call sqrt_sf_xy_ad(z,hwork,grd_loc%kbegin_loc,grd_loc%kend_loc)
+        call sqrt_sf_xy_ad(ig,z,hwork,grd_loc%kbegin_loc,grd_loc%kend_loc)
      end if
   end if
 
@@ -3989,6 +4147,7 @@ subroutine hybens_localization_setup
 !   12-05-2012  el akkraoui  hybrid beta parameters now vertically varying
 !   2012-10-16  wu - only call setup_ens_wgt if necessary
 !   2014-05-22  wu  modification to allow vertically varying localization scales in regional
+!   2022-09-15  yokota - add scale/variable/time-dependent localization
 !
 !   input argument list:
 !
@@ -4009,6 +4168,8 @@ subroutine hybens_localization_setup
    use hybrid_ensemble_parameters, only: readin_beta,beta_s,beta_e,beta_s0,beta_e0,sqrt_beta_s,sqrt_beta_e
    use hybrid_ensemble_parameters, only: readin_localization,create_hybens_localization_parameters, &
                                          vvlocal,s_ens_h,s_ens_hv,s_ens_v,s_ens_vv
+   use hybrid_ensemble_parameters, only: ntotensgrp,naensgrp,naensloc,ntlevs_ens,nsclgrp
+   use hybrid_ensemble_parameters, only: en_perts
    use gsi_io, only: verbose
 
    implicit none
@@ -4016,8 +4177,17 @@ subroutine hybens_localization_setup
    integer(i_kind),parameter   :: lunin = 47
    character(len=40),parameter :: fname = 'hybens_info'
    integer(i_kind) :: k,msig,istat,nz,kl
+   integer(i_kind) :: n,nk,m,ig,ic3,ic2,ipic,istatus,ii
+   integer(i_kind) :: ipc3d(nc3d),ipc2d(nc2d)
    logical         :: lexist,print_verbose
-   real(r_kind),allocatable:: s_ens_h_gu_x(:),s_ens_h_gu_y(:)
+   real(r_kind),allocatable:: s_ens_h_gu_x(:,:),s_ens_h_gu_y(:,:)
+   logical :: l_read_success
+   type(gsi_bundle) :: a_en(n_ens)
+   type(gsi_grid)  :: grid_ens
+   real(r_kind), pointer :: values(:) => NULL()
+   character(len=*),parameter::myname_=myname//'*hybens_localization_setup'
+
+   l_read_success=.false.
    print_verbose=.false. .and. mype == 0
    if(verbose .and. mype == 0)print_verbose=.true.
 
@@ -4044,10 +4214,25 @@ subroutine hybens_localization_setup
          endif
          if(mype==0) write(6,'(" LOCALIZATION, BETA_S, BETA_E VERTICAL PROFILES FOLLOW")')
          do k = 1,grd_ens%nsig
-            read(lunin,101) s_ens_hv(k), s_ens_vv(k), beta_s(k), beta_e(k)
-            if(mype==0) write(6,101) s_ens_hv(k), s_ens_vv(k), beta_s(k), beta_e(k)
+            read(lunin,101) s_ens_hv(k,1), s_ens_vv(k,1), beta_s(k), beta_e(k)
+            if(mype==0) write(6,101) s_ens_hv(k,1), s_ens_vv(k,1), beta_s(k), beta_e(k)
          enddo
+         do ig=2,naensloc
+            do k = 1,grd_ens%nsig
+               read(lunin,101,end=300) s_ens_hv(k,ig),s_ens_vv(k,ig)
+            enddo
+         enddo
+         l_read_success=.true.
          close(lunin)
+300      continue
+         if(.not.l_read_success) then
+            do ig=2,naensloc
+               do k = 1,grd_ens%nsig
+                  s_ens_hv(k,ig)=s_ens_hv(k,1)
+                  s_ens_vv(k,ig)=s_ens_vv(k,1)
+               enddo
+            enddo
+         endif
 
       else
 
@@ -4060,8 +4245,8 @@ subroutine hybens_localization_setup
          vvlocal = .true.
          nz = msig
          kl = grd_loc%kend_alloc-grd_loc%kbegin_loc+1
-         if(.not.allocated(s_ens_h_gu_x)) allocate(s_ens_h_gu_x(grd_loc%nsig*n_ens))
-         if(.not.allocated(s_ens_h_gu_y)) allocate(s_ens_h_gu_y(grd_loc%nsig*n_ens))
+         if(.not.allocated(s_ens_h_gu_x)) allocate(s_ens_h_gu_x(grd_loc%nsig*n_ens,naensloc))
+         if(.not.allocated(s_ens_h_gu_y)) allocate(s_ens_h_gu_y(grd_loc%nsig*n_ens,naensloc))
       endif
 
    endif ! if ( readin_localization .or. readin_beta )
@@ -4091,10 +4276,12 @@ subroutine hybens_localization_setup
    if ( .not. readin_localization ) then ! assign all levels to same value, s_ens_h, s_ens_v
       nz = 1
       kl = 1
-      if(.not.allocated(s_ens_h_gu_x)) allocate(s_ens_h_gu_x(1))
-      if(.not.allocated(s_ens_h_gu_y)) allocate(s_ens_h_gu_y(1))
-      s_ens_hv = s_ens_h
-      s_ens_vv = s_ens_v
+      if(.not.allocated(s_ens_h_gu_x)) allocate(s_ens_h_gu_x(1,naensloc))
+      if(.not.allocated(s_ens_h_gu_y)) allocate(s_ens_h_gu_y(1,naensloc))
+      do ig=1,naensloc
+         s_ens_hv(:,ig) = s_ens_h(ig)
+         s_ens_vv(:,ig) = s_ens_v(ig)
+      enddo
    endif
 
    ! Set up localization filters
@@ -4103,11 +4290,19 @@ subroutine hybens_localization_setup
    call normal_new_factorization_rf_z
 
    if ( regional ) then ! convert s_ens_h from km to grid units.
-      call convert_km_to_grid_units(s_ens_h_gu_x,s_ens_h_gu_y,nz)
       if ( vvlocal ) then
-         call init_rf_x(s_ens_h_gu_x(grd_loc%kbegin_loc:grd_loc%kend_alloc),kl)
-         call init_rf_y(s_ens_h_gu_y(grd_loc%kbegin_loc:grd_loc%kend_alloc),kl)
+         call convert_km_to_grid_units(s_ens_h_gu_x(1:nz,:),s_ens_h_gu_y(1:nz,:),nz)
+         do n=2,n_ens
+            nk=(n-1)*nz
+            do k=1,nz
+               s_ens_h_gu_x(nk+k,:)=s_ens_h_gu_x(k,:)
+               s_ens_h_gu_y(nk+k,:)=s_ens_h_gu_y(k,:)
+            enddo
+         enddo
+         call init_rf_x(s_ens_h_gu_x(grd_loc%kbegin_loc:grd_loc%kend_alloc,:),kl)
+         call init_rf_y(s_ens_h_gu_y(grd_loc%kbegin_loc:grd_loc%kend_alloc,:),kl)
       else
+         call convert_km_to_grid_units(s_ens_h_gu_x,s_ens_h_gu_y,nz)
          call init_rf_x(s_ens_h_gu_x,kl)
          call init_rf_y(s_ens_h_gu_y,kl)
       endif
@@ -4117,6 +4312,80 @@ subroutine hybens_localization_setup
       call init_sf_xy(jcap_ens)
    endif
 
+   if(ntotensgrp>1) then
+      call gsi_bundlegetpointer(en_perts(1,1,1),cvars3d,ipc3d,istatus)
+      if(istatus/=0) then
+         write(6,*) myname_,': cannot find 3d pointers'
+         call stop2(999)
+      endif
+      call gsi_bundlegetpointer(en_perts(1,1,1),cvars2d,ipc2d,istatus)
+      if(istatus/=0) then
+         write(6,*) myname_,': cannot find 2d pointers'
+         call stop2(999)
+      endif
+      if(nsclgrp>1) then
+         call gsi_gridcreate(grid_ens,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig)
+         allocate(values(grd_ens%latlon11*grd_ens%nsig*n_ens))
+         do ig=1,nsclgrp-1
+            ii=0
+            do n=1,n_ens
+               a_en(n)%values => values(ii+1:ii+grd_ens%latlon11*grd_ens%nsig)
+               call gsi_bundleset(a_en(n),grid_ens,'Ensemble Bundle',istatus,names3d=(/'a_en'/),bundle_kind=r_kind)
+               if (istatus/=0) then
+                  write(6,*) myname_,': error alloc(ensemble bundle)'
+                  call stop2(999)
+               endif
+               ii=ii+grd_ens%latlon11*grd_ens%nsig
+            enddo
+            do m=1,ntlevs_ens
+               do n=1,n_ens
+                  en_perts(n,ig+1,m)%valuesr4=en_perts(n,ig,m)%valuesr4
+               enddo
+               do ic3=1,nc3d
+                  ipic=ipc3d(ic3)
+                  do n=1,n_ens
+                     do k=1,grd_ens%nsig
+                        a_en(n)%r3(1)%q(:,:,k)=en_perts(n,ig,m)%r3(ipic)%qr4(:,:,k)
+                     enddo
+                  enddo
+                  call bkgcov_a_en_new_factorization(naensgrp+ig,a_en)
+                  do n=1,n_ens
+                     do k=1,grd_ens%nsig
+                        en_perts(n,ig,m)%r3(ipic)%qr4(:,:,k)=a_en(n)%r3(1)%q(:,:,k)
+                     enddo
+                  enddo
+               enddo
+               do ic2=1,nc2d
+                  ipic=ipc2d(ic2)
+                  do n=1,n_ens
+                     do k=1,grd_ens%nsig
+                        a_en(n)%r3(1)%q(:,:,k)=en_perts(n,ig,m)%r2(ipic)%qr4(:,:)
+                     enddo
+                  enddo
+                  call bkgcov_a_en_new_factorization(naensgrp+ig,a_en)
+                  do n=1,n_ens
+                     en_perts(n,ig,m)%r2(ipic)%qr4(:,:)=a_en(n)%r3(1)%q(:,:,1)
+                  enddo
+               enddo
+               do n=1,n_ens
+                  en_perts(n,ig+1,m)%valuesr4=en_perts(n,ig+1,m)%valuesr4-en_perts(n,ig,m)%valuesr4
+               enddo
+            enddo
+            do n=1,n_ens
+               call gsi_bundleunset(a_en(n),istatus)
+            enddo
+         enddo
+         deallocate(values)
+      endif
+      do ig=nsclgrp+1,ntotensgrp
+         do m=1,ntlevs_ens
+            do n=1,n_ens
+               en_perts(n,ig,m)%valuesr4=en_perts(n,ig-nsclgrp,m)%valuesr4
+            enddo
+         enddo
+      enddo
+   endif
+
    !!!!!!!! setup beta_s, beta_e!!!!!!!!!!!!
    ! vertical variation of static and ensemble weights
 
@@ -4138,12 +4407,14 @@ subroutine hybens_localization_setup
 
    ! write out final values for s_ens_hv, s_ens_vv, beta_s, beta_e
    if ( print_verbose ) then
-      write(6,*) 'HYBENS_LOCALIZATION_SETUP: s_ens_hv,s_ens_vv,beta_s,beta_e'
+      write(6,*) 'HYBENS_LOCALIZATION_SETUP: s_ens_hv(:,1),s_ens_vv(:,1),beta_s,beta_e'
       do k=1,grd_ens%nsig
-         write(6,101) s_ens_hv(k), s_ens_vv(k), beta_s(k), beta_e(k)
+         write(6,101) s_ens_hv(k,1), s_ens_vv(k,1), beta_s(k), beta_e(k)
       enddo
    endif
 
+   call setup_ensgrp2aensgrp
+
    return
 
 end subroutine hybens_localization_setup
@@ -4174,16 +4445,17 @@ subroutine convert_km_to_grid_units(s_ens_h_gu_x,s_ens_h_gu_y,nz)
 !$$$
 
   use kinds, only: r_kind,i_kind
-  use hybrid_ensemble_parameters, only: grd_loc,n_ens,s_ens_hv
+  use hybrid_ensemble_parameters, only: s_ens_hv
   use hybrid_ensemble_parameters, only: region_dx_ens,region_dy_ens
+  use hybrid_ensemble_parameters, only: naensloc
   use gsi_io, only: verbose
   implicit none
 
   integer(i_kind) ,intent(in   ) ::nz
-  real(r_kind),intent(  out) ::s_ens_h_gu_x(nz*n_ens),s_ens_h_gu_y(nz*n_ens)
+  real(r_kind),intent(  out) ::s_ens_h_gu_x(nz,naensloc),s_ens_h_gu_y(nz,naensloc)
   logical :: print_verbose
   real(r_kind) dxmax,dymax
-  integer(i_kind) k,n,nk
+  integer(i_kind) k
 
   print_verbose=.false.
   if(verbose) print_verbose=.true.
@@ -4197,25 +4469,14 @@ subroutine convert_km_to_grid_units(s_ens_h_gu_x,s_ens_h_gu_y,nz)
   end if
 
   do k=1,nz
-     s_ens_h_gu_x(k)=s_ens_hv(k)/(.001_r_kind*dxmax)
-     s_ens_h_gu_y(k)=s_ens_hv(k)/(.001_r_kind*dymax)
-     if(print_verbose) write(6,*)' in convert_km_to_grid_units,s_ens_h,s_ens_h_gu_x,y=', &
-                    s_ens_hv(k),s_ens_h_gu_x(k),s_ens_h_gu_y(k)
-
+     s_ens_h_gu_x(k,:)=s_ens_hv(k,:)/(.001_r_kind*dxmax)
+     s_ens_h_gu_y(k,:)=s_ens_hv(k,:)/(.001_r_kind*dymax)
+     if(print_verbose) write(6,*)' in convert_km_to_grid_units,s_ens_h(k,1),s_ens_h_gu_x,y(k,1)=', &
+                    s_ens_hv(k,1),s_ens_h_gu_x(k,1),s_ens_h_gu_y(k,1)
   enddo
 
-  if(nz>1)then
-  do n=2,n_ens
-     nk=(n-1)*grd_loc%nsig
-     do k=1,grd_loc%nsig
-        s_ens_h_gu_x(nk+k)=s_ens_h_gu_x(k)
-        s_ens_h_gu_y(nk+k)=s_ens_h_gu_y(k)
-     enddo
-  enddo
-  endif
   return
 
-
 end subroutine convert_km_to_grid_units
 
 subroutine grads1(f,nvert,mype,fname)
@@ -5146,4 +5407,116 @@ subroutine setup_pwgt
 
 end subroutine setup_pwgt
 
+subroutine setup_ensgrp2aensgrp
+!$$$  subprogram documentation block
+!                .      .    .                                       .
+! subprogram:    set a matrix of (naensgrp,naensgrp)
+!
+! program history log:
+!   2022-09-15  yokota  - add scale/variable/time-dependent localization
+!
+!   input argument list:
+!
+!   output argument list:
+!
+! remarks:
+!    need to reconcile grid in gsi_bundle w/ grid_ens/grid_anl
+!
+! attributes:
+!   language: f90
+!   machine:  ibm RS/6000 SP
+!
+!$$$ end documentation block
+  use constants, only: zero,one
+  use hybrid_ensemble_parameters, only: l_timloc_opt,i_ensloccov4tim,i_ensloccov4var,i_ensloccov4scl
+  use hybrid_ensemble_parameters, only: ensloccov4tim,ensloccov4var,ensloccov4scl
+  use hybrid_ensemble_parameters, only: ntotensgrp,naensgrp,ntlevs_ens,nsclgrp,ngvarloc
+  use hybrid_ensemble_parameters, only: ensgrp2aensgrp
+  use hybrid_ensemble_parameters, only: idaen2d,idaen3d
+  use hybrid_ensemble_parameters, only: alphacvarsclgrpmat
+  implicit none
+
+  integer (i_kind):: i,j
+  integer (i_kind):: ig,ibin,ic,iscl,itim1,itim2,igvar1,igvar2,iscl1,iscl2,ivargrp
+  integer (i_kind):: ntimloc,interval4aens
+
+  if (l_timloc_opt) then
+     ntimloc=ntlevs_ens
+     interval4aens=ntotensgrp
+  else
+     ntimloc=1
+     interval4aens=0
+  endif
+  if(naensgrp/=ntimloc*ngvarloc*nsclgrp) then
+     write(6,*)'setup_ensgrp2aensgrp: wrong naensgrp'
+     call stop2(666)
+  endif
+  if(ntotensgrp/=ngvarloc*nsclgrp) then
+     write(6,*)'setup_ensgrp2aensgrp: wrong ntotensgrp'
+     call stop2(666)
+  endif
+  ensgrp2aensgrp=-999
+  do ibin=1,ntlevs_ens
+     do ic=1,nc3d+nc2d
+        if(ngvarloc>1) then
+           if(ic<=nc3d) ivargrp=idaen3d(ic)
+           if(ic> nc3d) ivargrp=idaen2d(ic-nc3d)
+        else
+           ivargrp=1
+        endif
+        do iscl=1,nsclgrp
+           ig=(ivargrp-1)*nsclgrp+iscl
+           ensgrp2aensgrp(ig,ic,ibin)=(ibin-1)*interval4aens+(ivargrp-1)*nsclgrp+iscl
+        enddo
+     enddo
+  enddo
+
+  if (i_ensloccov4tim==0) then
+     ensloccov4tim=one
+  elseif (i_ensloccov4tim==1)then
+     ensloccov4tim=zero
+     ensloccov4tim(1)=one
+  else
+     write(6,*)'setup_ensgrp2aensgrp: wrong i_ensloccov4tim'
+     call stop2(666)
+  endif
+  if (i_ensloccov4var==0) then
+     ensloccov4var=one
+  elseif (i_ensloccov4var==1)then
+     ensloccov4var=zero
+     ensloccov4var(1)=one
+  else
+     write(6,*)'setup_ensgrp2aensgrp: wrong i_ensloccov4var'
+     call stop2(666)
+  endif
+  if (i_ensloccov4scl==0) then
+     ensloccov4scl=one
+  elseif (i_ensloccov4scl==1)then
+     ensloccov4scl=zero
+     ensloccov4scl(1)=one
+  else
+     write(6,*)'setup_ensgrp2aensgrp: wrong i_ensloccov4scl'
+     call stop2(666)
+  endif
+
+  do itim2=1,ntimloc
+     do itim1=1,ntimloc
+        do igvar1=1,ngvarloc
+           do igvar2=1,ngvarloc
+              do iscl1=1,nsclgrp
+                 do iscl2=1,nsclgrp
+                    i=(itim1-1)*interval4aens+(igvar1-1)*nsclgrp+iscl1
+                    j=(itim2-1)*interval4aens+(igvar2-1)*nsclgrp+iscl2
+                    alphacvarsclgrpmat(i,j)=ensloccov4tim(abs(itim1-itim2)+1) &
+                                           *ensloccov4var(abs(igvar1-igvar2)+1) &
+                                           *ensloccov4scl(abs(iscl1-iscl2)+1) !first for ttime covariance
+                 enddo
+              enddo
+           enddo
+        enddo
+     enddo
+  enddo
+
+end subroutine setup_ensgrp2aensgrp
+
 end module hybrid_ensemble_isotropic
diff --git a/src/gsi/hybrid_ensemble_parameters.f90 b/src/gsi/hybrid_ensemble_parameters.f90
index 9ca7770c41..7b1c963764 100644
--- a/src/gsi/hybrid_ensemble_parameters.f90
+++ b/src/gsi/hybrid_ensemble_parameters.f90
@@ -90,8 +90,19 @@ module hybrid_ensemble_parameters
 !      beta_e0 - default weight given to ensemble background error covariance
 !                (if .not. readin_beta). if beta_e0<0, then it is set to
 !                1.-beta_s0 (this is the default)
-!      s_ens_h:    horizontal localization correlation length (units of km), default = 2828.0
-!      s_ens_v:    vertical localization correlation length (grid units), default = 30.0
+!      s_ens_h: horizontal localization correlation length of Gaussian exp(-0.5*(r/L)**2)
+!               (units of km), default = 2828.0
+!      s_ens_v: vertical localization correlation length of Gaussian exp(-0.5*(r/L)**2)
+!               (grid units if s_ens_v>=0, or units of ln(p) if s_ens_v<0), default = 30.0
+!                  in scale/variable/time-dependent localization (SDL/VDL/TDL),
+!                  localization length for i-th scale, j-th variable, and k-th time is
+!                     s_ens_[hv]( i + nsclgrp*(j-1) + nsclgrp*ngvarloc*(k-1) )
+!                        in SDL(nsclgrp>1),         i = 1(largest scale)  .. nsclgrp(smallest scale)
+!                        in VDL(ngvarloc=2),        j = 1(itracer<=10)    .. 2(itracer>=11)
+!                        in TDL(l_timloc_opt=true), k = 1(first time bin) .. ntlevs_ens(last time bin)
+!                  in SDL, scale separation length for i-th scale is also set here as
+!                     s_ens_[hv]( naensgrp+i ) - naensgrp is the total number of localization lengths for SDL/VDL/TDL
+!                        in applying SDL only horizontally, set s_ens_v(naensgrp+i)=0.0
 !      generate_ens:  if .true., generate ensemble perturbations internally as random samples of background B.
 !                       (used primarily for testing/debugging)
 !                     if .false., read external ensemble perturbations
@@ -118,6 +129,19 @@ module hybrid_ensemble_parameters
 !      ensemble_path: path to ensemble members; default './'
 !      ens_fast_read: read ensemble in parallel; default '.false.'
 !      sst_staticB:   if .true. (default) uses only static part of B error covariance for SST
+!      nsclgrp:         number of scale-dependent localization lengths
+!      l_timloc_opt:    if true, then turn on time-dependent localization
+!      ngvarloc:        number of variable-dependent localization lengths
+!      naensloc:        total number of spatial localization lengths and scale separation lengths (should be naensgrp+nsclgrp-1)
+!      i_ensloccov4tim: flag of cross-temporal localization
+!                         =0: cross-temporal covariance is retained
+!                         =1: cross-temporal covariance is zero
+!      i_ensloccov4var: flag of cross-variable localization
+!                         =0: cross-variable covariance is retained
+!                         =1: cross-variable covariance is zero
+!      i_ensloccov4scl: flag of cross-scale localization
+!                         =0: cross-scale covariance is retained
+!                         =1: cross-scale covariance is zero
 !=====================================================================================================
 !
 !
@@ -151,6 +175,7 @@ module hybrid_ensemble_parameters
 !   2015-01-22  Hu      - add flag i_en_perts_io to control reading ensemble perturbation.
 !   2015-02-11  Hu      - add flag l_ens_in_diff_time to force GSI hybrid use ensembles not available at analysis time
 !   2015-09-18  todling - add sst_staticB to control use of ensemble SST error covariance 
+!   2022-09-15  yokota  - add scale/variable/time-dependent localization
 !
 ! subroutines included:
 
@@ -174,8 +199,10 @@ module hybrid_ensemble_parameters
 !                                   beta_s(:) = beta_s0     , vertically varying weights given to B ; 
 !                                   beta_e(:) = 1 - beta_s0 , vertically varying weights given to A.
 !                            If (readin_beta) then beta_s and beta_e are read from a file and beta_s0 is not used.
-!   def s_ens_h             - homogeneous isotropic horizontal ensemble localization scale (km)
-!   def s_ens_v             - vertical localization scale (grid units for now)
+!   def s_ens_h             - horizontal localization correlation length of Gaussian exp(-0.5*(r/L)**2)
+!                             (units of km), default = 2828.0
+!   def s_ens_v             - vertical localization correlation length of Gaussian exp(-0.5*(r/L)**2)
+!                             (grid units if s_ens_v>=0, or units of ln(p) if s_ens_v<0), default = 30.0
 !   def readin_localization - flag to read (.true.)external localization information file
 !   def eqspace_ensgrid     - if .true., then ensemble grid is equal spaced, staggered 1/2 grid unit off
 !                               poles.  if .false., then Gaussian grid assumed for ensemble (global only)
@@ -292,12 +319,20 @@ module hybrid_ensemble_parameters
   public :: en_perts,ps_bar
   public :: region_lat_ens,region_lon_ens
   public :: region_dx_ens,region_dy_ens
+  public :: naensgrp,ntotensgrp,nsclgrp,naensloc,ngvarloc
+  public :: ensgrp2aensgrp
+  public :: ensloccov4tim,ensloccov4var,ensloccov4scl
+  public :: alphacvarsclgrpmat
+  public :: l_timloc_opt
+  public :: i_ensloccov4tim,i_ensloccov4var,i_ensloccov4scl
+  public :: idaen3d,idaen2d
   public :: ens_fast_read
   public :: l_both_fv3sar_gfs_ens 
   public :: sst_staticB
   public :: limqens
 
   logical l_hyb_ens,uv_hyb_ens,q_hyb_ens,oz_univ_static,sst_staticB
+  logical l_timloc_opt
   logical aniso_a_en
   logical full_ensemble,pwgtflg
   logical generate_ens
@@ -317,11 +352,16 @@ module hybrid_ensemble_parameters
   integer(i_kind) i_en_perts_io
   integer(i_kind) n_ens,nlon_ens,nlat_ens,jcap_ens,jcap_ens_test
   integer(i_kind) n_ens_gfs,n_ens_fv3sar
-  real(r_kind) beta_s0,beta_e0,s_ens_h,s_ens_v,grid_ratio_ens
+  real(r_kind) beta_s0,beta_e0,grid_ratio_ens
+  integer(i_kind),parameter::max_naensloc=20
+  integer(i_kind),parameter::max_nvars=100
+  real(r_kind) s_ens_h(max_naensloc)
+  real(r_kind) s_ens_v(max_naensloc)
   type(sub2grid_info),save :: grd_ens,grd_loc,grd_sploc,grd_anl,grd_e1,grd_a1
   type(spec_vars),save :: sp_ens,sp_loc
   type(egrid2agrid_parm),save :: p_e2a,p_sploc2ens
-  real(r_kind),allocatable,dimension(:) :: s_ens_hv,s_ens_vv
+  real(r_kind),allocatable,dimension(:,:) :: s_ens_vv
+  real(r_kind),allocatable,dimension(:,:) :: s_ens_hv
   real(r_kind),allocatable,dimension(:) :: sqrt_beta_s,sqrt_beta_e
   real(r_kind),allocatable,dimension(:) :: beta_s,beta_e
   real(r_kind),allocatable,dimension(:,:,:) :: pwgt
@@ -336,6 +376,18 @@ module hybrid_ensemble_parameters
   integer(i_kind) regional_ensemble_option
   integer(i_kind) fv3sar_ensemble_opt 
   character(len=512),save :: ensemble_path
+  real(r_kind),allocatable,dimension(:,:) :: alphacvarsclgrpmat
+  integer(i_kind),allocatable,dimension(:,:,:) :: ensgrp2aensgrp
+  real(r_kind),allocatable,dimension(:) :: ensloccov4tim,ensloccov4var,ensloccov4scl
+  integer(i_kind) :: nsclgrp=1
+  integer(i_kind) :: naensgrp=1
+  integer(i_kind) :: ntotensgrp=1
+  integer(i_kind) :: naensloc=1
+  integer(i_kind) :: ngvarloc=1
+  integer(i_kind) :: i_ensloccov4tim=0
+  integer(i_kind) :: i_ensloccov4var=0
+  integer(i_kind) :: i_ensloccov4scl=0
+  integer(i_kind),allocatable,dimension(:) :: idaen3d,idaen2d
 
 ! following is for storage of ensemble perturbations:
 
@@ -343,7 +395,7 @@ module hybrid_ensemble_parameters
 !   def nelen               - length of one ensemble perturbation vector
 
   integer(i_kind) nelen
-  type(gsi_bundle),save,allocatable :: en_perts(:,:)
+  type(gsi_bundle),save,allocatable :: en_perts(:,:,:)
   real(r_single),dimension(:,:,:),allocatable:: ps_bar
   real(r_single):: limqens
 
@@ -384,6 +436,7 @@ subroutine init_hybrid_ensemble_parameters
   implicit none
 
   l_hyb_ens=.false.
+  l_timloc_opt=.false.
   full_ensemble=.false.
   pwgtflg=.false.
   uv_hyb_ens=.false.
@@ -435,15 +488,23 @@ subroutine create_hybens_localization_parameters
   use constants, only: zero
   implicit none
   
-  allocate( s_ens_hv(grd_ens%nsig),s_ens_vv(grd_ens%nsig) )
+  allocate( s_ens_hv(grd_ens%nsig,naensloc),s_ens_vv(grd_ens%nsig,naensloc) )
   allocate( beta_s(grd_ens%nsig),beta_e(grd_ens%nsig))
   allocate( sqrt_beta_s(grd_ens%nsig),sqrt_beta_e(grd_ens%nsig) )
   allocate( pwgt(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig) )
+  allocate( alphacvarsclgrpmat(naensgrp,naensgrp) )
+  allocate( ensgrp2aensgrp(ntotensgrp,max_nvars,ntlevs_ens) )
+  allocate( ensloccov4tim(ntlevs_ens),ensloccov4var(ngvarloc),ensloccov4scl(nsclgrp) )
   beta_s  =one
   beta_e  =zero
   sqrt_beta_s=one
   sqrt_beta_e=zero
   pwgt=zero
+  alphacvarsclgrpmat=one
+  ensgrp2aensgrp=1
+  ensloccov4tim=one
+  ensloccov4var=one
+  ensloccov4scl=one
   
 end subroutine create_hybens_localization_parameters
 
@@ -453,6 +514,8 @@ subroutine destroy_hybens_localization_parameters
   deallocate(s_ens_vv,s_ens_hv) 
   deallocate(beta_s,beta_e)
   deallocate(sqrt_beta_s,sqrt_beta_e,pwgt)
+  deallocate(alphacvarsclgrpmat)
+  deallocate(ensgrp2aensgrp)
 
 end subroutine destroy_hybens_localization_parameters
 
diff --git a/src/gsi/jfunc.f90 b/src/gsi/jfunc.f90
index d65b4c8fd3..1b92ad2e94 100644
--- a/src/gsi/jfunc.f90
+++ b/src/gsi/jfunc.f90
@@ -738,6 +738,7 @@ subroutine set_pointer
     use gsi_4dvar, only: nsubwin, lsqrtb
     use bias_predictors, only: setup_predictors
     use hybrid_ensemble_parameters, only: l_hyb_ens,n_ens,generate_ens,grd_ens,nval_lenz_en
+    use hybrid_ensemble_parameters, only: naensgrp
     implicit none
 
     integer(i_kind) n_ensz,nval_lenz_tot,nval_lenz_enz
@@ -748,8 +749,8 @@ subroutine set_pointer
     nval_levs=max(0,nc3d)*nsig+max(0,nc2d)
     nval_len=nval_levs*latlon11
     if(l_hyb_ens) then
-       nval_len=nval_len+n_ens*nsig*grd_ens%latlon11
-       nval_levs_ens=nval_levs+n_ens*nsig
+       nval_len=nval_len+naensgrp*n_ens*nsig*grd_ens%latlon11
+       nval_levs_ens=nval_levs+naensgrp*n_ens*nsig
     end if
     nsclen=npred*jpch_rad
     npclen=npredp*npcptype
diff --git a/src/gsi/stub_get_pseudo_ensperts.f90 b/src/gsi/stub_get_pseudo_ensperts.f90
index 7bababc0c3..ca39b4b6e2 100644
--- a/src/gsi/stub_get_pseudo_ensperts.f90
+++ b/src/gsi/stub_get_pseudo_ensperts.f90
@@ -32,7 +32,7 @@ subroutine get_pseudo_ensperts_dummy(this,en_perts,nelen)
     use kinds, only: i_kind
     implicit none
     class(get_pseudo_ensperts_class), intent(inout) :: this
-    type(gsi_bundle),allocatable, intent(in   ) :: en_perts(:,:)
+    type(gsi_bundle),allocatable, intent(in   ) :: en_perts(:,:,:)
     integer(i_kind),                   intent(in   ) :: nelen
   
     write(6,*)'get_pseudo_ensperts:  ***WARNING*** dummy call ... does nothing!'
diff --git a/src/gsi/stub_get_wrf_mass_ensperts.f90 b/src/gsi/stub_get_wrf_mass_ensperts.f90
index 3089cb9336..bf98f3106e 100644
--- a/src/gsi/stub_get_wrf_mass_ensperts.f90
+++ b/src/gsi/stub_get_wrf_mass_ensperts.f90
@@ -38,7 +38,7 @@ subroutine ens_spread_dualres_regional_dummy(this,mype,en_perts,nelen,en_bar)
     implicit none
     class(get_wrf_mass_ensperts_class), intent(inout) :: this
     integer(i_kind),intent(in):: mype
-    type(gsi_bundle),allocatable, intent(in   ) :: en_perts(:,:)
+    type(gsi_bundle),allocatable, intent(in   ) :: en_perts(:,:,:)
     integer(i_kind), intent(in   ):: nelen
     type(gsi_bundle),OPTIONAL,intent(in):: en_bar
     write(6,*)'get_wrf_mass_ensperts:  ***WARNING*** dummy call ... does nothing!'
@@ -50,7 +50,7 @@ subroutine get_wrf_mass_ensperts_dummy(this,en_perts,nelen,ps_bar)
     use gsi_bundlemod, only: gsi_bundle
     implicit none
     class(get_wrf_mass_ensperts_class), intent(inout) :: this
-    type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:)
+    type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:,:)
     integer(i_kind), intent(in   ):: nelen
     real(r_single),dimension(:,:,:),allocatable:: ps_bar
     write(6,*)'get_wrf_mass_ensperts:  ***WARNING*** dummy call ... does nothing!'
diff --git a/src/gsi/stub_get_wrf_nmm_ensperts.f90 b/src/gsi/stub_get_wrf_nmm_ensperts.f90
index 14ab869ffd..db27125d70 100644
--- a/src/gsi/stub_get_wrf_nmm_ensperts.f90
+++ b/src/gsi/stub_get_wrf_nmm_ensperts.f90
@@ -13,7 +13,7 @@ subroutine get_wrf_nmm_ensperts_dummy(this,en_perts,nelen,region_lat_ens,region_
     
     implicit none
     class(get_wrf_nmm_ensperts_class), intent(inout) :: this
-    type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:)
+    type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:,:)
     integer(i_kind), intent(in   ):: nelen
     real(r_kind),allocatable, intent(inout):: region_lat_ens(:,:),region_lon_ens(:,:)
     real(r_single),dimension(:,:,:),allocatable, intent(inout):: ps_bar