C***********************************************************************
C***********************************************************************
C Beginning of User-Specified Subroutines for FUNDY5
C Beginning of User-Specified Subroutines for FUNDY5
C***********************************************************************
C***********************************************************************
      SUBROUTINE BC(WW,U,V,NNDIM,III)
C Fixed global variables
      COMPLEX U(NNDIM),V(NNDIM)
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================Beginning-of-user-specified-instructions===============
C
C User-specified local assignments
C
      COMPLEX EYE,SPEED,VNORM
      character*72 bcsfile,header,meshname
      bcsfile='g2sm2.s2c'
C
C Beginning of executable code
C
      IF(III.EQ.0) GO TO 3
      IF(III-2)4,5,6
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C This section assigns the frequency of the forcing in SEC**-1 to WW.
C A return statement must follow the WW=___  statement.
C
    3 CONTINUE
      print*,'  Subroutine BC'
      open(unit=11,file=bcsfile,status='old')
      read(11,101)meshname
      read(11,101)header
      read(11,*)WW
      rewind(11)
      close(11)
      write(*,*)'    frequency = ',WW
      RETURN
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Elevation BC's are assigned in this section.  One statement of the
C   form U(I)=CVALUE must appear for each set depth node, where I is the
C   horizontal node number and CVALUE is the complex elevation for node 
C   I.  
C A return statement must follow the last elevation BC assignment.
C
    4 CONTINUE
      EYE=(0.,1.)
      FAC=-ACOS(-1.)/180.
C
C Read the data file of elevation BC's
C
      print*,'  Subroutine BC'
      write(*,*)'    Specifying BCS from ', bcsfile
 101  format(a)
      open(unit=11,file=bcsfile,status='old')
      read(11,101)header
      read(11,101)header
      read(11,*)WW
      do 10 i=1,nndim
         read(11,*,end=11)nodenum,elamp,elpha
         U(nodenum) = elamp*CEXP(EYE*FAC*elpha)
 10   continue
 11   continue
      rewind(1)
      close(1)
      RETURN
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Nonzero normal velocity BC's are assigned in this section.  One 
C   statement of the form U(I)=CVALUE must appear for each set nonzero
C   normal velocity node, where I is the horizontal node number and 
C   CVALUE is the complex value of 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
      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      U(I)=CXVEL
C      V(I)=CYVEL
C   I is the horizontal node number, (CXVEL,CYVEL) are the complex
C   values of the (X,Y) velocities for node I.  
C A return statement must follow the last velocity BC assignment.
C
    6 CONTINUE
      RETURN
C
C End of routine
C
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE ATMOS(ATMX,ATMY,NNDIM)
C Fixed global variables
      COMPLEX ATMX(NNDIM),ATMY(NNDIM)
C-----------------------------------------------------------------------
C This subroutine must be modified by the user to assign proper values
C   to arrays ATMX and ATMY, the right-sides of the conservative X- and
C   Y-momentum equations, respectively.  (i.e. H*PSI)  Two statements of
C   the following form must appear for each node where the atmospheric 
C   forcing is nonzero:
C      ATMX(I)=XVALUE
C      ATMY(I)=YVALUE
C   I is the node number and XVALUE,YVALUE are the complex values of the
C   forcing.  ATMX,ATMY are initialized with the value of zero by the
C   main routine before this routine is called.  Therefore, if no value
C   is assigned for a node in this routine, the default value of zero 
C   is used.
C A return statement must follow the last assignment.
C================Beginning-of-user-specified-instructions===============
C
C User-specified local assignments
C
      COMPLEX STRESS
C
C Beginning of Executable Code
C
C
C Specify the wind stress divided by the reference density
C
      smag=0.000
      sangt=137.0
C
C   Set complex stress
C
      STRESS=CMPLX(smag,0.0) 
C
C Write stress being applied to the CRT
C
      write(*,*)'  Subroutine ATMOS'
      write(*,*)'    Wind stress = ',smag
      write(*,*)'    Wind angle =  ',sangt
C
C Compute x and y components of the stress
C
      sangx=90.-sangt
      sx=cos(sangx*3.14159265/180.)
      sy=sin(sangx*3.14159265/180.)
C
C Load ATMX and ATMY arrays
C
      DO 10 I=1,NNDIM
         ATMX(I)=STRESS*sx
         ATMY(I)=STRESS*sy
 10   continue
C
C End of routine
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE VERTGRID5(NNDIM,NEDIM,NNVDIM,NN,NE,X,Y,IN,
     &NNV,Z,AK,ENZ,RX,RY,HRBARELX,HRBARELY,HDOWN)
C Fixed global arrays
      INTEGER IN(3,NEDIM)
      REAL X(NNDIM),Y(NNDIM),HDOWN(NNDIM),AK(NNDIM)
      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 NNV (the actual number of nodes in the vertical) and to the
