cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c User assigned psudeu-distance function.  This is used to select data points
c used in the objective analysis.  
      function pdist(x,y,h,t)
      real x,y,h,t,pdist
      pdist=x*x + y*y 
c      pdist=x*x + y*y + h*h*10000.
c      pdist=sqrt(x*x+y*y+h*h*10000)  ! RHE 04/12/2004
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      
      subroutine modelcov(c,x1,y1,z1,h1,t1,x2,y2,z2,h2,t2)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c User assigned covariance function.  This returns the covariance between 
c two points for the residual field.c=E[(u1-mu1)(u2-mu2)]
c c (output) covariance between (x1,y1,z1,h1,t1)and(x2,y2,z2,h2,t2) 
c x1 (x2) - x coordinate of the first(second) point
c y1 (y2) - y coordinate of the first(second) point
c z1 (z2) - z coordinate of the first(second) point
c h1 (h2) - bathymetric depth at x1,y1 (x2,y2)
c t1 (t2) - time(yearday) of the first(second) point
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      real c,x1,y1,z1,t1,x2,y2,z2,t2,s,h1,h2
      real v,lambda_x,lambda_z,lambda_h,lambda_t
      v=3.

c     lambda_x=1.e-10
c     lambda_z=0.0025
c     lambda_h=0.0001
c     lambda_t=0.0002778
     
c     s=sqrt(  ((x1-x2)**2 + (y1-y2)**2 )*lambda_x
c    & + ((z1-z2)**2)*lambda_z
c    & + ((h1-h2)**2)*lambda_h
c    & + ((t1-t2)**2)*lambda_t  )

      lambda_x=50000.
      lambda_z=10.        ! 20 m
      lambda_h=100.
      lambda_t=60.
      s=sqrt(  ( (x1-x2)**2 + (y1-y2)**2 )/(lambda_x**2)
     & + ((z1-z2)**2)/(lambda_z**2)
     & + ((h1-h2)**2)/(lambda_h**2)   ! &
     & + ((t1-t2)**2)/(lambda_t**2)  )

      c=v*exp(-s)
      
      return
      end 
      
      subroutine obsercov(c,x1,y1,z1,h1,s1,i1,x2,y2,z2,h2,s2,i2)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c User assigned covariance function.  This returns the observational error 
c covariance between two points,(x1,y1,z1,h1,t1)and(x2,y2,z2,h2,t2).  i1 
c and i2 are the instrument numdber for the first and second obsevation.
c c (output) observartional error covariance between 
c           (x1,y1,z1,h1,t1)and(x2,y2,z2,h2,t2) 
c x1 (x2) - x coordinate of the first(second) measurement
c y1 (y2) - y coordinate of the first(second) measurement
c z1 (z2) - z coordinate of the first(second) measurement
c h1 (h2) - bathymetric depth at x1,y1 (x2,y2)
c t1 (t2) - time(yearday) of the first(second) measurement
c s1 (s2) - standart error for first(second) measurement
c i1 (i2) - instrument number for first(second) measurement
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      real c,x1,y1,z1,s1,t1,x2,y2,z2,s2,t2,s,h1,h2
      integer i1,i2
      c=0.	
      if (i1.eq.i2) then
c       s=sqrt( (x1-x2)**2 + (y1-y2)**2 )/100.0 + abs(z1-z2)/1.0
       s=sqrt( (x1-x2)**2 + (y1-y2)**2 )/100.0 + abs(z1-z2)/10.0
        c=s1*s2*exp(-(s))
      endif
      return
      end
      

      subroutine oaci_out(nn,nnv,meshname,z,p,s,v,nr,spp) 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  subroutine oaci_out(nn,nnv,meshname,z,p,s,v,nr,spp) 
c  Write output from oaci interpolation
c 
c
c  nn - number of horizontal nodes in mesh
c  nnv - number of sigma layers inmesh
c  meshname - yep
c  z(n,j) z coordinate of sigma layer j at node n
c  p(n,j) estimate of perturbation to field at  sigma layer j of node n
c  s(n,j) prior estimate of scalar field at  sigma layer j of node n
c  v(n,j) posterior error variance at  sigma layer j of node n
c  spp(n,j) posterior estimate of scalar field at  sigma layer j of node n
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      implicit none
      include 'oaci.dim'
C Mesh
      real z(NNdim,nnvdim),s(nndim,nnvdim),p(nndim,nnvdim),
     &v(nndim,nnvdim),spp(nndim,nnvdim),nr(nndim,nnvdim)
      integer nn,nnv
      character*72 meshname,filename,header   
c unit 4 is oaci.inp      
      
      read(4,'(A)') filename
      if (filename(1:2).NE.'NO')then
        header='interpolated data-climatology mifit'
        call WRITES3R(nndim,nnvdim,filename,meshname,header,z,p,
     &                  nn,nnv)
      endif

      read(4,'(A)') filename
      if (filename(1:2).NE.'NO')then
        header='posterior variance'
        call WRITES3R(nndim,nnvdim,filename,meshname,header,z,v,
     &                  nn,nnv)
      endif
      
      read(4,'(A)') filename
      if (filename(1:2).NE.'NO')then
        header='analysis field'
        call WRITES3R(nndim,nnvdim,filename,meshname,header,z,spp,
     &                  nn,nnv)
      endif
      
      return
      end
