C***********************************************************************
C bcs subroutine for the SACO 
C***********************************************************************
      SUBROUTINE BCS2(KD,SECKD,ITER,DEGLAT,UBAR,VBAR,III)
C Include file for parameter statements
      INCLUDE '../CASCO.DIM'
      include '../startstop.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.

	real    bt1, bt2 
        integer i, p1, p2
	integer nbc,nodeii(nndim)

	real    bc_space(nrdim)
	character*120 filename
	character*80  header
        save    bc_space
 	real    dseconds
 	integer dday
	integer nramp
	real    rampup
	integer ii
	real    zres
	real    pi 
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 of bcs for saco' 
       write(*,*)'Start time ', kd_saco(1), sec_saco(1)
       write(*,*)'Final time ', kd_saco(nsteps), sec_saco(nsteps)

       write(*,*)'Enter name of residual bcs file'
       read(*,'(a)')filename 
       write(*,*)filename 
       open(1,file=filename,status='old')
       read(1,'(a)') header
       write(*,'(a)') header
       read(1,'(a)') header
       write(*,'(a)') header
       do i = 1, nr
	  read(1,*)ii, zres
          bc_space(i) = zres
	  nodeii(i)=ii
       enddo
       do i = 1, nr
	  write(*,*)i, lrvs(i), bc_space(i)
       enddo
       close(1) 
       
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  Create test bcs. 
C  use spatial structure from above and ramp up over first 4 days
C 
       call gday_diff(kd,seckd,kd_saco(1),sec_saco(1), dday,dseconds)
       nramp = 3
       if(dday.lt.nramp) then 
         rampup = (dday + dseconds/(24.*3600))/ float(nramp) 
       else
	  rampup = 1.0 
       endif
       write(*,*)'BCS2: ', iter, kd,seckd, rampup

      DO I=1,NR 
c         NBC = LRVS(I)
c         UBAR(NBC) = rampup*bc_space(i) 
         NBC = nodeii(i)
         UBAR(NBC) = rampup*bc_space(i) 
c         write(*,*)i,nbc,nodeii(i),ubar(nbc)
c	 if(i.eq.2) write(*,*)nbc,rampup, bc_space(i), ubar(nbc)
      ENDDO

      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
