C***********************************************************************
C***********************************************************************
C23456789012345678901234567890123456789012345678901234567890123456789012 
C
C Get the unexplained elevation 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 6 Nov 2002.  CGH.  based on getvelobs.f 
C              The removing the data outside the mesh seems to be 
C              disabled.  Why? 

      subroutine getelevobs  

      include 'CASCO.DIM'
      include 'observes.h' 
      include 'mesh.h' 
      include 'startstop.h' 
      include 'misc.h' 
C
C Local arrays for sorting
C
      Integer indx(nelevobsdim)
      Real    temp(nelevobsdim)
C
C local variables
      real Tobsgmt0(nelevobsdim)
      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(*,*)'GETELEVOBS: '
      write(*,*)'Enter unexplained elevation file (*.uxe)'
      write(*,*)'format: t x y elev *'
      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 elevation observations provided.'
	  n_eobs   = 0
	  return
      endif 
      lun = 49 
      open(unit=lun, file=filename, status='old')

      call getobs_4c(lun,n_eobs, x_eobs, y_eobs, Tobsgmt0,
     &      eobs,nelevobsdim,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(*,*)'**GETELEVOBS: Fatal Error**'
	 write(*,*)'Mesh name from *uxe 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) n_eobs,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,N_EOBS
            CALL BASIS2D(NNDIM,NEDIM,IN,X,Y,NN,NE,X_eobs(I),Y_eobs(I),
     &                NCLOSEST,NELEM,BASISXY)
            IF(NELEM.NE.0)THEN
               NOBSIN=NOBSIN+1
               X_eobs(NOBSIN) = X_eobs(i)
               Y_eobs(NOBSIN) = Y_eobs(i)
               Tobsgmt0(NOBSIN) = Tobsgmt0(i)
               Eobs(NOBSIN) = Eobs(i)
            ENDIF
         ENDDO
         WRITE(*,'(I6,A,I6,3A)')NOBSIN,' of the original ',N_EOBS,
     &   ' O*D entries are in the domain.'
         N_EOBS=NOBSIN
         WRITE(*,*)NOBSIN,N_EOBS

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(n_eobs, tobsgmt0, indx)
C
      Call SortDRL(      x_eobs, n_eobs, indx, temp)
      Call SortDRL(      y_eobs, n_eobs, indx, temp)
      Call SortDRL(  tobsgmt0, n_eobs, indx, temp)
      Call SortDRL(      eobs, n_eobs, indx, temp)
C
C The sorting is complete.
C



C
C Convert UTC0 to Gregorian time 
         do i = 1, n_eobs
            CALL UTC0_GDAY(IYEAR,Tobsgmt0(i),kd_eobs(i),sec_eobs(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,n_eobs
c           write(*,*)'obs i ',kd_eobs(i),sec_eobs(i),Tobsgmt0(i), iyear
	    ncheck1 = igday_compare(kd0,sec0,kd_eobs(i),sec_eobs(i))
	    ncheck2 = igday_compare(kdl,secl,kd_eobs(i),sec_eobs(i))
	    keep = .false. 
	    if((ncheck1.eq.-1).and.(ncheck2.eq.1)) keep = .true.

	    if(keep.eqv..true.) then 
c           write(*,*)'keep ',kd_eobs(i),sec_eobs(i),Tobsgmt0(i), iyear
               NOBSIN=NOBSIN+1
               X_eobs(NOBSIN) = X_eobs(i)
               Y_eobs(NOBSIN) = Y_eobs(i)
               kd_eobs(NOBSIN)  = kd_eobs(i)
               sec_eobs(NOBSIN) = sec_eobs(i)
               eobs(NOBSIN)   = eobs(i)
            endif
         ENDDO
         WRITE(*,'(I6,A,I6,3A)')NOBSIN,' of the original ',N_EOBS,
     &   ' O*D entries are in the correct time window.'
         N_EOBS=NOBSIN
         WRITE(*,*)NOBSIN,N_EOBS
 	 if(n_eobs.lt.1) then 
            write(*,*)'**FATAL ERROR: in getelevobs**'
 	    write(*,*)'No elevation observations are in the correct',
     &      ' space/time window.  I will stop now.'
 	   stop
	 endif

      RETURN
      END

