subroutine greg_sqrt(xens,yo,H,ngp,mobs,RRs,nens)

c      implicit none

	include 'mpif.h'

        integer   MASTER, FROM_MASTER, FROM_WORKER
        parameter (MASTER = 0)
        parameter (FROM_MASTER = 1)
        parameter (FROM_WORKER = 2)

        integer   taskid, ierr, numtasks, numworkers, snumworkers
        integer   source, dest, mtype, num, avenum, extra, offset
        integer   status(MPI_STATUS_SIZE)

! Arguments
      integer, intent(in)    :: nens, mobs, ngp
      real,    intent(inout) :: xens(ngp,nens)
      real,    intent(in)    :: yo(mobs), RRs(mobs), H(mobs,ngp)

! Local Variables
      integer :: rad, xob(mobs), yob(mobs), ind, nx, ny
      integer :: k, j, i, r2, kk, jj, ko, jo, kj, g
      real    :: PHT(ngp), HPHT, Ks(ngp), Khat(ngp), alpha
      real    :: xp(ngp), xa(ngp), zp(ngp,nens), R, Rs, boost
      real, parameter :: eps = 1.e-6 ! corresponds to verr = 0.001
! Filter Stuff
      integer :: d, fexp
      real    :: dr, fcoef, filt

        call MPI_COMM_RANK(MPI_COMM_WORLD,taskid,ierr)
        call MPI_COMM_SIZE(MPI_COMM_WORLD,numtasks,ierr)
        

      if (taskid.eq.MASTER) then


	write(6,*) 'a'

      nx = 33
      ny = 65

      Rs = RRs(1)
      R = sqrt(Rs)

      rad = 20
      rad = 100
      r2 = rad*rad

      fexp = 2
      fcoef = 0.05
      boost = 1.05
      boost =1.0

      zp = xens

!*** Transform 2D state vector into 1D vector....
!  do j = 1, nens
!     do k = 1, ny
!        zp(nx*(k-1)+1:k*nx,j) = xens(:,k,j)
!     end do
!  end do

	write(6,*) 'b'

!*** Find the initial ensemble mean
      do j = 1, ngp
         xp(j) = sum(zp(j,:))/float(nens) 
      end do

!*** Apply inflation factor to initial ensemble
      do k = 1, nens
         zp(:,k) = boost*(zp(:,k) - xp) + xp
      end do

	write(6,*) 'c'

!*** Find the xob and yob arrays from H
      do j = 1, ngp
         do k = 1, mobs
            if ( H(k,j) == 1. ) then
               xob(k) = mod(j-1,nx) + 1
               yob(k) = (j-1)/nx + 1
            end if
         end do
      end do


	write(6,*) 'd'

!*** Now process each observation sequentially abiding by cut-off radius
      do j = 1, mobs
         ind = nx*( yob(j) - 1 ) + xob(j)

	write(6,*) 'e'

!*** Find PH' first
         PHT = 0.0
         do jj = yob(j)-rad, yob(j)+rad 
            do kk = xob(j)-rad, xob(j)+rad 
               jo = jj
               ko = kk

	write(6,*) 'f', j, jo, ko
           
!*** Point is within block of radius, but it may not be within the 
!      basin boundaries
               if ( ko>0 .and. ko<=nx .and. jo>0 .and. jo<=ny ) then

!*** Since we've sequestered a square of side 2*rad and the
!      cut-off radius assumes a circle, we need to check to
!      make sure the point we're considering is actually 
!      within the circle.

	write(6,*) 'g'

                  d = (ko - xob(j))**2 + (jo - yob(j))**2   

	write(6,*) 'h'

                  if ( d <= r2 ) then
                 
                     dr = sqrt( float( d ) ) 
           
!*** The element of interest in the 1D vector according to addresses 
!      kk and jj is:
                     kj = nx*(jo-1) + ko

	write(6,*) 'i'

!*** Evaluate the filter coefficient based on distance from center d
                filt = 1.0 - exp( -fcoef*( (float(rad) - dr)**fexp ) )
		filt = 1.0

	write(6,*) 'j'

!*** Now contribute to PHT sum
                do g = 1, nens
                   PHT(kj) = PHT(kj) + filt*(zp(kj,g) - xp(kj))*  
     &                  (zp(ind,g) - xp(ind))
                end do

	write(6,*) 'k'
              
                  end if
               end if
            end do
         end do
         PHT = PHT/float(nens - 1)

!*** Now find HPH' from PH'.  Because of cut-off radius, there is a
!      (good) chance that HPH' will be zero.
         HPHT = PHT(ind)

!*** Evaluate Ks
         Ks = PHT/( HPHT + Rs )

!*** Update all effected elements in the mean
         xa = xp + Ks*( yo(j) - xp(ind) )

!*** Now update all ensemble members as perturbations about mean
         alpha = 1.0/( 1.0 + sqrt( Rs/( HPHT + Rs ) ) ) 
         Khat = alpha*Ks
            
         do g = 1, nens
            zp(:,g) = ((zp(:,g) - xp) - Khat*( zp(ind,g) - 
     &           xp(ind) )) + xa
         end do

!*** Use analysis ensemble as the background for the next observation
         xp = xa

      end do

      print*, 'EnSRF:: Done Mobs Loop'

      xens = zp

      endif

      return
  
      end subroutine greg_sqrt