c
c  Program detide_ADCIRC_elev.f
c    uses SUBROUTINE GDAY from nmlpaks
c    input: tides_only_bc.out  
c           fort.61            
c
c   subtracts the tides from the wind_and_tide solution
c
c    outputs: forecast detided elevations bcs for bank150
c             to a file
c
c  Cisco Werner and Rick Luettich (9 April 1999)
c  Last revised:
c  Cisco, Dennis and Craig (14 April 1999)
c  Cisco (23 April 1999)
c  Dennis (5 May 1999)
c
      character*80 title,header
c
c  arrays dimensioned exactly to # of bank150 bcs
c
      parameter(nmax=367)
      integer nnums(nmax)
      real tide(nmax),windandtide(nmax)
      character*80 infile,outfile1
c
      iyear=1999
      nbcsbank150 = 116

      open(unit=87,
     .     file='/usr/people/mcgillic/globec/MESH/bank150.elev.ncd',
     .     status='old')

      do 87 i=1,nbcsbank150
         read(87,*)nnums(i)
  87   continue

c
c  input file with wind plus tides elev values on bank150 bc nodes
c
      write(*,*)' enter name of fort.61 filename'
      read(*,80)infile
      open(2,file=infile,status='old')
      read(2,80)title
      read(2,*)nhoursw,nbcs,dt,nout,ndum

      if (nbcs .ne. nbcsbank150) then
         write(*,*)'Wrong number of bcs',nbcs,nbcsbank150
         stop
       endif
c
c
c  from wind+tides file read:
c    sec1Marw   = seconds since 1 March 1999, 0000Z
c    niter1Marw = number of iterations since 1 March 1999, 0000Z
c
      read(2,*)seckd,niter1Marw
      id=1
      im=3
      iyr=1999
      call gday(id,im,iyr,KD0)
      ndays=int(seckd/86400.)
      nsec=seckd-86400*ndays
      kd=kd0+ndays
      call dmy(id,im,iyr,kd)
      write(*,*)' start of forecast file:'
      write(*,*)' year    = ',iyr
      write(*,*)' month   = ',im
      write(*,*)' day     = ',id
      write(*,*)' seconds = ',nsec
      write(*,*)' '
c
c input file with tides only elevation values on bank150 bc nodes
c
      open(1,file='tides_only_bc.out',status='old')
      read(1,80)title
      read(1,*)nhours,nbcs,dt,nout,ndum
c
c  nhours = number of hourly output sets of elev bcs
c  nbcs   = number of boundary nodes on bank150
c  dt     = time interval between each set of elev bcs (3600 secs)
c  nout   = number of ADCIRC time steps between outputs (120)
c  ndum   = dummy 
c
c  open output file with detided elev values on bank150 bc nodes
c     and write header information, etc.
c
      write(*,*)' enter name of output bank150 detided .cbc filename'
      read(*,80)outfile1


c
c   align the tides_only and the tides_plus_winds files
c     
      do i=1,10000
        read(1,*)sec1Mart,niter1Mart
        if(sec1Mart.eq.seckd)then
           write(*,*)'sec1Mart and seckd= ',sec1Mart,seckd
           go to 10
        end if
        do ibcs=1,nbcs
	   read(1,*)l,tide(l)
        end do
      end do
  10  backspace(1)
      backspace(2)
c
c  detide the elevation data
c

        nhoursw = 0

 69     read(1,*)sec1Mart,niter1Mart

        nhoursw = nhoursw + 1
        write(*,*)nhoursw

c
c  sec1Mar   = seconds since 1 March 1999, 0000Z
c  niter1Mar = number of iterations since 1 March 1999, 0000Z
c
        read(2,*,end=70)sec1Mar,niter1Mar
	checkt=sec1Mar-sec1Mart
	if(checkt.ne.0.)then
	   write(*,*)' start times in the files dont line up -- '
           write(*,*)sec1Mar,sec1Mart
	   write(*,*)' program halted'
	   stop
	end if
	do ibcs=1,nbcs
	   read(1,*)l,tide(l)
	   read(2,*)l,windandtide(l)
           write(33,101)nhoursw,ibcs,nnums(ibcs),windandtide(l)-tide(l)
	end do

      goto 69

  70  continue
c
      open(3,file=outfile1,status='unknown')
      write(3,*)'bank150'
      header = 'detided ADCIRC cbc file'
      write(3,*)header
      write(3,99)id,im,iyr,nsec   ! FEW changed from float(nsec) to nsec
      write(3,100)nhoursw-1,dt
      write(3,100)nbcsbank150

      rewind(33)

 71   read(33,101,end=72)it,ibcs,nnums2,windandtide2
      write(3,101)it,ibcs,nnums2,windandtide2
      goto 71

 72   continue
