C***********************************************************************
C***********************************************************************
C Beginning of User-Specified Subroutines for TRUXTON6
C Beginning of User-Specified Subroutines for TRUXTON6
C***********************************************************************
C***********************************************************************
C USER IS RESPONSIBLE FOR 3 SUBROUTINES: VERTGRIDF6(FUNDY6 STANDARD)
C                                        OUTPUTT6  (TRUXTON6 STANDARD)
C                                        TRUXOUT6  (TRUXTON6 STANDARD)
C***********************************************************************
C***********************************************************************
      SUBROUTINE VERTGRIDF6(NN,NE,X,Y,HDOWN,IN,NLEV,ZL,NNV,
     &Z,AK,ENZ,RX,RY,HRBARELX,HRBARELY)
C Include file for parameter statements
      INCLUDE 'FUNDY6.DIM'
C Fixed global arrays
      INTEGER IN(3,NEDIM)
      REAL X(NNDIM),Y(NNDIM),HDOWN(NNDIM),AK(NNDIM),ZL(NLEVDIM)
      REAL Z(NNDIM,NNVDIM),ENZ(NNDIM,NNVDIM)
      COMPLEX RX(NNDIM,NNVDIM),RY(NNDIM,NNVDIM)
      COMPLEX HRBARELX(NEDIM),HRBARELY(NEDIM)
C-----------------------------------------------------------------------
C This subroutine must be modified by the user to assign proper values
C   to the arrays:
C          Z(I,J) - Vertical coordinate of 3-D nodes
C          AK(I) - Linear bottom stress coefficient, k
C          ENZ(I,J) - Vertical viscosity, Nz
C          RX(I,J) - RHS of x-momentum equation
C          RY(I,J) - RHS of y-momentum equation
C          HRBARELX(L) - RHS of elevation equation
C          HRBARELX(L) - RHS of elevation equation
C
C Arrays passed to this subroutine to facilitate computations:
C          X(I),Y(I) - Horizontal nodal coordinates
C          HDOWN(I) - Bathymetric depth at nodes
C          IN(3,L) - LABLIB STANDARD triangular element incidence list
C
C Index range:  I=1,NN, J=1,NNV
C          I=1,NN, J=1,NNV, L=1,NE
C          NN - number of horizontal nodes
C          NE - number of horizontal 2-D triangular elements
C          NNV - number of vertical nodes
C
C Z and J are positive upward:
C          Bottom:  Z=-HDOWN; J=1
C          Surface: Z= 0    ; J=NNV
C-----------------------------------------------------------------------
C Four standard subroutines are available to assign values to these
C   arrays:  Subroutines UNIGRID5, SINEGRID5, RHOGRID5, and BAROCLINIC5.
C
CCCCC
C
C   Subroutine UNIGRID5 assigns uniform ENZ, uniform AK, uniform
C     vertical spacing, and zero baroclinic forcing.  Assignment of NNV,
C     EN0 and AK0 is required before this subroutine is called.
C
C      CALL UNIGRID5(NNDIM,NEDIM,NNVDIM,NN,NE,NNV,EN0,AK0,
C     *      Z,AK,ENZ,RX,RY,HRBARELX,HRBARELY,HDOWN)
C
CCCCC
C
C   Subroutine SINEGRID5 assigns uniform ENZ, uniform AK, sinusiodal
C     vertical spacing, and zero baroclinic forcing.  It requires the
C     assignment of NNV, AKO, and ENO before it is called.
C
C      CALL SINEGRID5(NNDIM,NEDIM,NNVDIM,NN,NE,NNV,EN0,AK0,DZBL,
C     *      Z,AK,ENZ,RX,RY,HRBARELX,HRBARELY,HDOWN)
C
CCCCC
C
C   Subroutine RHOGRID5 assigns uniform ENZ, uniform AK, and reads NNV,
C     Z, baroclinic forcing from a LABLIB Standard .rho file identified
C     by the user in RHOGRID5.  Assignment of EN0 and AK0 is required 
C     before this subroutine is called.
C
C      CALL RHOGRID5(NNDIM,NEDIM,NNVDIM,NN,NE,NNV,EN0,AK0,
C     *      Z,AK,ENZ,RX,RY,HRBARELX,HRBARELY,HDOWN)
C
CCCCC
C
C   Subroutine BAROCLINIC5 computes the baroclinic forcing terms for a
C     given vertical structure (Z) and a LABLIB Standard .lst file 
C     identified by the user in BAROCLINIC5.  Assignment of Z is 
C     required before this subroutine is called
C
C      CALL BAROCLINIC5(NN,NE,NNV,X,Y,IN,Z,RX,RY,HRBARELX,HRBARELY)
C
C================Beginning-of-user-specified-instructions===============
C
C User-specified local assignments
C
      CHARACTER*72 FILEAK,FILEEN
