C***********************************************************************
C***********************************************************************
C23456789012345678901234567890123456789012345678901234567890123456789012 
C
C Get the unexplained velocity data (errors) which are used to force moody. 
C The errors are defined as Observations - Model 
C
C Manipulate the data: 
C   1) remove data which is outside the mesh.
C   2) sort the data by time (in progress).
C   3) compute the Gregorian time for use by casco. 
C   4) remove data which is outside the simulation time window. 
C 
C 15 Jan 1999  CGH and DRL     
C 11 Nov 1999  CGH check mesh name from *uxv with specified meshname 
C              in misc.h 
C 6 Nov 2002.  CGH. Change name as part of scheme to add elevation
C              Need to check for none and NONE as file names.  
C              The removing the data outside the mesh seems to be 
C              disabled.  Why? 


      subroutine getvelobs

      include 'CASCO.DIM'
      include 'observes.h' 
      include 'mesh.h' 
      include 'startstop.h' 
      include 'misc.h' 
C
C
C Local arrays for sorting
C
      Integer indx(nobsdim)
      Real    temp(nobsdim)
C
C local variables
      real Tobsgmt0(nobsdim)
      character*120 filename

      integer lun
      integer iyear
      integer i 
      integer nobsin, nclosest, nelem 
      real    basisxy(3)
      logical keep 
      integer ncheck1, ncheck2
      character*72 nmesh2 
      integer iend1, iend2 

C functions
      integer igday_compare

      write(*,*)'GETVELOBS: '
      write(*,*)'Enter unexplained velocity file (*.uxv)'
      write(*,*)'format: t x y u v zmin zmax zbot *'
      write(*,*)'Enter "none" if none is desired ' 
      read(*,'(a)')filename
      if(filename(1:4).eq.'none'.or.filename(1:4).eq.'NONE') then
	  write(*,*)'No velocity observations provided.'
          nobsin = 0
	  nobs   = 0
	  return
      endif 
      lun = 49 
      open(unit=lun, file=filename, status='old')

      call getobs_8c(lun,Nobs, Xobs, Yobs, Tobsgmt0,
     &      Uobs,Vobs,Zmaxobs,Zminobs,Zbottomobs,nobsdim,iyear,nmesh2)

      iend1 = index(meshname,' ') - 1
      iend2 = index(nmesh2,' ') - 1
      write(*,*)meshname(1:iend1)
      write(*,*)nmesh2(1:iend2)
      if(meshname(1:iend1).ne.nmesh2(1:iend2)) then 
         write(*,*)'**GETVELOBS: Fatal Error**'
	 write(*,*)'Mesh name from *uxv file is not the same as the '
	 write(*,*)'mesh name given to the forward model.'
	 write(*,*)meshname(1:iend1)
	 write(*,*)nmesh2(1:iend2)
	 STOP 
      endif
      write(2,21) nobs,filename
      write(*,*)nmesh2(1:5)
21    format('GetObs read',i5,' observations from file ',(a))
      write(2,*) ' '
C
C 
C   Discard observations which lie outside the boundary of the finite
C   element mesh
         NOBSIN=0
         DO I=1,NOBS
            CALL BASIS2D(NNDIM,NEDIM,IN,X,Y,NN,NE,Xobs(I),Yobs(I),
     &                NCLOSEST,NELEM,BASISXY)
            IF(NELEM.NE.0)THEN
               NOBSIN=NOBSIN+1
               Xobs(NOBSIN) = Xobs(i)
               Yobs(NOBSIN) = Yobs(i)
               Tobsgmt0(NOBSIN) = Tobsgmt0(i)
               Uobs(NOBSIN) = Uobs(i)
               Vobs(NOBSIN) = Vobs(i)
               Zmaxobs(NOBSIN) = Zmaxobs(i)
               Zminobs(NOBSIN) = Zminobs(i)
               Zbottomobs(NOBSIN) = Zbottomobs(i)
            ENDIF
         ENDDO
         WRITE(*,'(I6,A,I6,3A)')NOBSIN,' of the original ',NOBS,
     &   ' O*D entries are in the domain.'
         NOBS=NOBSIN
         WRITE(*,*)NOBSIN,NOBS

