      PROGRAM NCOM_TIDEBODY
      IMPLICIT NONE 
!DAN======================================================================  
c
C  NCOM_tidebody - Usage:  
C            NCOM_tidebody tidcon wstart wend whrinc itest jtest [grid.a] 
C                 outputs tidal body forcing (m) every whrinc hours: 
C                 from wind day wstart 
C                 to wind day    wend
C
C                 tidcon 1 digit per constituent (Q1K2P1N2O1K1S2M2), 0=off,1=on
C                 itest & jtest are grid points for a trace:
C                 if  1 <= itest <=idm  and 1 <= jtest <= jdm then a file:
c                      NCOM_bodytide_trace.txt  will be written which prints 
c                        the values of the sum and all 8 (filtered by tidcon)
c                        principal tidal components for each time step.
c                 else
c                      NCOM_bodytide_trace.txt   just contains some setup
c                               diagnostic information
c                 endif 
C
C                 grid.a is a NCOM grid file, [default regional.grid.a].
C                 Note that the corresponding grid.b must also exist.
C
C                 idm,jdm are taken from grid.a
c
c          Output file:  NCOM_tideforce.a is written (unformatted)
c    Standard Output: Can generate a NCOM_tideforce.b  file (formatted)
C
C  this version for "serial" Unix systems.
C
C  Dan Moore  QinetiQ NA - September 2011
C
!DAN============================================================================
      REAL*4,  ALLOCATABLE :: F(:,:,:),PLAT(:,:),PLON(:,:),AMASK(:,:),
     &                        PAD(:)
      REAL*4               :: xmin,xmax
      INTEGER IOS,IForce_File_Number,i,j,k,ipt
      INTEGER      IARGC,itest,jtest
      INTEGER      NARG,n2pad
      CHARACTER*240 CARG,CFILE
      CHARACTER*2   TideMode(8)
      CHARACTER*24  Tides
      DATA TideMode/'M2','S2','K1','O1','N2','P1','K2','Q1'/
C
      INTEGER       IDM,JDM,NPAD,TIDCON,NRECL,TIDCON1,n2drec
      REAL*8        wstart,wstop,whrinc,TT
      CHARACTER*6   CVARIN
      CHARACTER*240 CFILEB,CFILEG
      CHARACTER*80  CLINE
      LOGICAL tide_on(8),Print_Trace
C
C     READ ARGUMENTS.
C
      NARG = IARGC()
C
      IF     (NARG.EQ.6) THEN
        CALL GETARG(1,CARG)
        READ(CARG,*) tidcon
        CALL GETARG(2,CARG)
        READ(CARG,*) wstart
        CALL GETARG(3,CARG)
        READ(CARG,*) wstop
        CALL GETARG(4,CARG)
        READ(CARG,*) whrinc
        CALL GETARG(5,CARG)
        READ(CARG,*) itest
        CALL GETARG(6,CARG)
        READ(CARG,*) jtest
        CFILEG = 'regional.grid.a'
      ELSEIF     (NARG.EQ.7) THEN
        CALL GETARG(1,CARG)
        READ(CARG,*) tidcon
        CALL GETARG(2,CARG)
        READ(CARG,*) wstart
        CALL GETARG(3,CARG)
        READ(CARG,*) wstop
        CALL GETARG(4,CARG)
        READ(CARG,*) whrinc
        CALL GETARG(5,CARG)
        READ(CARG,*) itest
        CALL GETARG(6,CARG)
        READ(CARG,*) jtest
        CALL GETARG(7,CFILEG)
      ELSE
        WRITE(6,*) 
     +   'Usage:  NCOM_tidebody tidcon wstart wend whrinc  [grid.a] '
        CALL EXIT(1)
      ENDIF
      OPEN(UNIT=66,FILE='NCOM_tidebody_trace.txt',FORM='FORMATTED',
     &     ACTION='WRITE')
      WRITE(66,*)'Argument List processed:'
      WRITE(66,*)' tidcon = ',tidcon
      WRITE(66,*)' wstart = ',wstart
      WRITE(66,*)' wstop  = ',wstop
      WRITE(66,*)' whrinc = ',whrinc
      WRITE(66,*)' itest = ',itest
      WRITE(66,*)' jtest = ',jtest
      WRITE(66,'(a,a)')'Regional Data File = ',TRIM(CFILEG)
C
C     GET IDM,JDM FROM regional.grid.b.
C
      CFILEB = CFILEG(1:LEN_TRIM(CFILEG)-1) // 'b'

      WRITE(66,'(a,a)')' Grid data file = ',TRIM(CFILEB)
C
      OPEN(UNIT=11,FILE=CFILEB,FORM='FORMATTED',
     &     STATUS='OLD',ACTION='READ')
C
      READ( 11,*) IDM,CVARIN
      IF (CVARIN.NE.'idm   ') THEN
        WRITE(6,*) 'NCOM_tidelat: bad header file ',
     &             CFILEB(1:LEN_TRIM(CFILEB))
        WRITE(66,*) 'NCOM_tidelat: bad header file ',
     &             CFILEB(1:LEN_TRIM(CFILEB))
        CALL EXIT(2)
      ENDIF
      READ( 11,*) JDM,CVARIN
      IF (CVARIN.NE.'jdm   ') THEN
        WRITE(6,*) 'NCOM_tidelat: bad header file ',
     &             CFILEB(1:LEN_TRIM(CFILEB))
        WRITE(66,*) 'NCOM_tidelat: bad header file ',
     &             CFILEB(1:LEN_TRIM(CFILEB))
        CALL EXIT(2)
      ENDIF
      READ(11,'(a80)')CLINE
C
      CLOSE(UNIT=11)
      write(66,*)' IDM  = ',IDM
      write(66,*)' JDM  = ',JDM
  116  format (
     & i5,4x,'''idm   '' = longitudinal array size'/
     & i5,4x,'''jdm   '' = latitudinal  array size'/
     & a70)
       write(6 ,116)IDM,JDM,trim(CLINE)
       write(66,116)IDM,JDM,trim(CLINE)

C
      ALLOCATE( F(IDM,JDM,9), STAT=IOS )
      IF     (IOS.NE.0) THEN
        WRITE(6,*) 'Error in NCOM_tidebody: could not allocate ',
     +             IDM*JDM,' words for F'
        WRITE(66,*) 'Error in NCOM_tidebody: could not allocate ',
     +             IDM*JDM,' words for F'
        CALL EXIT(2)
      ENDIF
      ALLOCATE( PLAT(IDM,JDM), STAT=IOS )
      IF     (IOS.NE.0) THEN
        WRITE(6,*) 'Error in NCOM_tidebody: could not allocate ',
     +             IDM*JDM,' words for PLAT'
        WRITE(66,*) 'Error in NCOM_tidebody: could not allocate ',
     +             IDM*JDM,' words for PLAT'
        CALL EXIT(2)
      ENDIF

      ALLOCATE( PLON(IDM,JDM), STAT=IOS )
      IF     (IOS.NE.0) THEN
        WRITE(6,*) 'Error in NCOM_tidebody: could not allocate ',
     +             IDM*JDM,' words for PLON'
        WRITE(66,*) 'Error in NCOM_tidebody: could not allocate ',
     +             IDM*JDM,' words for PLON'
        CALL EXIT(2)
      ENDIF

      ALLOCATE( AMASK(IDM,JDM), STAT=IOS )
      IF     (IOS.NE.0) THEN
        WRITE(6,*) 'Error in NCOM_tidebody: could not allocate ',
     +             IDM*JDM,' words for AMASK'
        WRITE(66,*) 'Error in NCOM_tidebody: could not allocate ',
     +             IDM*JDM,' words for AMASK'
        CALL EXIT(2)
      ENDIF
C----------------------------------------------------------------
C
C     INPUT PLAT, PLON ARRAYS.
C
      
C
#ifdef CRAY
#ifdef t3e
      IF     (MOD(NRECL,4096).EQ.0) THEN
        WRITE(CASN,8000) NRECL/4096
 8000   FORMAT('-F cachea:',I4.4,':1:0')
        IU8 = 11
        CALL ASNUNIT(IU8,CASN,IOS8)
        IF     (IOS8.NE.0) THEN
          write(6,*) 'Error: can''t asnunit 11'
          write(6,*) 'ios  = ',ios8
          write(6,*) 'casn = ',casn
          write(66,*) 'Error: can''t asnunit 11'
          write(66,*) 'ios  = ',ios8
          write(66,*) 'casn = ',casn
          CALL EXIT(5)
        ENDIF
        IU8 = 21
        CALL ASNUNIT(IU8,CASN,IOS8)
        IF     (IOS8.NE.0) THEN
          write(6,*) 'Error: can''t asnunit 21'
          write(6,*) 'ios  = ',ios8
          write(6,*) 'casn = ',casn
          write(66,*) 'Error: can''t asnunit 21'
          write(66,*) 'ios  = ',ios8
          write(66,*) 'casn = ',casn
          CALL EXIT(5)
        ENDIF
      ENDIF