c
c Set parameters for vertical grid
      NNV=21
      AK0=0.004
      EN0=0.04
      DZBL=2.5     
      CALL SINEGRID5(NNDIM,NEDIM,NNVDIM,NN,NE,NNV,EN0,AK0,DZBL,
     *      Z,AK,ENZ,RX,RY,HRBARELX,HRBARELY,HDOWN)
c
c Overwrite friction parameters with values from specified files
      WRITE(*,'(2A)')'Enter the name of the file containing the time-',
     &'average values of the linearized bottom stress coefficient:'
      READ(*,'(A)')FILEAK
      WRITE(*,'(2A)')'Enter the name of the file containing the time-',
     &'average values of the vertical eddy viscosity:'
      READ(*,'(A)')FILEEN
      OPEN(71,FILE=FILEAK,STATUS='OLD')
      READ(71,'(A)')HEADER
      READ(71,'(A)')HEADER
      DO I=1,NN
         READ(71,*)II,AK(I)
      ENDDO
      IF(II.NE.NN)THEN
         WRITE(*,*)'ERROR READING ',FILEAK
         STOP
      ENDIF
      OPEN(72,FILE=FILEEN,STATUS='OLD')
      READ(72,'(A)')HEADER
      READ(72,'(A)')HEADER
      READ(72,*)NNV
      DO I=1,NN
         DO J=1,NNV
            READ(72,*)II,ZIJ,ENZ(I,J)
         ENDDO
      ENDDO
      IF(II.NE.NN)THEN
         WRITE(*,*)'ERROR READING ',FILEEN
         STOP
      ENDIF
C
C End of routine
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE OUTPUTT6(MESHNAME,NN,NE,X,Y,IN,NNV,Z,
     &ZETA,UBAR,VBAR,U,V,W,FREQ,IYEAR,YDGMT0,LRVS,NR)
C Include file for parameter statements
      INCLUDE 'FUNDY6.DIM'
C Fixed global arrays
      REAL X(NNDIM),Y(NNDIM),Z(NNDIM,NNVDIM)
      INTEGER IN(3,NEDIM),LRVS(NNDIM)
      COMPLEX ZETA(NNDIM),UBAR(NNDIM),VBAR(NNDIM),
     &U(NNDIM,NNVDIM),V(NNDIM,NNVDIM),W(NNDIM,NNVDIM)
      CHARACTER*72 MESHNAME
C-----------------------------------------------------------------------
C User-defined output routine
C
C Arrays and constants:
C          X(I),Y(I) - Horizontal nodal coordinates
C          IN(3,L) - LABLIB STANDARD triangular element incidence list
C          Z(I,J) - Vertical coordinate of 3-D nodes
C          ZETA(I) - Free surface elevation at nodes
C          UBAR(I) - Vertically averaged velocity in x-direction
C          VBAR(I) - Vertically averaged velocity in y-direction
C          U(I,J),V(I,J),W(I,J) - x,y,and z components of velocity
C          FREQ - Frequency in rad/sec
C
C Index range:  I=1,NN, J=1,NNV
C          NN - number of horizontal nodes
C          NE - number of horizontal 2-D triangular elements
C          NNV - number of vertical nodes
C
C Z and J are positive upward:
C          Bottom:  Z=-HDOWN; J=1
C          Surface: Z= 0    ; J=NNV
C-----------------------------------------------------------------------
C A standard routine 'SUBROUTINE DUMP5' writes all computed results to a
C   user-specified file.  A single line in this subroutine suffices:
C
C      CALL DUMP5(X,Y,Z,ZETA,UBAR,VBAR,U,V,W,
C     &NN,NNV,NNDIM,NNVDIM)
C
C DUMP takes time and produces a large file.
C================Beginning-of-user-specified-instructions===============
C
C User-specified local assignments
C
      character*72 filename,begin,meshname
      begin='truxton6_simulation/'
      indxb=index(begin,' ')-1
      indxm=index(meshname,' ')-1
