diff --git a/src/EnKF/gfs/src/getsfcensmeanp.fd/getsfcensmeanp.f90 b/src/EnKF/gfs/src/getsfcensmeanp.fd/getsfcensmeanp.f90
index 2589a086..e0efd0a4 100644
--- a/src/EnKF/gfs/src/getsfcensmeanp.fd/getsfcensmeanp.f90
+++ b/src/EnKF/gfs/src/getsfcensmeanp.fd/getsfcensmeanp.f90
@@ -42,14 +42,15 @@ program getsfcensmeanp
   logical:: nemsio, sfcio, ncio
 
   type(Dataset) :: dset,dseto
-  type(Dimension) :: londim,latdim
+  type(Dimension) :: londim,latdim,levdim
   real(4), allocatable, dimension(:,:) :: values_2d, values_2d_avg
+  real(4), allocatable, dimension(:,:,:):: values_3d, values_3d_avg
 
   character*500 filenamein,filenameout,datapath,fileprefix
   character*3 charnanal
   integer lunin,lunout,iret,nanals,k
   integer mype,mype1,npe,orig_group, new_group, new_comm
-  integer nrec, lonb, latb, n, npts, nvar
+  integer nrec, lonb, latb, levs, n, npts, nvar
   integer,dimension(7):: idate
   integer,dimension(:),allocatable:: new_group_members
   real(8) rnanals
@@ -308,7 +309,9 @@ program getsfcensmeanp
      elseif (ncio) then
         londim = get_dim(dset,'grid_xt'); lonb = londim%len
         latdim = get_dim(dset,'grid_yt'); latb = latdim%len
+        levdim = get_dim(dset,'pfull');   levs = levdim%len
         allocate(values_2d_avg(lonb,latb))
+        allocate(values_3d_avg(lonb,latb,levs))
         if (mype == 0) then
            dseto = create_dataset(filenameout, dset, copy_vardata=.true.)
            print *,'opened netcdf file ',trim(filenameout)
@@ -328,16 +331,30 @@ program getsfcensmeanp
                trim(dset%variables(nvar)%name) == 'orog') then
               cycle
            endif
-           call read_vardata(dset,trim(dset%variables(nvar)%name),values_2d)
-           call mpi_allreduce(values_2d,values_2d_avg,lonb*latb,mpi_real4,mpi_sum,new_comm,iret)
-           values_2d_avg = values_2d_avg * rnanals
-           if (mype == 0) then
-             print *,'writing ens mean ',trim(dset%variables(nvar)%name)
-             call write_vardata(dseto,trim(dset%variables(nvar)%name),values_2d_avg)
+           if (dset%variables(nvar)%ndims == 3) then
+              call read_vardata(dset,trim(dset%variables(nvar)%name),values_2d)
+              call mpi_allreduce(values_2d,values_2d_avg,lonb*latb,mpi_real4,mpi_sum,new_comm,iret)
+              values_2d_avg = values_2d_avg * rnanals
+              if (mype == 0) then
+                 print *,'writing ens mean ',trim(dset%variables(nvar)%name)
+                 call write_vardata(dseto,trim(dset%variables(nvar)%name),values_2d_avg)
+              endif
+           elseif (dset%variables(nvar)%ndims == 4) then
+              call read_vardata(dset,trim(dset%variables(nvar)%name),values_3d)
+              call mpi_allreduce(values_3d,values_3d_avg,lonb*latb*levs,mpi_real4,mpi_sum,new_comm,iret)
+              values_3d_avg = values_3d_avg * rnanals
+              if (mype == 0) then
+                 print *,'writing ens mean ',trim(dset%variables(nvar)%name)
+                 call write_vardata(dseto,trim(dset%variables(nvar)%name),values_3d_avg)
+              endif
+           else
+              write(6,*)'***ERROR*** invalid ndims= ',dset%variables(nvar)%ndims,&
+                   ' for variable ',trim(dset%variables(nvar)%name)
+              call MPI_Abort(MPI_COMM_WORLD,98,iret)
            endif
         enddo
         if (mype == 0) call close_dataset(dseto)
-        deallocate(values_2d,values_2d_avg)
+        deallocate(values_2d,values_2d_avg,values_3d,values_3d_avg)
      endif
 ! Jump here if more mpi processors than files to process
   else