C***********************************************************************
C***********************************************************************
      PROGRAM INTEGO3D
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
      INCLUDE 'TRUXTON6ARRAYS.DIM'
      PARAMETER(NO3DDIM=10*NNDIM,MDIM=6)
C
C O3D Data Arrays
      REAL ALON(NO3DDIM),ALAT(NO3DDIM),YDOBS(NO3DDIM),
     &UOBS(NO3DDIM),VOBS(NO3DDIM),
     &DMINOBS(NO3DDIM),DMAXOBS(NO3DDIM),DOBS(NO3DDIM)
      REAL XOBS(NO3DDIM),YOBS(NO3DDIM)
      REAL ZMAXOBS(NO3DDIM),ZMINOBS(NO3DDIM),ZBOTOBS(NO3DDIM)
      INTEGER NFLAG(NO3DDIM)
C
C FEM Mesh Arrays
      REAL X(NNDIM),Y(NNDIM),BATHY(NNDIM)
      INTEGER IN(3,NEDIM)
C
C Tidal Constituent Arrays
      REAL Z(NNDIM,NNVDIM)
      COMPLEX U(NNDIM,NNVDIM,MDIM),V(NNDIM,NNVDIM,MDIM),
     &W(NNDIM,NNVDIM,MDIM)
      COMPLEX UBAR(NNDIM,MDIM),VBAR(NNDIM,MDIM)
      REAL FREQ(MDIM)
      COMPLEX UBARK,VBARK,EYE
C
C   as required by IOS timing routines
      REAL F(MDIM),VU(MDIM),SIG(MDIM)
      CHARACTER*5 KON(MDIM),CDUM
      CHARACTER*72 KONFILE(MDIM)
C
C Interpolation Arrays
      REAL BAS(3)
      REAL UBAROBS(NO3DDIM),VBAROBS(NO3DDIM),DEPTH(NO3DDIM)
C
C Character Assignments
      CHARACTER*72 PATHMESH,FILEO3D,FILEOM3D,FILEMM3D,HEADER
      CHARACTER*1 ANS
C
C 1.Write user information
      WRITE(*,'(A)')'This program samples model output to compare'//
     &' with 3-D (t,long,lat) observational data sets'//
     &' (i.e. o3d observational file format).'//
     &'  The program is currently able to perform these comparisons'//
     &' for the following observational data sets:'
*      WRITE(*,'(A)')'  Vertically averaged horizontal currents'//
*     &' (i.e. v2c model output file format).'
      WRITE(*,'(A)')'  Vertically averaged horizontal currents'//
     &' (i.e. vel model output file format).'
      WRITE(*,'(A)')' '
      WRITE(*,'(A)')'This program currently uses SUBROUTINE MERCGOM2'//
     &' to perform the mercator projection'//
     &' (i.e. gb2 family of NML meshes).'
C
C 2.Settings:
C
C   constants
      EYE=CMPLX(0.0,1.0)
      TWOPI=2.0*ACOS(-1.0)
      DEGLAT=43.50
C   constituents
      WRITE(*,'(A)')' '
      WRITE(*,'(A)')' '
      WRITE(*,'(A)')'Data Entry:'
      WRITE(*,'(A)')'-----------'
*      WRITE(*,'(A)')'Enter the number of tidal constituents for'//
*     &' which you have v2c model output.'
      WRITE(*,'(A)')'Enter the number of tidal constituents for'//
     &' which you have vel model output.'
      READ(*,*)M
      WRITE(*,'(A,I2,A)')'Enter the constituent name (IOS convention)'//
     &' and the name of the corresponding .v2c file for the ',M,
     &' constituents (formatted read:ch*5,ch*72).'
      WRITE(*,'(A)')'KON..FILENAME......'
      DO K=1,M
         READ(*,'(A,A)')KON(K),KONFILE(K)
      ENDDO
C   filenames
      WRITE(*,'(A)')'Enter the path for your .nod, .ele, and .bat'//
     &' files (i.e. g2s => g2s.nod, g2s.ele, and g2s.bat):'
      READ(*,'(A)')PATHMESH
      WRITE(*,'(A)')'Enter the name for input .o3d file:'
      READ(*,'(A)')FILEO3D
      WRITE(*,'(A)')'Enter names for output .m3d files:'
      WRITE(*,'(A)')'1.Echo of .o3d file with (x,y) .vs. (long,lat):'
      READ(*,'(A)')FILEOM3D
      WRITE(*,'(A)')'2.Model .o3d file:'
      READ(*,'(A)')FILEMM3D
      WRITE(*,'(A)')'-----------'
