module get_fv3_regional_ensperts_mod
use abstract_get_fv3_regional_ensperts_mod,only: abstract_get_fv3_regional_ensperts_class
  use kinds, only : i_kind
  use general_sub2grid_mod, only: sub2grid_info
  use constants, only:max_varname_length
  type, extends(abstract_get_fv3_regional_ensperts_class) :: get_fv3_regional_ensperts_class
  contains
    procedure, pass(this) :: get_fv3_regional_ensperts => get_fv3_regional_ensperts_run
    procedure, pass(this) :: ens_spread_dualres_regional => ens_spread_dualres_regional_fv3_regional
    procedure, pass(this) :: general_read_fv3_regional
  end type get_fv3_regional_ensperts_class
    type(sub2grid_info):: grd_fv3lam_ens_dynvar_io_nouv,grd_fv3lam_ens_tracer_io_nouv,grd_fv3lam_ens_uv
    character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_ens_io_dynmetvars3d_nouv 
                               ! copy of cvars3d excluding uv 3-d fields   
    character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_ens_io_tracermetvars3d_nouv
                               ! copy of cvars3d excluding uv 3-d fields   
    character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_ens_io_dynmetvars2d_nouv
    character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_ens_io_tracermetvars2d_nouv
contains
  
  subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar)
  !$$$  subprogram documentation block
  !                .      .    .                                       .
  ! subprogram:    get_fv3_regional_ensperts  read arw model ensemble members
  !   prgmmr: Ting            org: EMC/NCEP            date: 2018-12-13
  !
  ! abstract: read ensemble members from the fv3 regional (fv3_SAR)
  ! model,following Wanshu's programs to read those background files 
  !
  !
  ! program history log:
  !   2011-08-31  todling - revisit en_perts (single-prec) in light of extended bundle
  !   2019-04-22  CAPS(C. Tong) - add direct reflectivity DA option
  !
  !   2021-08-10  lei     - modify for fv3-lam ensemble spread output
  !   2021-11-01  lei     - modify for fv3-lam parallel IO
  !   input argument list:
  !
  !   output argument list:
  !
  ! attributes:
  !   language: f90
  !   machine: 
  !
  !$$$ end documentation block
  
     use kinds, only: r_kind,i_kind,r_single
     use constants, only: zero,one,half,zero_single,rd_over_cp,one_tenth
     use mpimod, only: mpi_comm_world,ierror,mype
     use hybrid_ensemble_parameters, only: n_ens,grd_ens
     use hybrid_ensemble_parameters, only: ntlevs_ens,ensemble_path
     use control_vectors, only: cvars2d,cvars3d,nc2d,nc3d
     use gsi_bundlemod, only: gsi_bundlecreate
     use gsi_bundlemod, only: gsi_grid
     use gsi_bundlemod, only: gsi_bundle
     use gsi_bundlemod, only: gsi_bundlegetpointer
     use gsi_bundlemod, only: gsi_bundledestroy
     use gsi_bundlemod, only: gsi_gridcreate
     use gsi_4dvar, only: ens_fhrlevs
     use gsi_rfv3io_mod, only: type_fv3regfilenameg
     use hybrid_ensemble_parameters, only: write_ens_sprd
     use directDA_radaruse_mod, only: l_use_cvpqx, cvpqx_pval, cld_nt_updt
     use directDA_radaruse_mod, only: l_use_dbz_directDA
     use directDA_radaruse_mod, only: l_cvpnr
     use general_sub2grid_mod, only: sub2grid_info,general_sub2grid_create_info
     use gridmod,only: regional
     use gsi_rfv3io_mod, only: fv3lam_io_dynmetvars3d_nouv,fv3lam_io_tracermetvars3d_nouv 
     use gsi_rfv3io_mod, only: fv3lam_io_dynmetvars2d_nouv,fv3lam_io_tracermetvars2d_nouv 
    

     implicit none
     class(get_fv3_regional_ensperts_class), intent(inout) :: this
     type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:)
     integer(i_kind), intent(in   ):: nelen
     real(r_single),dimension(:,:,:),allocatable,intent(inout):: ps_bar
 
     real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig):: u,v,tv,oz,rh
     real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2):: ps
     real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig)::w,ql,qi,qr,qg,qs,qnr
 
     real(r_single),pointer,dimension(:,:,:):: w3 =>NULL()
     real(r_single),pointer,dimension(:,:):: w2 =>NULL()
     real(r_kind),pointer,dimension(:,:,:):: x3 =>NULL()
     real(r_kind),pointer,dimension(:,:):: x2 =>NULL()
     type(gsi_bundle),allocatable,dimension(:):: en_bar
     type(gsi_grid):: grid_ens
     real(r_kind):: bar_norm,sig_norm,kapr,kap1

     character(len=64),dimension(:,:),allocatable:: names
     character(len=64),dimension(:,:),allocatable:: uvnames
     integer(i_kind),dimension(:,:),allocatable:: lnames
     integer(i_kind),dimension(:,:),allocatable:: uvlnames
 
     integer(i_kind):: i,j,k,n,mm1,istatus
     integer(i_kind):: ndynvario2d,ntracerio2d
     integer(r_kind):: ndynvario3d,ntracerio3d
     integer(i_kind):: inner_vars,numfields
     integer(i_kind):: ilev,ic2,ic3
     integer(i_kind):: m

     
     character(255) ensfilenam_str
     type(type_fv3regfilenameg)::fv3_filename 

     call gsi_gridcreate(grid_ens,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig)
     ! Allocate bundle to hold mean of ensemble members
     allocate(en_bar(ntlevs_ens))

