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, ijuli, isteps 
        integer nens            ! ensemble size
        integer iters           ! number of iterations
        integer n               ! state size, calculated from
	integer nx, ny, nz      ! physical dimension sizes, and
        integer nf              ! number of observed variables
        integer mobs, isai
	integer ii, jj, zob, xob, yob
#include "hadley3.h"
	integer iobsloc(mobs)   ! obs locations, in terms of 
                                ! dimensions in the state vector.
                                ! Indicated with the matrix H in KF 
                                ! literature. 
                                ! In practical terms, the numbers
                                ! the locations of the words in the 
                                ! pickup file
        integer  indxa
	real*8    ytrue(n)      ! state vector
        real*8    yobs(mobs)    ! current observations of state 
        real*8    ytmp(n)       ! tmp state vector
	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.0005		! obs error variance - 1/2 mm/s
                                ! only good for u and v 
                                ! need to vary this for
                                ! different variables

c*** 	Here is where you read iobsloc and set the corresponding
c     variances.  
	call ReadObsLoc(3,mobs,iobsloc)
	write (*,*) 'iobsloc set',mobs, iobsloc(mobs), nens
	var_obs=eps**2
	do i=1,mobs
c     why the if/then?
           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*** get initial observations -- assumes that model has
c    already spun up (how long?)

        write (*,*) 'getting Truth'	   
        call ReadPickup(0,n,nx,ny,nz,ytrue)
        write (*,*) 'got Truth'

c***  fantasy observations
        write (*,*) 'generating observations'	   
	call random2(pert,n)
        yobsfull = ytrue+eps*pert
c*** apply observation mask
	do i = 1,mobs
	   indxa = iobsloc(i)
	   yobs(i) = yobsfull(indxa)
	end do

        write (*,*) 'generating ensemble'	   
        do i=1,nens
           call random2(pert,n)           
           do j=1,n
              xens(j,i)=yobsfull(j)+eps*pert(j)
              ytmp(j)=xens(j,i)
           enddo
           call  WritePickUp(i,n,nx,ny,nz,ytmp)
        enddo     

c***  main iteration loop
 	do i=1,iters

c***  step truth forward
	   write (*,*) 'Stepping Truth'
	   call Model(0)
c     call ReadPickup(0,n,nx,ny,nz,ytmp)  
c     call  WritePickUp(0,n,nx,ny,nz,ytmp)
	   write (*,*) 'Stepped Truth'

           write (*,*) 'getting Truth'	   
           call ReadPickup(0,n,nx,ny,nz,ytrue)
           write (*,*) 'got Truth'
           
c***  fantasy observations
           write (*,*) 'generating observations'	   
           call random2(pert,n)
           yobsfull = ytrue+eps*pert
c***  apply observation mask
           do j = 1,mobs
              indxa = iobsloc(j)
              yobs(j) = yobsfull(indxa)
           end do

c***  step ensemble forward
	   write (*,*) 'Stepping Ensemble'
	   do j=1,nens
	      call Model(j)
	      call  ReadPickUp(j,n,nx,ny,nz,ytmp)
	      do k=1,n
		 xens(k,j)=ytmp(k)
	      enddo
	   enddo	   
	   write (*,*) 'Stepped Ensemble'

           
           call EnSRF3d(xens,yobs,iobsloc,n,mobs,R,nens,nx,ny,nz,mask)	   

	   do j=0,nens
	      do k=1,n
		 ytmp(k)=xens(k,j)
	      enddo
	      call  WritePickUp(j,n,nx,ny,nz,ytmp)
	   enddo	   


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(yobs,yobsfull,n)
           analerr=distsub(yobs,ave,n)
           write(6,*) "\tOBSERR\t", "\tANALERR"
           write(6,*) obserr, analerr
           write(2,*) iters, obserr, analerr
           write(3,'(250f18.10)') (ave(i), i=1,n)

           call flush()
           
	enddo
	
	return
	end