C   open files
      OPEN(1,FILE=FILEO3D,STATUS='OLD')
      OPEN(2,FILE=FILEOM3D,STATUS='UNKNOWN')
      REWIND(2)
      OPEN(3,FILE=FILEMM3D,STATUS='UNKNOWN')
      REWIND(3)
C
C 2.Load FEM data:
C   Mesh
      WRITE(*,'(A)')' '
      WRITE(*,'(A)')' '
      WRITE(*,'(A)')'File Reading:'
      WRITE(*,'(A)')'------------'
      CALL LOADMESHGOM(NNDIM,NEDIM,PATHMESH,X,Y,IN,BATHY,NN,NE)
C   Vertically integrated currents
*      CALL LOADV2CM(NNDIM,NN,MDIM,M,KON,KONFILE,UBAR,VBAR,FREQ)
C   3-D currents
      CALL LOADVELM(NNDIM,NN,NNVDIM,NNV,MDIM,M,KON,KONFILE,Z,U,V,W,FREQ)
C
C 3.Load O3D data
C
C   read and discard header
      HEADER=' '
      DO WHILE(HEADER(1:4).NE.'XXXX')
         READ(1,'(A)')HEADER
         WRITE(2,'(A)')HEADER
         WRITE(3,'(A)')HEADER
      ENDDO
C
C   read comment line, year, and number of columns
      READ(1,'(A)')HEADER
      WRITE(2,'(A)')PATHMESH
      WRITE(3,'(A)')PATHMESH
      READ(1,'(A)')HEADER
      WRITE(2,'(A)')HEADER
      WRITE(3,'(100A)')(KON(K),K=1,M),' - ',HEADER
      READ(1,*)IYEAR
      WRITE(2,'(I4)')IYEAR
      WRITE(3,'(I4)')IYEAR
      READ(1,*)NCOL
      WRITE(2,'(I4)')8
      WRITE(3,'(I4)')8
C
C   now read data
      DO I=1,NO3DDIM
         READ(1,*,END=11)YDOBS(I),ALON(I),ALAT(I),UOBS(I),VOBS(I),
     &   DMINOBS(I),DMAXOBS(I),DOBS(I)
         NOBS=I
         CALL MERCGOM2(1,ALON(I),ALAT(I),XOBS(I),YOBS(I))
      ENDDO
 11   CLOSE(1)
      WRITE(*,*)NOBS,' O3D entries read from ',FILEO3D
      WRITE(*,'(A)')'------------'
C
C 4.Set-up for tidal synthesis
C
C   Use IOS files to load F(K) and VU(K) for the beginning of the 
C   day on which the cruise began:
      CALL YDAY0_GDAY(IYEAR,YDOBS(1),KD0,SECKD0)
      KH0=KD0*24
      CALL OPNVUF(KH0,CDUM,DEGLAT,RDUM,RDUM,RDUM)
      CLOSE(8)
      DO K=1,M
         CALL VUF(KH0,KON(K),DEGLAT,F(K),VU(K),SIG(K))
      ENDDO
C
C   Echo results to CRT
      WRITE(*,'(A)')' '
      WRITE(*,'(A)')' '
      WRITE(*,'(A)')'Parameters sent to IOS routines:'
      WRITE(*,'(A)')'--------------------------------'
      WRITE(*,'(A,I10)')  'Year:                             ',IYEAR
      WRITE(*,'(A,F10.4)')'Reference Time (yearday_gmt0):    ',YDOBS(1)
      WRITE(*,'(A,I10)')  'Reference Time (Gregorian Hours): ',KH0
      WRITE(*,'(A,F10.4)')'Degrees Latitude:                 ',DEGLAT
      WRITE(*,'(A)')'--------------------------------'
      WRITE(*,'(A)')' '
      WRITE(*,'(A)')' '
      WRITE(*,'(A)')'Tidal parameters loaded from IOS routines:'
      WRITE(*,'(A)')'Constituent      F          VU          SIG'
      WRITE(*,'(A)')'----------- ---------- -----------  ----------'
      DO K=1,M
         WRITE(*,'(4X,A,1X,3e12.4)')KON(K),F(K),VU(K),SIG(K)
      ENDDO
      WRITE(*,'(A)')'----------- ---------- -----------  ----------'