!clt setup varnames for IO
     ndynvario2d=0
     ntracerio2d=0
     ndynvario3d=size(fv3lam_io_dynmetvars3d_nouv)
     ntracerio3d=size(fv3lam_io_tracermetvars3d_nouv)
     if (allocated(fv3lam_io_dynmetvars2d_nouv))then
       ndynvario2d=size(fv3lam_io_dynmetvars2d_nouv)
     endif
     if (allocated(fv3lam_io_tracermetvars2d_nouv)) then
       ntracerio2d=size(fv3lam_io_tracermetvars2d_nouv)
     endif
     allocate(fv3lam_ens_io_dynmetvars3d_nouv(ndynvario3d))
     allocate(fv3lam_ens_io_tracermetvars3d_nouv(ndynvario3d))
     fv3lam_ens_io_dynmetvars3d_nouv=fv3lam_io_dynmetvars3d_nouv
     fv3lam_ens_io_tracermetvars3d_nouv=fv3lam_io_tracermetvars3d_nouv
     if (ndynvario2d > 0 ) then
       allocate(fv3lam_ens_io_dynmetvars2d_nouv(ndynvario2d))
       fv3lam_ens_io_dynmetvars2d_nouv=fv3lam_io_dynmetvars2d_nouv
     endif
     if (ntracerio2d > 0 ) then
       allocate(fv3lam_ens_io_tracermetvars2d_nouv(ntracerio2d))
       fv3lam_ens_io_tracermetvars2d_nouv =fv3lam_io_tracermetvars3d_nouv 
     endif


     inner_vars=1
     numfields=inner_vars*(ndynvario3d*grd_ens%nsig+ndynvario2d)  
     allocate(lnames(1,numfields),names(1,numfields))
     ilev=1
     do i=1,ndynvario3d
       do k=1,grd_ens%nsig
         lnames(1,ilev)=k
         names(1,ilev)=fv3lam_ens_io_dynmetvars3d_nouv(i)
         ilev=ilev+1
       enddo
     enddo
     do i=1,ndynvario2d
       lnames(1,ilev)=1
       names(1,ilev)=fv3lam_ens_io_dynmetvars2d_nouv(i)
       ilev=ilev+1
     enddo
   

     call general_sub2grid_create_info(grd_fv3lam_ens_dynvar_io_nouv,inner_vars,grd_ens%nlat,&
          grd_ens%nlon,grd_ens%nsig,numfields,regional,names=names,lnames=lnames)

     inner_vars=1
     numfields=inner_vars*(ntracerio3d*grd_ens%nsig+ntracerio2d)  
     deallocate(lnames,names)
     allocate(lnames(1,numfields),names(1,numfields))
     ilev=1
     do i=1,ntracerio3d
      do k=1,grd_ens%nsig
        lnames(1,ilev)=k
        names(1,ilev)=fv3lam_ens_io_tracermetvars3d_nouv(i)
        ilev=ilev+1
      enddo
    enddo
    do i=1,ntracerio2d
       lnames(1,ilev)=1
       names(1,ilev)=fv3lam_ens_io_tracermetvars2d_nouv(i)
       ilev=ilev+1
    enddo
   

    call general_sub2grid_create_info(grd_fv3lam_ens_tracer_io_nouv,inner_vars,grd_ens%nlat,&
         grd_ens%nlon,grd_ens%nsig,numfields,regional,names=names,lnames=lnames)





    numfields=grd_ens%nsig
    inner_vars=2
    allocate(uvlnames(inner_vars,numfields),uvnames(inner_vars,numfields))
    ilev=1
    do k=1,grd_ens%nsig
      uvlnames(1,ilev)=k
      uvlnames(2,ilev)=k
      uvnames(1,ilev)='u'
      uvnames(2,ilev)='v'
      ilev=ilev+1
    enddo
    call general_sub2grid_create_info(grd_fv3lam_ens_uv,inner_vars,grd_ens%nlat,&
         grd_ens%nlon,grd_ens%nsig,numfields,regional,names=uvnames,lnames=uvlnames)

    deallocate(lnames,names,uvlnames,uvnames)


    do m=1,ntlevs_ens
       call gsi_bundlecreate(en_bar(m),grid_ens,'ensemble',istatus,names2d=cvars2d,names3d=cvars3d,bundle_kind=r_kind)
       if(istatus/=0) then
          write(6,*)' get_fv3_regional_ensperts_netcdf: trouble creating en_bar bundle'
          call stop2(9991)
       endif
    enddo ! for m 

     ! print info message for dirZDA
    if(mype==0)then
       if (l_use_cvpqx) then
         if ( cvpqx_pval == 0._r_kind ) then        ! CVlogq
             write(6,*) 'general_read_fv3_regional_dirZDA: convert qr/qs/qg to log transform.'
         else if ( cvpqx_pval > 0._r_kind ) then   ! CVpq
             write(6,*) 'general_read_fv3_regional_dirZDA: ',     &
                        'reset minimum of qr/qs/qg to specified values before power transform.' 
             write(6,*) 'general_read_fv3_regional_dirZDA: convert qr/qs/qg with power transform.'
         end if
       end if
       if ( cld_nt_updt > 0 .and. l_cvpnr) then
         write(6,*) 'general_read_fv3_regional_dirZDA: ',     &
                    'reset minimum of qnr to a specified value before power transform'
         write(6,*) 'general_read_fv3_regional_dirZDA: convert qnr with power transform .'
       end if
    end if


    do m=1,ntlevs_ens



 !
 ! INITIALIZE ENSEMBLE MEAN ACCUMULATORS
        en_bar(m)%values=zero
 
        do n=1,n_ens
           en_perts(n,m)%valuesr4 = zero
        enddo
 
        mm1=mype+1
        kap1=rd_over_cp+one
        kapr=one/rd_over_cp
 !
 ! LOOP OVER ENSEMBLE MEMBERS 
        do n=1,n_ens
           write(ensfilenam_str,22) trim(adjustl(ensemble_path)),ens_fhrlevs(m),n
