C***********************************************************************
C***********************************************************************
      PROGRAM INTEGM3D
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
      INCLUDE 'FUNDY5xARRAYS.DIM'
      PARAMETER(NCOLDIM=13)
C
C Argument Declaration
C
C   FEM Mesh Arrays
      REAL X(NNDIM),Y(NNDIM),BATHY(NNDIM)
      REAL ZI(NNVDIM),UI(NNVDIM),VI(NNVDIM)
      INTEGER IN(3,NEDIM)
      CHARACTER*72 PATHMESH
C
C   Observational data
      REAL AOBS1(NOBSDIM,NCOLDIM),AOBS2(NOBSDIM,NCOLDIM),
     &AOBS3(NOBSDIM,NCOLDIM)
      REAL UBAR1(NOBSDIM),VBAR1(NOBSDIM)
      CHARACTER*72 OBSIN1,OBSIN2,OBSOUT1,OBSOUTE
C
C Enter filenames
      WRITE(*,'(A)')'Enter name for input 21 vertical node .m3d file:'
      READ(*,'(A)')OBSIN1
      WRITE(*,'(A)')'Enter name for input partial integral .m3d file:'
      READ(*,'(A)')OBSIN2
      WRITE(*,'(A)')'Enter name for output partial integral .m3d file:'
      READ(*,'(A)')OBSOUT1
      WRITE(*,'(A)')'Enter the name output error .m3d file:'
      READ(*,'(A)')OBSOUTE

C
C Load observational data
      CALL READOBS(NOBSDIM,NCOLDIM,OBSIN1,IYEAR1,NOBS1,NCOL1,AOBS1)
      CALL READOBS(NOBSDIM,NCOLDIM,OBSIN2,IYEAR2,NOBS2,NCOL2,AOBS2)
C
C
      NNV=21
      NOBS3=0
      NTMAX=0
      NBMIN=NNV
      NCOL3=8
      DO I1=1,NOBS1,NNV
         DO J=1,NNV
            ZI(J)=AOBS1(I1+J-1,4)
            UI(J)=AOBS1(I1+J-1,5)
            VI(J)=AOBS1(I1+J-1,6)
         ENDDO
         YM3D=AOBS1(I1,1)
         I2=0
 10      I2=I2+1
         YO3D=AOBS2(I2,1)
         IF(ABS(YO3D-YM3D).LT.0.0001)THEN
            NOBS3=NOBS3+1
C
            UBAR1(NOBS3)=AOBS2(I2,4)
            VBAR1(NOBS3)=AOBS2(I2,5)
            HO3D=AOBS2(I2,8)
            HM3D=ZI(NNV)-ZI(1)
            ZTOP=ZI(NNV)-AOBS2(I2,6)/HO3D*HM3D
            ZBOT=ZI(NNV)-AOBS2(I2,7)/HO3D*HM3D
            ZBAT=ZI(1)
C
            UINT=0.0
            VINT=0.0
            GO TO 15
C Total Integral
 14         DO J=2,NNV
               DZ=ZI(J)-ZI(J-1)
               UINT=UINT+0.5*(UI(J-1)+UI(J))*DZ
               VINT=VINT+0.5*(VI(J-1)+VI(J))*DZ
            ENDDO
            UBAR=UINT/HM3D
            VBAR=VINT/HM3D
*            GO TO 16
C Partial Integral
 15         CALL BASIS1D1(NNV,ZI,ZBOT,B2B,N2B,B1B,N1B)
            CALL BASIS1D1(NNV,ZI,ZTOP,B2T,N2T,B1T,N1T)
C  Correct for BASIS1D1 structure
            IF(N2B.EQ.1)THEN
               N2B=2
               B2B=0.0
               N1B=1
               B1B=1.0
            ENDIF
            IF(N1T.EQ.NNV)THEN
               N2T=NNV
               B2T=1.0
               N1T=NNV-1
               B1T=0.0
            ENDIF
            NTMAX=MAX(NTMAX,N2T)
            NBMIN=MIN(NBMIN,N1B)