#else
      CALL ASNUNIT(11,'-F syscall -N ieee',IOS)
      IF     (IOS.NE.0) THEN
        write(6,*) 'Error: can''t asnunit 11'
        write(6,*) 'ios = ',ios
        write(66,*) 'Error: can''t asnunit 11'
        write(66,*) 'ios = ',ios
        CALL EXIT(5)
      ENDIF
      CALL ASNUNIT(21,'-F syscall -N ieee',IOS)
      IF     (IOS.NE.0) THEN
        write(6,*) 'Error: can''t asnunit 21'
        write(6,*) 'ios = ',ios
        write(66,*) 'Error: can''t asnunit 21'
        write(66,*) 'ios = ',ios
        CALL EXIT(5)
      ENDIF
#endif
#endif
 
      n2drec=((IDM*JDM+4095)/4096)*4096
      NPAD=n2drec-IDM*JDM
      if(NPAD.ne.0)allocate(PAD(NPAD))
c
      NRECL=n2drec
C
      OPEN(UNIT=11, FILE=CFILEG, FORM='UNFORMATTED', STATUS='OLD',
     +         ACCESS='DIRECT', RECL=NRECL, IOSTAT=IOS)
      IF     (IOS.NE.0) THEN
        write(6,*) 'Error: can''t open ',TRIM(CFILEG)
        write(6,*) 'ios   = ',ios
        write(6,*) 'nrecl = ',nrecl
        write(66,*) 'Error: can''t open ',TRIM(CFILEG)
        write(66,*) 'ios   = ',ios
        write(66,*) 'nrecl = ',nrecl
        CALL EXIT(3)
      ENDIF
C
      READ(11,REC=1,IOSTAT=IOS) PLON
#ifdef ENDIAN_IO
      CALL ENDIAN_SWAP(PLON,IDM*JDM)
#endif
      IF     (IOS.NE.0) THEN
        WRITE(6,*) 'can''t read ',TRIM(CFILEG)
        WRITE(66,*) 'can''t read ',TRIM(CFILEG)
        CALL EXIT(4)
      ENDIF
c
      READ(11,REC=2,IOSTAT=IOS) PLAT
#ifdef ENDIAN_IO
      CALL ENDIAN_SWAP(PLAT,IDM*JDM)
#endif
      IF     (IOS.NE.0) THEN
        WRITE(6,*) 'can''t read ',TRIM(CFILEG)
        WRITE(66,*) 'can''t read ',TRIM(CFILEG)
        CALL EXIT(4)
      ENDIF
C
      CLOSE(UNIT=11)
      write(66,*)'Array PLon read'
      write(66,*)'Array PLat read'
C
C
C
C     OUTPUT FILE.
C
      CFILE='NCOM_tide_force.a'
      OPEN(UNIT=21, FILE=CFILE, FORM='UNFORMATTED', STATUS='NEW',
     +         ACCESS='DIRECT', RECL=NRECL, IOSTAT=IOS)
      IF     (IOS.NE.0) THEN
        write(6,*) 'Error: can''t open ',TRIM(CFILE)
        write(6,*) 'ios   = ',ios
        write(6,*) 'nrecl = ',nrecl
        write(66,*) 'Error: can''t open ',TRIM(CFILE)
        write(66,*) 'ios   = ',ios
        write(66,*) 'nrecl = ',nrecl
        CALL EXIT(3)
      ENDIF
C----------------------------------------------------------------
          tidcon1 = tidcon
          do i =1,8
            tide_on(i) = mod(tidcon1,10) .eq. 1
            tidcon1    =     tidcon1/10  ! shift by one decimal digit
          enddo
        
      TIDES='                        ' 
      ipt=1
      do i=1,8
        if(tide_on(i))then
           TIDES(ipt:ipt+1)=TideMode(i)
           ipt=ipt+3
        endif
      end do
c      WRITE(6,'(a,a)')'Tidal Modes included: ',trim(TIDES)
      WRITE(66,'(a,a)')'Tidal Modes included: ',trim(TIDES)
      TT=wstart
      IForce_File_Number=0
      F=0.0
      Print_Trace=(itest.ge.1).and.(itest.le.idm).and.
     &            (jtest.ge.1).and.(jtest.le.jdm)
   10 CALL NCOM_TIDE_FORCE(F,PLAT,PLON,AMASK,TT,IDM,JDM,tide_on)
      if(Print_Trace)
     &write(66,'(a,i4,a,i4,a,3F12.2,9G15.5)')'TIDE_FORCE(',itest,',',
     &jtest,') for T,lat,lon = ',TT,plat(itest,jtest),
     &PLON(itest,jtest),F(itest,jtest,9),(F(itest,jtest,k),k=1,8)
c      do k=1,9
      IForce_File_Number=IForce_File_Number+1
        IF(NPAD.EQ.0) THEN
          WRITE(21,REC=IForce_File_Number,IOSTAT=IOS) F(:,:,9)
        ELSE
          WRITE(21,REC=IForce_File_Number,IOSTAT=IOS) F(:,:,9),PAD
        ENDIF
c      end do
c
c Need special code to only sample over Ocean points?
c Need  Ip array?
c      
      xmin=minval(F(:,:,9))
      xmax=maxval(F(:,:,9))
      write(6,'(F10.3,a,2g15.5)')TT,'  Time: tide force: min,max = '
     & ,xmin,xmax

      TT=TT+whrinc/24.d0
      IF(TT.LE.wstop)GO TO 10

     
      CLOSE(21)
      CLOSE(66)

      CALL EXIT(0)
 5000 FORMAT(I4)
      END  
      SUBROUTINE NCOM_TIDE_FORCE(F,PLAT,PLON,AMASK,TT,
     &IDM,JDM,tide_on)
      IMPLICIT NONE
      REAL*4       F(IDM,JDM,9),PLAT(IDM,JDM),PLON(IDM,JDM),
     &             AMASK(IDM,JDM),dayfrac,ramp,timesec
      REAL*8       TT,TT0      
      CHARACTER*240 CFILEG
      INTEGER      IDM,JDM,tidcon,itest,jtest,icall,idate
      INTEGER      ihh,isec,imm,itime,iruno,iec(8)
      INTEGER, allocatable :: is(:),ie(:)
      LOGICAL tide_on(8)
      SAVE TT0,icall,is,ie
      common/par2o/ iruno( 4,1)
      DATA icall/0/
C
C     MOST OF WORK IS DONE HERE.
C
#ifdef sun
      INTEGER      IR_ISNAN
C
#endif
      CHARACTER*18 CASN
      INTEGER      I,J,K,IOS,NRECL
#ifdef CRAY
      INTEGER*8    IU8,IOS8
#endif

      if(icall.eq.0)then
         allocate(is(1:JDM))
         allocate(ie(1:JDM))


         TT0=TT
c
c     Convert HYCOM start day-time to NCOM Start time format
c
       write(6,*)'HYCOM Start Date:',TT
       call dayceni(0.0,TT,idate,dayfrac)
       write(6,*)'NCOM equivalent idate1:',idate
       write(6,*)'NCOM dayfrac          :',dayfrac
       ihh = int(dayfrac*24.)
       imm = int(dayfrac*1440.) - ihh*60
       isec = nint(dayfrac*86400.) - imm*60 - ihh*3600
       itime = ihh*1000000 + imm*10000 + isec*100
       write(6,*)'ihh,imm,isec,iteme=',ihh,imm,isec,itime
       iruno(1,1)=idate
       iruno(2,1)=itime
      endif
      icall=icall+1
      timesec=(TT-TT0)*DBLE(24*3600)
      ramp=1.0
       do j=1,jdm
         call tidepot(timesec,ramp,idm,jdm,idm,jdm,iec,is,ie,j,
     &          amask,plon,plat,F(1:idm,j,9),F(1:idm,j,3),
     & F(1:idm,j,4),F(1:idm,j,6),F(1:idm,j,8),F(1:idm,j,7),
     & F(1:idm,j,1),F(1:idm,j,5),F(1:idm,j,2))
       end do
       do j=1,jdm
          do i=1,idm
             F(i,j,9)=0
             do k=1,8
                if(tide_on(k))then
                  F(i,j,9)=F(i,j,9)+F(i,j,k)
                else
                  F(i,j,k)=0
                endif
              end do
          end do
        end do


C

#ifdef ENDIAN_IO
      CALL ENDIAN_SWAP(F,IDM*JDM)
#endif
      RETURN
      END
      module NCOM_Tidepot_Mod
      implicit none
!#include "../../include/PARAM.h"c
c  Define halo width (nm0).
      integer    nmh,nm0
      parameter (nmh=0, nm0=1-nmh)
c
c  Define maximum horizontal array dimensions for total grid (nmxa) and
c  for a single processor (nmx), and maximum vertical dimension (lmx).
      integer    nmxa,nmx,lmx
      parameter (nmxa=504, nmx=nmxa, lmx=385)
c
c  Define maximum number of scalar (nrmx) and turbulence (nqmx) fields.
      integer    nrmx,nqmx
      parameter (nrmx=2, nqmx=2)
