c***  ensemble kalman filter subroutine

	subroutine EnKF(A,z,H,ns,no,varmes,nrsamp)

c***  A is the matrix of ensemble forecasts, z is the observation,
c***  H is the observation operator, ns is the size of the model
c***  state, no is the number of observation locations, varmes is
c***  the observational uncertainty, and nrsamp is the ensemble size.

	integer   ns, no, nrsamp
	integer   ijim

	real, intent(inout)         :: A(ns,nrsamp)
	real, intent(in)            :: z(no)

   	integer i,j,m,mm,info,idum,ibawk
        real    rcond
	real    H(no,ns)
	real    vave(ns), vvar(ns), varmes(no)
	real    temp

   	real, allocatable    :: K(:,:)
   	real, allocatable    :: RR(:,:)
	real, allocatable    :: RRR(:,:)
   	real, allocatable    :: w(:)
   	real, allocatable    :: tmp(:,:)
   	real, allocatable    :: d(:)
	real, allocatable    :: pd(:,:)
	real, allocatable    :: transA(:,:)
	real, allocatable    :: transtmp(:,:)
	real, allocatable    :: tmp2(:,:)


   	allocate (K(no,ns), RR(no,no), RRR(no,no))
	allocate (transA(nrsamp,ns),transtmp(nrsamp,no))
   	allocate (w(ns), tmp(no,nrsamp))
   	allocate (d(no), pd(1,ns), tmp2(1,ns))

c***  initialising the observational error covariance matrix and
c***  the projector from model space to observational space.
   	 RR=0.0

c***  define obs error covariance matrix as a diagonal matrix
   	 do m=1,no
	  RR(m,m)=varmes(m)
   	 enddo

c***  calculate ensemble mean
   	 call ranmean2(A,w,ns,nrsamp)

c***  remove ensemble mean from the ensemble members
   	 do i=1,nrsamp
   	  A(:,i)=A(:,i)-w(:)
   	 enddo

c***  tmp takes on the product of the observation projector and 
c***  the ensemble anomalies.
   	 tmp=matmul(H,A)

c***  K takes on the value of H(A-ABAR)(A-ABAR)^T
	 K=matmul(tmp,transpose(A))

c***  K is what is called B in eqn (6)
   	 K=(1.0/(float(nrsamp)-1.0))*K

c***  restore A by adding back ensemble mean value.
   	 do i=1,nrsamp
   	  A(:,i)=A(:,i)+w(:)
   	 enddo  

c***  RR is obs error cov matrix, and the forecast error cov is added to
c***  it, making it the first half of the right hand side of equation (8).
c***  if one was to wipe out spurious long distance correlations, it would
c***  happen here.
   	RR=RR+(1./(float(nrsamp)-1.))*matmul(tmp,transpose(tmp))

c***  standard solver for coefficients beta
   	 call spoco(RR,no,no,rcond,w,info)

c***  correct each member of the ensemble
   	 do j=1,nrsamp

   	  call random2(d,no)
	  d=z+sqrt(varmes)*d-matmul(H,A(:,j))
   	  call sposl(RR,no,no,d)
   	  A(:,j)=A(:,j)+matmul(d,K) 

   	 enddo  		!A has now been corrected by the data.

   	deallocate (K, RR, RRR)
	deallocate (transA, transtmp)
   	deallocate (w, tmp)
   	deallocate (d,pd,tmp2)

   	return
   	end