
      implicit none
      include 'oaci.dim'
C Mesh
      real hdown(NNdim), x(NNdim), y(NNdim)
      real z(NNdim,nnvdim),s(nndim,nnvdim),p(nndim,nnvdim),
     &v(nndim,nnvdim),spp(nndim,nnvdim)
      integer nn,ne,nnv,in(3,nedim)
      character*72 meshname,icsname,obsname,filename,header   
c
c Observations
c              
      real tobs(nobsdim),xobs(nobsdim),yobs(nobsdim),zobs(nobsdim),
     &sobs(nobsdim),stdobs(nobsdim),hdobs(nobsdim),r(nobsdim),t
      integer nobs,inst(nobsdim)
      real nr(nndim,nnvdim)
      integer k,j,n

      open(unit=2,file='oaci.out', status='unknown')
      open(unit=4,file='oaci.inp', status='unknown')

C
C Get the mesh
C
      read(4,'(A)') meshname
      call loadmeshgom(nndim,nedim,meshname,x,y,in,hdown,nn,ne)      
C
C load 3d climatological field
C
      read(4,'(A)') icsname
      call LOADS3R(nndim,nnvdim,nn,nnv,icsname,s,z)
c
c load obsercation with insturment noise terms from mir3d file
c      
      read(4,'(A)') obsname
      call loadmir3d(nobsdim,nobs,obsname,tobs,xobs,yobs,zobs,sobs,
     &stdobs,inst)
      write(*,*) 'read ', nobs, ' observations from ',obsname
c
c compute difference residual field, r, and remove obs outside mesh.
c      
      write(*,*)'computing residual'
      call compute_dmm(nobs,tobs,xobs,yobs,zobs,sobs,inst,r,s,z,x,
     &y,in,nn,ne,nnv)
      write(*,*) nobs, ' of the origonal observations in domain'
c
c write residual field, r to residual.dat
c
      open(7,file='residual.dat')
      do k=1,nobs
       write(7,*)tobs(k),xobs(k),yobs(k),zobs(k),r(k),stdobs(k),inst(k)
      enddo
      close(7)
c
c intepolate hdown to observation points
c
      call interphdown2obs(nobs,xobs,yobs,hdobs,x,y,hdown,in,nn,ne)

      read(4,*)t
c
c carry out oa back to initial 3d mesh nodes
c
      call oa2mesh(nn,nnv,x,y,z,hdown,t,nobs,tobs,xobs,yobs,zobs,hdobs,
     &             stdobs,inst,r,p,v,nr)
c
c add interpolated residual field to background estimate
c
      do n=1,nn
       do j=1,nnv
        spp(n,j)=s(n,j)+p(n,j)
       enddo
      enddo
      
      call oaci_out(nn,nnv,meshname,z,p,s,v,nr,spp)

      write(*,*) 'Finished.'
      close(2)
      close(4)
      STOP
      END
C      
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