22  format(a,'fv3SAR',i2.2,'_ens_mem',i3.3)
 ! DEFINE INPUT FILE NAME
           fv3_filename%grid_spec=trim(ensfilenam_str)//'-fv3_grid_spec' !exmaple thinktobe
           fv3_filename%ak_bk=trim(ensfilenam_str)//'-fv3_akbk'
           fv3_filename%dynvars=trim(ensfilenam_str)//'-fv3_dynvars'
           fv3_filename%tracers=trim(ensfilenam_str)//"-fv3_tracer"
           fv3_filename%sfcdata=trim(ensfilenam_str)//"-fv3_sfcdata"
           fv3_filename%couplerres=trim(ensfilenam_str)//"-coupler.res"
 ! 
 ! READ ENEMBLE MEMBERS DATA
           if (mype == 0) write(6,'(a,a)') &
               'CALL READ_FV3_REGIONAL_ENSPERTS FOR ENS DATA with the filename str : ',trim(ensfilenam_str)
           if (.not.l_use_dbz_directDA ) then ! Read additional hydrometers and w for dirZDA
              call this%general_read_fv3_regional(fv3_filename,ps,u,v,tv,rh,oz) 
           else
              call this%general_read_fv3_regional(fv3_filename,ps,u,v,tv,rh,oz,   &
                                                  g_ql=ql,g_qi=qi,g_qr=qr,g_qs=qs,g_qg=qg,g_qnr=qnr,g_w=w)
           end if
 
 ! SAVE ENSEMBLE MEMBER DATA IN COLUMN VECTOR
           do ic3=1,nc3d
 
              call gsi_bundlegetpointer(en_perts(n,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)
              end if
              call gsi_bundlegetpointer(en_bar(m),trim(cvars3d(ic3)),x3,istatus)
              if(istatus/=0) then
                 write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for en_bar'
                 call stop2(9993)
              end if

              select case (trim(cvars3d(ic3)))
 
                 case('sf','SF')
    
                    do k=1,grd_ens%nsig
                       do i=1,grd_ens%lon2
                          do j=1,grd_ens%lat2
                             w3(j,i,k) = u(j,i,k)
                             x3(j,i,k)=x3(j,i,k)+u(j,i,k)
                          end do
                       end do
                    end do
 
                 case('vp','VP')
 
                    do k=1,grd_ens%nsig
                       do i=1,grd_ens%lon2
                          do j=1,grd_ens%lat2
                             w3(j,i,k) = v(j,i,k)
                             x3(j,i,k)=x3(j,i,k)+v(j,i,k)
                          end do
                       end do
                    end do
 
                 case('t','T')
 
                    do k=1,grd_ens%nsig
                       do i=1,grd_ens%lon2
                          do j=1,grd_ens%lat2
                             w3(j,i,k) = tv(j,i,k)
                             x3(j,i,k)=x3(j,i,k)+tv(j,i,k)
                          end do
                       end do
                    end do
 
                 case('q','Q')
 
                    do k=1,grd_ens%nsig
                       do i=1,grd_ens%lon2
                          do j=1,grd_ens%lat2
                             w3(j,i,k) = rh(j,i,k)
                             x3(j,i,k)=x3(j,i,k)+rh(j,i,k)
                          end do
                       end do
                    end do
 
                 case('oz','OZ')
 
                    do k=1,grd_ens%nsig
                       do i=1,grd_ens%lon2
                          do j=1,grd_ens%lat2
                             w3(j,i,k) = oz(j,i,k)
                             x3(j,i,k)=x3(j,i,k)+oz(j,i,k)
                          end do
                       end do
                    end do
 
! save additional ensemble varaible data for direct reflectivity DA

                 case('ql','QL')

                    do k=1,grd_ens%nsig
                       do i=1,grd_ens%lon2
                          do j=1,grd_ens%lat2
                             w3(j,i,k) = ql(j,i,k)
                             x3(j,i,k)=x3(j,i,k)+ql(j,i,k)
                          end do
                       end do
                    end do

                 case('qi','QI')

                    do k=1,grd_ens%nsig
                       do i=1,grd_ens%lon2
                          do j=1,grd_ens%lat2
                             w3(j,i,k) = qi(j,i,k)
                             x3(j,i,k)=x3(j,i,k)+qi(j,i,k)
                          end do
                       end do
                    end do

                 case('qr','QR')

                    do k=1,grd_ens%nsig
                       do i=1,grd_ens%lon2
                          do j=1,grd_ens%lat2
                             w3(j,i,k) = qr(j,i,k)
                             x3(j,i,k)=x3(j,i,k)+qr(j,i,k)
                          end do
                       end do
                    end do

                 case('qs','QS')

                    do k=1,grd_ens%nsig
                       do i=1,grd_ens%lon2
                          do j=1,grd_ens%lat2
                             w3(j,i,k) = qs(j,i,k)
                             x3(j,i,k)=x3(j,i,k)+qs(j,i,k)
                          end do
                       end do
                    end do

                 case('qg','QG')

                    do k=1,grd_ens%nsig
                       do i=1,grd_ens%lon2
                          do j=1,grd_ens%lat2
                             w3(j,i,k) = qg(j,i,k)
                             x3(j,i,k)=x3(j,i,k)+qg(j,i,k)
                          end do
                       end do
                    end do

                 case('qnr','QNR')

                      do k=1,grd_ens%nsig
                         do i=1,grd_ens%lon2
                            do j=1,grd_ens%lat2
                               if ( l_use_dbz_directDA ) then ! direct reflectivity DA
                                  if ( cld_nt_updt > 0 ) then ! Update Nc 
                                     w3(j,i,k) = qnr(j,i,k)
                                     x3(j,i,k)=x3(j,i,k)+qnr(j,i,k)
                                  end if
                               else     ! .not. l_use_dbz_directDA
                                  w3(j,i,k) = qnr(j,i,k)
                                  x3(j,i,k)=x3(j,i,k)+qnr(j,i,k)
                               end if
                            end do
                         end do
                      end do

                 case('w','W')
                    do k=1,grd_ens%nsig
                       do i=1,grd_ens%lon2
                          do j=1,grd_ens%lat2
                             w3(j,i,k) = w(j,i,k)
                             x3(j,i,k)=x3(j,i,k)+w(j,i,k)
                          end do
                       end do
                    end do
              end select

       

           end do
 
           do ic2=1,nc2d
    
              call gsi_bundlegetpointer(en_perts(n,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)
              end if
              call gsi_bundlegetpointer(en_bar(m),trim(cvars2d(ic2)),x2,istatus)
              if(istatus/=0) then
                 write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for en_bar'
                 call stop2(9995)
              end if
 
              select case (trim(cvars2d(ic2)))
 
                 case('ps','PS')
 
                    do i=1,grd_ens%lon2
                       do j=1,grd_ens%lat2
                          w2(j,i) = ps(j,i)
                          x2(j,i)=x2(j,i)+ps(j,i)
                       end do
                    end do
 
                 case('sst','SST')
 ! IGNORE SST IN HYBRID for now
 
                    do i=1,grd_ens%lon2
                       do j=1,grd_ens%lat2
                          w2(j,i) = zero
                          x2(j,i)=zero
                       end do
                    end do
 
              end select
           end do
        enddo 
 !
 ! CALCULATE ENSEMBLE MEAN
        bar_norm = one/float(n_ens)
        en_bar(m)%values=en_bar(m)%values*bar_norm
 
 ! Copy pbar to module array.  ps_bar may be needed for vertical localization
 ! in terms of scale heights/normalized p/p
        do ic2=1,nc2d
  
           if(trim(cvars2d(ic2)) == 'ps'.or.trim(cvars2d(ic2)) == 'PS') then
 
              call gsi_bundlegetpointer(en_bar(m),trim(cvars2d(ic2)),x2,istatus)
              if(istatus/=0) then
                 write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for en_bar to get ps_bar'
                 call stop2(9996)
              end if
  
              do i=1,grd_ens%lon2
                 do j=1,grd_ens%lat2
                    ps_bar(j,i,1)=x2(j,i)
                 end do
              end do
              exit
           end if
        end do
 
        call mpi_barrier(mpi_comm_world,ierror)
 !
 !
 ! CONVERT ENSEMBLE MEMBERS TO ENSEMBLE PERTURBATIONS
        sig_norm=sqrt(one/max(one,n_ens-one))
 
        do n=1,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
           end do
        end do

    enddo ! it 4d loop
 ! CALCULATE ENSEMBLE SPREAD
    write_ens_sprd=.true.
    if(write_ens_sprd ) then
        call this%ens_spread_dualres_regional(mype,en_perts,nelen)
        call mpi_barrier(mpi_comm_world,ierror) ! do we need this mpi_barrier here? 
    endif
    do m=1,ntlevs_ens
      call gsi_bundledestroy(en_bar(m),istatus)
      if(istatus/=0) then
        write(6,*)' in get_fv3_regional_ensperts_netcdf: trouble destroying en_bar bundle'
               call stop2(9997)
      endif
    end do

    deallocate(en_bar)
  !
  
   return

30 write(6,*) 'get_fv3_regional_ensperts_netcdf: open filelist failed '
   call stop2(555)
20 write(6,*) 'get_fv3_regional_ensperts_netcdf: read WRF-ARW ens failed ',n
   call stop2(555)

  end subroutine get_fv3_regional_ensperts_run
  
  subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g_rh,g_oz, &
                                        g_ql,g_qi,g_qr,g_qs,g_qg,g_qnr,g_w)
  !$$$  subprogram documentation block
  !     first compied from general_read_arw_regional           .      .    .                                       .
  ! subprogram:    general_read_fv3_regional  read fv3sar model ensemble members
  !   prgmmr: Ting             org: emc/ncep            date: 2018
  !
  ! abstract: read ensemble members from the fv3sar model in "restart" or "cold start"  netcdf format
  !           for use with hybrid ensemble option. 
  !
  ! program history log:
  !   2018-  Ting      - intial versions  
  !
  !   input argument list:
  !
  !   output argument list:
  !
  ! attributes:
  !   language: f90
  !   machine:  ibm RS/6000 SP
  !
  !$$$ end documentation block
  
    use netcdf, only: nf90_nowrite
    use netcdf, only: nf90_open,nf90_close
    use netcdf, only: nf90_inq_dimid,nf90_inquire_dimension
    use netcdf, only: nf90_inq_varid,nf90_inquire_variable,nf90_get_var
    use kinds, only: r_kind,r_single,i_kind
    use gridmod, only: eta1_ll,eta2_ll
    use constants, only: zero,one,fv,zero_single,one_tenth,h300
    use hybrid_ensemble_parameters, only: grd_ens,q_hyb_ens
    use hybrid_ensemble_parameters, only: fv3sar_ensemble_opt 

    use mpimod, only: mpi_comm_world,mpi_rtype
    use netcdf_mod, only: nc_check
    use gsi_rfv3io_mod,only: type_fv3regfilenameg
    use gsi_rfv3io_mod,only:n2d 
    use constants, only: half,zero
    use gsi_rfv3io_mod, only: gsi_fv3ncdf_read 
    use gsi_rfv3io_mod, only: gsi_fv3ncdf_read_v1
    use gsi_rfv3io_mod, only: gsi_fv3ncdf_readuv
    use gsi_rfv3io_mod, only: gsi_fv3ncdf_readuv_v1
    use gsi_rfv3io_mod, only: gsi_fv3ncdf2d_read_v1
    use directDA_radaruse_mod, only: l_use_dbz_directDA
    use gsi_bundlemod, only: gsi_bundle
    use gsi_bundlemod, only: gsi_gridcreate
    use gsi_bundlemod, only: gsi_grid
    use gsi_bundlemod, only: gsi_bundlecreate,gsi_bundledestroy
    use gsi_bundlemod, only: gsi_bundlegetvar
    use hybrid_ensemble_parameters, only: grd_ens
    use directDA_radaruse_mod, only: l_use_cvpqx, cvpqx_pval, cld_nt_updt
    use directDA_radaruse_mod, only: l_cvpnr, cvpnr_pval



    implicit none
!
! Declare passed variables
    class(get_fv3_regional_ensperts_class), intent(inout) :: this
    type (type_fv3regfilenameg)                  , intent (in)   :: fv3_filenameginput
    real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),intent(out)::g_u,g_v,g_tv,g_rh,g_oz
    real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),optional,intent(out)::g_ql,g_qi,g_qr
    real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),optional,intent(out)::g_qs,g_qg,g_qnr,g_w

    real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2),intent(out):: g_ps
    real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig+1) ::g_prsi 
    real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig) ::g_prsl ,g_tsen,g_q,g_delp