C
C   Allow user to overwrite F(K)=1 and SIG(K)=FREQ(K), in case they
C   are hoping to perfectly invert model data:
      WRITE(*,'(A)')' '
      WRITE(*,'(A)')' '
      WRITE(*,'(A)')'Do you desire to overwrite F(K)=1 and SIG(K) values
     & with values from your .v2c files?'
      READ(*,'(A)')ANS
      IF(ANS.EQ.'Y'.OR.ANS.EQ.'y')THEN
         DO K=1,M
            F(K)=1.0
            SIG(K)=FREQ(K)/2.0/ACOS(-1.0)*3600.0
         ENDDO
         WRITE(*,'(A)')' '
         WRITE(*,'(A)')' '
         WRITE(*,'(A)')'Revised tidal parameters:'
         WRITE(*,'(A)')'Constituent      F          VU          SIG'
         WRITE(*,'(A)')'----------- ---------- -----------  ----------'
         DO K=1,M
            WRITE(*,'(4X,A,1X,3e12.4)')KON(K),F(K),VU(K),SIG(K)
         ENDDO
         WRITE(*,'(A)')'----------- ---------- -----------  ----------'
      ELSE
         WRITE(*,'(A)')'Tidal parameters not being revised.'
      ENDIF
C
C 5.Interpolate FEM data to O3D locations
      WRITE(*,'(A)')' '
      WRITE(*,'(A)')' '
      write(*,'(A)')'Interpolating O3D locations:'
      write(*,'(A)')'----------------------------'
      NOBSIN=0
      DO 20 I=1,NOBS
C
C   Determine basis functions at O3D location
         IF(MOD(I,100).EQ.0)WRITE(*,*)I
         CALL BASIS2D(NNDIM,NEDIM,IN,X,Y,NN,NE,XOBS(I),YOBS(I),
     &                NCLOSEST,NELEM,BAS)
C
C   For observations outside of the domain, perform no interpolations
         IF(NELEM.EQ.0)THEN
            NFLAG(I)=0
            GO TO 20
         ENDIF
C
C   For observations inside the domain, set flag, counters, time since
C   the beginning of the day on which the cruise began, and
C   loop over frequencies to reconstruct time-domain signal
         NFLAG(I)=1
         NOBSIN=NOBSIN+1
         CALL YDAY0_GDAY(IYEAR,YDOBS(I),KD,SECKD)
         TUVSET=(KD-KD0)*24.0*3600.0+SECKD
         N1=IN(1,NELEM)
         N2=IN(2,NELEM)
         N3=IN(3,NELEM)
*         CALL INTERPV2CM(NNDIM,MDIM,UBAR,VBAR,M,F,VU,SIG,TUVSET,
*     &   N1,N2,N3,BAS,UBAROBS(I),VBAROBS(I))
         CALL INTERPVELM(MDIM,NNV,Z,U,V,W,M,F,VU,SIG,
     &   DMINOBS(I),DMAXOBS(I),DOBS(I),TUVSET,N1,N2,N3,BAS,
     &   UBAROBS(I),VBAROBS(I),ZMAXOBS(I),ZMINOBS(I),ZBOTOBS(I))
 20   CONTINUE
      write(*,'(A)')' '
      CALL CSTRING(PATHMESH,IS,IE)
      WRITE(*,'(I6,A,I6,3A)')NOBSIN,' of the original ',NOBS,
     &' O3D entries are in the ',PATHMESH(IS:IE),' domain.'
      write(*,'(A)')'----------------------------'
C
C Output Results
      write(*,'(A)')' '
      write(*,'(A)')' '
      write(*,'(A)')'Writing Output Files:'
      write(*,'(A)')'---------------------'
      write(*,'(A)')FILEOM3D
      write(*,'(A)')FILEMM3D
      DO I=1,NOBS
         IF(NFLAG(I).EQ.1)THEN
            WRITE(2,'(F10.4,2F10.0,2F10.4,3F10.2)')
     &      YDOBS(I),XOBS(I),YOBS(I),UOBS(I),VOBS(I),
     &      ZMAXOBS(I),ZMINOBS(I),ZBOTOBS(I)
