C***********************************************************************
      SUBROUTINE NUMERICAL_UV_INTEG(KD,SECKD,ITER,DELT,IN,XNOD,YNOD,
     &BATHY,NN,NNV,NE,UBAR,VBAR,Z,U,V,W,RHO,T,S,Q2,Q2L,Cd,ENZM,MESH,
     &OBSFILE,nf1,nf2)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
C
C Argument Declaration
C
C   FEM Arrays
      INTEGER IN(3,NEDIM)
      REAL XNOD(NNDIM),YNOD(NNDIM),UBAR(NNDIM),VBAR(NNDIM)
      REAL Z(NNDIM,NNVDIM),U(NNDIM,NNVDIM),V(NNDIM,NNVDIM),
     &W(NNDIM,NNVDIM),RHO(NNDIM,NNVDIM),T(NNDIM,NNVDIM),S(NNDIM,NNVDIM),
     &Q2(NNDIM,NNVDIM),Q2L(NNDIM,NNVDIM)
      CHARACTER*72 MESH,OBSFILE
C
C   Observational data
      REAL AOBS(NOBSDIM,NCOLDIM)
      INTEGER KDOBS(NOBSDIM)
      REAL SECOBS(NOBSDIM)
      REAL XOBS(NOBSDIM),YOBS(NOBSDIM)
      REAL UADCP(NOBSDIM),VADCP(NOBSDIM)
      REAL UOBS(NOBSDIM),VOBS(NOBSDIM)
      REAL ZMAXOBS(NOBSDIM),ZMINOBS(NOBSDIM),ZBOTTOMOBS(NOBSDIM)
C
C   Dependent variables and horizontal basis functions for
C   observation locations
      REAL ZZOBS(NOBSDIM,NNVDIM),
     &UZOBS(NOBSDIM,NNVDIM),VZOBS(NOBSDIM,NNVDIM),WZOBS(NOBSDIM,NNVDIM),
     &ROBS(NOBSDIM,NNVDIM),TOBS(NOBSDIM,NNVDIM),SOBS(NOBSDIM,NNVDIM)
      REAL BASISOBS(5,NOBSDIM)
      REAL BXY(3)
      CHARACTER*72 MODFILE
C
C Save necessary arguments
      SAVE XOBS,YOBS,UADCP,VADCP
      SAVE ZZOBS,UZOBS,VZOBS,WZOBS,BASISOBS,IOBSDELT,IOBS,NOBSDONE

C
C When ITER=0
      IF(ITER.EQ.0)THEN
C
C   Read Observations File
         CALL CROPSTRING(OBSFILE,IS,IB,IE)
         WRITE(*,*)'SUBROUTINE NUMERICAL_UV_INTEG: ',OBSFILE(IS:IE)
         CALL READOBS(NOBSDIM,NCOLDIM,OBSFILE,IYEAR,NOBS,NCOL,AOBS)
         DO I=1,NOBS
            XOBS(I)=AOBS(I,2)
            YOBS(I)=AOBS(I,3)
         ENDDO

C   Discard observations which lie outside the boundary of the finite
C   element mesh and determine the following for observations which lie
C   within the boundary of the finite element mesh:
C     -Gregrorian time
C     -horizontal basis functions
C     -Set XOBS,YOBS,ZMINOBS,ZMAXOBS,ZBOTTOMOBS
         NOBSIN=0
         DO I=1,NOBS
            CALL BASIS2D(NNDIM,NEDIM,IN,XNOD,YNOD,NN,NE,AOBS(I,2),
     &                AOBS(I,3),NCLOSEST,NELEM,BXY)
	    IF (NELEM.EQ.0) THEN
	      WRITE(2,*)I,AOBS(I,2),AOBS(I,3)
	    ENDIF
            IF(NELEM.NE.0)THEN
               NOBSIN=NOBSIN+1
               CALL GMT0_GDAY(IYEAR,AOBS(I,1),
     &            KDOBS(NOBSIN),SECOBS(NOBSIN))
               BASISOBS(1,NOBSIN)=NELEM
               BASISOBS(2,NOBSIN)=BXY(1)
               BASISOBS(3,NOBSIN)=BXY(2)
               BASISOBS(4,NOBSIN)=BXY(3)
               XOBS(NOBSIN)=AOBS(I,2)
               YOBS(NOBSIN)=AOBS(I,3)
               UADCP(NOBSIN)=AOBS(I,4)
               VADCP(NOBSIN)=AOBS(I,5)
               ZMINOBS(NOBSIN)=AOBS(I,6)
               ZMAXOBS(NOBSIN)=AOBS(I,7)
               ZBOTTOMOBS(NOBSIN)=AOBS(I,8)
            ENDIF
         ENDDO
         WRITE(*,'(1X,I6,A,I6,3A)')NOBSIN,' of the original ',NOBS,
     &   ' O*D entries are in the domain.'
         NOBS=NOBSIN
C
C   Initialize IOBS and NOBSDONE for subsequent calls
         IOBS=1
         NOBSDONE=0
C
C   Open Results Files and Write Headers
         MODFILE=OBSFILE(IS:IE-4)//'.Q5.m3d'
         OPEN(nf1,FILE=MODFILE,STATUS='UNKNOWN')
         REWIND(nf1)
         WRITE(nf1,'(A)')'XXXX'
         WRITE(nf1,'(A)')MESH
         WRITE(nf1,'(A)')'Q5 results at locs in '//OBSFILE(IB:IE)
         WRITE(nf1,'(I4)')IYEAR
         WRITE(nf1,'(I4)')10
         MODFILE=OBSFILE(IS:IE-4)//'-Q5.m3d'
         OPEN(nf2,FILE=MODFILE,STATUS='UNKNOWN')
         REWIND(nf2)
         WRITE(nf2,'(A)')'XXXX'
         WRITE(nf2,'(A)')MESH
         WRITE(nf2,'(A)')'Q5 results at locs in '//OBSFILE(IB:IE)
         WRITE(nf2,'(I4)')IYEAR
         WRITE(nf2,'(I4)')8
      ENDIF

C
C For all ITER
C
C   Return if all observations has been passed
      IF(NOBSDONE.EQ.NOBS)RETURN
C
C   Set values if the model timestep just completed spans the
C   the observation time
      DO II=IOBS-IOBSDELT,IOBS-1
         CALL GDAY_GMT0(KDOBS(II),SECOBS(II),IYEAR,YDAYGMT0)
         NELEM=INT(BASISOBS(1,II))
         BXY(1)=BASISOBS(2,II)
         BXY(2)=BASISOBS(3,II)
         BXY(3)=BASISOBS(4,II)
         DO J=1,NNV
            CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,Z,ZB2J)
            CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,U,UB2J)
            CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,V,VB2J)
            CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,W,WB2J)
            CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,RHO,RB2J)
            CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,T,TB2J)
            CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,S,SB2J)
            ZZOBS(II,J)=ZB2J-(ZB2J-ZZOBS(II,J))*BASISOBS(5,II)
            UZOBS(II,J)=UB2J-(UB2J-UZOBS(II,J))*BASISOBS(5,II)
            VZOBS(II,J)=VB2J-(VB2J-VZOBS(II,J))*BASISOBS(5,II)
            WZOBS(II,J)=WB2J-(WB2J-WZOBS(II,J))*BASISOBS(5,II)
            ROBS(II,J)=RB2J-(RB2J-ROBS(II,J))*BASISOBS(5,II)
            TOBS(II,J)=TB2J-(TB2J-TOBS(II,J))*BASISOBS(5,II)
            SOBS(II,J)=SB2J-(SB2J-SOBS(II,J))*BASISOBS(5,II)
            WRITE(nf1,'(F10.4,2F10.0,7E14.6)')YDAYGMT0,XOBS(II),
     &      YOBS(II),ZZOBS(II,J),UZOBS(II,J),VZOBS(II,J),WZOBS(II,J),
     &      ROBS(II,J),TOBS(II,J),SOBS(II,J)
         ENDDO
         CALL UV_VERT_INTEG(II,NNV,ZZOBS,UZOBS,VZOBS,
     &   ZMAXOBS,ZMINOBS,ZBOTTOMOBS,UOBS,VOBS,ZZMIN,ZZMAX)
         WRITE(nf2,'(F10.4,2F10.0,2E14.6,3f10.4)')YDAYGMT0,
     &   XOBS(II),YOBS(II),UADCP(II)-UOBS(II),VADCP(II)-VOBS(II),
     &   ZZMAX,ZZMIN,ZZOBS(II,1)
      ENDDO
      NOBSDONE=IOBS-1