1001     continue 
C
C
C
C Sort the observations by time, in ascending order.
C  First, index the time array.  INDEXX is a Recipes routine which
C     computes the ordering array indx(i).
C  Then, sort all the 1-D arrays.
C
      Call INDEXX(nobs, tobsgmt0, indx)
C
      Call SortDRL(      xobs, nobs, indx, temp)
      Call SortDRL(      yobs, nobs, indx, temp)
      Call SortDRL(  tobsgmt0, nobs, indx, temp)
      Call SortDRL(      uobs, nobs, indx, temp)
      Call SortDRL(      vobs, nobs, indx, temp)
      Call SortDRL(   zmaxobs, nobs, indx, temp)
      Call SortDRL(   zminobs, nobs, indx, temp)
      Call SortDRL(zbottomobs, nobs, indx, temp)
C
C The sorting is complete.
C



C
C Convert UTC0 to Gregorian time 
         do i = 1, nobs
            CALL UTC0_GDAY(IYEAR,Tobsgmt0(i),kdobs(i),secobs(i))
         enddo

C Now eliminate all observations which do not within the
C integration time of the model 
C Defined by KD0, SEC0, KDL, SECL

         write(*,*)'start ',KD0, SEC0
         write(*,*)'stop  ',KDL, SECL
         nobsin = 0 
         DO I=1,NOBS
c           write(*,*)'obs i ',kdobs(i),secobs(i),Tobsgmt0(i), iyear
	    ncheck1 = igday_compare(kd0,sec0,kdobs(i),secobs(i))
	    ncheck2 = igday_compare(kdl,secl,kdobs(i),secobs(i))
	    keep = .false. 
	    if((ncheck1.eq.-1).and.(ncheck2.eq.1)) keep = .true.

	    if(keep.eqv..true.) then 
c           write(*,*)'keep ',kdobs(i),secobs(i),Tobsgmt0(i), iyear
               NOBSIN=NOBSIN+1
               Xobs(NOBSIN) = Xobs(i)
               Yobs(NOBSIN) = Yobs(i)
               kdobs(NOBSIN)  = kdobs(i)
               secobs(NOBSIN) = secobs(i)
               Uobs(NOBSIN)   = Uobs(i)
               Vobs(NOBSIN)   = Vobs(i)
               Zmaxobs(NOBSIN) = Zmaxobs(i)
               Zminobs(NOBSIN) = Zminobs(i)
               Zbottomobs(NOBSIN) = Zbottomobs(i)
            endif
         ENDDO
         WRITE(*,'(I6,A,I6,3A)')NOBSIN,' of the original ',NOBS,
     &   ' O*D entries are in the correct time window.'
         NOBS=NOBSIN
         WRITE(*,*)NOBSIN,NOBS
 	 if(nobs.lt.1) then 
            write(*,*)'**FATAL ERROR: in getvelobs**'
 	    write(*,*)'No velocity observations are in the correct',
     &      ' space/time window.  I will stop now.'
 	   stop
	 endif


C write the unexplained velocity observations to the error variables
C moved to error_init
c     do i = 1, nobs
c uerror(i) = Uobs(i)
c verror(i) = Vobs(i)
c     enddo
C

      RETURN
      END


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      Subroutine SortDRL(a, nobs, indx, temp)
C
C arrays passed:
C
      real    a(nobs), temp(nobs)
      integer indx(nobs)
      integer nobs
c
      integer i

      do i=1,nobs
         temp(i)=a(i)
      enddo
c
      do i=1,nobs
         a(i)=temp(indx(i))
      enddo
c
      return
      end