c
c  Define maximum number of open bndy pts for total grid (nobmxt),
c  max number of open bndy pts for a single tile (nobmx), and max
c  number of tidal constituents (ntcmx).
      integer nobmxt,nobmx,ntcmx
      parameter (nobmxt=8000,nobmx=8000,ntcmx=8)
c
c  Define maximum number of rivers (nrivmx).
      integer nrivmx
      parameter (nrivmx=2000)
c
c  Define maximum number of nested grids (mxgrdso).
      integer    mxgrdso
      parameter (mxgrdso=7)
c
c  Define maximum number of individual model grid points at which
c  output data can be saved.  Note that this is limited by the 
c  available output unit numbers (see subroutine outpt).
c  Currently, unit numbers 61 to 98 are reserved for this.
      integer nsavmx
      parameter (nsavmx=40)
c
c  Define named integer constants (used by subroutine xctilr).
      integer    ityphs,itypss,itypus,itypvs
      integer    ityphv,itypsv,itypuv,itypvv
      parameter (ityphs= 1, itypss= 2, itypus= 3, itypvs= 4)
      parameter (ityphv=11, itypsv=12, itypuv=13, itypvv=14)
c
c  Define named real constants
c
      real       zero,zp5,one,two
      parameter (zero=0.0, zp5=0.5, one=1.0, two=2.0)
      private

      type t_tidenest
        real, pointer :: p(:,:,:)
      end type t_tidenest
      public :: t_tidenest

      integer, save :: nc(mxgrdso), itspec(ntcmx,mxgrdso)
      real   , save :: tidefq(ntcmx,mxgrdso)
      real   , save :: tideph(ntcmx,mxgrdso), amp(ntcmx,mxgrdso)
      type(t_tidenest), save ::
     &  tplat(mxgrdso), coslon(mxgrdso), sinlon(mxgrdso)

      public :: nc, itspec, tidefq, tideph, amp
      public :: tplat, coslon, sinlon

      end module NCOM_Tidepot_Mod
c=======================================================================
      subroutine tidepot_init(nt,mt,n,m,is,ie,j,amsk,elon,alat,ep)
c-----------------------------------------------------------------------
c  subroutine to allocate memory and compute long-lat-dependent
c  tidal potential functions.  Must be called only once per nest.
c  created 2007-06-11, Tim Campbell, NRL. (adapted from Paul Martin)
c  last updated 2007-06-11, Tim Campbell, NRL.
c-----------------------------------------------------------------------
c       idate1 = 1st date of form YYYYMMDD (1901 <= YYYY <= 2099).
c       itime1 = 1st time of form HHMMSSCC (CC is hundredths of sec).
c       idate2 = 2nd date of form YYYYMMDD (1901 <= YYYY <= 2099).
c       itime2 = 2nd time of form HHMMSSCC (CC is hundredths of sec).
c
c  returned arguments:
c       iday   = number of days.
c       ihr    = number of hours.
c       min    = number of minutes.
c       isec   = number of seconds.
c       ihsec  = number of hundredths of sec.
c
      use NCOM_Tidepot_Mod
c
      implicit none
!#include "../../include/PARAM.h"
!#include "../../include/COMMON.h"
c

c
c  Define named real constants
      real       zero,zp5,one,two
      parameter (zero=0.0, zp5=0.5, one=1.0, two=2.0)
      
      integer nesto,ntcmx
      PARAMETER(nesto=1,ntcmx=8)
           
      integer io_unit_offset, istdo_unit, iruno
      common/par2o/ iruno( 4,1)
      data istdo_unit/6/
      
c  Define halo width (nm0).
      integer    nmh,nm0
      parameter (nmh=0, nm0=1-nmh)
c
c  Define maximum horizontal array dimensions for total grid (nmxa) and
c  for a single processor (nmx), and maximum vertical dimension (lmx).
      integer    nmxa,nmx,lmx
      parameter (nmxa=504, nmx=nmxa, lmx=385)
c
c  Declare passed variables.
      integer nt,mt,n,m,j
      integer is(nm0:m+nmh),ie(nm0:m+nmh)
      real    amsk(nm0:n+nmh,nm0:m+nmh)
      real    elon(nm0:n+nmh,nm0:m+nmh)
      real    alat(nm0:n+nmh,nm0:m+nmh),ep(nm0:n)
c
c  Declare local temporary variables.
      character*5 tidecn(ntcmx)
      integer i,jj,k,ic,idate,itime,iyear,month,iday,ihr
      real    tidenf(ntcmx)
      real    pi,raddeg,alatave,a,b
c
c  define/set some parameters.
c      idate=iruno(1,nesto)
c      itime=iruno(2,nesto)
      idate=20110101
      itime=18000000
      pi=3.1415926535
      raddeg=pi/180.
c
c  allocate storage for spatial dependent functions
      if (.not.associated(tplat(nesto)%p)) then
        allocate(  tplat(nesto)%p(nm0:n+nmh,nm0:m+nmh,0:2) )
        allocate( coslon(nesto)%p(nm0:n+nmh,nm0:m+nmh,1:2) )
        allocate( sinlon(nesto)%p(nm0:n+nmh,nm0:m+nmh,1:2) )
      else
        call xcstop('Error in tidepot_init: only call once per nest')
      endif
c
c  define tidal constituents for which tidal potential is to be calculated.
c  note that the tidal constituents for which the tidal potential is
c  calculated may be different from the tidal constituents for which the
c  open bc are calculated.
c  get list of tidal constituents from input file.
c      call rw_tpcn(1,nesto,nc(nesto),tidecn)
      nc(1)=8
      tidecn(1)='K1'
      tidecn(2)='O1'
      tidecn(3)='P1'
      tidecn(4)='Q1'
      tidecn(5)='K2'
      tidecn(6)='M2'
      tidecn(7)='N2'
      tidecn(8)='S2'
c
c  define tidal species (i.e., declinational k=0, diurnal k=1,
c  semi-diurnal k=2):
      do ic=1,nc(nesto)
        if (tidecn(ic)(1:2) .eq. 'MF' .or.
     &      tidecn(ic)(1:2) .eq. 'MM') then
          itspec(ic,nesto)=0
        else if (tidecn(ic)(1:2) .eq. 'K1' .or. 
     &           tidecn(ic)(1:2) .eq. 'O1' .or. 
     &           tidecn(ic)(1:2) .eq. 'P1' .or. 
     &           tidecn(ic)(1:2) .eq. 'Q1') then
          itspec(ic,nesto)=1
        else if (tidecn(ic)(1:2) .eq. 'K2' .or. 
     &           tidecn(ic)(1:2) .eq. 'M2' .or. 
     &           tidecn(ic)(1:2) .eq. 'N2' .or. 
     &           tidecn(ic)(1:2) .eq. 'S2') then
          itspec(ic,nesto)=2
        else
c          if (mnproc.eq.1) then
            write(istdo_unit,'(a,i2,a,i3,a,a)')
     &        'tidepot_init:  nesto=',nesto,'  ic =',ic,
     &        '  tidecn = ',tidecn(ic)
c            call zhflsh(istdo_unit)
c          endif
         call xcstop('Error in tidepot_init: invalid tidal constituent')
        endif
      enddo
c
c  get tidal frequencies and tidal node factors for date and hr.
c  note that these are only computed to nearest hr.!!!!!!!!
c  Hence, run must be started on the hour (need to fix this).!!!!!!!!!
c      call idt2ymd(idate,iyear,month,iday)
      iyear=2011
      month=1
      iday=1
      ihr=itime/1000000
c      if (mnproc .eq. 1) then
        write(istdo_unit,'(/a,i5,a,i3,a,i3,a,i3)')
     &    'tidepot_init:  year=',iyear,
     &    '  month=',month,'  day=',iday,'  hr=',ihr
c        call zhflsh(istdo_unit)
c      endif
c      call xceget(alatave,alat,n,m, (nt+1)/2, (mt+1)/2)
      alatave=alat((mt+1)/2,(nt+1)/2)
      call tide_fac(ihr,iday,month,iyear,alatave,nc(nesto),0,
     &  tidecn,tidefq(1,nesto),tidenf,tideph(1,nesto))
c
c  if the year is greater than 1901, use tidal amplitude and phase
c  adjustments appropriate for the current date and time that were
c  computed above.  Otherwise, use equilibrium tidal forcing, i.e.,
c  with no adjustments.
c  check that run starts on hour if amp and phase are adjusted.
      if (idate .lt. 19020101) then
        do ic=1,nc(nesto)
          tidenf(ic)=one
          tideph(ic,nesto)=zero
        enddo
      else
        if (itime-1000000*ihr .ne. 0) then
c          if (mnproc.eq.1) then
            write(istdo_unit,'(2(a,i8))')
     &        'tidepot_init:  idate = ',idate,'  itime =',itime
c            call zhflsh(istdo_unit)
c          endif
          call xcstop('Error in tidepot_init:  run not started on hour')
        endif
      endif
c
c  define tidal amplitude for each constituent.
c  correct amplitude for earth-tide (0.69) and node factor (tidenf).
      do ic=1,nc(nesto)
        call tc_amp(tidecn(ic),amp(ic,nesto))
        amp(ic,nesto)=0.69*tidenf(ic)*amp(ic,nesto)
      enddo
