C***********************************************************************
C first cut of bcs subroutine for the SACO portion of SACO_MOODY
C***********************************************************************
      SUBROUTINE BCS2(KD,SECKD,ITER,DEGLAT,UBAR,VBAR,III)
C Include file for parameter statements
      INCLUDE 'CASCO.DIM'
      include 'startstop.h'
      include 'deltabcs.h' 
      include 'sacotime.h'
      include 'mesh.h'

C Need to include mesh.h to get NR and LRVS.
C should simply pass nr and lrvs to BCS2 
C
C DEGLAT could be removed 

C Fixed global variables
      INTEGER ITER,III
      INTEGER KD
      REAL    SECKD
      REAL    DEGLAT
      REAL UBAR(NNDIM),VBAR(NNDIM)
C--
C-----------------------------------------------------------------------
C This subroutine must be modified by the user to assign proper
C   Boundary Conditions.  There are four sections which require
C   modification, beginning at statement labels 3,4,5, and 6.  
C Each section must begin with the appropriately labelled continue 
C   statement, and end with a return statement.
C This subroutine is called once at the beginning of the simulation 
C   (with ITER=0 and III=0) for initialization purposes; then it is
C   called as necessary during each time step to set the desired
C   horizontal boundary conditions.
C================Beginning-of-user-specified-instructions===============
C
C User-specified local assignments
C Gregorian time at the intervals of the BCS from MOODY.
c       integer kdbc(NTBCDIM)
c	real    secbc(NTBCDIM)
c	save    kdbc, secbc

	real    bt1, bt2 
        integer i, p1, p2
	integer nbc
c	real    dseconds
c	integer dday
C
C  Declarations for multiple tidal constituents:  
C
C Beginning of executable code
C
      IF(III.EQ.0)GO TO 3
      IF(III-2)4,5,6
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C This section may be used for initialization purposes.
C
    3 CONTINUE

C We will build the time arrays that go with the input boundary conditions.
C It would be better to build the times on the fly but roundoff gets
C in the way of my simple-minded calculation.  But I know that there
C must be a better way.
C
C      
       write(*,*)'BCS2: initialize for synthesis of moody bcs' 
       write(*,*)'DELTTBC, ntbc ', DELTTBC, ntbc
c      do i = 1, ntbcdim
c        kdbc(i) = 0
c	 secbc(i) = 0.0
c       enddo

c       kdbc(1) = kd0
c       secbc(1) = sec0
c      if(ntbc.lt.1) then 
c	  WRITE(*,*)'bcs1: bad initialization. ntbc <1',ntbc
c	  WRITE(2,*)'bcs1: bad initialization. ntbc <1',ntbc
c	  stop
c       endif
c       do i = 2, ntbc
c	  kdbc(i) = kdbc(i-1)
c	  secbc(i) = secbc(i-1) + DELTTBC
c	  write(*,*)'before ',i, kdbc(i), secbc(i)
c	  if(secbc(i).ge.86400) call up_date(kdbc(i),secbc(i))
c	  write(*,*)'after ',i, kdbc(i), secbc(i)
c       enddo
c       write(*,*)'BCS1 startup', kd0,sec0 
       write(*,*)'First BCS time ',kdbc(1),secbc(1)
c       write(*,*)'Second BCS time ',kdbc(2),secbc(2)
       write(*,*)'Final BCS time ',kdbc(ntbc),secbc(ntbc), ntbc

C 
      RETURN
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Surface Elevation BC's are assigned in this section.  One statement of
C   the form UBAR(I)=ELEV must appear for each specified surface
C   elevation node, where I is the horizontal node number and ELEV is
C   the elevation for node I.  
C A return statement must follow the last elevation BC assignment.
C
    4 CONTINUE
c
c  Synthesize the bcs at saco time (kd, seckd) from the time 
C  series of boundary conditions from MOODY which are given at 
C  DELTTBC intervals starting at KD0, SEC0. 

C  find interval which contains the current time
c      write(*,*)'BCS2, kd, seckd,iter', kd,seckd,iter
c     dday = kd - kd0 
c     dseconds = dday*24*3600 + (seckd - sec0)
c     dday = kd - kdbc(1) 
c     dseconds = dday*24*3600 + (seckd - secbc(1))
c     p1 = int(dseconds/DELTTBC)  + 1
c     if(p1.le.0) then 
c	 write(*,*)'problems in bcs1', kd,kdbc(1),seckd,secbc(1),p1
c	 write(*,*)dday,dseconds,DELTTBC
c	 stop 
c      endif
c      p2 = p1+1 
c      write(*,*)p1,p2,dseconds, dseconds/DELTTBC

C  now calculate the time domain basis functions

c     dseconds = (kd - kdbc(p1))*24.*3600. + seckd - secbc(p1)
c     BT2 = dseconds/DELTTBC 

      if(iter.eq.0) then 
	 p1 = 1
	 p2 = 1
	 BT2 = 0
	 BT1 = 0
      else
         p2 = pkj(iter)
         p1 = p2 - 1
      endif
      if(p1.lt.0)  p1  = 1 

      BT2 = kjweight(iter)
      BT1 =  1.0 - BT2 
c     write(*,*)'BCS2: ',dseconds,BT1, BT2
c     write(*,*)'BCS2: ',p1,p2,BT1, BT2

      DO I=1,NR 
         NBC = LRVS(I)
         UBAR(NBC)=0.0
         UBAR(NBC) = BT1*DBC(I,P1) + BT2*DBC(I,P2)
C	 if(i.eq.2) write(*,*)nbc,DBC(I,P1), DBC(I,P2),ubar(nbc)
      ENDDO
c     i = 21 
c     nbc = lrvs(i) 
c     write(53,*)nbc,kd,seckd,
c    $       (kd-kdbc(1))*24. + (seckd-secbc(1))/3600., ubar(nbc)
c      write(*,*)'leaving bcs2'
      RETURN
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Nonzero normal velocity BC's are assigned in this section.  One 
C   statement of the form UBAR(I)=VNORM must appear for each set nonzero
C   normal velocity node, where I is the horizontal node number and 
C   VNORM is the normal velocity for node I.  
C   Positive normal velocity indicates flow out of the system.  
C A return statement must follow the last normal velocity BC assignment.
C
    5 CONTINUE
      WRITE(*,*)'SUBROUTINE BC4 called in error with III=2'
      STOP
      RETURN
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Nonzero velocity BC's are assigned in this section.  Two statements
C   of the following form must appear for each nonzero velocity node:
C      UBAR(I)=XVEL
C      VBAR(I)=YVEL
C   I is the horizontal node number, (XVEL,YVEL) are the
C   values of the (X,Y) velocities for node I.  
C A return statement must follow the last velocity BC assignment.
C
    6 CONTINUE
      WRITE(*,*)'SUBROUTINE BC4 called in error with III=3'
      STOP
      RETURN
C
C End of routine
C
      END