c
c Write tidal elevation bcs.s2c file
c
      write(*,*)'    Enter name of the elevation .bcs.s2c results file'
      write(*,*)'      (entering NO => no file will be created)'
      read(*,100)filename         
      write(*,100)filename
      if(filename.NE.'NO')then        
         OPEN (UNIT=7,FILE=filename,STATUS='NEW')
         write(7,100)meshname(1:indxm)
         write(7,101)IYEAR,YDGMT0,begin(1:indxb)//'tidal elevation BCS'
         write(7,8000)FREQ
         write(7,8002)(lrvs(i),cabs(zeta(lrvs(i))),
     &                 phaselagd(zeta(lrvs(i))),i=1,nr)
         close(7)
      endif
c
c Write tidal elevation .s2c file
c
      write(*,*)'    Enter name of the elevation .s2c results file'
      write(*,*)'      (entering NO => no file will be created)'
      read(*,100)filename         
      write(*,100)filename
      if(filename.NE.'NO')then        
         OPEN (UNIT=7,FILE=filename,STATUS='NEW')
         write(7,100)meshname(1:indxm)
         write(7,101)IYEAR,YDGMT0,begin(1:indxb)//'tidal elevations'
         write(7,8000)FREQ
         write(7,8002)(i,cabs(zeta(i)),phaselagd(zeta(i)),I=1,nn)
         close(7)
      endif
c
c Write depth average tidal velocity .v2c file
c
      write(*,*)'    Enter name of the average velocity .v2c results fil
     &e'
      write(*,*)'      (entering NO => no file will be created)'
      read(*,100)filename        
      write(*,100)filename
      if(filename.NE.'NO')then        
         open (unit=7,file=filename,status='new')
         write(7,100)meshname(1:indxm)
         write(7,101)IYEAR,YDGMT0,begin(1:indxb)//'tidal velocities'
         write(7,8000)FREQ
         write(7,8004)(I,cabs(UBAR(I)),phaselagd(UBAR(I)),
     &cabs(VBAR(I)),phaselagd(VBAR(I)),I=1,nn)
         close(7)
      endif
c
c Format Statements
c
 100  format(a)
 101  format(i4,1x,f8.4,1x,a)
 8000 format(1e14.6)
 8002 format(i5,2e14.6)
 8004 format(i5,4e14.6)
C
C End of routine
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      Subroutine TruxOut6(NOBS,Xobs,Yobs,Tobs,Tobsgmt0,IYEAR,
     & Ubarobs,Vbarobs,Zmaxobs,Zminobs,Zbotob,
     & Nf, Nf1, Nf2, freq, infer, amprat, 
     & NBC,ZCamp,LRVS,CCov, 
     & Ubarls, Vbarls, Upartial, Vpartial, 
     & Uerr, Verr, RMSV,rcond,irank,maxrank)
C     
      Include 'TRUXTON6ARRAYS.DIM' 
C
      Real Xobs(*),Yobs(*),Tobs(*),Tobsgmt0(*),
     & Ubarobs(*),Vbarobs(*),Zmaxobs(*),Zminobs(*),Zbotob(*),
     & freq(*), infer(*), amprat(*), 
     & LRVS(*), UBARls(*), VBARls(*), 
     & Upartial(nobsdim, nfdim), Vpartial(nobsdim, nfdim), 
     & Uerr(*), Verr(*)
C     
      Complex ZCamp(nbcdim,nfdim)
      Complex CCov(NBCdim*NFdim,NBCdim*NFdim)