c
c  print out tidal parameters.
c      if (mnproc .eq. 1) then
        write(istdo_unit,'(/a,a,i4,a,i2,a,i2,a,i2,a,f5.0)')
     &    'tidepot_init:  parameters for tidal potential for',
     &    ' year=',iyear,'  month=',month,'  day=',iday,'  hour=',ihr,
     &    '  mean lat=',alatave
        write(istdo_unit,'(2a)') ' index constit species  period (h)',
     &    '     amp (m)     nodefac phase (deg)'
        do ic=1,nc(nesto)
          a=two*pi/(3600.0*tidefq(ic,nesto))
          write(istdo_unit,'(i6,3x,a5,i8,f12.4,f12.4,f12.4,f12.2)')
     &      ic,tidecn(ic),itspec(ic,nesto),a,amp(ic,nesto),
     &      tidenf(ic),tideph(ic,nesto)
        enddo
c        call zhflsh(istdo_unit)
c      endif
c
c  convert tidal phase from deg to radians.
      do ic=1,nc(nesto)
        tideph(ic,nesto)=tideph(ic,nesto)*raddeg
      enddo
c
c  define latitude- and longitude-dependent functions.
c  these depend on the tidal "species", i.e., declinational (k=0),
c  diurnal (k=1), or semi-diurnal (k=2).
      a=two*raddeg
      do jj=nm0,m+nmh
        do i=nm0,n+nmh
           tplat(nesto)%p(i,jj,0)=3.0*sin(raddeg*alat(i,jj))-one
           tplat(nesto)%p(i,jj,1)=sin(a*alat(i,jj))
           tplat(nesto)%p(i,jj,2)=cos(raddeg*alat(i,jj))**2
          coslon(nesto)%p(i,jj,1)=cos(elon(i,jj)*raddeg)
          sinlon(nesto)%p(i,jj,1)=sin(elon(i,jj)*raddeg)
          coslon(nesto)%p(i,jj,2)=cos(elon(i,jj)*a)
          sinlon(nesto)%p(i,jj,2)=sin(elon(i,jj)*a)
        enddo
      enddo
c      write(6,*)'191:tplat(1)%p(1,1,1)=',tplat(1)%p(1,1,1)
c      write(6,*)'191:coslon(1)%p(1,1,1)=',coslon(1)%p(1,1,1)
c      write(6,*)'191:sinlon(1)%p(1,1,1)=',sinlon(1)%p(1,1,1)
c      write(6,*)'191:tplat(1)%p(1,1,2)=',tplat(1)%p(1,1,2)
c      write(6,*)'191:coslon(1)%p(1,1,2)=',coslon(1)%p(1,1,2)
c      write(6,*)'191:sinlon(1)%p(1,1,2)=',sinlon(1)%p(1,1,2)
c
      return
      end subroutine tidepot_init

      subroutine tidepot(times,ramp,nt,mt,n,m,iec,is,ie,j,amsk,elon,
     &  alat,ep,ep1,ep2,ep3,ep4,ep5,ep6,ep7,ep8)   
c-----------------------------------------------------------------------
c  subroutine to calculate tidal potential (ep).
c  set up only for K1 O1 P1 Q1 K2 M2 N2 S2 MF MM constituents.!!!!
c  set up for max of nesto = 7.!!!!
c       times = model time in seconds.
c       ramp  = ramp to turn on forcing slowly.
c       nt,mt = global grid dimensions.
c       n,m   = local grid dimensions (i.e., on local processor).
c       iec   = integer flags denoting tile bndys as interior/exterior.
c       is,ie = shrink-wrap indexes.
c       j     = index of x-z section to be calculated.
c       amsk  = land-sea mask at elevation pts.
c       elon  = longitude in degrees E.
c       alat  = latitude in degrees N.
c       ep    = returned values of tidal potental (m).
c  created 2003-08-18, Paul J Martin, NRL.
c  modified 2007-05-29 to compute tidal potential for more than 1 grid.
c  adapted into ncom_tide_mod module 2007-06-11, Tim Campbell, NRL.
c  last updated 2007-06-11, Tim Campbell, NRL.
c-----------------------------------------------------------------------
c
      use NCOM_Tidepot_Mod
c
      implicit none
!#include "../../include/PARAM.h"
!#include "../../include/COMMON.h"
c
c  Define named real constants
      real       zero,zp5,one,two
      parameter (zero=0.0, zp5=0.5, one=1.0, two=2.0)
      
      integer nesto,ntcmx,mxgrdso
      PARAMETER(nesto=1,ntcmx=8,mxgrdso=1)
      
      integer io_unit_offset, istdo_unit
      data istdo_unit/6/
c  Define halo width (nm0).
      integer    nmh,nm0
      parameter (nmh=0, nm0=1-nmh)
c
c  Define maximum horizontal array dimensions for total grid (nmxa) and
c  for a single processor (nmx), and maximum vertical dimension (lmx).
      integer    nmxa,nmx,lmx
      parameter (nmxa=504, nmx=nmxa, lmx=385)
c
c  Declare passed variables.
      integer nt,mt,n,m,iec(8),j
      integer is(nm0:m+nmh),ie(nm0:m+nmh)
      real    times,ramp,epq
      real    amsk(nm0:n+nmh,nm0:m+nmh)
      real    elon(nm0:n+nmh,nm0:m+nmh)
      real    alat(nm0:n+nmh,nm0:m+nmh),ep(nm0:n)
      real ep1(nm0:n),ep2(nm0:n),ep3(nm0:n),ep4(nm0:n)
      real ep5(nm0:n),ep6(nm0:n),ep7(nm0:n),ep8(nm0:n)
c
c  declare local temporary variables.
      integer i,ic,k
      real    a,b
c
c  declare local saved variables.
      integer, save :: init(mxgrdso)
c
c  handle initialization
      data init/mxgrdso*0/
c      write(6,*)'tidepot:259,init(1),j=',init(1),j
      if (init(nesto) .ne. 0) go to 050
      init(nesto)=1
c      write(6,*)'About to call tidepot_init'
      call tidepot_init(nt,mt,n,m,is,ie,j,amsk,elon,alat,ep)
c      write(6,*)'Back from tidepot_init'
050   continue
c
c  apply at bndy pts??  Yes, I think so.
c
c  initialize ep.
c
!     do i=is(j)-1,ie(j)+iec(2)
      do i=1,n 
        do k=1,9   
         ep(i)=0.0
        enddo
      enddo
c      write(6,*)'Finished initializing htide(1:idm,j)'
      
c
c  sum over tidal constituents.
c  apply ramp to increase tidal forcing over a period of time.
c  for the declinational (fortnightly) tides (MF, MM), there is no
c  longitude dependence (i.e., coslon=1 and sinlon=0), so use the
c  simplified calculation for these.
      do ic=1,nc(nesto)
c          write(6,*)'tidepot:284,ic,amp(1,1),tidefq(1,1),times,',
c     &  'tideph(1,1),ramp='
c          write(6,*)ic,amp(1,1),tidefq(1,1),times,tideph(1,1),ramp
        a=ramp*amp(ic,nesto)
     &   *cos(tidefq(ic,nesto)*times+tideph(ic,nesto))
        b=ramp*amp(ic,nesto)
     &   *sin(tidefq(ic,nesto)*times+tideph(ic,nesto))
        k=itspec(ic,nesto)
c        write(6,*)"292:k,ic,nest0,tfq=",k,ic,nesto,tidefq(ic,nesto)
        if (tidefq(ic,nesto) .lt. 1.0e-5) then
c            write(6,*)'300::'
          do i=1,n  
c          do i=is(j)-1,ie(j)+iec(2)
            ep(i)=ep(i)+tplat(nesto)%p(i,j,k)*a
          enddo
        else
c          do i=is(j)-1,ie(j)+iec(2)
c      write(6,*)'307:n,a,b,ep(1),j,k=',n,a,b,ep(1),j,k
c      write(6,*)'307:n,a,b,ep(1),j,k=',n,a,b,ep(1),j,k
c      write(6,*)'j,k,tplat(1)%p(1,j,k)=',j,k,tplat(1)%p(1,j,k)
c      write(6,*)'j,k,tplat(1)%p(1,j,k)=',j,k,tplat(1)%p(1,j,k)
          do i=1,n
c             write(6,*)'310,i,ep(i)=',i,ep(i) 
c             write(6,*)'310,i,ep(i)=',i,ep(i) 
            epq=tplat(nesto)%p(i,j,k)
     &           *(coslon(nesto)%p(i,j,k)*a-sinlon(nesto)%p(i,j,k)*b)
            ep(i)=ep(i)+epq
            if(ic.eq.1)ep1(i)=epq
            if(ic.eq.2)ep2(i)=epq
            if(ic.eq.3)ep3(i)=epq
            if(ic.eq.4)ep4(i)=epq
            if(ic.eq.5)ep5(i)=epq
            if(ic.eq.6)ep6(i)=epq
            if(ic.eq.7)ep7(i)=epq
            if(ic.eq.8)ep8(i)=epq
          enddo
        endif
      enddo