C
C   Return if the next timestep does not span the temporal location
C   of the next observational location
      IOBSDELT=0
 2    DTSEC=(KDOBS(IOBS)-KD)*24.0*3600.0+SECOBS(IOBS)-SECKD
      IF(DTSEC.LT.0.0)RETURN
      IF(DTSEC.GE.DELT)RETURN
C
C   Evaluate temporal basis function at this time level for this
C   observation, perform horizontal interpolation, and increment
C   counter
      BASISOBS(5,IOBS)=1.0-DTSEC/DELT
      NELEM=INT(BASISOBS(1,IOBS))
      BXY(1)=BASISOBS(2,IOBS)
      BXY(2)=BASISOBS(3,IOBS)
      BXY(3)=BASISOBS(4,IOBS)
      DO J=1,NNV
       CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,Z,ZZOBS(IOBS,J))
       CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,U,UZOBS(IOBS,J))
       CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,V,VZOBS(IOBS,J))
       CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,W,WZOBS(IOBS,J))
       CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,RHO,ROBS(IOBS,J))
       CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,T,TOBS(IOBS,J))
       CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,S,SOBS(IOBS,J))
      ENDDO
      IOBSDELT=IOBSDELT+1
      IOBS=IOBS+1
      GO TO 2
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE UV_VERT_INTEG(II,NNV,ZZOBS,UZOBS,VZOBS,
     &ZMAXOBS,ZMINOBS,ZBOTTOMOBS,UOBS,VOBS,ZZMIN,ZZMAX)
C-----------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
C
C   Observational data
      REAL ZMAXOBS(NOBSDIM),ZMINOBS(NOBSDIM),ZBOTTOMOBS(NOBSDIM)
C
C Dependent variables for observation locations
      REAL ZZOBS(NOBSDIM,NNVDIM),
     &UZOBS(NOBSDIM,NNVDIM),VZOBS(NOBSDIM,NNVDIM)
      REAL UOBS(NOBSDIM),VOBS(NOBSDIM)
      REAL ZZOBSII(NNVDIM)
C
C Create 1-D array of vertical node locations
      DO J=1,NNV
         ZZOBSII(J)=ZZOBS(II,J)
      ENDDO
C
C Determine observational range over which averaging should occur
      ZZMIN=ZZOBSII(1)+(ZBOTTOMOBS(II)-ZMAXOBS(II))*
     &(ZZOBSII(NNV)-ZZOBSII(1))/ZBOTTOMOBS(II)
      ZZMAX=ZZOBSII(1)+(ZBOTTOMOBS(II)-ZMINOBS(II))*
     &(ZZOBSII(NNV)-ZZOBSII(1))/ZBOTTOMOBS(II)
C
C Evaluate vertical basis functions and dependent variables at
C extremes of averaging interval
      CALL BASIS1D1(NNV,ZZOBSII,ZZMIN,B2,N2,B1,N1)
      UZMINOBS=B1*UZOBS(II,N1)+B2*UZOBS(II,N2)
      VZMINOBS=B1*VZOBS(II,N1)+B2*VZOBS(II,N2)
      N2MIN=N2
      CALL BASIS1D1(NNV,ZZOBSII,ZZMAX,B2,N2,B1,N1)
      UZMAXOBS=B1*UZOBS(II,N1)+B2*UZOBS(II,N2)
      VZMAXOBS=B1*VZOBS(II,N1)+B2*VZOBS(II,N2)
      N1MAX=N1
C
C Perform vertical averaging
      IF(N2MIN.GT.N1MAX)THEN
         UOBS(II)=0.5*(UZMINOBS+UZMAXOBS)
         VOBS(II)=0.5*(VZMINOBS+VZMAXOBS)
      ELSE
         DZ=ZZOBSII(N2MIN)-ZZMIN
         UOBS(II)=0.5*(UZMINOBS+UZOBS(II,N2MIN))*DZ
         VOBS(II)=0.5*(VZMINOBS+VZOBS(II,N2MIN))*DZ
         DO J=N2MIN,N1MAX-1
            DZ=(ZZOBSII(J+1)-ZZOBSII(J))
            UOBS(II)=UOBS(II)+0.5*(UZOBS(II,J)+UZOBS(II,J+1))*DZ
            VOBS(II)=VOBS(II)+0.5*(VZOBS(II,J)+VZOBS(II,J+1))*DZ
         ENDDO
         DZ=ZZMAX-ZZOBSII(N1MAX)
         DZTOTAL=ZZMAX-ZZMIN
         UOBS(II)=(UOBS(II)+0.5*(UZOBS(II,N1MAX)+UZMAXOBS)*DZ)/DZTOTAL
         VOBS(II)=(VOBS(II)+0.5*(VZOBS(II,N1MAX)+VZMAXOBS)*DZ)/DZTOTAL
      ENDIF
      RETURN
      END

C***********************************************************************
C***********************************************************************
      SUBROUTINE CRUISE_AVERAGE(KD,SECKD,ITER,DELT,IN,XNOD,YNOD,NN,NNV,
     &NE,UBAR,VBAR,Z,U,V,Cd,ENZM,MESH,OBSFILE,nf)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
C
C Argument Declaration
C
C   Global Assignments
      INTEGER IN(3,NEDIM)
      REAL XNOD(NNDIM),YNOD(NNDIM),UBAR(NNDIM),VBAR(NNDIM)
      REAL Z(NNDIM,NNVDIM),U(NNDIM,NNVDIM),V(NNDIM,NNVDIM),
     &ENZM(NNDIM,NNVDIM)
      CHARACTER*72 MESH,OBSFILE
C
C   Local assignments
C
C    Observational data
      CHARACTER*72 OBSFILEL,MODFILE
      REAL AOBS(NOBSDIM,NCOLDIM)
C
C    Interpolation arrays, summation arrays, and nodal values of ENZM
      REAL BASISATXY(3)
      REAL AK(NNDIM),AKSUM(NNDIM),USUM(NNDIM),VSUM(NNDIM)
      REAL ZSUM(NNDIM,NNVDIM),ENZMSUM(NNDIM,NNVDIM),ENZMIJ(NNDIM,NNVDIM)
C
C Save necessary arguments
      SAVE OBSFILEL,NFRICT,ITERSTART,ITEREND
C
C When ITER=0
      IF(ITER.EQ.0)THEN
C
C   Read Observations File
         CALL CROPSTRING(OBSFILE,IS,IB,IE)
         WRITE(*,*)'SUBROUTINE CRUISE_AVERAGE: ',OBSFILE(IS:IE)
         CALL READOBS(NOBSDIM,NCOLDIM,OBSFILE,IYEAR,NOBS,NCOL,AOBS)
         OBSFILEL=OBSFILE
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,XNOD,YNOD,NN,NE,AOBS(I,2),
     &                AOBS(I,3),NCLOSEST,NELEM,BASISXY)
            IF(NELEM.NE.0)THEN
               NOBSIN=NOBSIN+1
               DO J=1,NCOL
                  AOBS(NOBSIN,J)=AOBS(I,J)
               ENDDO
            ENDIF
         ENDDO
         WRITE(*,'(1X,I6,A,I6,3A)')NOBSIN,' of the original ',NOBS,
     &   ' O*D entries are in the domain.'
         NOBS=NOBSIN