*            WRITE(3,'(F10.4,2F10.0,2E12.4,I5)')
CDJM        WRITE(3,*)
             WRITE(3,'(F10.4,2F10.0,2E12.4,3f8.2)')
     &      YDOBS(I),XOBS(I),YOBS(I),UBAROBS(I),VBAROBS(I),
     &      ZMAXOBS(I),ZMINOBS(I),ZBOTOBS(I)
         ENDIF
      ENDDO
      write(*,'(A)')'---------------------'
      CLOSE(1)
      CLOSE(2)
      STOP
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE YDAY0_GDAY(IYEAR,YDAY,KD,SECKD)
C-----------------------------------------------------------------------
      CALL GDAY(01,01,IYEAR,KDNEWYEAR)
      KD=KDNEWYEAR+INT(AINT(YDAY))
      SECKD=24.0*3600.0*(YDAY-AINT(YDAY))
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE YDAY1_GDAY(IYEAR,YDAY,KD,SECKD)
C-----------------------------------------------------------------------
      CALL GDAY(01,01,IYEAR,KDNEWYEAR)
      KD=KDNEWYEAR+INT(AINT(YDAY))-1
      SECKD=24.0*3600.0*(YDAY-AINT(YDAY))
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE LOADV2CM(NNDIM,NN,MDIM,M,KON,KONFILE,U,V,FREQ)
C-----------------------------------------------------------------------
C purpose: This routine loads multiple .v2c files
C
C index range: i=1,nn
C
C inputs:  nndim - parameter used to dimension node arrays in the
C            calling routine
C          nn - number of nodes for which data should exist
C          filename - name of file to be read
C
C outputs: u(i,mdim) - x-component of the complex vector at node i 
C          v(i,mdim) - y-component of the complex vector at node i
C          freq(mdim) - frequency in radians/sec
C 
C file format: See CEN 3/29/93 Memorandum
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           1992 JUL 30 - LOADV2C
C           1994 APR 25 - Revised LOADV2C
C           1997 JAN 20 - LOADV2CM
C-----------------------------------------------------------------------
C
C ARGUMENTS
      complex u(nndim,mdim),v(nndim,mdim)
      real freq(mdim)
      character*5 kon(mdim)
      character*72 filename,header,konfile(mdim)
C
C START OF EXECUTABLE CODE
C
C Check for oversize of arrays
      nflag=0
      if(nn.GT.nndim)then
	 write(*,*)'Dimensioning problem in SUBROUTINE LOADV2C:'
         write(*,*)'  nn exceeds nndim => rectify and recompile'
         nflag=1
      endif
      if(nflag.EQ.1)stop
C
C Open file
      do k=1,m
         filename=konfile(k)
         open(63,file=filename,status='old')
C
C Read headers
         read(63,'(A)')header
         read(63,'(A)')header
         read(63,*)freq(k)
C
C Read data
         fac=2.0*ATAN2(1.0,0.0)/180.0
         do i=1,nndim
            read(63,*,end=11)node,uamp,upha,vamp,vpha
            u(i,k)=CMPLX(uamp*COS(-upha*fac),uamp*SIN(-upha*fac))
            v(i,k)=CMPLX(vamp*COS(-vpha*fac),vamp*SIN(-vpha*fac))
            ii=i
         enddo
 11      close(63)
C
C Message to standard out
         if(ii.NE.nn)then
         write(*,*)'There is a mismatch between the number of horizontal
     & nodes in your .v2c file and nn in SUBROUTINE LOADV2C.'
         STOP
         else
         write(*,*)ii,' complex 2-D vector values loaded from ',filename
         endif
      enddo
C
C END OF ROUTINE
      return
      end
C***********************************************************************
C***********************************************************************
      SUBROUTINE INTERPV2CM(NNDIM,MDIM,UBAR,VBAR,M,F,VU,SIG,TUVSET,
     &N1,N2,N3,BAS,UBARTOT,VBARTOT)
C-----------------------------------------------------------------------
C
C Model Tidal Constituent Arrays
      COMPLEX UBAR(NNDIM,MDIM),VBAR(NNDIM,MDIM)
      REAL F(MDIM),VU(MDIM),SIG(MDIM)
      COMPLEX UBARK,VBARK,EYE
C
C Interpolation Arrays:  BAS(1)@N1,BAS(2)@N2,BAS(3)@N3,
      REAL BAS(3)
C
C Assignments
      EYE=CMPLX(0.0,1.0)
      TWOPI=2.0*ACOS(-1.0)