c
      return
      end subroutine tidepot

      subroutine tidepot_final
c-----------------------------------------------------------------------
c  subroutine to free memory allocated for tidal potential forcing
c  *** called only once and end of run ***
c  created 2007-06-11, Tim Campbell, NRL.
c  last updated 2007-06-11, Tim Campbell, NRL.
c-----------------------------------------------------------------------
c
      use NCOM_Tidepot_Mod
c
      implicit none
!#include "../../include/PARAM.h"
!#include "../../include/COMMON.h"
c
c  Define named real constants
      real       zero,zp5,one,two
      parameter (zero=0.0, zp5=0.5, one=1.0, two=2.0)
      
      integer nesto,ntcmx
      PARAMETER(nesto=1,ntcmx=8)
      
      integer io_unit_offset, istdo_unit
      data istdo_unit/6/
c  Define halo width (nm0).
      integer    nmh,nm0
      parameter (nmh=0, nm0=1-nmh)
c
c  Define maximum horizontal array dimensions for total grid (nmxa) and
c  for a single processor (nmx), and maximum vertical dimension (lmx).
      integer    nmxa,nmx,lmx
      parameter (nmxa=504, nmx=nmxa, lmx=385)
c
c  declare local temporary variables.
      integer nn
c
c      do nn=1,nnesto
      do nn=1,nesto
        if (associated(tplat(nn)%p)) then
          deallocate(tplat(nn)%p)
          nullify(tplat(nn)%p)
        endif
        if (associated(coslon(nn)%p)) then
          deallocate(coslon(nn)%p)
          nullify(coslon(nn)%p)
        endif
        if (associated(sinlon(nn)%p)) then
          deallocate(sinlon(nn)%p)
          nullify(sinlon(nn)%p)
        endif
      enddo
c
      return
      end subroutine tidepot_final
C==========================================================
!     
! File:   tide_fac.for
! Author: DanMoore
!
! Created on September 16, 2011, 4:53 PM
!
      subroutine tide_fac(ihh,idd,imm,iyear,xlat,ntides,iprint,
     &  kon2,freqx2,fx2,vud2)
      implicit none
c!include "../../include/PARAM.h"
c
      integer ihh,idd,imm,iyear,ntides,iprint
      character*5 kon2(ntides)
      real    xlat,freqx2(ntides),fx2(ntides),vud2(ntides)
c
c  subroutine to calculate tidal data needed for predicting the tides
c  for a particular time period.
c  inputs arguments:
c       ihh    = hour of day (0 - 23, gmt).
c       idd    = day of month (1 - 31).
c       imm    = month (1 - 12).
c       iyear  = year (4-digit integer).
c       xlat   = mean latitude of region on degrees n.
c       ntides = number of tidal constituents to be obtained.
c       kon2   = array names of tidal constituents.  these must
c                be 5 characters long and in caps (see data file
c                model_tide/ios_tidetbl for list of valid names).
c       iprint = flag to print output (if iprint=1).
c  output arguments:
c       freqx2 = tidal frequency for each constituent in radians per sec.
c       fx2    = nodal factor for each tidal constituent.  this is a
c                small (usually less than 10%) correction to the tidal
c                amplitude.
c       vud2   = tidal phase for each constituent for input date.
c  note that the input date is gmt.
c
c  a particular tidal constituent can be calculated as
c
c       tide = fx2*amp*cos( freqx2*(t - t0) - phase + vud2)
c
c  where amp and phase are the equilibrium amplitude and phase for the 
c  tidal constituent at a particular location, t is the time, and t0
c  is the time for which the tidal data were calculated (i.e., the
c  input date to tide_fac).
c
c' created 5-7-97 from Cheryl Ann Blain's program tide_fac.
c
c'...............Cheryl Ann Blain's program tide_fac....................
c
c       program tide_fac.f
c
c       description:  this program computes tidal constituent   
c                     frequencies, nodal factors, and equilibrium
c                     arguments for a specific date and latitude.
c 
c       required input files:
c                 ios tidal constituent database    (unit =  8)
c                 user defined analysis values      (unit =  9)
c                    date: ihr,iday,imonth,iyear
c                    latitude (degrees)
c                    number of constituents to be analyzed
c                    name of each constituent (5 characters, uppercase)
c
c       output file:
c                 table of tidal constituent information (unit = 12)
c
c       written by:   Dr. Cheryl Ann Blain
c                     Naval Research Laboratory
c                     NRL Code 7322
c                     Stennis Space Center, MS 39529
c                     blain@nrlssc.navy.mil  (601)688-5450
c
c       date:  october 28, 1996
c
!include "../../include/COMMON.h"
c
c   local variables
c
      integer i, kd, kh, mtd, istdo_unit
      real    cyc, five, fk, freqk, freqx, fx, pi, twopi, vud, vuk, vux
      character*5 konk,kon(170),konx
c
      real   degree
c
c  Define named real constants
      real       zero,zp5,one,two
      parameter (zero=0.0, zp5=0.5, one=1.0, two=2.0)
      
      degree(cyc)=(cyc-ifix(cyc)-sign(zp5,cyc)+zp5)*360.

      mtd=9
      pi=3.1415926536
      twopi=2.*pi
c
c*   date for analysis: ihh,idd,imm,iyy,icc
c*
c*     ihh,idd,imm,iyear = hour, day, month, and year, of
c*                           the date for analysis
c*
c*   central latitude for computation
c*     a latitude is needed for the purpose of calculating
c*     satellite to main constituent amplitude ratios.
c*
c*     xlat= latitude in degrees
c*
c*
c*   number of constituents and names of those constituents
c*     146 constituents are possible and listed in the
c*     data file "ios_tidetbl"
c*
      do i=1,ntides
         kon(i)=kon2(i)
      enddo

      five=5.0
      if(abs(xlat).lt.five) xlat=sign(five,xlat)
      if(idd)1080,1070,1080
1070  continue
      call xcstop('Error in tide_fac:  xlat cannot be zero.')
      stop
c
c***********************************************************************
c* convert the date into a Gregorian day number
c*
1080  call gday(idd,imm,iyear,kd)
      kh=kd*24+ihh
      kd=kh/24
c
c***********************************************************************
c* astronomical arguments, satellites to the main constituents, and
c* shallow water constituents are read on device 199. (file: ios_tidetbl)
c
      call opnvuf(kh,konk,xlat,fk,vuk,freqk)
c
c***********************************************************************
c* v (= astronomical phase argument), u & f (nodal modulation
c* phase and amplitude corrections) are calculated for all mtot
c* constituents.
c
      call setvuf(kh,konk,xlat,fk,vuk,freqk)
c
c*********************************************************************
c* output the frequency, nodal factor and equilibruim argument
c*  for each constituent
c
cpjm  change output unit number from 12 to 6.
cpjm  print output only if iprint=1.
c
      if (iprint .eq. 1) then
c      if (mnproc.eq.1) then
      istdo_unit=6
      write(istdo_unit,1084)imm, idd, iyear, ihh
1084  format(1x,'DATE (M-D,YR HR): ',2x,i2,'-',i2,',',1x,i4,2x,i2)
      write(istdo_unit,*)
      write(istdo_unit,*)'LATITUDE (deg): ',xlat
      write(istdo_unit,*)
      write(istdo_unit,*)
     &'   TIDE         FREQUENCY (rad/s)      NOD. FAC.    EQU. ARG.'
      write(istdo_unit,*)
     &'   ____       ____________________     _________    _________'
      write(istdo_unit,*)
c      call zhflsh(istdo_unit)
c      endif
      endif
      do i=1,ntides
         konx=kon(i)
         call vuf(kd,konx,xlat,fx,vux,freqx)
c*
c* convert frequency from cycles/hr to radians/sec
         freqx=freqx*twopi/3600.
c*
c* convert equilibruim argument from cycles to degrees
         vud=degree(vux)
         if (iprint .eq. 1) then
c           if (mnproc.eq.1) then
           write(istdo_unit,1082)konx(1:2),freqx,fx,vud
1082       format(5x,a2,5x,e20.12,5x,f9.6,5x,f7.2)
!           write(istdo_unit,1082)konx,freqx,fx,vud
!1082       format(5x,a5,5x,e20.15,5x,f8.6,5x,f7.2)
c           call zhflsh(istdo_unit)
c           endif
         endif

cpjm  define output values for each tidal constituent.
         freqx2(i)=freqx
         fx2(i)=fx
         vud2(i)=vud

      enddo
cpjm  put a return here to convert to subroutine.
      return
c     end of tide_fac.
      end
c==========================================================
      subroutine opnvuf(kh,konk,xlat,fk,vuk,freqk)
      implicit none
!!#include "../../include/PARAM.h"
c
      character*5 konk
      integer kh
      real    xlat,fk,vuk,freqk,one,zero