C
C USER-BUILT Subroutine for Truxton1 Output
C
C Arrays passed:  
C
C--------------------------------------------------------------
C Observations:
C   NOBS = # of Velocity Observations
C   (Xobs,Yobs,Tobs,Tobsgmt0)_i = Space-time coords of observ'n i
C   (Ubarobs,Vbarobs)_i         = Observed velocity at i
C   (Zmaxobs,Zminobs)           = Averaging depth interval for (U,V)_i
C    Zbotobs_i                  = Bathymetric depth observed at i
C--------------------------------------------------------------
C Spectrum of tidal constituents
C    Nf   = number of constituents = Nf1+Nf2
C    Nf1  = number of primary constituents
c    Nf2  = number of additional constituents to be inferred
C    (freq, infer, amprat)_i, i=i, Nf: 
c       frequency is the radian frequency, sec**(-1)
c       infer_i is the inference frequency number
c            infer_i  =  i    for i <= Nf1
c            infer_i <= Nf1   for i >  Nf1
c       amprat_i is the complex amplitude ratio for an inferred const.
C            ZCamp_i  = amprat_i * ZCamp_(infer_i) for  i >  Nf1
C            ZCamp_i  = (1.,0.)                    for  i <= Nf1
C------------------------------------------------------------
C Optimal BC's:
C   NBC      = Number of Boudary Nodes
C   ZCamp_ij = COMPLEX BC amplitude:  Node i,    frequency j
c                                   I=1, NBC  //   j=1,Nf
C   LRVS_i   = Global node number pointer for the ith BC
C   CCov_ij  = COMPLEX covariance of the BC estimates 
C     The indices are numbered as follows:
C       i,j  = 1, NBC            :BC's for all nodes at frequency 1
C       i,j  = NBC+1, 2NBC       :BC's for all nodes at frequency 2
C       i,j  = (k-1)NBC+1, k*NBC :BC's for all nodes at frequency k     

C     CCov is normalized by the mean square velocity error, RMSV**2
c     se(i)  = sqrt(cov(i,i))is the standard error for the ith BC
c     SE(i)  = RMSV*se(i)    is the Standard Error for the ith BC
c   Zero Covariance is reported for Inference frequencies.  
C--------------------------------------------------------------
C Modeled Velocity:
C   (Ubarls, Vbarls)_i      = Model result at observation point i
C   (Upartial, Vpartial)_ij = Model result at i due to frequency j
c          i=i, nobs; j=i, nf1
C------------------------------------------------------------
C Errors (Model - Observation)
C   RMSV = Root Mean Square Velocity Error 
C   (Uerr, Verr)_i = Velocity error at observation point i 
C     i=1,Nobs 
C------------------------------------------------------------
C Warning: Unit (2) is the echo file for the run and is open 
C   during execution of TruxOut.  
C------------------------------------------------------------
C================Beginning-of-user-specified-instructions===============
C
C User-specified local assignments
C
      character*72 filename
c
c Output RMSV to standard out
c
      write(*,*)' '
      write(*,*)'   Root Mean Square Velocity Error =',RMSV
      write(*,*)'   Condition Number                =',rcond
      write(*,*)'   Actual Rank                     =',irank
      write(*,*)'   Maximum Possible Rank           =',maxrank
      write(*,*)' '
c
c Write output .dat file
c
      write(*,*)'    Enter name for the .dat file'
      write(*,*)'      (entering NO => no file will be created)'
      read(*,'(A)')filename         
      write(*,'(A)')filename
      if(filename.NE.'NO')then        
         OPEN (UNIT=7,FILE=filename,STATUS='NEW')
         DO I=1,NOBS
            write(7,'(8e14.6)')TOBSGMT0(I),XOBS(I),YOBS(I),
     &            UBARLS(I),VBARLS(I),
     &            UERR(I),VERR(I),TOBS(I)
         enddo
         close(7)
      endif
C
C End of routine
C
      Return
      End
C***********************************************************************
C***********************************************************************
C End of User-Specified Subroutines for TRUXTON6
C End of User-Specified Subroutines for TRUXTON6
C***********************************************************************
C**********************************************************************