C  Now perform vertical integration
            UIB=UI(N2B)*B2B+UI(N1B)*B1B
            VIB=VI(N2B)*B2B+VI(N1B)*B1B
            DZ=ZI(N2B)-ZBOT
            UINT=0.5*(UIB+UI(N2B))*DZ
            VINT=0.5*(VIB+VI(N2B))*DZ
            DO J=N2B+1,N1T
               DZ=ZI(J)-ZI(J-1)
               UINT=UINT+0.5*(UI(J-1)+UI(J))*DZ
               VINT=VINT+0.5*(VI(J-1)+VI(J))*DZ
            ENDDO
            UIT=UI(N2T)*B2T+UI(N1T)*B1T
            VIT=VI(N2T)*B2T+VI(N1T)*B1T
            DZ=ZTOP-ZI(N1T)
            UINT=UINT+0.5*(UI(N1T)+UIT)*DZ
            VINT=VINT+0.5*(VI(N1T)+VIT)*DZ
            UBAR=UINT/(ZTOP-ZBOT)
            VBAR=VINT/(ZTOP-ZBOT)
C Write Results to AOBS3 Array
 16         AOBS3(NOBS3,1)=YM3D
            AOBS3(NOBS3,2)=AOBS1(I1,2)
            AOBS3(NOBS3,3)=AOBS1(I1,3)
            AOBS3(NOBS3,4)=UBAR
            AOBS3(NOBS3,5)=VBAR
            AOBS3(NOBS3,6)=ZTOP
            AOBS3(NOBS3,7)=ZBOT
            AOBS3(NOBS3,8)=ZBAT
         ELSE IF(YO3D.LT.YM3D)THEN
            GO TO 10
         ELSE
         ENDIF
      ENDDO
      WRITE(*,*)'Number of observations which survived:',NOBS3
      WRITE(*,*)'Highest vertical node used:           ',NTMAX
      WRITE(*,*)'Lowest vertical node used:            ',NBMIN
C
C Write observational data
 20   CALL WRITEOBS(NOBSDIM,NCOLDIM,OBSOUT1,IYEAR2,NOBS3,NCOL3,AOBS3)
      DO I3=1,NOBS3
         AOBS3(I3,4)=UBAR1(I3)-AOBS3(I3,4)
         AOBS3(I3,5)=VBAR1(I3)-AOBS3(I3,5)
      ENDDO   
      CALL WRITEOBS(NOBSDIM,NCOLDIM,OBSOUTE,IYEAR2,NOBS3,NCOL3,AOBS3)
      STOP
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE READOBS(NOBSDIM,NCOLDIM,OBSFILE,IYEAR,NOBS,NCOL,AOBS)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C Argument Declaration
      REAL AOBS(NOBSDIM,NCOLDIM)
      CHARACTER*72 OBSFILE,HEADER
C
C Open file
      OPEN(1,FILE=OBSFILE,STATUS='OLD')
C
C Read and discard header
      HEADER=' '
      DO WHILE(HEADER(1:4).NE.'XXXX')
         READ(1,'(A)')HEADER
      ENDDO
C
C Read comment lines, year, and number of columns
      READ(1,'(A)')HEADER
      READ(1,'(A)')HEADER
      READ(1,*)IYEAR
      READ(1,*)NCOL
      IF(NCOL.GT.NCOLDIM)THEN
         WRITE(*,*)'Number of columns exceeds NCOLDIM'
         WRITE(*,*)'NCOL=',NCOL,'             NCOLDIM=',NCOLDIM
         STOP
      ENDIF
C
C Read data
      DO I=1,NOBSDIM
         READ(1,*,END=11)(AOBS(I,J),J=1,NCOL)
         NOBS=I
      ENDDO
      WRITE(*,*)'Error reading observational data'
      STOP
 11   CLOSE(1)
      WRITE(*,*)NOBS,' x',NCOL,' O*D entries read from ',OBSFILE
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE WRITEOBS(NOBSDIM,NCOLDIM,OBSFILE,IYEAR,NOBS,NCOL,AOBS)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C Argument Declaration
      REAL AOBS(NOBSDIM,NCOLDIM)
      CHARACTER*72 OBSFILE,HEADER
C
C Open file
      OPEN(1,FILE=OBSFILE,STATUS='UNKNOWN')
      REWIND(1)
C
C Write Header Lines
      WRITE(1,'(A)')'XXXX'
      WRITE(1,'(A)')'comment'
      WRITE(1,'(A)')'comment'
      WRITE(1,'(I4)')IYEAR
      WRITE(1,'(I4)')NCOL
C
C Write data
      DO I=1,NOBS
         WRITE(1,*)(AOBS(I,J),J=1,NCOL)
      ENDDO
      CLOSE(1)
      WRITE(*,*)NOBS,' x',NCOL,' O*D entries written to ',OBSFILE
      RETURN
      END
C***********************************************************************
C***********************************************************************