c
c opnvuf - reads in kontab and data and calls setvuf
c
c  saved common variables for vuf, opnvuf and setvuf.
      character*5   konco,kontab
      common/vufc5/ konco(320),kontab(170)
      save  /vufc5/
      integer       ii,jj,kk,ll,mm,nn,ldel,mdel,ndel,ir,nj,
     &              ntidal,ntotal
      common/vufi4/ ii(50),jj(50),kk(50),ll(50),mm(50),nn(50),
     &              ldel(205),mdel(205),ndel(205),ir(205),nj(170),
     &              ntidal, ntotal
      save  /vufi4/
      real          freq,ee,ph,semi,coef,f,vu
      common/vufr4/ freq(170),ee(205),ph(205),semi(50),coef(320),
     &              f(170),vu(170)
      save  /vufr4/
c
!!!#include "../../include/CAF.h"
c
c     local variables
      integer j, j1, j4, jbase, jl, k, k1
      character*13 cenv

      character*5  kblank
      data kblank/'     '/
      one=1
      zero=0

cpjm  change unit number of input file of tidal data from 8 to 199.

c      write(cenv,'(''NCOM_'',a,''_'',i1,a)') 'TIDES',0,'D'
c      call zhopne(199,cenv,'formatted','old',0)
       OPEN(199,FILE="IOS_tidetbl.D",STATUS="old",ACTION="read") 
c
c***********************************************************************
c*  here the main constituents and their doodson numbers are read in
c*  format (6x,a5,1x,6i3,f5.2,i4). the values are respectively
c*     kon    = constituent name
c*  ii,jj,kk,ll,mm,nn = the six doodson numbers
c*     semi   = phase correction
c*     nj     = the number of satellites for this constituent.
c*  the end of all main constituents is denoted by a blank card.
c

      jbase=0
      do 90 k=1,50
c      read(unit=199,fmt=60 CAF_A) kontab(k),ii(k),jj(k),kk(k),ll(k),
      read(199,60) kontab(k),ii(k),jj(k),kk(k),ll(k),mm(k),nn(k),
     &semi(k),nj(k)
c      write(6,*)"k,kontab(k),ii(k),jj(k),kk(k),ll(k),mm(k),nn(k),"
c     &,"semi(k),nj(k)="
c      write(6,*)k,kontab(k),ii(k),jj(k),kk(k),ll(k),mm(k),nn(k),semi(k),
c     & nj(k)
   60 format(6x,a5,1x,6i3,f5.2,i4)
      if(kontab(k).eq.kblank) go to 100
70    j1=jbase+1
      if(nj(k).ge.1) go to 75
      nj(k)=1
      jl=j1
      ph(j1)=0.
      ee(j1)=0.
      ldel(j1)=0
      mdel(j1)=0
      ndel(j1)=0
      ir(j1)=0
      go to 90
   75 jl=jbase+nj(k)
c
c***********************************************************************
c*  if nj>0, information on the satellite constituents is read , three
c*  satellites per card, in the format (11x,3(3i3,f4.2,f7.4,1x,i1,1x)).
c*  for each satellite the values read are
c*     ldel,mdel,ndel = the changes in the last three doodson numbers
c*                      from those of the main constituent.
c*     ph     = the phase correction
c*     ee     = the amplitude ratio of the satellite tidal potential to
c*              that of the main constituent.
c*     ir     = 1 if the amplitude ratio has to be multiplied by the
c*                latitude correction factor for diurnal constituents
c*              2 if the amplitude ratio has to be multiplied by the
c*                latitude correction factor for semi-diurnal consti-
c*                tuents.
c*              otherwise if no correction is required to the amplitude
c*                ratio.
c

c      read(unit=199,fmt=80 CAF_A)
      read(unit=199,fmt=80)
     &   (ldel(j),mdel(j),ndel(j),ph(j),ee(j),ir(j),j=j1,jl)
   80 format((11x,3(3i3,f4.2,f7.4,1x,i1,1x)))

90    jbase=jl
100   ntidal=k-1        !for us this will be 45

c
c***********************************************************************
c*  the shallow water constituents and the main constituents from which
c*  they are derived are read in here with the format
c*  (6x,a5,i1,2x,4(f5.2,a5,5x)). the values are respectively
c*     kon    = name  of the shallow water constituent
c*     nj     = number of main constituents from which it is derived.
c*     coef,konco = combination number and name of these main
c*                  constituents.
c*  the end of these constituents is denoted by a blank card.
c
      jbase=0
      k1=ntidal+1
      do 160 k=k1,1000
      j1=jbase+1
      j4=j1+3

c      read(unit=199,fmt=130 CAF_A)
      read(unit=199,fmt=130)
     &   kontab(k),nj(k),(coef(j),konco(j),j=j1,j4)
  130 format(6x,a5,i1,2x,4(f5.2,a5,5x))
      if(kontab(k).eq.kblank) go to 170
160   jbase=jbase+nj(k)
170   ntotal=k-1        ! for us this will be 146
c      call zhclos(199)
CLOSE(199)
      call setvuf(kh,konk,xlat,fk,vuk,freqk)
      return
c     end of opnvuf.
      end

      subroutine setvuf(kh,konk,xlat,fk,vuk,freqk)   
      implicit none
!!#include "../../include/PARAM.h"
c
      character*5 konk
      integer kh
      real    xlat,fk,vuk,freqk,zero,one
c
c setvuf - evaluates f and vu for all constituents in kontab
c
c***********************************************************************
c*  ntidal is the number of main constituents
c*  ntotal is the number of constituents (main + shallow water)
c*  for  the given time kh, the table of f and v+u values is
c*  calculated for all the constituents.
c*     f is the nodal modulation adjustment factor for amplitude
c*     u is the nodal modulation adjustment factor for phase
c*     v is the astronomical argument adjustment for phase.
c
!!#include "../../include/COMMON.h"
c
c  saved common variables for vuf, opnvuf and setvuf.
      character*5   konco,kontab
      common/vufc5/ konco(320),kontab(170)
      save  /vufc5/
      integer       ii,jj,kk,ll,mm,nn,ldel,mdel,ndel,ir,nj,
     &              ntidal,ntotal,istdo_unit
      common/vufi4/ ii(50),jj(50),kk(50),ll(50),mm(50),nn(50),
     &              ldel(205),mdel(205),ndel(205),ir(205),nj(170),
     &              ntidal, ntotal
      save  /vufi4/
      real          freq,ee,ph,semi,coef,f,vu
      common/vufr4/ freq(170),ee(205),ph(205),semi(50),coef(320),
     &              f(170),vu(170)
      save  /vufr4/
c
c     local variables
      integer iflag, int24, intdys, inthh, iuu, iv, j, j1, jbase,
     &        jl, k, k1, kd, km1, l
      real    pi, rr, slat, sumc, sums, twopi, uu, uudbl, v, vdbl
      real*8  d1,h,pp,s,p,enp,dh,dpp,ds,dp,dnp,hh,tau,dtau
      
      istdo_unit=6
      zero=0
      one=1

      pi=3.1415926536
      twopi=2.*pi

      slat=sin(pi*xlat/180.)

c
c***********************************************************************
c*  the astronomical arguments are calculated by linear approximation
c*  at the mid point of the analysis period.
c*     ds,dh,dp,dnp,dpp are their respective rates of change over a 365
c*     day period as of 000 et 1/1/1976.
       
      d1=kh/24.d0
      call gday(31,12,1899,kd)
      d1=d1-kd-0.5d0
      call astr(d1,h,pp,s,p,enp,dh,dpp,ds,dp,dnp)
      int24=24
      intdys=kh/int24
      inthh=kh-intdys*int24
      hh=inthh
      tau=hh/24.d0+h-s
      dtau=365.d0+dh-ds


c
c***********************************************************************
c*  only the fractional part of a solar day need be retained for compu-
c*  ting the lunar time tau.
c

      jbase=0
      do 210 k=1,ntidal

        freq(k)=(ii(k)*dtau+jj(k)*ds+kk(k)*dh+ll(k)*dp+mm(k)*dnp+
     1          nn(k)*dpp)/(365.*24.)
!      write(6,*)'setvuf Line 217,k,freq(k)=',k,freq(k)
        vdbl=ii(k)*tau+jj(k)*s+kk(k)*h+ll(k)*p+
     &          mm(k)*enp+nn(k)*pp+semi(k)

        iv=vdbl
        iv=(iv/2)*2
        v=vdbl-iv
        j1=jbase+1
        jl=jbase+nj(k)
        sumc=1.
        sums=0.

        do 200 j=j1,jl
c
c***********************************************************************
c*  here the satellite amplitude ratio adjustment for latitude is made
c

          rr=ee(j)
          l=ir(j)+1
          go to (901,902,903),l
 902      rr=ee(j)*0.36309*(1.-5.*slat*slat)/slat
          go to 901
 903      rr=ee(j)*2.59808*slat
 901      continue
          uudbl=ldel(j)*p+mdel(j)*enp+ndel(j)*pp+ph(j)
          iuu=uudbl
          uu=uudbl-iuu
          sumc=sumc+rr*cos(uu*twopi)
          sums=sums+rr*sin(uu*twopi)

  200   continue

        f(k)=sqrt(sumc*sumc+sums*sums)
        vu(k)=v+atan2(sums,sumc)/twopi
      