C
C   Initialize iteration counters and arrays for averaging
         CALL GMT0_GDAY(IYEAR,AOBS(1,1),KDST,SECKDST)
         DELSEC=(KDST-KD)*24.0*3600.0+(SECKDST-SECKD)
         ITERSTART=NINT(DELSEC/DELT)
         CALL GMT0_GDAY(IYEAR,AOBS(NOBS,1),KDEND,SECKDEND)
         DELSEC=(KDEND-KD)*24.0*3600.0+(SECKDEND-SECKD)
         ITEREND=NINT(DELSEC/DELT)
         NFRICT=0
         DO I=1,NN
            AKSUM(I)=0.0
            USUM(I)=0.0
            VSUM(I)=0.0
         ENDDO
         DO J=1,NNV
            DO I=1,NN
               ZSUM(I,J)=0.0
               ENZMSUM(I,J)=0.0
            ENDDO
         ENDDO
         WRITE(*,*)' Mean ak, enzm, & vbar to be computed for interval:'
         CALL GDAY_GMT0(KD,SECKD+ITERSTART*DELT,IYEAR,YDGMT0)
         WRITE(*,*)' ITER1 IYEAR1 YDGMT01 =',ITERSTART,IYEAR,YDGMT0
         CALL GDAY_GMT0(KD,SECKD+ITEREND*DELT,IYEAR,YDGMT0)
         WRITE(*,*)' ITER2 IYEAR2 YDGMT02 =',ITEREND,IYEAR,YDGMT0
      ENDIF
C
C Return if outside desired temporal interval
      IF(ITER.LT.ITERSTART)RETURN
      IF(ITER.GT.ITEREND)RETURN
C
C Perform averaging for desired temporal interval
      NFRICT=NFRICT+1
      AKMIN=0.00001
      CALL QUADSTRESS(AK,Cd,AKMIN,NN,U,V,NNDIM)
      CALL NODALENZ(NNDIM,NNVDIM,NN,NNV,Z,ENZM,ENZMIJ)
      DO I=1,NN
        AKSUM(I)=AKSUM(I)+AK(I)
        USUM(I)=USUM(I)+UBAR(I)
        VSUM(I)=VSUM(I)+VBAR(I)
      ENDDO
      DO J=1,NNV
         DO I=1,NN
            ZSUM(I,J)=ZSUM(I,J)+Z(I,J)
            ENZMSUM(I,J)=ENZMSUM(I,J)+ENZMIJ(I,J)
         ENDDO
      ENDDO
C
C Output files if at the end of desired temporal interval
      IF(ITER.EQ.ITEREND)THEN
         CALL CROPSTRING(OBSFILEL,IS,IB,IE)
         MODFILE=OBSFILEL(IS:IE-4)//'.Q5.ak.s2r'
         OPEN(nf,FILE=MODFILE,STATUS='UNKNOWN')
         REWIND(nf)
         WRITE(nf,'(A)')MESH
         WRITE(nf,*)NFRICT,ITERSTART,ITEREND
         DO I=1,NN
           AKSUM(I)=AKSUM(I)/REAL(NFRICT)
           WRITE(nf,*)I,AKSUM(I)
         ENDDO
         CLOSE(nf)
         MODFILE=OBSFILEL(IS:IE-4)//'.Q5.vbar.v2r'
         OPEN(nf,FILE=MODFILE,STATUS='UNKNOWN')
         REWIND(nf)
         WRITE(nf,'(A)')MESH
         WRITE(nf,*)NFRICT,ITERSTART,ITEREND
         DO I=1,NN
           USUM(I)=USUM(I)/REAL(NFRICT)
           VSUM(I)=VSUM(I)/REAL(NFRICT)
           WRITE(nf,*)I,USUM(I),VSUM(I)
         ENDDO
         CLOSE(nf)
         MODFILE=OBSFILEL(IS:IE-4)//'.Q5.en.s3r'
         OPEN(nf,FILE=MODFILE,STATUS='UNKNOWN')
         REWIND(nf)
         WRITE(nf,'(A)')MESH
         WRITE(nf,*)NFRICT,ITERSTART,ITEREND
         WRITE(nf,*)NNV
         DO I=1,NN
            DO J=1,NNV
              ZSUM(I,J)=ZSUM(I,J)/REAL(NFRICT)
              ENZMSUM(I,J)=ENZMSUM(I,J)/REAL(NFRICT)
              WRITE(nf,*)I,ZSUM(I,J),ENZMSUM(I,J)
            ENDDO
         ENDDO
         CLOSE(nf)
      ENDIF
      RETURN
      END

      SUBROUTINE NUMERICAL_MOORINGS(NN,NNV,NE,XNOD,YNOD,Z,IN,U,V,W,RHO,
     &T,S,ENZM,ENZH,ENZQ,Q2,Q2L,KD,SECKD,NTZ,MOORFILE,nf1)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C Include file for parameter statements
      INCLUDE 'DCMS.DIM'
C
C Global arrays
      real XNOD(nndim),YNOD(nndim)
      real Z(nndim,nnvdim)
      real U(nndim,nnvdim),V(nndim,nnvdim),W(nndim,nnvdim)
      real RHO(nndim,nnvdim),T(nndim,nnvdim),S(nndim,nnvdim)
      real ENZM(nndim,nnvdim),ENZH(nndim,nnvdim),ENZQ(nndim,nnvdim)
      real Q2(nndim,nnvdim),Q2L(nndim,nnvdim)
      integer IN(3,nedim)
      character*72 MOORFILE,INTERNAL,CM,IMFILE
C
C Local declarations
      PARAMETER(MMDIM=9)
      REAL XM(9),YM(9)
      real basis(3)
      character*72 filetzq4d
C
C When NTZ=0, read MOORFILE and create MOOR*_tz.q4d files
C  MOORING 1
      IF(NTZ.EQ.0)THEN
         CALL CSTRING(MOORFILE,IS,IE)
         WRITE(*,*)'SUBROUTINE NUMERICAL_MOORINGS: ',MOORFILE(IS:IE)
         OPEN(nf1,FILE=MOORFILE(IS:IE),STATUS='OLD')
         DO IM=1,MMDIM
            IF(IM.EQ.MMDIM+1)THEN
               WRITE(*,*)'MAXIMUM NUMBER OF MOORINGS EXCEEDED:',MMDIM
               STOP
            ENDIF
            READ(nf1,*,END=10)II,XM(IM),YM(IM)
            IMOOR=IM
            WRITE(INTERNAL,'(i10)')IM
            READ(INTERNAL,'(A)')CM
            CALL CSTRING(CM,IS1,IE1)
            IMFILE=MOORFILE(IS:IE-4)//CM(IS1:IE1)//'.Q5.q4d'
            WRITE(*,*)IMFILE
            OPEN(nf1+IM,FILE=IMFILE,STATUS='UNKNOWN')
            REWIND(nf1+IM)
         ENDDO
 10   ENDIF
C
C Loop over moorings, writing results to respective files
      CALL GDAY_GMT0(kd,seckd,iyear,ydgmt0)
      DO IM=1,IMOOR
         CALL GDAY_GMT0(kd,seckd,iyear,ydgmt0)
         CALL BASIS2D(nndim,nedim,in,xnod,ynod,nn,ne,
     &   XM(IM),YM(IM),nclosest,nelem,basis)
      do j=1,nnv
         CALL POIVAL3D(nndim,nnvdim,nedim,in,j,nelem,basis,z,zpt)
         CALL POIVAL3D(nndim,nnvdim,nedim,in,j,nelem,basis,u,upt)
         CALL POIVAL3D(nndim,nnvdim,nedim,in,j,nelem,basis,v,vpt)
         CALL POIVAL3D(nndim,nnvdim,nedim,in,j,nelem,basis,w,wpt)
         CALL POIVAL3D(nndim,nnvdim,nedim,in,j,nelem,basis,rho,rhopt)
         CALL POIVAL3D(nndim,nnvdim,nedim,in,j,nelem,basis,t,tpt)
         CALL POIVAL3D(nndim,nnvdim,nedim,in,j,nelem,basis,s,spt)
*         CALL POIVAL3D(nndim,nnvdim,nedim,in,j,nelem,basis,q2,q2pt)
*         CALL POIVAL3D(nndim,nnvdim,nedim,in,j,nelem,basis,q2l,q2lpt)
*         CALL POIVAL3D(nndim,nnvdim,nedim,in,j,nelem,basis,enzm,enzmpt)
*         CALL POIVAL3D(nndim,nnvdim,nedim,in,j,nelem,basis,enzh,enzhpt)
*         CALL POIVAL3D(nndim,nnvdim,nedim,in,j,nelem,basis,enzq,enzqpt)
*         write(nf1+IM,'(f9.5,12e13.5)')ydgmt0,
*     &   zpt,upt,vpt,wpt,rhopt,tpt,spt,q2pt,q2lpt,enzmpt,enzhpt,enzqpt
         write(nf1+IM,'(f9.5,7e13.5)')ydgmt0,
     &   zpt,upt,vpt,wpt,rhopt,tpt,spt
      enddo
      ENDDO