C
C Perform Interpolation
      UBARTOT=0.0
      VBARTOT=0.0
      DO K=1,M
         UBARK=+BAS(1)*UBAR(N1,K)+BAS(2)*UBAR(N2,K)+BAS(3)*UBAR(N3,K)
         VBARK=+BAS(1)*VBAR(N1,K)+BAS(2)*VBAR(N2,K)+BAS(3)*VBAR(N3,K)
         ANGLEK=TWOPI*(SIG(K)*TUVSET/3600.0+VU(K))
         UBARTOT=UBARTOT+F(K)*REAL(UBARK*CEXP(EYE*ANGLEK))
         VBARTOT=VBARTOT+F(K)*REAL(VBARK*CEXP(EYE*ANGLEK))
      ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE LOADVELM(NNDIM,NN,NNVDIM,NNV,MDIM,M,KON,KONFILE,
     &Z,U,V,W,FREQ)
C-----------------------------------------------------------------------
C purpose: This routine loads multiple .vel files
C
C index range: i=1,nn,j=1,nnv,k=1,m
C
C inputs:  nndim - parameter used to dimension node arrays in the
C            calling routine
C          nn - number of nodes for which data should exist
C          filename - name of file to be read
C
C outputs: u(i,j,k) - x-component of the complex velocity
C          v(i,j,k) - y-component of the complex velocity
C          w(i,j,k) - z-component of the complex velocity
C          freq(k) - frequency in radians/sec
C 
C file format: See CEN 3/29/93 Memorandum
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           1992 JUL 30 - LOADV2C
C           1994 APR 25 - Revised LOADV2C
C           1997 JAN 20 - LOADV2CM
C           1997 MAR 17 - LOADVELM
C-----------------------------------------------------------------------
C
C ARGUMENTS
      real z(nndim,nnvdim)
      complex u(nndim,nnvdim,mdim),v(nndim,nnvdim,mdim),
     &w(nndim,nnvdim,mdim)
      real freq(mdim)
      character*5 kon(mdim)
      character*72 filename,header,konfile(mdim)
      complex eye
      eye=cmplx(0.0,1.0)
C
C START OF EXECUTABLE CODE
C
C Check for oversize of arrays
      nflag=0
      if(nn.GT.nndim)then
	 write(*,*)'Dimensioning problem in SUBROUTINE LOADV2C:'
         write(*,*)'  nn exceeds nndim => rectify and recompile'
         nflag=1
      endif
      if(nflag.EQ.1)stop
C
C Open file
      do k=1,m
         filename=konfile(k)
         open(63,file=filename,status='old')
C
C Read headers
         read(63,'(A)')header
         read(63,'(A)')header
         read(63,*)nnv
         if(k.EQ.1)then
            nnv1=nnv
         else
            if(nnv.NE.nnv1)then
               write(*,*)'nnv.NE.nnv1:',nnv,nnv1
               stop
            endif
         endif
         read(63,*)nfreq
         read(63,*)freq(k)
C
C Read data
         do i=1,nndim
            do j=1,nnv
               read(63,*,end=11)node,zval,ua,up,va,vp,wa,wp
               if(k.EQ.1)z(i,j)=zval
               u(i,j,k)=ua*CEXP(-EYE*up)
               v(i,j,k)=va*CEXP(-EYE*vp)
               w(i,j,k)=wa*CEXP(-EYE*wp)
               ii=i
            enddo
         enddo
 11      close(63)
C
C Message to standard out
         if(ii.NE.nn)then
         write(*,*)'There is a mismatch between the number of horizontal
     & nodes in your .vel file and nn in SUBROUTINE LOADVELM.'
         STOP
         else
         write(*,*)ii,'x',nnv,' complex 2-D vector values loaded from'
     &             ,filename
         endif
      enddo
C
C END OF ROUTINE
      return
      end
C***********************************************************************
C***********************************************************************
      SUBROUTINE INTERPVELM(MDIM,NNV,Z,U,V,W,M,F,VU,SIG,DTOP,DBOT,DTOT,
     &TUVSET,N1,N2,N3,BAS,UBAR,VBAR,ZTOP,ZBOT,ZTOT)
C-----------------------------------------------------------------------
      INCLUDE 'TRUXTON6ARRAYS.DIM'
C
C Model Tidal Constituent Arrays
      REAL Z(NNDIM,NNVDIM)
      COMPLEX U(NNDIM,NNVDIM,MDIM),V(NNDIM,NNVDIM,MDIM),
     &W(NNDIM,NNVDIM,MDIM)
      REAL F(MDIM),VU(MDIM),SIG(MDIM)
      COMPLEX UK,VK,WK,EYE
C
C Interpolation Arrays:  BAS(1)@N1,BAS(2)@N2,BAS(3)@N3,
      REAL BAS(3)