210   jbase=jl

c
c***********************************************************************
c*  here f and v+u of the shallow water constituents are computed from
c*  the values of the main constituent from which they are derived.
c

      jbase=0
      k1=ntidal+1
      if(k1.gt.ntotal) return
      do 270 k=k1,ntotal
        f(k)=one
        vu(k)=zero
        freq(k)=zero
        iflag=1
      go to 268
 268  j1=jbase+1
      jl=jbase+nj(k)
      do 260 j=j1,jl
        km1=k-1
        do 240 l=1,km1
          if(kontab(l).eq.konco(j)) go to 250
 240    continue
c        if (mnproc.eq.1) then
        write(istdo_unit,241)konco(j)
 241    format('   SETVUF STOP.',a5)
c        call zhflsh(istdo_unit)
c        endif
        call xcstop('Error in setvuf:  kontab ne konco')
        stop
 250  f(k)=f(k)*f(l)**abs(coef(j))
      vu(k)=vu(k)+coef(j)*vu(l)
      if (iflag.eq.1) freq(k)=freq(k)+coef(j)*freq(l)
 260   continue
 270  jbase=jl
c      write(6,*)'opnvuf:Line 289 end of setvuf'
      return
c     end of setvuf.
      end
      
            subroutine vuf(kh,konk,xlat,fk,vuk,freqk)
      implicit none
!#include "../../include/PARAM.h"
c
c vuf - finds appropriate f,vu and sig for a specified constituent
c
      integer kh
      real    xlat,fk,vuk,freqk
      character*5 konk
c
!#include "../../include/COMMON.h"
c
c      subroutine vuf(kd0,konk,xlat,fk,vuk,freqk)
c      implicit none
c  obtained from Cheryl Ann Blain, 5-7-97.  (only slightly modified.)
c
c in opnvuf - reads in kontab and data and calls setvuf
c in setvuf - evaluates f and vu for all constituents in kontab
c in vuf - finds appropriate f,vu and sig for a specified constituent
c
c' taken from mike foreman's ios subroutine tide1.f
c  alterations by wendy gentleman,   july 1995
c  1) one or more records for each main constituent containing the
c     following information in format (6x,a5,1x,6i3,f5.2,i4)
c      kon    = constituent name
c   ii,jj,kk,ll,mm,nn = the six doodson numbers
c      semi   = phase correction
c      nj     = the number of satellites for this constituent.
c   if nj>0, information on the satellite constituents is read , three
c   satellites per card, in the format (11x,3(3i3,f4.2,f7.4,1x,i1,1x)).
c   for each satellite the values read are
c      ldel,mdel,ndel = the changes in the last three doodson numbers
c                       from those of the main constituent.
c      ph     = the phase correction
c      ee     = the amplitude ratio of the satellite tidal potential to
c               that of the main constituent.
c      ir     = 1 if the amplitude ratio has to be multiplied by the
c                 latitude correction factor for diurnal constituents
c               2 if the amplitude ratio has to be multiplied by the
c                 latitude correction factor for semi-diurnal consti-
c                 tuents.
c               otherwise if no correction is required to the amplitude
c                 ratio.
c
c  2) one blank record
c
c  3) one record for each shallow water constituent specifying in format
c     (11x,3(3i3,f4.2,f7.4,1x,i1,1x)), the following information
c      kon    = name  of the shallow water constituent
c      nj     = number of main constituents from which it is derived.
c      coef,konco = combination number and name of these main
c                   constituents.
c
c  4) one blank record
c
c=======================================================================
c
c*  mtot           is the number of constituents included in the data
c*                 file. at present this is 146.
c*  m              is the number of constituents included in the
c*                 analysis. at present it will have a maximum of 69 +
c                  the number of shallow water constituents specially
c*                 designated for inclusion (not in the standard 69
c*                 constituent package) in the analysis.
c*  kontab(i)      is the array containing all the constituent names as
c*                 they are read in from the data file. it should have
c*                 the minimum dimension mtot.
c*  kon(i)         is the array containing the names of all the
c*                 constituents that are to be included in the least
c*                 squares analysis to the tidal record. they are chosen
c*                 from the kontab list via the rayleigh criterion. the
c*                 array should have dimension at least m.
c*  freq(i)        is the array of frequencies (cycles/ hr.) correspond-
c*                 ing to the constituents contained in kontab. it
c*                 should have dimension at least mtot.
c*  sig(i)         is the array of frequencies corresponding to the
c*                 constituents contained in kon. it should have minimum
c*                 dimension m.
c*  v              astronomical phase argument
c*  u & f          nodal modulation &phase and amplitude corrections
c
c right now these (v, u & f) are calculated mtot constituents
c
c minimal dimension requirements
c ==============================
c the dimension of
c   kon and nj should be at least equal to the total number of possible
c   constituents plus one (presently 146+1), the dimension of ii, jj,
c   kk, ll, mm, nn, and semi should be at least equal to the number of
c   main constituents plus one (presently 45+1), the dimension of ee,
c   ldel, mdel, ndel, ir, and ph should be at least equal to the total
c   number of satellites to all the main constituents plus the number
c   of constituents with no satellites (presently 162+8 if there are no
c   third order satellites for both n2 and l2, and 167+8 if all third
c   order satellites are included for both n2 and l2),
c   and the dimension of konco, and coeff should be at least equal to
c   the sum for all shallow water constituents of the number of main
c   constituents from which each is derived,plus four (presently 251+4).
c=======================================================================
c
c  saved common variables for vuf, opnvuf and setvuf.
      character*5   konco,kontab,kkkkk
      common/vufc5/ konco(320),kontab(170)
      save  /vufc5/
      integer       ii,jj,kk,ll,mm,nn,ldel,mdel,ndel,ir,nj,
     &              ntidal,ntotal
      common/vufi4/ ii(50),jj(50),kk(50),ll(50),mm(50),nn(50),
     &              ldel(205),mdel(205),ndel(205),ir(205),nj(170),
     &              ntidal, ntotal
      save  /vufi4/
      real          freq,ee,ph,semi,coef,f,vu
      common/vufr4/ freq(170),ee(205),ph(205),semi(50),coef(320),
     &              f(170),vu(170)
      save  /vufr4/
c
c     local variables
      integer indx, l

c determines the index for the specified constituents
c      write(6,*)'opnvuf:line 413 konk = ',konk

        do  209 l=1,ntotal
            kkkkk=kontab(l)
          if(konk(1:2).eq.kkkkk(1:2)) then
            indx=l
            go to 211 
           end if
 209    continue
 
c        call xcstop('Error in VUF:  invalid tidal constituent name')
         WRITE(6,*)'Error in VUF:  invalid tidal constituent name'
        stop
 211    continue

        fk=f(indx)
        vuk=vu(indx)
        freqk=freq(indx)

cc        if (mnproc.eq.1) then
cc        write(6,*) konk(1:2),vuk*360.,fk,freqk
cc        call zhflsh(istdo_unit)
cc        endif
!98     format(' ',a5,5x,2f12.5)
cc98     format(' ',a2,5x,2f12.5)
      return
      end
C============================================================

      subroutine tc_amp(tidecn,amp)
c  subroutine to return amplitude for equilibrium tide for specified
c  tidal constituent.
c      tidecn = input tidal constituent (K1 O1 P1 Q1 K2 M2 N2 S2 MF MM).
c      amp    = returned amplitude in m.
c  note:  the returned equilibrium tidal amplitudes need to be corrected
c  for the "earth tide" (by multiplying by factor of ~ 0.69) and, if
c  simulating a particular time period, corrected for that particular
c  time period by muliplying by a "node factor" (this is generally a
c  small < 10% correction).  These corrections are NOT done here.
c  created 2003-08-18, Paul J Martin, NRL.
c
      implicit none
!#include "../../include/PARAM.h"
!#include "../../include/COMMON.h"
c
c  declare passed variables.
      character*(*) tidecn
      real    amp
c
c  declare local temporary variables.
      integer ic
c
c  declare local saved variables.
      integer nc
      parameter(nc=10)
      character*2 cons(nc)
      real      amps(nc)
      save      cons,amps
c
      data cons/'K1','O1','P1','Q1','K2','M2','N2','S2','MF','MM'/
      data amps/0.141565, 0.100514, 0.046843, 0.019256, 0.030704,
     &          0.242334, 0.046398, 0.112841, 0.041742, 0.022026/
c
      do ic=1,nc
        if (tidecn(1:2) .eq. cons(ic)(1:2)) then
          amp=amps(ic)
          return
        endif
      enddo
c      if (mnproc.eq.1) then
        write(6,'(a,i3,2a)')
     &    'tc_amp: ic =',ic,'  constituent = ',tidecn
c        call zhflsh(istdo_unit)
c      endif
      call xcstop('Error in tc_amp:  invalid tidal constituent')
      end
c=====================================================================
      subroutine astr(d1,h,pp,s,p,np,dh,dpp,ds,dp,dnp)
      implicit none