C
C END OF ROUTINE
      return
      end

C DJM
C added idriftw to argument list
      SUBROUTINE NUMERICAL_DRIFTERS(KD,SECKD,ITER,DELT,
     &NN,NE,NNV,Z,U,V,W,Q2,Q2L,RHO,T,S,ENZM,ENZH,ENZQ,INDFILE,GR2PATH,
     &idriftw)
C End DJM

C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
      PARAMETER(MAXDROG=5000)
      CHARACTER*72 INDFILE,GR2PATH
      CHARACTER*72 HEADER
C
C Argument Declaration
      REAL Z(NNDIM,NNVDIM),U(NNDIM,NNVDIM),V(NNDIM,NNVDIM),
     &W(NNDIM,NNVDIM),Q2(NNDIM,NNVDIM),Q2L(NNDIM,NNVDIM),
     &RHO(NNDIM,NNVDIM),T(NNDIM,NNVDIM),S(NNDIM,NNVDIM),
     &ENZM(NNDIM,NNVDIM),ENZH(NNDIM,NNVDIM),ENZQ(NNDIM,NNVDIM)
      REAL WZERO(NNDIM,NNVDIM)
      INTEGER ITERL(MAXDROG),IYEARL(MAXDROG)
      REAL YDGMT0L(MAXDROG)
      SAVE KD0,SECKD0,ITERL,WZERO
C DJM
      integer idriftw
C End DJM

C
C If ITER=0; Perform a variety of initialization functions
C
      IF(ITER.EQ.0)THEN
C  initialize timing parameters;
         KD0=KD
         SECKD0=SECKD
C  read number of drogues and launch times from .ind file,
C  then set release iteration for each drogue.
         OPEN(1,FILE=INDFILE,STATUS='OLD')
         DO I=1,11
            READ(1,'(A)')HEADER
         ENDDO
         READ(1,*)NDROG
         IF(NDROG.GT.MAXDROG)THEN
            WRITE(*,*)'MAXIMUM NUMBER OF DROGUES EXCEEDED'
            WRITE(*,*)'INCREASE MAXDROG AND RECOMPILE'
            STOP
         ENDIF
         READ(1,'(A)')HEADER
         DO I=1,NDROG
            READ(1,*)XDROG,YDROG,ZDROG,IYEARL(I),YDGMT0L(I)
            CALL GMT0_GDAY(IYEARL(I),YDGMT0L(I),KDL,SECKDL)
            ITERL(I)=NINT(((KDL-KD)*24.0*3600.0+SECKDL-SECKD)/DELT)
            ITERL(I)=MAX(ITERL(I),0)
         ENDDO
         CLOSE(1)
C DJM
         if (idriftw .eq. 0) then

C  create array of vertical velocities WZERO(I,J)=0
C  (for fixed depth tracking)
         DO I=1,NN
            DO J=1,NNV
               WZERO(I,J)=0.0
            ENDDO
         ENDDO

         endif
C End DJM

C  echo information to standard out
         CALL CSTRING(INDFILE,IS,IE)
         WRITE(*,*)'SUBROUTINE NUMERICAL_DRIFTERS: ',INDFILE(IS:IE)
         WRITE(*,*)' Drifter launch info (DRIFTER# ITER IYEAR YDGMT0):'
         DO I=1,NDROG
            WRITE(*,*)I,ITERL(I),IYEARL(I),YDGMT0L(I)
         ENDDO
C
      ENDIF
C
C For all ITER, Track numerical drifters using DROG3DDT
C
      ELAPSEC=(KD-KD0)*24.0*3600.0+SECKD-SECKD0

C DJM
      if (idriftw .eq. 0) then
         CALL DROG3DDT(NN,NE,NNV,ELAPSEC,ELAPSEC-DELT,ITER,ITERL,
     &                    Z,U,V,WZERO,Q2,Q2L,RHO,T,S,ENZM,ENZH,ENZQ,
     &                    GR2PATH,INDFILE)
        else
         CALL DROG3DDT(NN,NE,NNV,ELAPSEC,ELAPSEC-DELT,ITER,ITERL,
     &                    Z,U,V,W,Q2,Q2L,RHO,T,S,ENZM,ENZH,ENZQ,
     &                    GR2PATH,INDFILE)
       endif
C End DJM

      RETURN
      END
***********************************************************************
C***********************************************************************
      SUBROUTINE NUMERICAL_VPR(KD,SECKD,ITER,DELT,IN,XNOD,YNOD,BATHY,
     &NN,NNV,NE,UBAR,VBAR,Z,U,V,W,RHO,T,S,Q2,Q2L,Cd,ENZM,MESH,
     &OBSFILE,nf)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
C
C Argument Declaration
C
C   FEM Arrays
      INTEGER IN(3,NEDIM)
      REAL XNOD(NNDIM),YNOD(NNDIM),UBAR(NNDIM),VBAR(NNDIM)
      REAL Z(NNDIM,NNVDIM),U(NNDIM,NNVDIM),V(NNDIM,NNVDIM),
     &W(NNDIM,NNVDIM),RHO(NNDIM,NNVDIM),T(NNDIM,NNVDIM),S(NNDIM,NNVDIM),
     &Q2(NNDIM,NNVDIM),Q2L(NNDIM,NNVDIM)
      CHARACTER*72 MESH,OBSFILE
C
C   Observational data
      REAL AOBS(NOBSDIM,NCOLDIM)
      INTEGER KDOBS(NOBSDIM)
      REAL SECOBS(NOBSDIM)
      REAL XOBS(NOBSDIM),YOBS(NOBSDIM),ZOBS(NOBSDIM)
C
C   Dependent variables and horizontal basis functions for
C   observation locations
      REAL ZZOBS(NOBSDIM,NNVDIM),
     &UZOBS(NOBSDIM,NNVDIM),VZOBS(NOBSDIM,NNVDIM),WZOBS(NOBSDIM,NNVDIM),
     &ROBS(NOBSDIM,NNVDIM),TOBS(NOBSDIM,NNVDIM),SOBS(NOBSDIM,NNVDIM)
      REAL BXY(3),BASISOBS(5,NOBSDIM),ZZOBSII(NNVDIM)
      CHARACTER*72 MODFILE
C
C Save necessary arguments
      SAVE XOBS,YOBS,ZOBS,KDOBS,SECKDOBS
      SAVE ZZOBS,UZOBS,VZOBS,WZOBS,BASISOBS,IOBSDELT,IOBS,NOBSDONE
C
C When ITER=0
      IF(ITER.EQ.0)THEN
C
C   Read Observations File
         CALL CROPSTRING(OBSFILE,IS,IB,IE)
         WRITE(*,*)'SUBROUTINE NUMERICAL_VPR: ',OBSFILE(IS:IE)
         CALL READOBS(NOBSDIM,NCOLDIM,OBSFILE,IYEAR,NOBS,NCOL,AOBS)
C
C   Discard observations which lie outside the boundary of the finite
C   element mesh and determine the following for observations which lie
C   within the boundary of the finite element mesh:
C     -Gregrorian time
C     -horizontal basis functions
C     -Set XOBS,YOBS,ZOBS
         NOBSIN=0
         DO I=1,NOBS
            CALL BASIS2D(NNDIM,NEDIM,IN,XNOD,YNOD,NN,NE,AOBS(I,2),
     &                AOBS(I,3),NCLOSEST,NELEM,BXY)
            IF(NELEM.NE.0)THEN
               NOBSIN=NOBSIN+1
               CALL GMT0_GDAY(IYEAR,AOBS(I,1),
     &            KDOBS(NOBSIN),SECOBS(NOBSIN))
               BASISOBS(1,NOBSIN)=NELEM
               BASISOBS(2,NOBSIN)=BXY(1)
               BASISOBS(3,NOBSIN)=BXY(2)
               BASISOBS(4,NOBSIN)=BXY(3)
               XOBS(NOBSIN)=AOBS(I,2)
               YOBS(NOBSIN)=AOBS(I,3)
               ZOBS(NOBSIN)=AOBS(I,4)
            ENDIF
         ENDDO
         WRITE(*,'(1X,I6,A,I6,3A)')NOBSIN,' of the original ',NOBS,
     &   ' O*D entries are in the domain.'
         NOBS=NOBSIN