C
C Local Assignments
      REAL ZOBS(NNVDIM),UOBS(NNVDIM),VOBS(NNVDIM),WOBS(NNVDIM)
      EYE=CMPLX(0.0,1.0)
      TWOPI=2.0*ACOS(-1.0)
C
C Interpolation 3-D variables to observed locations
      DO J=1,NNV
         ZOBS(J)=+BAS(1)*Z(N1,J)+BAS(2)*Z(N2,J)+BAS(3)*Z(N3,J)
         UJ=0.0
         VJ=0.0
         WJ=0.0
         DO K=1,M
            UK=+BAS(1)*U(N1,J,K)+BAS(2)*U(N2,J,K)+BAS(3)*U(N3,J,K)
            VK=+BAS(1)*V(N1,J,K)+BAS(2)*V(N2,J,K)+BAS(3)*V(N3,J,K)
            WK=+BAS(1)*W(N1,J,K)+BAS(2)*W(N2,J,K)+BAS(3)*W(N3,J,K)
            ANGLEK=TWOPI*(SIG(K)*TUVSET/3600.0+VU(K))
            UJ=UJ+F(K)*REAL(UK*CEXP(EYE*ANGLEK))
            VJ=VJ+F(K)*REAL(VK*CEXP(EYE*ANGLEK))
            WJ=WJ+F(K)*REAL(WK*CEXP(EYE*ANGLEK))
         ENDDO
         UOBS(J)=UJ
         VOBS(J)=VJ
         WOBS(J)=WJ
      ENDDO
C
C Map extremes of observed water column to model vertical grid
      HTOT=ZOBS(NNV)-ZOBS(1)
      ZTOT=-HTOT
      IF(DTOP.LE.0.0)THEN
         WRITE(*,*)'DTOP.LE.0:',DTOP
         WRITE(*,*)'SETTING DTOP=0.0'
         DTOP=0.0
      ENDIF
      IF(DBOT.LT.DTOP)THEN
         WRITE(*,*)'DBOT.LT.DTOP:',DBOT,DTOP
         WRITE(*,*)'SETTING DBOT=DTOP'
         DBOT=DTOP
      ENDIF
      IF(DBOT.GE.DTOT)THEN
         WRITE(*,*)'DBOT.GE.DTOT:',DBOT,DTOT
         WRITE(*,*)'SETTING DBOT=DTOT'
         DBOT=DTOT
      ENDIF
      ZBOT=ZOBS(NNV)-DBOT/DTOT*HTOT
      ZTOP=ZOBS(NNV)-DTOP/DTOT*HTOT
      CALL BASIS1D1(NNV,ZOBS,ZBOT,BV2,NV2,BV1,NV1)
      UBOT=UOBS(NV1)*BV1+UOBS(NV2)*BV2
      VBOT=VOBS(NV1)*BV1+VOBS(NV2)*BV2
      NBOTP=NV2
      CALL BASIS1D1(NNV,ZOBS,ZTOP,BV2,NV2,BV1,NV1)
      UTOP=UOBS(NV1)*BV1+UOBS(NV2)*BV2
      VTOP=VOBS(NV1)*BV1+VOBS(NV2)*BV2
      NTOPM=NV1
C
C Now vertically integrate
C  Initialize sums with bottom interval
      UINT=0.5*(UBOT+UOBS(NBOTP))*(ZOBS(NBOTP)-ZBOT)
      VINT=0.5*(VBOT+VOBS(NBOTP))*(ZOBS(NBOTP)-ZBOT)
C  Add intermediate levels
      DO J=NBOTP,NTOPM-1
         UINT=UINT+0.5*(UOBS(J)+UOBS(J+1))*(ZOBS(J+1)-ZOBS(J))
         VINT=VINT+0.5*(VOBS(J)+VOBS(J+1))*(ZOBS(J+1)-ZOBS(J))
      ENDDO
C  Add top interval      
      UINT=UINT+0.5*(UOBS(NTOPM)+UTOP)*(ZTOP-ZOBS(NTOPM))
      VINT=VINT+0.5*(VOBS(NTOPM)+VTOP)*(ZTOP-ZOBS(NTOPM))
C
C Now divide by range of depth to determine vertical average
      HOBS=ZTOP-ZBOT
      UBAR=UINT/HOBS
      VBAR=VINT/HOBS
      RETURN
      END
C***********************************************************************
C***********************************************************************