!!!#include "../../include/PARAM.h"
c
      real*8 d1,h,pp,s,p,np,dh,dpp,ds,dp,dnp
c
c  obtained from Cheryl Ann Blain, 5-7-97.  (not modified.)
c     this subroutine calculates the following five ephermides
c     of the sun and moon
c     h = mean longitude of the sum
c     pp = mean longitude of the solar perigee
c     s = mean longitude of the moon
c     p = mean longitude of the lunar perigee
c     np = negative of the longitude of the mean ascending node
c     and their rates of change.
c     units for the ephermides are cycles and for their derivatives
c     are cycles/365 days
c     the formulae for calculating this ephermides were taken from
c     pages 98 and 107 of the explanatory supplement to the
c     astronomical ephermeris and the american ephermis and
c     nautical almanac (1961)
c
      real*8 d2, f, f2
c
      d2=d1*1.d-4
      f=360.d0
      f2=f/365.d0
      h=279.696678d0+.9856473354d0*d1+.00002267d0*d2*d2
      pp=281.220833d0+.0000470684d0*d1+.0000339d0*d2*d2+
     1  .00000007d0*d2**3
      s=270.434164d0+13.1763965268d0*d1-.000085d0*d2*d2+
     1  .000000039d0*d2**3
      p=334.329556d0+.1114040803d0*d1-.0007739d0*d2*d2-
     1  .00000026d0*d2**3
      np=-259.183275d0+.0529539222d0*d1-.0001557d0*d2*d2-
     1  .00000005d0*d2**3
      h=h/f
      pp=pp/f
      s=s/f
      p=p/f
      np=np/f
      h=h-int(h)
      pp=pp-int(pp)
      s=s-int(s)
      p=p-int(p)
      np=np-int(np)
      dh=.9856473354d0+2.d-8*.00002267d0*d1
      dpp=.0000470684d0+2.d-8*.0000339d0*d1
     1  +3.d-12*.00000007d0*d1**2
      ds=13.1763965268d0-2.d-8*.000085d0*d1+
     1  3.d-12*.000000039d0*d1**2
      dp=.1114040803d0-2.d-8*.0007739d0*d1-
     1  3.d-12*.00000026d0*d1**2
      dnp=+.0529539222d0-2.d-8*.0001557d0*d1-
     1  3.d-12*.00000005d0*d1**2
      dh=dh/f2
      dpp=dpp/f2
      ds=ds/f2
      dp=dp/f2
      dnp=dnp/f2
      return
c     end of astr.
      end
C==========================================================
      subroutine dayceni(dcen,dcend,idate,timed)
c  subroutine to convert a fractional day-of-the-century, which is
c  defined as the number of days since 00Z, Jan 1, 1900, to a date and
c  time.  The year must lie between 1901 and 2099.
c  This subroutine is the inverse of subroutine daycen.
c       dcen  = number days since 00Z Jan 1, 1900 in single precision.
c       dcend = number days since 00Z Jan 1, 1900 in double precision.
c               note that either dcen or dcend must be set to zero,
c               only the non-zero value (dcen or dcend) is used.
c       idate = returned integer date of form YYYYMMDD.
c       timed = returned fractional time of day (0.0 to 1.0).
c  created 2-23-99, Paul J Martin, NRL.
c
      implicit none
!#include "../../include/PARAM.h"
!#include "../../include/COMMON.h"
      integer idate
      real    timed,dcen
      real*8  dcend,a,b
c
      integer iyear,month,iday,leap,idays,i,iyrest,idayr
c
      integer mon(12,2)
      save mon
      data mon/0,31,59,90,120,151,181,212,243,273,304,334,
     &         0,31,60,91,121,152,182,213,244,274,305,335/
c
c  note:  all years divisible by 4 are leap years except for century
c         years not divisible by 400, e.g., 1900 was NOT a leap year.
c
c  note:  use dcend (double precision) or if dcend is zero, use dcen
c         (single precision).
c
c  calc fractional time of day.
      if (dcen .lt. 365.0) then
        idays=int(dcend)
        a=float(idays)
        b=dcend-a
        timed=b
      else
        idays=int(dcen)
        timed=dcen-float(idays)
      endif
c
c  find year.
      iyrest=1900+(idays+0.125)/365.25
cc      write(istdo_unit,'(/a,i8,a,i8)')
cc     &  'dayceni:  idays=',idays,'  iyrest =',iyrest
      do i=iyrest+1,iyrest-1,-1
        idayr=idays-(i-1900)*365-(i-1901)/4
cc        write(istdo_unit,'(/a,i8,a,i12)')
cc     &    'dayceni:  i=',i,'  idayr=',idayr
        if (idayr .ge. 0) then
          iyear=i
          go to 210
        endif
      enddo
      stop 'Error in dayceni:  year not found'
210   continue
      if (iyear .lt. 1901) stop 'Error in dayceni:  year < 1901'
      if (iyear .gt. 2099) stop 'Error in dayceni:  year > 2099'
c
c  find month and day of month.
c  note that idayr here is actually the day of year - 1.
      if (mod(iyear,4) .eq. 0 .and. .not.
     &  (mod(iyear,100) .eq. 0 .and. .not. mod(iyear,400) .eq. 0)) then
        leap=2
      else
        leap=1
      endif
      do i=2,12
        if (mon(i,leap) .gt. idayr) then
          month=i-1
          go to 220
        endif
      enddo
      month=12
220   continue
      iday=idayr+1-mon(month,leap)
c
c  convert year, month, and day of month to idate (YYYYMMDD).
      call ymd2idt(iyear,month,iday,idate)
      return
      end

            subroutine ymd2idt(iyear,month,iday,idate)
c  subroutine to convert the year, month, and day to an integer
c  of the form YYYYMMDD.
c       iyear = output year (4-digit year).
c       month = output month (1 to 12).
c       iday  = output day of month (1 to 31).
c       idate = integer of form YYYYMMDD.
c  created 2-23-99, Paul J Martin, NRL.
c
      implicit none
      integer iyear,month,iday,idate
c
      idate=iyear*10000+month*100+iday
      return
      end
      
            subroutine idt2ymd(idate,iyear,month,iday)
c  subroutine to convert an integer of the form YYYYMMDD to year,
c  month, and day.
c       idate = integer of form YYYYMMDD.
c       iyear = output year (4-digit year).
c       month = output month (1 to 12).
c       iday  = output day of month (1 to 31).
c  created 2-23-99, Paul J Martin, NRL.
c
      implicit none
      integer idate,iyear,month,iday
      integer i
c
      i=idate
      iyear=i/10000
      i=i-iyear*10000
      month=i/100
      iday=i-month*100
      return
      end
      
            subroutine gday(idd,imm,iyear,kd)
      implicit none
!!!#include "../../include/PARAM.h"
c
      integer idd,imm,iyear,kd
c
c  obtained from Cheryl Ann Blain (5-7-97).  (not modified).
c
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
!!!#include "../../include/COMMON.h"
c
      integer icc, iyy,istdo_unit
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/
      
      istdo_unit=6
c!
c make iyy and icc variables
        icc=iyear/100
        iyy=iyear-icc*100

c!  test for invalid input:
      if(icc.lt.0)then
c        if (mnproc.eq.1) then
         write(istdo_unit,5000)icc
c        call zhflsh(istdo_unit)
c        endif
        call xcstop('Error in gday:  icc < 0')
        stop
      endif
      if(iyy.lt.0.or.iyy.gt.99)then
c        if (mnproc.eq.1) then
        write(istdo_unit,5010)iyy
c        call zhflsh(istdo_unit)
c        endif
        call xcstop('Error in gday:  iyy < 0 or > 99')
        stop
      endif
      if(imm.le.0.or.imm.gt.12)then
c        if (mnproc.eq.1) then
        write(istdo_unit,5020)imm
c        call zhflsh(istdo_unit)
c        endif
        call xcstop('Error in gday:  imm <= 0 or > 12')
        stop
      endif
      if(idd.le.0)then
c        if (mnproc.eq.1) then
        write(istdo_unit,5030)idd
c        call zhflsh(istdo_unit)
c        endif
        call xcstop('Error in gday:  ldd <= 0')
        stop
      endif
      if(imm.ne.2.and.idd.gt.ndm(imm))then
c        if (mnproc.eq.1) then
        write(istdo_unit,5030)idd
c        call zhflsh(istdo_unit)
c        endif
        call xcstop('Error in gday:  imm ne 2 and idd > ndm')
        stop
      endif
      if(imm.eq.2.and.idd.gt.29)then
c        if (mnproc.eq.1) then
        write(istdo_unit,5030)idd
c        call zhflsh(istdo_unit)
c        endif
        call xcstop('Error in gday:  imm = 2 and idd > 29')
        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
c        if (mnproc.eq.1) then
        write(istdo_unit,5030)idd
c        call zhflsh(istdo_unit)
c        endif
        call xcstop('Error in gday:  imm = 2 and idd > 28')
        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     end of gday.
      end
C======================================================================

       subroutine xcstop(string)
       character *(*) string
       write(6,*)trim(string)
       return
       end
      