C
C   Initialize IOBS and NOBSDONE for subsequent calls
         IOBS=1
         NOBSDONE=0
C
C   Open Results Files and Write Headers
         MODFILE=OBSFILE(IS:IE-4)//'.Q5.m4d'
         OPEN(nf,FILE=MODFILE,STATUS='UNKNOWN')
         REWIND(nf)
         WRITE(nf,'(A)')'XXXX'
         WRITE(nf,'(A)')MESH
         WRITE(nf,'(A)')'Q5 results at locs in '//OBSFILE(IB:IE)
         WRITE(nf,'(I4)')IYEAR
         WRITE(nf,'(I4)')10
      ENDIF
C
C For all ITER
C
C   Return if all observations has been passed
      IF(NOBSDONE.EQ.NOBS)RETURN
C
C   Set values if the model timestep just completed spans the
C   the observation time
      DO II=IOBS-IOBSDELT,IOBS-1
         CALL GDAY_GMT0(KDOBS(II),SECOBS(II),IYEAR,YDAYGMT0)
         NELEM=INT(BASISOBS(1,II))
         BXY(1)=BASISOBS(2,II)
         BXY(2)=BASISOBS(3,II)
         BXY(3)=BASISOBS(4,II)
         DO J=1,NNV
            CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,Z,ZB2J)
            CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,U,UB2J)
            CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,V,VB2J)
            CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,W,WB2J)
            CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,RHO,RB2J)
            CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,T,TB2J)
            CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,S,SB2J)
            ZZOBS(II,J)=ZB2J-(ZB2J-ZZOBS(II,J))*BASISOBS(5,II)
            ZZOBSII(J)=ZZOBS(II,J)
            UZOBS(II,J)=UB2J-(UB2J-UZOBS(II,J))*BASISOBS(5,II)
            VZOBS(II,J)=VB2J-(VB2J-VZOBS(II,J))*BASISOBS(5,II)
            WZOBS(II,J)=WB2J-(WB2J-WZOBS(II,J))*BASISOBS(5,II)
            ROBS(II,J)=RB2J-(RB2J-ROBS(II,J))*BASISOBS(5,II)
            TOBS(II,J)=TB2J-(TB2J-TOBS(II,J))*BASISOBS(5,II)
            SOBS(II,J)=SB2J-(SB2J-SOBS(II,J))*BASISOBS(5,II)
         ENDDO
C
C   Perform vertical interpolation to determine dependent variables at
C   desired vertical position
         CALL BASIS1D1(NNV,ZZOBSII,ZOBS(II),B2,N2,B1,N1)
         UMOD=B1*UZOBS(II,N1)+B2*UZOBS(II,N2)
         VMOD=B1*VZOBS(II,N1)+B2*VZOBS(II,N2)
         WMOD=B1*WZOBS(II,N1)+B2*WZOBS(II,N2)
         RMOD=B1*ROBS(II,N1)+B2*ROBS(II,N2)
         TMOD=B1*TOBS(II,N1)+B2*TOBS(II,N2)
         SMOD=B1*SOBS(II,N1)+B2*SOBS(II,N2)
         WRITE(nf,'(F10.4,2F10.0,F10.4,6E14.6)')YDAYGMT0,XOBS(II),
     &   YOBS(II),ZOBS(II),UMOD,VMOD,WMOD,RMOD,TMOD,SMOD
      ENDDO
      NOBSDONE=IOBS-1
C   Return if the next timestep does not span the temporal location
C   of the next observational location
      IOBSDELT=0
 2    DTSEC=(KDOBS(IOBS)-KD)*24.0*3600.0+SECOBS(IOBS)-SECKD
      IF(DTSEC.LT.0.0)RETURN
      IF(DTSEC.GE.DELT)RETURN
C
C   Evaluate temporal basis function at this time level for this
C   observation, perform horizontal interpolation, and increment
C   counter
      BASISOBS(5,IOBS)=1.0-DTSEC/DELT
      NELEM=INT(BASISOBS(1,IOBS))
      BXY(1)=BASISOBS(2,IOBS)
      BXY(2)=BASISOBS(3,IOBS)
      BXY(3)=BASISOBS(4,IOBS)
      DO J=1,NNV
       CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,Z,ZZOBS(IOBS,J))
       CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,U,UZOBS(IOBS,J))
       CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,V,VZOBS(IOBS,J))
       CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,W,WZOBS(IOBS,J))
       CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,RHO,ROBS(IOBS,J))
       CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,T,TOBS(IOBS,J))
       CALL POIVAL3D(NNDIM,NNVDIM,NEDIM,IN,J,NELEM,BXY,S,SOBS(IOBS,J))
      ENDDO
      IOBSDELT=IOBSDELT+1
      IOBS=IOBS+1
      GO TO 2
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE TAU_OF_TIME(KD,SECKD,ITER,ATMXT,ATMYT,
     &KDSTART,SECKDSTART,RAMPKD,OBSFILE,nf)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
C
C Argument Declaration
C
C   Observational data
      REAL AOBS(NOBSDIM,NCOLDIM)
      CHARACTER*72 OBSFILE
      CHARACTER*72 MODFILE
C
C Save necessary arguments
      SAVE NOBS,NCOL,AOBS,GMT0OBS1,KDOBS1,SECKDOBS1,NOBSINT,RHOREF
C
C When ITER=0
      IF(ITER.EQ.0)THEN
C
C  Read Observations File, determine Gregorian time at beginning of
C  observation file, initialize counter, and set reference density
         CALL CROPSTRING(OBSFILE,IS,IB,IE)
         CALL GDAY_GMT0(KDSTART,SECKDSTART,IYEAR,YDGMT0)
         WRITE(*,*)'SUBROUTINE TAUXY_OF_TIME: ',OBSFILE(IS:IE)
         WRITE(*,*)' Wind stress inforced beginning at:'
         WRITE(*,*)' IYEAR YDGMT0 =',IYEAR,YDGMT0
         WRITE(*,*)' Ramp up period (days) =',RAMPKD
         CALL READOBS(NOBSDIM,NCOLDIM,OBSFILE,IYEAR,NOBS,NCOL,AOBS)
         GMT0OBS1=AOBS(1,1)
         CALL GMT0_GDAY(IYEAR,GMT0OBS1,KDOBS1,SECKDOBS1)
         NOBSINT=0
         CALL EQSTATE2(10.0,32.0,SIGTREF)
         RHOREF=1000.0+SIGTREF
         PRINT*, RHOREF
C
C  Now:
C   AOBS(I,1) are the times corresponding to the wind stress values
C   AOBS(I,2) are the taux values at time AOBS(I,1)
C   AOBS(I,3) are the tauy values at time AOBS(I,1)
C
C  Open model output file
         MODFILE=OBSFILE(IS:IE-4)//'.Q5.m1d'
         OPEN(nf,FILE=MODFILE,STATUS='UNKNOWN')
         REWIND(nf)
         WRITE(nf,'(A)')'XXXX'
         WRITE(nf,'(A)')'null comment line'
         WRITE(nf,'(A)')'Q5 wind stress values from '//OBSFILE(IB:IE)
         WRITE(nf,'(I4)')IYEAR
         WRITE(nf,'(I4)')3
      ENDIF
C
C For all ITER
C
C  Determine current simulation time in reference frame of .o1d file
      GMT0NOW=GMT0OBS1+REAL(KD-KDOBS1)+(SECKD-SECKDOBS1)/24.0/3600.0
C
C  If simulation time preceeds first observational time,
C  utilize first observation
      IF(GMT0NOW.LE.AOBS(1,1))THEN
         TAUXT=AOBS(1,2)
         TAUYT=AOBS(1,3)
         GO TO 20
      ENDIF
C
C  If simulation time exceeds final observational time,
C  utilize last observation
      IF(GMT0NOW.GE.AOBS(NOBS,1))THEN
         TAUXT=AOBS(NOBS,2)
         TAUYT=AOBS(NOBS,3)
         GO TO 20
      ENDIF
