program hadley3

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:    Sai Ravela
c***  Last modified: Feb 2, 2003

        implicit none

	integer i, j, k, n, ijim, ijuli, isteps, nens, ics
	integer nx, ny, nz, nf, nnx, nny, nnz, mobs, isai
	integer ii, jj, zob, xob, yob
#include "hadley3.h"
	integer iobsloc(mobs), indxa
	real*8    ytrue(n),  yobs(mobs), ytmp(n)
	real*8    pert(n),  R(mobs), yobsfull(n), ave(n), var(n)
        real*8    xens(n,nens)
	real*8    var_obs, eps, dt, x
	real*8    obserr, analerr, distsub, distobs
        integer iret, system
	real*4    mask(ny,nx)
	
c***  open some output files
        open(unit=2,file='test.dat',status='unknown')
c        open(unit=3,file='fields.dat',status='unknown')

c***  initialisations
	eps=0.005		! obs error - 1/2 mm/s
	isteps=1	! number of steps per obs cycle

c*** 	Here is where you read iobsloc and corresponding
c	variances. You should know
c	mobs in advance...in hadley

	call ReadObs(3,mobs,iobsloc)
	write (*,*) 'iobsloc set',mobs, iobsloc(mobs), nens
	var_obs=eps**2
	do i=1,mobs
	 if (iobsloc(i) < nx*ny*nz*3) then
	    R(i)=var_obs
	 else
	    R(i) = 0.1    ! 0.9 deg other
	 end if
	enddo


c***  loop over initial conditions
 	do ijim=1,ics
	   
	   call ReadPickup(0,n,nx,ny,nz,ytrue)
	   write (*,*) 'Picked Truth'
	   do k = 1,mobs
	      indxa = iobsloc(k)
	      yobs(k) = ytrue(indxa)
	   end do
	   do isai=1,nens
	      call ReadPickUp(isai,n,nx,ny,nz,ytmp)
	      do k=1,n
		 xens(k,isai)=ytmp(k)
	      enddo
	   enddo     
	   write (*,*) 'Picked Ensemble'
	  call EnSRF3d(xens,yobs,iobsloc,n,mobs,R,nens,nx,ny,nz,mask)	   
c***    get ensemble mean and variance
	   call ranmean2(xens,ave,n,nens)
c***    calculate obs and anal err
	   analerr=distobs(yobs,ave,iobsloc,mobs,n)
	   write(6,*) ijim, analerr
	   write(2,*) ijim, analerr
	   analerr=distsub(ytrue,ave,n)
	   write(6,*) '>',ijim, analerr
	   write(2,*) '>',ijim, analerr
c       write(3,'(250f18.10)') (ave(i), i=1,n)
	   

c***    step truth forward
	   call Model(0)
	   call ReadPickup(0,n,nx,ny,nz,ytmp)  
	   call  WritePickUp(0,n,nx,ny,nz,ytmp)
	   write (*,*) 'Stepped Truth'
c***    step ensemble forward
	   do isai=1,nens
	      do k=1,n
		 ytmp(k)=xens(k,isai)
	      enddo
	      call  WritePickUp(isai,n,nx,ny,nz,ytmp)
	      call Model(isai)
	   enddo	   
	   write (*,*) 'Stepped Ensemble'
	   call flush()
	   
	enddo
	
	return
	end