!
! Declare local parameters
    real(r_kind),parameter:: r0_01 = 0.01_r_kind
    real(r_kind),parameter:: r10   = 10.0_r_kind
    real(r_kind),parameter:: r100  = 100.0_r_kind
   !
!   Declare local variables
    
    integer(i_kind):: i,j,k,kp
 
    integer(i_kind) iderivative

    
    logical ice

    character(len=24),parameter :: myname_ = 'general_read_fv3_regional'
    type(gsi_bundle) :: gsibundle_fv3lam_ens_dynvar_nouv
    type(gsi_bundle) :: gsibundle_fv3lam_ens_tracer_nouv
    type(gsi_grid):: grid_ens

    character(len=:),allocatable :: grid_spec !='fv3_grid_spec'            
    character(len=:),allocatable :: ak_bk     !='fv3_akbk'
    character(len=:),allocatable :: dynvars   !='fv3_dynvars'
    character(len=:),allocatable :: tracers   !='fv3_tracer'
    character(len=:),allocatable :: sfcdata   !='fv3_sfcdata'
    character(len=:),allocatable :: couplerres!='coupler.res'
    integer (i_kind) ier,istatus

    
    associate( this => this ) ! eliminates warning for unused dummy argument needed for binding
    end associate



    call gsi_gridcreate(grid_ens,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig)
    grid_spec=fv3_filenameginput%grid_spec
    ak_bk=fv3_filenameginput%ak_bk
    dynvars=fv3_filenameginput%dynvars
    tracers=fv3_filenameginput%tracers
    sfcdata=fv3_filenameginput%sfcdata
    couplerres=fv3_filenameginput%couplerres


     
     
    if (allocated(fv3lam_ens_io_dynmetvars2d_nouv) ) then   
       call gsi_bundlecreate(gsibundle_fv3lam_ens_dynvar_nouv,grid_ens,'gsibundle_fv3lam_ens_dynvar_nouv', istatus,&
                names2d=fv3lam_ens_io_dynmetvars2d_nouv,names3d=fv3lam_ens_io_dynmetvars3d_nouv)
    else
       call gsi_bundlecreate(gsibundle_fv3lam_ens_dynvar_nouv,grid_ens,'gsibundle_fv3lam_ens_dynvar_nouv',istatus, &
                names3d=fv3lam_ens_io_dynmetvars3d_nouv)
    endif

     
    if (allocated(fv3lam_ens_io_tracermetvars2d_nouv) ) then   
       call gsi_bundlecreate(gsibundle_fv3lam_ens_tracer_nouv,grid_ens,'gsibundle_fv3lam_ens_tracer_nouv', istatus,&
                names2d=fv3lam_ens_io_tracermetvars2d_nouv,names3d=fv3lam_ens_io_tracermetvars3d_nouv)
    else
       call gsi_bundlecreate(gsibundle_fv3lam_ens_tracer_nouv,grid_ens,'gsibundle_fv3lam_ens_tracer_nouv',istatus, &
                names3d=fv3lam_ens_io_tracermetvars3d_nouv)
    endif
     
  


     
    if(fv3sar_ensemble_opt == 0 ) then  
      call gsi_fv3ncdf_readuv(grd_fv3lam_ens_uv,g_u,g_v,fv3_filenameginput)
    else
      call gsi_fv3ncdf_readuv_v1(grd_fv3lam_ens_uv,g_u,g_v,fv3_filenameginput)
    endif
    if(fv3sar_ensemble_opt == 0) then
      call gsi_fv3ncdf_read(grd_fv3lam_ens_dynvar_io_nouv,gsibundle_fv3lam_ens_dynvar_nouv,&
                            fv3_filenameginput%dynvars,fv3_filenameginput)
      call gsi_fv3ncdf_read(grd_fv3lam_ens_tracer_io_nouv,gsibundle_fv3lam_ens_tracer_nouv,&
                            fv3_filenameginput%tracers,fv3_filenameginput)
    else
      call gsi_fv3ncdf_read_v1(grd_fv3lam_ens_dynvar_io_nouv,gsibundle_fv3lam_ens_dynvar_nouv,&
                               fv3_filenameginput%dynvars,fv3_filenameginput)
      call gsi_fv3ncdf_read_v1(grd_fv3lam_ens_tracer_io_nouv,gsibundle_fv3lam_ens_tracer_nouv,&
                               fv3_filenameginput%tracers,fv3_filenameginput)
    endif
    ier=0
    call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_dynvar_nouv, 'tsen' ,g_tsen ,istatus );ier=ier+istatus
    call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'q'  ,g_q ,istatus );ier=ier+istatus
    call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'oz'  ,g_oz ,istatus );ier=ier+istatus
    if (l_use_dbz_directDA) then
       call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'ql' ,g_ql ,istatus );ier=ier+istatus
       call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'qi' ,g_qi ,istatus );ier=ier+istatus
       call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'qr' ,g_qr ,istatus );ier=ier+istatus
       call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'qs' ,g_qs ,istatus );ier=ier+istatus
       call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'qg' ,g_qg ,istatus );ier=ier+istatus
       call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'qnr',g_qnr ,istatus );ier=ier+istatus
       call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_dynvar_nouv, 'w' , g_w ,istatus );ier=ier+istatus
    end if
    

    if (fv3sar_ensemble_opt == 0) then 
      call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_dynvar_nouv, 'delp'  ,g_delp ,istatus );ier=ier+istatus
      g_prsi(:,:,grd_ens%nsig+1)=eta1_ll(grd_ens%nsig+1) !thinkto be done , should use eta1_ll from ensemble grid
      do i=grd_ens%nsig,1,-1
         g_prsi(:,:,i)=g_delp(:,:,i)*0.001_r_kind+g_prsi(:,:,i+1)
      enddo
      g_ps(:,:)=g_prsi(:,:,1)
    else  ! for the ensemble processed frm CHGRES
      call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_dynvar_nouv, 'ps'  ,g_ps ,istatus );ier=ier+istatus
      g_ps=g_ps*0.001_r_kind
      do k=1,grd_ens%nsig+1
        g_prsi(:,:,k)=eta1_ll(k)+eta2_ll(k)*g_ps
      enddo

    endif
     