C   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 The first line of vertgrid following the dimensioning should assign
C   NNV (on return to the main program, a check of NNV relative to 
C   NNVDIM will be made).
C
C Then, there must be assignments to fill the aforementioned 
C   real (Z,AK,ENZ) and complex (RX,RY,HRBARELX,HRBARELY) arrays over 
C   the index range:
C          I=1,NN, J=1,NNV, L=1,NE
C          NN - number of horizontal nodes
C          NE - the number of horizontal 2-D triangular elements
C          NNV - actual 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 return statement terminates this subroutine.
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
      NNV=21
      AK0=0.005
      EN0=0.02
      DZBL=2.5
      CALL UNIGRID5(NNDIM,NEDIM,NNVDIM,NN,NE,NNV,EN0,AK0,
     *      Z,AK,ENZ,RX,RY,HRBARELX,HRBARELY,HDOWN)
C      CALL SINEGRID5(NNDIM,NEDIM,NNVDIM,NN,NE,NNV,EN0,AK0,DZBL,
C     *      Z,AK,ENZ,RX,RY,HRBARELX,HRBARELY,HDOWN)
      RETURN
C
C End of routine
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE OUTPUT5(X,Y,Z,ZETA,UBAR,VBAR,U,V,W,OMEGA,
     * NN,NNV,NNDIM,NNVDIM)
C Fixed global arrays
      DIMENSION X(NNDIM),Y(NNDIM),Z(NNDIM,NNVDIM)
      COMPLEX ZETA(NNDIM),UBAR(NNDIM),VBAR(NNDIM),
     *     U(NNDIM,NNVDIM),V(NNDIM,NNVDIM),W(NNDIM,NNVDIM)
