program driver

c***  This subroutine calculates analysis errors using greg's
c***  ensemble square root filter, or the traditional ensemble
c***  kalman filter.
c***
c***  Written by:    Jim Hansen
c***  Last modified: Nov 27, 2002

        implicit none

	integer i, j, k, n, ijim, ijuli, isteps, nens, ics
	integer nx, ny, nz, nf, nnx, nny, nnz, mobs
#include "parameters.h"
	integer iobsloc(mobs)
	real*8    y(n), dydx(n), yobs(mobs), ytmp(n)
	real*8    pert(n), H(mobs,n), R(mobs), yobsfull(n)
        real*8    xens(n,nens), par(n), var(n), ave(n)
	real*8    var_obs, eps, dt, x
	real*8    obserr, analerr, distsub
	real*8    bcup, bcdown, bcin, bcout

	common /boundary/ bcup,bcdown,bcin,bcout
	common /dim/      nnx,nny,nnz

        external derivsL953d

c***  open some output files
        open(unit=2,file='test.dat',status='unknown')
        open(unit=3,file='fields.dat',status='unknown')

c***  initialisations
	eps=0.2		! obs error variance
	isteps=1	! number of steps per obs cycle
	dt=0.05		! integration stepsize
	do i=1,n
	 par(i)=8.	! forcing magnitude
	enddo
	bcup=5.		! boundary conditions 
	bcdown=5.
	bcin=5.
	bcout=5.
	nnx=nx
	nny=ny
	nnz=nz

c***  specify observation locations
	do i=1,mobs
	 iobsloc(i)=i
	enddo
	do i=1,mobs
	 do j=1,n
	  H(i,j)=0.
	 enddo
	enddo
	do i=1,mobs
	 H(i,iobsloc(i))=1.
	enddo

c***  define obs variance
	var_obs=eps**2
	do i=1,mobs
	 R(i)=var_obs
	enddo

c***  make up some initial conditions
	call random2(pert,n)
	do i=1,n
	 y(i)=par(1)*pert(i)
	enddo

c***  remove transient
	do i=1,2**14	
         call stepit(x,y,dydx,par,dt,n,isteps,derivsL953d)
	enddo

c***  define an observation about truth
	call random2(pert,n)
        do i=1,n
         yobsfull(i)=y(i)+eps*pert(i)
        enddo

c***  generate an nens member ensemble of initial conditions.
        do i=1,nens
         call random2(pert,n)
         do j=1,n
          xens(j,i)=yobsfull(j)+eps*pert(j)
         enddo
        enddo   

c***  loop over initial conditions
 	do ijim=1,ics

c***  step truth forward
         call stepit(x,y,dydx,par,dt,n,isteps,derivsL953d)

c***  define an observation about truth   
	 call random2(pert,n)
         do k=1,n
          yobsfull(k)=y(k)+eps*pert(k)
         enddo
	 do k=1,mobs
	  yobs(k)=yobsfull(iobsloc(k))
	 enddo

c***  step ensemble forward
	 do ijuli=1,nens
	  do k=1,n
	   ytmp(k)=xens(k,ijuli)
	  enddo
	  call stepit(x,ytmp,dydx,par,dt,n,isteps,derivsL953d)
	  do k=1,n
	   xens(k,ijuli)=ytmp(k)
	  enddo
	 enddo

c***  call the filter
c	 call EnSRF3dO(xens,yobs,H,n,mobs,R,nens,nx,ny,nz)
	 call EnKF(xens,yobs,H,n,mobs,R,nens)

c***  get ensemble mean and variance
	 call ranmean2(xens,ave,n,nens)
	 call ranvar(xens,ave,n,nens,var)

c***  calculate obs and anal err
	 obserr=distsub(y,yobsfull,n)
	 analerr=distsub(y,ave,n)
	 write(6,*) "\tOBSERR\t", "\tANALERR"
	 write(6,*) obserr, analerr
	 write(2,*) ijim, obserr, analerr
	 write(3,'(250f18.10)') (ave(i), i=1,n)

	 call flush()

	enddo
	
	return
	end