C
C  Else find current observation interval and apply temporal basis
C  functions to determine values (note that this requires time increases
C  monotonically in observation file
 15   IF(GMT0NOW.GE.AOBS(NOBSINT+1,1))THEN
         NOBSINT=NOBSINT+1
         GO TO 15
      ENDIF
      B1=(AOBS(NOBSINT+1,1)-GMT0NOW)/(AOBS(NOBSINT+1,1)-AOBS(NOBSINT,1))
      TAUXT=AOBS(NOBSINT+1,2)-B1*(AOBS(NOBSINT+1,2)-AOBS(NOBSINT,2))
      TAUYT=AOBS(NOBSINT+1,3)-B1*(AOBS(NOBSINT+1,3)-AOBS(NOBSINT,3))
C
C  Set ATMXT,ATMYT using observational data and desired ramp-up
 20   DELTAKD=REAL(KD-KDSTART)+(SECKD-SECKDSTART)/24.0/3600.0
      IF(RAMPKD.EQ.0)THEN
         RAMPUP=1.0
      ELSE
         RAMPUP=DELTAKD/RAMPKD
         RAMPUP=AMAX1(RAMPUP,0.0)
         RAMPUP=AMIN1(RAMPUP,1.0)
      ENDIF
      TAUXT=RAMPUP*TAUXT
      TAUYT=RAMPUP*TAUYT
      ATMXT=TAUXT/RHOREF
      ATMYT=TAUYT/RHOREF
C
C  Write wind stress values to .m1d file
      WRITE(nf,'(F10.4,2F10.5)')GMT0NOW,TAUXT,TAUYT
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE QFLUX_OF_TIME(KD,SECKD,ITER,BTEMPT,
     &KDSTART,SECKDSTART,RAMPKD,OBSFILE,nf)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
C
C Argument Declaration
C
C   Observational data
      REAL AOBS(NOBSDIM,NCOLDIM)
      CHARACTER*72 OBSFILE
      CHARACTER*72 MODFILE
C
C Save necessary arguments
      SAVE NOBS,NCOL,AOBS,GMT0OBS1,KDOBS1,SECKDOBS1,NOBSINT,RHOREF
C
C When ITER=0
      IF(ITER.EQ.0)THEN
C
C  Read Observations File, determine Gregorian time at beginning of
C  observation file, initialize counter, and set reference density
         CALL CROPSTRING(OBSFILE,IS,IB,IE)
         CALL GDAY_GMT0(KDSTART,SECKDSTART,IYEAR,YDGMT0)
         WRITE(*,*)'SUBROUTINE QFLUX_OF_TIME: ',OBSFILE(IS:IE)
         WRITE(*,*)' Heat flux inforced beginning at:'
         WRITE(*,*)' IYEAR YDGMT0 =',IYEAR,YDGMT0
         WRITE(*,*)' Ramp up period (days) =',RAMPKD
         CALL READOBS(NOBSDIM,NCOLDIM,OBSFILE,IYEAR,NOBS,NCOL,AOBS)
         GMT0OBS1=AOBS(1,1)
         CALL GMT0_GDAY(IYEAR,GMT0OBS1,KDOBS1,SECKDOBS1)
         NOBSINT=0
         CALL EQSTATE2(10.0,32.0,SIGTREF)
         RHOREF=1000.0+SIGTREF
C
C  Now:
C   AOBS(I,1) are the times corresponding to the heat flux values
C   AOBS(I,2) are heat flux values (watts/m2) at times in AOBS(I,1)
C
C  Open model output file
         MODFILE=OBSFILE(IS:IE-4)//'.Q5.m1d'
         OPEN(nf,FILE=MODFILE,STATUS='UNKNOWN')
         REWIND(nf)
         WRITE(nf,'(A)')'XXXX'
         WRITE(nf,'(A)')'null comment line'
         WRITE(nf,'(A)')'Q5 heat flux values from '//OBSFILE(IB:IE)
         WRITE(nf,'(I4)')IYEAR
         WRITE(nf,'(I4)')2
      ENDIF
C
C For all ITER
C
C  Determine current simulation time in reference frame of .o1d file
      GMT0NOW=GMT0OBS1+REAL(KD-KDOBS1)+(SECKD-SECKDOBS1)/24.0/3600.0
C
C  If simulation time preceeds first observational time,
C  utilize first observation
      IF(GMT0NOW.LE.AOBS(1,1))THEN
         QFLUX=AOBS(1,2)
         GO TO 20
      ENDIF
C
C  If simulation time exceeds final observational time,
C  utilize last observation
      IF(GMT0NOW.GE.AOBS(NOBS,1))THEN
         QFLUX=AOBS(NOBS,2)
         GO TO 20
      ENDIF
C
C  Else find current observation interval and apply temporal basis
C  functions to determine values (note that this requires time increases
C  monotonically in observation file
 15   IF(GMT0NOW.GE.AOBS(NOBSINT+1,1))THEN
         NOBSINT=NOBSINT+1
         GO TO 15
      ENDIF
      B1=(AOBS(NOBSINT+1,1)-GMT0NOW)/(AOBS(NOBSINT+1,1)-AOBS(NOBSINT,1))
      QFLUX=AOBS(NOBSINT+1,2)-B1*(AOBS(NOBSINT+1,2)-AOBS(NOBSINT,2))
C
C  Implement desired ramp-up
 20   DELTAKD=REAL(KD-KDSTART)+(SECKD-SECKDSTART)/24.0/3600.0
      IF(RAMPKD.EQ.0)THEN
         RAMPUP=1.0
      ELSE
         RAMPUP=DELTAKD/RAMPKD
         RAMPUP=AMAX1(RAMPUP,0.0)
         RAMPUP=AMIN1(RAMPUP,1.0)
      ENDIF
      QFLUX=RAMPUP*QFLUX
      Cp=4186.0
      BTEMPT=QFLUX/(RHOREF*CP)
C
C  Write heat flux values to .m1d file
      WRITE(nf,'(F10.4,2F10.5)')GMT0NOW,QFLUX
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE OVERWRITEMASSVAR(TSFILEPTH,NN,NNV,Z,RHO,T,S,nf)
C Include file for parameter statements
      INCLUDE 'DCMS.DIM'
C Fixed global variables
      REAL Z(NNDIM,NNVDIM),
     &RHO(NNDIM,NNVDIM),T(NNDIM,NNVDIM),S(NNDIM,NNVDIM)
      CHARACTER*72 TSFILEPTH
      CHARACTER*72 TSFILE,HEADER
C Local variables
      REAL ZIC(NNVDIM),TSIC(NNVDIM)
c
c Set filenames
      CALL CSTRING(TSFILEPTH,IS,IE)
      WRITE(*,*)'SUBROUTINE OVERWRITEMASSVAR:'
c
c Set reference value of sigmat and density
      CALL EQSTATE2(10.0,32.0,SIGTREF)
      RHOREF=1000.0+SIGTREF
c
c   Load temperature
      TSFILE=TSFILEPTH(IS:IE)//'ICS.T.s3r'
      OPEN(nf,FILE=TSFILE,STATUS='OLD')
      READ(nf,'(A)')HEADER
      READ(nf,'(A)')HEADER
      READ(nf,*)NNVIC
      PRINT*,NNVIC
      IF(NNVIC.GT.NNVDIM)THEN
      WRITE(*,*)'THE NUMBER OF VERTICAL NODES IN THE FILE USED TO INITIA
     &LIZE THE TEMPERATURE FIELD EXCEEDS NNVDIM.  EXECUTION TERMINATED
     &IN SUBROUTINE OVERWRITEMASSVAR.'
      STOP
      ENDIF
      DO I=1,NN
         DO J=1,NNVIC
            READ(nf,*)II,ZIC(J),TSIC(J)
         ENDDO
         DO J=1,NNV
            CALL BASIS1D1(NNVIC,ZIC,Z(I,J),B2,N2,B1,N1)
            T(I,J)=B1*TSIC(N1)+B2*TSIC(N2)
         ENDDO
      ENDDO
      CLOSE(nf)
      WRITE(*,'(1X,I6,A,I6,A,A)')NN,' x ',NNVIC,' values loaded from ',
     &TSFILEPTH(IS:IE)//'ICS.T.s3r'
c   Load salinity and compute rho
      TSFILE=TSFILEPTH(IS:IE)//'ICS.S.s3r'
      OPEN(nf,FILE=TSFILE,STATUS='OLD')
      READ(nf,'(A)')HEADER
      READ(nf,'(A)')HEADER
      READ(nf,*)NNVIC
      IF(NNVIC.GT.NNVDIM)THEN
      WRITE(*,*)'THE NUMBER OF VERTICAL NODES IN THE FILE USED TO INITIA
     &LIZE THE SALINITY FIELD EXCEEDS NNVDIM.  EXECUTION TERMINATED IN S
     &UBROUTINE OVERWRITEMASSVAR.'
      STOP
      ENDIF
      DO I=1,NN
         DO J=1,NNVIC
            READ(nf,*)II,ZIC(J),TSIC(J)
         ENDDO
         DO J=1,NNV
            CALL BASIS1D1(NNVIC,ZIC,Z(I,J),B2,N2,B1,N1)
            S(I,J)=B1*TSIC(N1)+B2*TSIC(N2)
            CALL EQSTATE2(T(I,J),S(I,J),SIGTMID)
            RHO(I,J)=(SIGTMID-SIGTREF)/RHOREF
         ENDDO
      ENDDO
      CLOSE(nf)  
      WRITE(*,'(1X,I6,A,I6,A,A)')NN,' x ',NNVIC,' values loaded from ',
     &TSFILEPTH(IS:IE)//'ICS.S.s3r'
      RETURN
      END
***********************************************************************
C***********************************************************************
      SUBROUTINE QUODDY2VSQ5(nn,nnv,ne,xnod,ynod,z,in,u,v,w,rho,t,s,
     &enzm,enzh,enzq,q2,q2l,kd,seckd,
     &px1,py1,px2,py2,nseg,sq4file,elefile)
C-----------------------------------------------------------------------
C purpose: This subroutine interpolates density, velocity, and
C            turbulence parameter data from
C            a 3-D mesh to a 2-D vertical slice and writes LABLIB
C            Standard .sq4 and .ele files containing the results.
C
C restrictions: QUODDY style 3-D mesh
C               currently uses mercator projection mercgom2
C               currently uses vertical mesh generator VSIGMAMESH
C               requires the include file DCMS.DIM
C
C index range: i=1,nn, j=1,nnv, k=1,ne, m=0,nseg, n=1,nnv*(nseg+1)
C
C inputs from calling routine:
C         3-D mesh input data
C          nn - number of horizontal nodes
C          nnv - number of nodes in the vertical
C          ne - number of horizontal elements
C          x(i),y(i) - cartesian horizontal nodal coordinates
C          z(i,j) - 3-D vertical nodal coordinates
C          in(3,k) - incidence list
C          rho(i,j),enzm(i,j),q2(i,j),q2l(i,j) - 3-D nodal arrays
C          u(i,j),v(i,j),w(i,j) - 3-D nodal velocities (REAL DATA)
C         Vertical slice horizontal coordinates
C          (px1,py1) is the location of endpoint 1
C          (px2,py2) is the location of endpoint 2
C         Uniform spacing 2-D vertical mesh data
C          nseg - number of horizontal sections
C
C outputs:Vertical slice 2-D mesh data files
C          sq4file - 2-D vertical slice .sq4 file
C          elefile -  2-D vertical slice .ele file
C
C history: Written by Christopher E. Naimie
C          Dartmouth College
C          OUTPUTQTONRVa - FEB 24, 1994
C          QUODDY2VNRV - JUL 27, 1994
C          QUODDY2VRVT - SEP 27, 1994
C          QUODDY2VSQ5 - JUL 03, 1998
C-----------------------------------------------------------------------
C Include file for parameter statements
      INCLUDE 'DCMS.DIM'
C
C Global arrays
      real xnod(nndim),ynod(nndim)
      real z(nndim,nnvdim)
      real u(nndim,nnvdim),v(nndim,nnvdim),w(nndim,nnvdim)
      real rho(nndim,nnvdim),t(nndim,nnvdim),s(nndim,nnvdim)
      real q2(nndim,nnvdim),q2l(nndim,nnvdim)
      real enzm(nndim,nnvdim),enzh(nndim,nnvdim),enzq(nndim,nnvdim)
      integer in(3,nedim)
C
C Arrays filled by VSIGMAMESH or VUNIMESH
      integer incd(3,nedim)
      real px(0:nndim),py(0:nndim),ord(nndim),abs(nndim)
C
C Assignments for 2-D mesh => .sq4 file output
C    (ord,abs) are the (horizontal,vertical) coordinates of the nodal
C      locations:
C        ord is the distance from node 2
C        abs is the value of z interpolated from the quoddy4 grid
C    (ut,un,uz) is the right-handed coordinate system components of
C      velocity:  (left to right, into the page, vertical) for vertical
C      slices and (x,y,z) for horizontal slices.
C    (vt,vs,vq2,vq2l,venzm,venzh,venzq) are the other quoddy4 dependent
C      variables
      real ut(nndim),un(nndim),uz(nndim)
      real vrho(nndim),vt(nndim),vs(nndim)
      real vq2(nndim),vq2l(nndim)
      real venzm(nndim),venzh(nndim),venzq(nndim)
      real basis(3)
C
C Other assignments
      character*(*) sq4file,elefile
C
C START OF EXECUTABLE ROUTINE
C
C Check array sizes
      nstop=0
      if(nn.GT.nndim)then
         write(*,*)'Increase nndim in QUODDY2VSQ5 to at least ',nn
         nstop=1
      endif
      if(nnv.GT.nnvdim)then
         write(*,*)'Increase nnvdim in QUODDY2VSQ5 to at least ',nnv
         nstop=1
      endif
      if(ne.GT.nedim)then
         write(*,*)'Increase nedim in QUODDY2VSQ5 to at least ',ne
         nstop=1
      endif
      if(nstop.EQ.1)stop
C
C Create vertical mesh using VSIGMAMESH or VUNIMESH
          CALL VSIGMAMESH(nndim,nedim,nnvdim,nn,xnod,ynod,nnv,z,ne,in,
     &nseg,nu,px,py,abs,ord,nele,incd,unitx,unity,stot,bmax,
     &px1,py1,px2,py2)
C
C Interpolate nodal 3-D quantities to the 2-D vertical mesh nodes
      n=0
      do 20 m=nseg,0,-1
        CALL BASIS2D(nndim,nedim,in,xnod,ynod,nn,ne,px(m),py(m),
     &               nclosest,nelem,basis)
        do 20 j=1,nnv
        n=n+1
        CALL POIVAL3D(nndim,nnvdim,nedim,in,j,nelem,basis,u,upt)
        CALL POIVAL3D(nndim,nnvdim,nedim,in,j,nelem,basis,v,vpt)
        CALL POIVAL3D(nndim,nnvdim,nedim,in,j,nelem,basis,w,uz(n))
        CALL POIVAL3D(nndim,nnvdim,nedim,in,j,nelem,basis,rho,vrho(n))
        CALL POIVAL3D(nndim,nnvdim,nedim,in,j,nelem,basis,t,vt(n))
        CALL POIVAL3D(nndim,nnvdim,nedim,in,j,nelem,basis,s,vs(n))
        CALL POIVAL3D(nndim,nnvdim,nedim,in,j,nelem,basis,q2,vq2(n))
        CALL POIVAL3D(nndim,nnvdim,nedim,in,j,nelem,basis,q2l,vq2l(n))
        CALL POIVAL3D(nndim,nnvdim,nedim,in,j,nelem,basis,enzm,venzm(n))
        CALL POIVAL3D(nndim,nnvdim,nedim,in,j,nelem,basis,enzh,venzh(n))
        CALL POIVAL3D(nndim,nnvdim,nedim,in,j,nelem,basis,enzq,venzq(n))
C
C Compute along and cross bank velocities
        ut(n)=upt*unity-vpt*unitx
        un(n)=upt*unitx+vpt*unity
 20   continue
C
C Check array sizes
      nstop=0
      if(n.GT.nndim)then
         write(*,*)'Increase nndim in QUODDY2VSQ5 to at least ',n
         nstop=1
      endif
      if(nele.GT.nedim)then
         write(*,*)'Increase nedim in QUODDY2VSQ5 to at least ',nele
         nstop=1
      endif
      if(nstop.EQ.1)stop
C
C Write .sq4 and .ele files for vertical slice
      open(21,file=sq4file,status='unknown')
      rewind(21)
      write(21,'(13e13.5)')(abs(n),ord(n),
     &ut(n),un(n),uz(n),vrho(n),vt(n),vs(n),vq2(n),vq2l(n),
     &venzm(n),venzh(n),venzq(n),n=1,nu)
      close(21)
      open(22,file=elefile,status='unknown')
      rewind(22)
      write(22,'(4i7)')(n,incd(1,n),incd(2,n),incd(3,n),n=1,nele)
      close(22)
C
C END OF ROUTINE
      return
      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
      CHARACTER*72 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
C
C Read data
      DO I=1,NOBSDIM
         READ(1,*,END=11)(AOBS(I,J),J=1,NCOL)
         NOBS=I
      ENDDO
 11   CLOSE(1)
      WRITE(*,'(1X,I6,A,I6,A)')NOBS,' x',NCOL,' O*D entries read from fi
     &le'
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE NODALENZ(NNDIM,NNVDIM,NN,NNV,Z,ENZELEM,ENZNODAL)
C-----------------------------------------------------------------------
C Convert the elemental based values of ENZ to nodal based values using
C ENZNODAL(I,J)=(DZ(I,J-1)*ENZELEM(I,J-1)+DZ(I,J)*ENZELEM(I,J))
C              /(DZ(I,J-1)+DZ(I,J))
C Recall that the vertical nodal index ranges from J=1,NNV and that
C the vertical elemental index ranges from J=1,NNV-1 (i.e. NEV=NNV-1).
C-----------------------------------------------------------------------
      REAL Z(NNDIM,NNVDIM),ENZELEM(NNDIM,NNVDIM),ENZNODAL(NNDIM,NNVDIM)
      NEV=NNV-1
      DO I=1,NN
         ENZNODAL(I,1)=ENZELEM(I,1)
         DO J=2,NEV
            DZ1=Z(I,J)-Z(I,J-1)
            DZ2=Z(I,J+1)-Z(I,J)
            DZTOT=DZ1+DZ2
            ENZNODAL(I,J)=(DZ1*ENZELEM(I,J-1)+DZ2*ENZELEM(I,J))/DZTOT
         ENDDO
         ENZNODAL(I,NNV)=ENZELEM(I,NEV)
      ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE WRITES2RV2RIOS(nndim,nn,nnv,Z,UBAR,VBAR,iter,mesh,
     &kd,seckd,deglat)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C Global declarations
      real Z(nndim,*),UBAR(*),VBAR(*)
      character*(*) mesh
C
C Local declarations
      character*72 in123,citer,filename
C
C Determine day, month, and year from kd
      CALL DMY(iday,imonth,iyear,kd)
C
C Write result files
      write(in123,'(i10)')iter
      read(in123,'(a)')citer
      CALL CSTRING(citer,istart,iend)
C
      filename='ZETA_ITER_'//citer(istart:iend)//'.s2r'
      write(*,*)'writing results file: ',filename
      open(19,file=filename,status='unknown')
      rewind(19)
      write(19,102)mesh,iday,imonth,iyear,seckd,deglat
 102  format(a,/,i2,1x,i2,1x,i4,1x,f7.1,1x,f10.6,
     &' (day,month,year,sec,deglat): free surface elevations')
      write(19,'(i5,1e12.4)')(i,z(i,nnv),i=1,nn)
      close(19)
C
      filename='VBAR_ITER_'//citer(istart:iend)//'.v2r'
      write(*,*)'writing results file: ',filename
      open(20,file=filename,status='unknown')
      rewind(20)
      write(20,201)mesh,iday,imonth,iyear,seckd,deglat
 201  format(a,/,i2,1x,i2,1x,i4,1x,f7.1,1x,f10.6,
     &' (day,month,year,sec,deglat): vertically average velocities')
      write(20,'(i5,2e12.4)')(i,UBAR(i),VBAR(i),i=1,nn)
      close(20)
      return
      end
C***********************************************************************
C***********************************************************************
      SUBROUTINE WRITES2RV2R(nndim,nn,nnv,Z,U,V,iter,mesh)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C Global declarations
      real Z(nndim,*),U(*),V(*)
      character*(*) mesh
C
C Local declarations
      character*72 in123,citer
C
C Write result files
      write(in123,'(i10)')iter
      read(in123,'(a)')citer
      CALL CSTRING(citer,istart,iend)
C
      write(*,*)'writing results file:  reszeta',
     &citer(istart:iend),'.s2r'
      open(19,file='reszeta'//citer(istart:iend)//'.s2r',
     &status='unknown')
      rewind(19)
      write(19,102)mesh,citer(istart:iend)
 102  format(a,/,
     &'Residual surface elevations - Iteration # ',a)
      write(19,'(i5,1e12.4)')(i,z(i,nnv),i=1,nn)
      close(19)
C
      write(*,*)'writing results file:  resvbar',
     &citer(istart:iend),'.v2r'
      open(20,file='resvbar'//citer(istart:iend)//'.v2r',
     &status='unknown')
      rewind(20)
      write(20,201)mesh,citer(istart:iend)
 201  format(a,/,
     &'Residual average velocities - Iteration # ',a)
      write(20,'(i5,2e12.4)')(i,U(i),V(i),i=1,nn)
      close(20)
      return
      end
C***********************************************************************
C***********************************************************************
      SUBROUTINE WRITETOPBOTV2R(nndim,nn,nnv,Z,U,V,iter,mesh)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C Global declarations
      real Z(nndim,*),U(nndim,*),V(nndim,*)
      character*(*) mesh
C
C Local declarations
      character*72 in123,citer
C
C Write result files
      write(in123,'(i10)')iter
      read(in123,'(a)')citer
      CALL CSTRING(citer,istart,iend)
C
      write(*,*)'writing results file:  resvtop',
     &citer(istart:iend),'.v2r'
      open(19,file='resvtop'//citer(istart:iend)//'.v2r',
     &status='unknown')
      rewind(19)
      write(19,102)mesh,citer(istart:iend)
 102  format(a,/,
     &'Residual surface velocities - Iteration # ',a)
      write(19,'(i5,2e12.4)')(i,U(i,nnv),V(i,nnv),i=1,nn)
      close(19)
C
      write(*,*)'writing results file:  resvbot',
     &citer(istart:iend),'.v2r'
      open(20,file='resvbot'//citer(istart:iend)//'.v2r',
     &status='unknown')
      rewind(20)
      write(20,201)mesh,citer(istart:iend)
 201  format(a,/,
     &'Residual bottom velocities - Iteration # ',a)
      write(20,'(i5,2e12.4)')(i,U(i,1),V(i,1),i=1,nn)
      close(20)
      return
      end
C***********************************************************************
C***********************************************************************

C***********************************************************************
C***********************************************************************
        SUBROUTINE EQSTATE2(T,S,R)
C-----------------------------------------------------------------------
C purpose: This subroutine computes the perturbation density (R; in
C          sigmat units) as a function of two variables; temperature
C          (T; 0 to 40 degrees C) and salinity (S; 0 to 42 psu).
C
C reference: UNESCO equation of state as published in Gill(1982)
C
C Contributed by Douglas T. Morgan
C-----------------------------------------------------------------------
C
C...Ensure T and S are within Range
C
      IF((T.GE.0).AND.(T.LE.40))GO TO 10
      WRITE(2,'(///A)')'EXECUTATION TERMINATED IN EQSTATE2 AS T LIES OUT
     &SIDE PERMISSIBLE RANGE (0 to 40 degrees C).'
      STOP
 10   IF((S.GE.0).AND.(S.LE.42))GO TO 20
      WRITE(2,'(///A)')'EXECUTATION TERMINATED IN EQSTATE2 AS S LIES OUT
     &SIDE PERMISSIBLE RANGE (0 to 42 psu).'
      STOP
C
C...Density of pure water
C
 20   R=-0.157406
     &  +6.793952E-2*T-9.095290E-3*T**2.0
     &  +1.001685E-4*T**3.0-1.120083E-6*T**4.0
     &  +6.536332E-9*T**5.0
C
C...Density at one standard atmosphere
C
      R=R+S*(0.824493-4.0899E-3*T+7.6438E-5*T**2.0
     &          -8.2467E-7*T**3.0+5.3875E-9*T**4.0)
     &   +S**1.5*(-5.72466E-3+1.0227E-4*T-1.6546E-6*T**2.0)
     &   +4.8314E-4*S**2.0
      RETURN
      END