C-----------------------------------------------------------------------
C User-defined output routine
C
C Arrays and constants:
C          X(I),Y(I) - cartesian nodal coordinates
C          Z(I,J) - Vertical coordinate of 3-D nodes
C          ZETA(I) - the 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
C Index range:  I=1,NN, J=1,NNV
C          NN - number of horizontal nodes
C          NNV - actual 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-----------------------------------------------------------------------
C A standard routine 'SUBROUTINE DUMP' writes all computed results to a 
C   user-specified file.  A single line in this subroutine suffices:
C
C      CALL DUMP(X,Y,Z,ZETA,UBAR,VBAR,U,V,W,P,RHO,
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 arrays
C
      character*72 filename,meshname,optset
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C BEGINNING OF EDITING SECTION
C
C A. This section is used to specify what set of standard output files 
C     (see NML Memo DRL 11/11/91) the user desires.  As discussed in 
C     NML MEMO DRL 11/11/91 each type of data file has a standard set
C     of header lines. 
C
C      line 1: the geometric meshname.  This is obtained through the use
C               of the common block BLK1.  
C      line 2: the user's description of the file.  This output routine 
C               is designed to allow the user to specify the beginning 
C               of this line.  The code automatically adds to this 
C               beginning a phrase indicative of the results being 
C               written.  
C      other header lines:  certian types of files contain other header 
C               lines.  These are automatically accounted for by this 
C               routine.
C
C
C B.  Regardless of the output files desired, the user must specify the 
C       CHARACTER STRING begin and dimension it accordingly
C
C C. Two sets of user standard output files are possible:
C     1. optset='real'     =>  the user selects which of the following 
C                                real output files are desired: 
C                                  .vel Eulerian 3-D velocities
C                                  .v2r Eulerian average velocities
C                                  .v2r Eulerian surface velocities
C                                  .v2r Eulerian bottom velocities
C                                  .s2r tidal elevations
C     2. optset='complex'  =>  the user selects which of the following
C                                complex output files are desired: 
C                                  .vel 3-D velocities
C                                  .v2c average velocities
C                                  .s2c tidal elevations
C
C
C
      character*30 begin
      meshname='g2s.fromgom2'
      begin='g2s.fromgom2 barotropic tides/'
      optset='complex'
C
C END OF EDITING SECTION
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      write(*,*)'  Subroutine OUTPUT'
C
C Write real results files if optset='real'
C 
      if(optset.EQ.'real')then     
      write(*,*)'    Enter name of 3-D velocity .vel results file'
      write(*,*)'      (entering NO => no file will be created)'
      read(*,102)filename
      write(*,102)filename
      if(filename.NE.'NO')then        
         open (unit=7,file=filename,status='new')
         write(7,102)meshname
         write(7,102)begin//' 3-D velocity field'
         write(7,*)nnv
         write(7,*)'1'
         write(7,*)OMEGA
         write(7,8007)((i,z(i,j),
     &      CABS(u(i,j))*COS(phaselag(u(i,j))),0.0,
     &      CABS(v(i,j))*COS(phaselag(v(i,j))),0.0,
     &      CABS(w(i,j))*COS(phaselag(w(i,j))),0.0,j=1,nnv),i=1,nn)
         close(7)
      endif
C
      write(*,*)'    Enter name of the average velocity .v2r results fil
     &e'
      write(*,*)'      (entering NO => no file will be created)'
      read(*,102)filename        
      write(*,102)filename
      if(filename.NE.'NO')then        
         open (unit=8,file=filename,status='new')
         write(8,102)meshname
         write(8,102)begin//' average Eulerian velocities'
         write(8,8002)(I,REAL(UBAR(I)),REAL(VBAR(I)),i=1,nn)
         close(8)
      endif
C
      write(*,*)'    Enter name of the surface velocity .v2r results fil
     &e'
      write(*,*)'      (entering NO => no file will be created)'
      read(*,102)filename        
      write(*,102)filename
      if(filename.NE.'NO')then        
         open (unit=10,file=filename,status='new')
         write(10,102)meshname
         write(10,102)begin//' surface Eulerian velocities'
         write(10,8002)(I,REAL(U(I,NNV)),REAL(V(I,NNV)),i=1,nn)
         close(10)
      endif
C
      write(*,*)'    Enter name of the bottom velocity .v2r results file
     &'
      write(*,*)'      (entering NO => no file will be created)'
      read(*,102)filename        
      write(*,102)filename
      if(filename.NE.'NO')then        
         open (unit=9,file=filename,status='new')
         write(9,102)meshname
         write(9,102)begin//' bottom Eulerian velocities'
         write(9,8002)(I,REAL(U(I,1)),REAL(V(I,1)),i=1,nn)
         close(9)
      endif
C
      write(*,*)'    Enter name of elevation .s2r results file'
      write(*,*)'      (entering NO => no file will be created)'
      read(*,102)filename         
      write(*,102)filename
      if(filename.NE.'NO')then        
         OPEN (UNIT=19,FILE=filename,STATUS='NEW')
         write(19,102)meshname
         write(19,102)begin//' time-independent tidal elevations'
         write(19,8001)(i,cabs(zeta(i))*COS(phaselag(zeta(i))),i=1,nn)
         close(19)
      endif
C
C Write complex results files if optset='complex'
C
      else if(optset.EQ.'complex')then
      write(*,*)'    Enter name of 3-D velocity .vel results file'
      write(*,*)'      (entering NO => no file will be created)'
      read(*,102)filename
      write(*,102)filename
      if(filename.NE.'NO')then        
         open (unit=7,file=filename,status='new')
         write(7,102)meshname
         write(7,102)begin//' 3-D velocity field'
         write(7,*)nnv
         write(7,*)'1'
         write(7,*)OMEGA
         write(7,8007)((i,z(i,j),
     &                  CABS(u(i,j)),phaselag(u(i,j)),
     &                  CABS(v(i,j)),phaselag(v(i,j)),
     &                  CABS(w(i,j)),phaselag(w(i,j)),j=1,nnv),i=1,nn)
         close(7)
      endif
C
      write(*,*)'    Enter name of the average velocity .v2c results fil
     &e'
      write(*,*)'      (entering NO => no file will be created)'
      read(*,102)filename        
      write(*,102)filename
      if(filename.NE.'NO')then        
         open (unit=8,file=filename,status='new')
         write(8,102)meshname
         write(8,102)begin//' tidal velocities'
         write(8,*)OMEGA
         write(8,8004)(I,cabs(UBAR(I)),phaselagd(UBAR(I)),
     &cabs(VBAR(I)),phaselagd(VBAR(I)),I=1,nn)
         close(8)
      endif
C
      write(*,*)'    Enter name of the elevation .s2c results file'
      write(*,*)'      (entering NO => no file will be created)'
      read(*,102)filename         
      write(*,102)filename
      if(filename.NE.'NO')then        
         OPEN (UNIT=19,FILE=filename,STATUS='NEW')
         write(19,102)meshname
         write(19,102)begin//' tidal elevations'
         write(19,*)OMEGA
         write(19,8002)(i,cabs(zeta(i)),phaselagd(zeta(i)),I=1,nn)
         close(19)
      endif
C
C
C
      else
         write(*,*)'optset=',optset,' is not an available option.'
      endif
C
C FORMAT STATEMENTS
C
 102  format(a)
 8001 format(i5,1e12.4)
 8002 format(i5,2e12.4)
 8004 format(i5,4e12.4)
 8007 format(i5,7e12.4)
C
C End of routine
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
C End of User-Specified Subroutines for FUNDY5
C End of User-Specified Subroutines for FUNDY5
C***********************************************************************
C***********************************************************************