c
      close(1)
      close(2)
      close(3)
  80  format(a)
  99  format(5x,i2,1x,i2,1x,i4,3x,i6)
 100  format(5x,i8,1x,f8.3)
 101  format(3(1x,i6),1x,e16.8)
      stop
      end
C*********************************************************************
C*********************************************************************
      SUBROUTINE GDAY(IDD,IMM,IYear,KD)
C!
C!  GIVEN DAY,MONTH,(EACH 2 DIGITS) and year (four digits), GDAY RETURNS
C!  THE DAY#, KD BASED ON THE GREGORIAN CALENDAR.
C!  THE GREGORIAN CALENDAR, CURRENTLY 'UNIVERSALLY' IN USE WAS
C!  INITIATED IN EUROPE IN THE SIXTEENTH CENTURY. NOTE THAT GDAY
C!  IS VALID ONLY FOR GREGORIAN CALENDAR DATES.
C
C   KD=1 CORRESPONDS TO JANUARY 1, 0000
C	
c 	Note that the Gregorian reform of the Julian calendar 
c	omitted 10 days in 1582 in order to restore the date
c	of the vernal equinox to March 21 (the day after
c	Oct 4, 1582 became Oct 15, 1582), and revised the leap 
c	year rule so that centurial years not divisible by 400
c	were not leap years.
c
C   THIS ROUTINE WAS WRITTEN BY EUGENE NEUFELD, AT IOS, IN JUNE 1990.
C
      INTEGER NDP(13)
      INTEGER NDM(12)
	
      DATA NDP/0,31,59,90,120,151,181,212,243,273,304,334,365/
      DATA NDM/31,28,31,30,31,30,31,31,30,31,30,31/
C!
	lp=6
c make iyy and icc variables
        icc=iyear/100
        iyy=iyear-icc*100

C!  TEST FOR INVALID INPUT:
      IF(ICC.LT.0)THEN
	 WRITE(LP,5000)ICC
	 STOP
      ENDIF
      IF(IYY.LT.0.OR.IYY.GT.99)THEN
	 WRITE(LP,5010)IYY
	 STOP
      ENDIF
      IF(IMM.LE.0.OR.IMM.GT.12)THEN
	 WRITE(LP,5020)IMM
	 STOP
      ENDIF
      IF(IDD.LE.0)THEN
	 WRITE(LP,5030)IDD
	 STOP
      ENDIF
      IF(IMM.NE.2.AND.IDD.GT.NDM(IMM))THEN
	 WRITE(LP,5030)IDD
	 STOP
      ENDIF
      IF(IMM.EQ.2.AND.IDD.GT.29)THEN
	 WRITE(LP,5030)IDD
	 STOP
      ENDIF
      IF(IMM.EQ.2.AND.IDD.GT.28.AND.((IYY/4)*4-IYY.NE.0.OR.(IYY.EQ.0.AND
     .    .(ICC/4)*4-ICC.NE.0)))THEN
	 WRITE(LP,5030)IDD
	 STOP
      ENDIF
5000  FORMAT(' INPUT ERROR. ICC = ',I7)
5010  FORMAT(' INPUT ERROR. IYY = ',I7)
5020  FORMAT(' INPUT ERROR. IMM = ',I7)
5030  FORMAT(' INPUT ERROR. IDD = ',I7)
C!
C!  CALCULATE DAY# OF LAST DAY OF LAST CENTURY:
      KD = ICC*36524 + (ICC+3)/4
C!
C!  CALCULATE DAY# OF LAST DAY OF LAST YEAR:
      KD = KD + IYY*365 + (IYY+3)/4
C!
C!  ADJUST FOR CENTURY RULE:
C!  (VIZ. NO LEAP-YEARS ON CENTURYS EXCEPT WHEN THE 2-DIGIT
C!  CENTURY IS DIVISIBLE BY 4.)
      IF(IYY.GT.0.AND.(ICC-(ICC/4)*4).NE.0) KD=KD-1
C!  KD NOW TRULY REPRESENTS THE DAY# OF THE LAST DAY OF LAST YEAR.
C!
C!  CALCULATE DAY# OF LAST DAY OF LAST MONTH:
      KD = KD + NDP(IMM)
C!
C!  ADJUST FOR LEAP YEARS:
      IF(IMM.GT.2.AND.((IYY/4)*4-IYY).EQ.0.AND.((IYY.NE.0).OR.
     .   (((ICC/4)*4-ICC).EQ.0)))   KD=KD+1
C!  KD NOW TRULY REPRESENTS THE DAY# OF THE LAST DAY OF THE LAST
C!  MONTH.
C!
C!  CALCULATE THE CURRENT DAY#:
      KD = KD + IDD
      RETURN