!!  tsen2tv  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    do k=1,grd_ens%nsig
       do j=1,grd_ens%lon2
          do i=1,grd_ens%lat2
             g_tv(i,j,k)=g_tsen(i,j,k)*(one+fv*g_q(i,j,k))
          enddo
       enddo
    enddo
    if (.not.q_hyb_ens) then
      ice=.true.
      iderivative=0
      do k=1,grd_ens%nsig
        kp=k+1
        do j=1,grd_ens%lon2
          do i=1,grd_ens%lat2
            g_prsl(i,j,k)=(g_prsi(i,j,k)+g_prsi(i,j,kp))*half
          end do
        end do
      end do
      call genqsat(g_rh,g_tsen(1,1,1),g_prsl(1,1,1),grd_ens%lat2,grd_ens%lon2,grd_ens%nsig,ice,iderivative)
      do k=1,grd_ens%nsig
        do j=1,grd_ens%lon2
          do i=1,grd_ens%lat2
            g_rh(i,j,k) = g_q(i,j,k)/g_rh(i,j,k)
          end do
        end do
      end do
    else
        do k=1,grd_ens%nsig
          do j=1,grd_ens%lon2
            do i=1,grd_ens%lat2
              g_rh(i,j,k) = g_q(i,j,k)
            end do
           end do
         end do
    end if


! CV transform
    do k=1,grd_ens%nsig
       do i=1,grd_ens%lon2
          do j=1,grd_ens%lat2
             if (l_use_cvpqx) then
                ! Qr/qr
                if (g_qr(j,i,k) <= 1.0E-5_r_kind) then
                    g_qr(j,i,k) = 1.0E-5_r_kind
                end if
                if (cvpqx_pval > 0.0_r_kind ) then !CVpq
                    g_qr(j,i,k)=((g_qr(j,i,k)**cvpqx_pval)-1)/cvpqx_pval
                else  ! CVlogq
                    g_qr(j,i,k) = log(g_qr(j,i,k))
                end if
                ! Qs/qs
                if (g_qs(j,i,k) <= 1.0E-5_r_kind) then
                    g_qs(j,i,k) = 1.0E-5_r_kind
                end if
                if (cvpqx_pval > 0.0_r_kind ) then !CVpq
                    g_qs(j,i,k)=((g_qs(j,i,k)**cvpqx_pval)-1)/cvpqx_pval
                else  ! CVlogq
                    g_qs(j,i,k) = log(g_qs(j,i,k))
                end if
                ! Qg/qg
                if (g_qg(j,i,k) <= 1.0E-5_r_kind) then
                   g_qg(j,i,k) = 1.0E-5_r_kind
                end if
                if (cvpqx_pval > 0.0_r_kind ) then !CVpq
                   g_qg(j,i,k)=((g_qg(j,i,k)**cvpqx_pval)-1)/cvpqx_pval
                else  ! CVlogq
                   g_qg(j,i,k) = log(g_qg(j,i,k))
                end if
             end if
             if ( cld_nt_updt > 0 .and. l_cvpnr) then ! CVpnr
                ! Qnr/qnr
                if (g_qnr(j,i,k) < one) then
                   g_qnr(j,i,k) = one
                end if
                g_qnr(j,i,k)=((g_qnr(j,i,k)**cvpnr_pval)-1)/cvpnr_pval
             end if
          enddo
       enddo
    enddo
    call gsi_bundledestroy(gsibundle_fv3lam_ens_dynvar_nouv)
    call gsi_bundledestroy(gsibundle_fv3lam_ens_tracer_nouv)


  return       
  end subroutine general_read_fv3_regional

  subroutine ens_spread_dualres_regional_fv3_regional(this,mype,en_perts,nelen)
  !$$$  subprogram documentation block
  !                .      .    .                                       .
  ! subprogram:    ens_spread_dualres_regional
  !   prgmmr: mizzi            org: ncar/mmm            date: 2010-08-11
  !
  ! abstract:
  !
  !
  ! program history log:
  !   2010-08-11  parrish, initial documentation
  !   2011-04-05  parrish - add pseudo-bundle capability
  !   2011-08-31  todling - revisit en_perts (single-prec) in light of extended bundle
  !   2021-08-10  lei     - modify for fv3-lam ensemble option
  !   input argument list:
  !      mype  - current processor number
  !
  !   output argument list:
  !
  ! attributes:
  !   language: f90
  !   machine:  ibm RS/6000 SP
  !
  !$$$ end documentation block
  !
    use kinds, only: r_single,r_kind,i_kind
    use hybrid_ensemble_parameters, only: n_ens,grd_ens,grd_anl
    use general_sub2grid_mod, only: sub2grid_info,general_sub2grid_create_info,general_sube2suba
    use constants, only:  zero,two,half,one
    use control_vectors, only: cvars2d,cvars3d
    use gsi_bundlemod, only: gsi_bundlecreate
    use gsi_bundlemod, only: gsi_grid
    use gsi_bundlemod, only: gsi_bundle
    use gsi_bundlemod, only: gsi_bundlegetpointer
    use gsi_bundlemod, only: gsi_bundledestroy
    use gsi_bundlemod, only: gsi_gridcreate
    use write_fv3_enspread_mod, only: write_fv3_enspread
    implicit none

    class(get_fv3_regional_ensperts_class), intent(inout) :: this
    integer(i_kind),intent(in):: mype
    type(gsi_bundle),allocatable, intent(in   ) :: en_perts(:,:)
    integer(i_kind), intent(in   ):: nelen
  
    type(gsi_bundle):: sube,suba
    type(gsi_grid):: grid_ens,grid_anl
    real(r_kind) sp_norm

    integer(i_kind) i,n
    integer(i_kind) istatus
    character(255) enspreadfilename 

    associate( this => this ) ! eliminates warning for unused dummy argument needed for binding
    end associate
          if(mype==0) write(6,*) "Enter ens_spread_dualres_regional_fv3_regional"  
          enspreadfilename="fv3_enspread_regll"
  !      create simple regular grid
          call gsi_gridcreate(grid_anl,grd_anl%lat2,grd_anl%lon2,grd_anl%nsig)
          call gsi_gridcreate(grid_ens,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig)
  
  !      create two internal bundles, one on analysis grid and one on ensemble grid
  
         call gsi_bundlecreate (suba,grid_anl,'ensemble work',istatus, &
                                   names2d=cvars2d,names3d=cvars3d,bundle_kind=r_kind)
         if(istatus/=0) then
            write(6,*)' in ens_spread_dualres_regional: trouble creating bundle_anl bundle'
            call stop2(9998)
         endif
         call gsi_bundlecreate (sube,grid_ens,'ensemble work ens',istatus, &
                                   names2d=cvars2d,names3d=cvars3d,bundle_kind=r_kind)
         if(istatus/=0) then
            write(6,*)' ens_spread_dualres_regional: trouble creating bundle_ens bundle'
            call stop2(9999)
         endif
  
         sp_norm=(one/float(n_ens))
  
         sube%values=zero
  !
  
         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))
            end do
         end do
     
         do i=1,nelen
           sube%values(i) = sqrt(sube%values(i))
         end do
  
         call write_fv3_enspread(grd_ens,enspreadfilename,sube,1)  !"1" as a placeholder
  
    return
  end subroutine ens_spread_dualres_regional_fv3_regional
end module get_fv3_regional_ensperts_mod