C!
C!
      ENTRY DMY(IDD,IMM,IYear,KD)
C!
C!  GIVEN THE (GREGORIAN) DAY#, KD, AS CALCULATED ABOVE IN THIS ROUTINE,
C!  ENTRY DMY RETURNS THE (GREGORIAN) DAY, MONTH, YEAR AND CENTURY.
C!
C!  TEST FOR VALID INPUT:
      IF(KD.LE.0) WRITE(LP,5040)KD
5040  FORMAT(' KD = ',I7,'  INVALID INPUT. DMY STOP.')
C!
C!  SAVE KD
      KKD=KD
C!  CALCULATE ICC AND SUBTRACT THE NUMBER OF DAYS REPRESENTED BY ICC
C!  FROM KKD
C!  JFH IS THE NUMBER OF 400 YEAR INTERVALS UP TO KKD
C!  JCC IS THE NUMBER OF ADDITIONAL CENTURIES UP TO KKD
      JFH = KKD/146097
      KKD = KKD - JFH*146097
      IF(KKD.LT.36525)THEN
	 JCC = 0
      ELSE
	 KKD = KKD - 36525
	 JCC = 1 + KKD/36524
	 KKD = KKD - (JCC-1)*36524
      END IF
      ICC = 4*JFH + JCC
      IF(KKD.EQ.0)THEN
	 ICC = ICC-1
	 IYY = 99
	 IMM = 12
	 IDD = 31
c	 RETURN
	go to 110
      ENDIF
C!
C!  CALCULATE IYY. JFY IS THE NUMBER OF FOUR YEAR INTERVALS IN THE
C!  CURRENT CENTURY. THE FIRST FOUR YEAR INTERVAL IS SHORT (1460 DAYS
C!  RATHER THAN 1461)IF THE CURRENT CENTURY IS NOT DIVISIBLE BY 4, AND
C!  IN THIS CASE JCC.NE.0 AS CALCULATED ABOVE.
C!
C!  CALCULATE JFY:
      JFY = 0
      IF(JCC.EQ.0)GOTO 10
      IF(KKD.LT.1460)GOTO 10
      JFY = 1
      KKD = KKD - 1460
10    KK = KKD/1461
      JFY = JFY + KK
      KKD = KKD - KK*1461
C!
C!  CALCULATE JYY, THE REMAINING YEARS OF THE CURRENT CENTURY UP TO THE
C!  CURRENT DAY:
      JYY = 0
C!  THE NEXT YEAR IS NOT A LEAP YEAR IF JFY=0 AND JCC.NE.0.
      IF(JFY.EQ.0.AND.JCC.NE.0)GOTO 20
      IF(KKD.LT.366)GOTO 30
      JYY = 1
      KKD = KKD - 366
20    JYYY = KKD/365
      JYY = JYY + JYYY
      KKD = KKD - JYYY*365
30    IYY = 4*JFY + JYY
      IF(KKD.EQ.0) THEN
	 IYY=IYY-1
	 IMM=12
	 IDD=31
c	 RETURN
	go to 110
      END IF
C!
C!  SET L=1 IF WE HAVE A LEAP YEAR.
      L=0
      IF(IYY-(IYY/4)*4.NE.0)GOTO 40
      IF(IYY.EQ.0.AND.(ICC-(ICC/4)*4).NE.0)GOTO 40
      L=1
C!
C!  CALCULATE IMM AND IDD
40    IF(KKD.GT.31) GOTO 50
      IMM=1
      IDD=KKD
c      RETURN
       go to 110
C!
50    IF(KKD.GT.59)GOTO 60
      IMM = 2
      IDD = KKD-31
c      RETURN
	go to 110
C!
60    IF(KKD.GT.60)GOTO 70
      IF(L.EQ.0)GOTO 70
      IMM = 2
      IDD = 29
c      RETURN
	go to 110
C!
70    IF(L.EQ.1) KKD=KKD-1
      DO 80 I=4,13
	 IF(KKD.GT.NDP(I))GOTO 80
	 IMM = I-1
	 IDD = KKD - NDP(I-1)
c	 RETURN
	go to 110
C!
80    CONTINUE

 110 	iyear=icc*100+iyy
	return

90    WRITE(LP,5050)
5050  FORMAT(' ERROR IN DMY.')
      STOP
      END

	Subroutine Up_date(kd,ssec)
C 
C this subroutine updates the current gregorian day number kd
C and resets the time, ssec, accordingly
C it inherently assumes that the variable is updated every day!

      integer kd
      real ssec
	
      kd=kd+aint(ssec/8.64E4)
      ssec=amod(ssec,8.64E4)
 
      return
      end
