C***********************************************************************
C***********************************************************************
C Beginning of Main Program for QUODDY5 (RELEASE 1.1)
C Beginning of Main Program for QUODDY5 (RELEASE 1.1)
C***********************************************************************
C***********************************************************************
      PROGRAM QUODDY5
C-----------------------------------------------------------------------
C Developed at Dartmouth College with the help of many colleagues from
C a number of other institutions.
C-----------------------------------------------------------------------
C README file: quoddy5_1.0_readme
C
C Please see the README file for an overview of the differences between 
C this version of QUODDY5 (i.e. RELEASE 1.1) and the release
C on which this revision is based (i.e. QUODDY5 RELEASE 1.0).
C
C The README file also contains information regarding the software
C development team as well as contact information for the individuals to
C whom it is suggested that questions initially be directed.
C-----------------------------------------------------------------------
C Arrays and Constants:
C
C  Index range: I=1,NN; J=1,NNV; JE=1,NEV; JL=1,NLEV; 
C               K1=1,NBE; K2=1,NGEO; L=1,NE;
C               M=1,NR; NB1=1,NHBW+1; NB2=1,2*NHBW+1
C
C   NN - the number of horizontal nodes
C   NNV - the number of nodes in the vertical
C   NEV - the number of elements in the vertical (NEV=NNV-1)
C   NLEV - the number of level surfaces for the baroclinic
C          pressure calculation
C   NBE - the number of boundary elements
C   NE - the number of horizontal elements
C   NHBW - the half bandwidth of elevation LHS matrix:
C             THETA=0.0 => gravity waves are handled explicitly 
C                       => NHBW=0
C             0.0<THETA<1.0 => gravity waves are handled semi-implicitly 
C                           => NHBW=half bandwidth of horizontal 
C                              finite element mesh
C   NR - the number of elevation boundary conditions
C
C +++++++++++++++++
C + INPUT SECTION +
C +++++++++++++++++
C
C STANDARD IN:
C
C   FILinq - name of the .inq file being used to conduct the simulation
C
C SUBROUTINE INPUTQ5:
C
C   FILicq - name of hot start .icq5 file
C            (set to 'COLD-START' if a cold start is desired)
C   IDD0,IMM0,IYEAR0,SEC0 - timing parameters for the beginning
C                           of the simulation (only set set if 
C                           a cold start simulation is being
C                           conducted)
C   IDDL,IMML,IYEARL,SECL - timing parameters for the end
C                           of the simulation
C   DEGLAT - degrees latitude to be used for coriolis terms (f-plane)
C   DelT - time step
C   NNV - number of vertical nodes
C
C      Horizontal Mesh Information:
C
C   MESHNAME - name of the horizontal mesh
C   NN - number of horizontal nodes
C   NE - the number of horizontal 2-D triangular elements
C   X(I),Y(I) - cartesian nodal coordinates
C   IN(3,L) - LABLIB STANDARD triangular element incidence list
C   Hdown(I) - Bathymetric depth at nodes
C   SUMx(I) - cosine value for boundary nodes where V.n is 
C           specified
C   SUMy(I) - sine value for boundary nodes where V.n is 
C           specified
C   DSbdry(I) - surface integral information
C   INBE(2,K1) - LABLIB STANDARD boundary element incidence list
C   NBE - number of boundary elements
C   INGEO(2,K2) - LABLIB STANDARD geostrophic boundary element
C           incidence list
C   NGEO - number of geostrophic boundary elements
C   NORM(I) - boundary node numbers where normal velocity is set
C           I =1,NVN       are V.n .NE. 0.0    (NCD=3)
C           I =KST,NNdim   are V.n=0.0         (NCD=1,2,6)
C   LRVS(I) - node numbers of set depth and set velocity nodes
C           I=1,NR         are set depth nodes (NCD=5,6,9)
C           I=NR+1,NR+NPCN are V=0.0           (NCD=10)
C           I=LST,NNdim    are V .NE. 0.0      (NCD=8,11)
C
C       Vertical Level Mesh Information:
C
C   NLEV - number of level surfaces for baroclinic pressure
C          calculation
C   ZL(I,JL) - vertical position of level surfaces for baroclinic
C        pressure calculation
C
C MAIN PROGRAM (may be overwritten individually in SUBROUTINE SETQ5):
C
C   T0 - Reference temperature
C   S0 - Reference salinity
C   G - Gravitational acceleration
C   OMEGA  - Radian frequency (rad/sec)
C   NONLIN - integer flag to indicate whether or not nonlinear
C            advection is to be included in the momentum equation:
C              NONLIN=0 => no advection in momentum equation
C              NONLIN=1 => use advection in momemtum equation
C   NLBS - integer flag to indicate linearized bottom stress
C          coefficient (aK) parameterizaton in horizontal momentum
C          equations (bottom stress=aK*bottom velocity):
C              NLBS=0 => apply constant aK=aKMIN
C              NLBS=1 => apply quadratic stress condition
C                        => aK=Cd*bottom speed
C   Cd - bottom slip coefficient; used when NLBS=1
C   ISMAG - integer flag to indicate type of horizontal 
C           viscosity/diffusivity implemented within horizontal 
C           diffusion operator:
C              ISMAG=0 => apply uniform viscosity/diffusivity 
C                         (equal to AH)
C              ISMAG=1 => viscosity/diffusivity based on the horizontal 
C                         shear and grid scale (i.e. similar in spirit  
C                         to Smagorinsky(1963)), with multiplicative
C                         constant set to AH.
C   AH - real data utilized for horizontal viscosity/diffusivity closure
C        (see ISMAG for implementation)
C   AHMIN - minimum horizontal viscosity/diffusivity (only implemented 
C           when ISMAG=1)
C   PRESSURE - Character string to indicate pressure computation:
C               'BAROTROPIC' or 'BAROCLINIC'
C   MASSVAR - Character string to indicate treatment of mass variables:
C             'TWO-PROG' => Temperature and salinity transport equations
C                           are solved every time-step and 
C                           SUBROUTINE EQUSTATE2 is used 
C                           to update RHO(T,S) after each time-step
C             'TWO-DIAG' => Initial values of temperature and salinity 
C                           are used by SUBROUTINE EQUSTATE2 to set
C                           RHO(T,S), which is subsequently held
C                           constant throughout the simulation
C             'ONE-PROG' => Temperature transport equation is solved
C                           every time-step and SUBROUTINE EQUSTATE1 is 
C                           used to update RHO(T) after each time-step
C             'ONE-DIAG' => Initial values of temperature are used by
C                           SUBROUTINE EQUSTATE1 to set RHO(T), 
C                           which is subsequently held
C                           constant throughout the simulation
C             'ZERO'=> No mass variables are implemented => temperature,
C                      salinity, and RHO are homogeneous and fixed
C   IHHBC - Binary switch to indicate BC type for temperature and
C           salinity at horizontal open boundaries:
C             1 => Fixed at the initial conditions
C             2 => No flux BC
C   CLOSURE - Character string to indicate turbulence closure:
C             'CONSTANT' => vertical eddy viscosity/diffusivities set
C                           to respective mins: EKMMIN, EKHMIN, EKQMIN
C             'MY25' => vertical eddy viscosity/diffusivities computed
C                       using Mellor and Yamada(1984,1982) turbulence
C                       closure with adjustments of Galperin et. al.
C                       (1988) and Blumberg et. al. (1992)
C   Q20 - Reference value of Q2 used to initialization
C   Q2L0 -  Reference value of Q2L used to initialization
C   ZLOGBOT - Height of the bottom node in the computational domain
C             above the true bottom
C   ZLOGTOP - Height of the top node in the computational domain
C             below the true surface
C   icqADVDIF - Binary switch to indicate treatment of advection and
C               horizontal diffusion in the Q2 and Q2L transport
C               equations:
C               0 => No advection and horizontal diffusion
C               1 => Include advection and horizontal diffusion
C   IQ2TBC  - Binary switch to indicate BC type for the Q2 transport
C             equation at the top of the water column:
C             1 => Type 1 BC given in Mellor and Yamada(1982)
C             2 => No flux BC
C   IQ2BBC  - Binary switch to indicate BC type for the Q2 transport
C             equation at the bottom of the water column:
C             1 => Type 1 BC given in Mellor and Yamada(1982)
C             2 => No flux BC
C   IQ2LTBC - Binary switch to indicate BC type for the Q2L transport
C             equation at the top of the water column:
C             1 => Type 1 BC given in Naimie(1995)[Ph.D. Thesis]
C             2 => No flux BC
C   IQ2LBBC - Binary switch to indicate BC type for the Q2L transport
C             equation at the bottom of the water column:
C             1 => Type 1 BC given in Naimie(1995)[Ph.D. Thesis]
C             2 => No flux BC
cJAS \/
C   IUSTARQ2 - Switch to determine usage of friction velocity
C              in the surface BC for Q2 and Q2L
C              0 => Use windstress in surface bcs (canonical setting)
C              1 => Use friction velocity in surface bcs
C              2 => q2l forcing mult by 3.
cJas /\
C   Tau0 - constant for weight applied to the continuity equation
C          when it is added to the vertically averaged momentum
C          equation to obtain the wave equation (typically .0002)
C   EPSN - implicit time weighting factor for vertical diffusion
C          of momentum:
C              EPSN=0.0 => explicit (at time level K)
C              EPSN=0.5 => centered (at time level K+1/2)
C              EPSN=1.0 => implicit (at time level K+1)
C   EPSH - implicit time weighting factor for vertical diffusion
C          of temperature and salinity:
C              EPSH=0.0 => explicit (at time level K)
C              EPSH=0.5 => centered (at time level K+1/2)
C              EPSH=1.0 => implicit (at time level K+1)
C   EPSQ - implicit time weighting factor for vertical diffusion
C          of Q2 and Q2L:
C              EPSQ=0.0 => explicit (at time level K)
C              EPSQ=0.5 => centered (at time level K+1/2)
C              EPSQ=1.0 => implicit (at time level K+1)
C   THETA - implicit gravity wave factor for the wave equation:
C              THETA=0.0 => explicit (at time level K-1)
C              THETA=1.0 => implicit (at time level K)
C   aKMIN - minimum linearized bottom stress coefficient
C   EKMmin - minimum vertical eddy diffusivity for momentum.
C   EKHmin - minimum vertical eddy diffusivity for heat and salt.
C   EKQmin - minimum vertical eddy diffusivity for turbulence.
C   Q2min,Q2Lmin,ELLmin - Minimum values of Q2, Q2L, and the mixing
C                         length.  Used only when CLOSURE='MY25'
C
C SUBROUTINE SETQ5: Allows the user to overwrite any of the following
C                   default settings:
C
C     T0,S0,G,OMEGA,NONLIN,NLBS,Cd,ISMAG,AH,
C     PRESSURE,MASSVAR,IHHBC,CLOSURE,
C     Q20,Q2L0,ZLOGBOT,ZLOGTOP,IQADVDIF,IQ2TBC,IQ2BBC,IQ2LTBC,IQ2LBBC,
C     IUSTARQ2,Tau0,EPSN,EPSH,EPSQ,THETA,aKMIN,Q2min,Q2Lmin,ELLmin,
cJAS  ^^^^^^^^
cJAS  EKMMIN,EKHMIN,EKQMIN)
C
C MAIN PROGRAM:
C
C   COR - Coriolis parameter (f-plane)
C   NEV - Number of vertical elements (NEV=NNV-1)
C   NHBW - Half bandwidth of elevation LHS matrix (Note: NHBW=0 when
C            gravity waves are handled explicitly [i.e. THETA=0]
C            => NBW=1)
C   NBW - Minimum allowable dimension size of NBWE
C 
C ++++++++++++++++++++++++++
C + INITIALIZATION SECTION +
C ++++++++++++++++++++++++++
C
C   IDD0,IMM0,IYEAR0,SEC0 - timing parameters for the beginning
C                           of the simulation, as set in the
C                           .icq5 hot-start file (the input
C                           values are only overwritten
C                           if a hot-start is begin conducted)
C
C   *OLD arrays are initialized with dependent variable data at one
C   time step prior to the beginning of the simulation  
C
C   ZETAOLD(I) - sea surface elevation
C   Uold(I),Vold(I) - x,y components of vertically averaged
C               horizontal velocity
C   Zold(I,J) - 3-D nodal grid
C
C   *MID arrays are initialized with dependent variable data at the
C   beginning of the simulation  
C
C   ZETAMID(I) - sea surface elevation
C   Umid(I),Vmid(I) - x,y components of vertically averaged
C                     horizontal velocity
C   Zmid(I,J) - 3-D nodal grid
C   UZmid(I,J),VZmid(I,J),WZmid(I,J) - x,y,z components of 
C                                      velocity
C   Q2mid(I,J),Q2Lmid(I,J) - twice the turbulent kinetic energy
C                            and the master length scale times
C                            twice the turbulent kinetic energy
C   Tmid(I,J),Smid(I,J) - temperature and salinity
C
C       Other Information:
C
C   BCFTR(I)  - Real binary multiplicative factor which deletes
C               the Coriolis term in momentum equation at solid
C               boundaries
C           BCFTR(I)=0.0 for land boundary nodes
C           BCFTR(I)=1.0 for all other nodes
C   INCONV(I) - Convection integer.
C           NONLIN=1 => INCONV(I)=1 for interior nodes
C                       INCONV(I)=0 for boundary nodes
C           NONLIN=0 => INCONV(I)=0 for all nodes
C
C ++++++++++++++++++++++++++++++++++++++++
C + STATIONARY ARRAY CALCULATION SECTION +
C ++++++++++++++++++++++++++++++++++++++++
C
C SUBROUTINE STATIONARYQ5:
C
C   AR(L) - elemental areas
C   DX(3,L) - element based array of DELTAx for each local node
C   DY(3,L) - element based array of DELTAy for each local node
C   NFTR - sparse factor
C   IQ(I),JQ(I*NFTR) - SPRSPAK STANDARD pointer arrays
C   PPx(I*NFTR) - SPRSPAK STANDARD derivative matrix  <Pi*dP/dxj>
C   PPy(I*NFTR) - SPRSPAK STANDARD derivative matrix  <Pi*dP/dyj>
C   AGPGP(I*NFTR) - SPRSPAK STANDARD Laplacian matrix
C                  AGPGP=-<AH*grad(phi_i) . grad(phi_j)>
C                  (note: only computed if ISMAG=0 => AH=constant)
C   SV(I) - stationary diagonal mass matrix
C   SH(NB1,I) - transpose of the stationary depth matrix
C   DIRICH(NB1,M) - matrix used for the assymetric terms resulting
C    from the elevation boundary conditions.
C    (note:  The current implementation renders SH symmetric and
C     DIRICH assymetric, which facilitates the storage of SH in
C     NB1=NHBW+1 [instead of NB2=2*NHBW+1] columns.  To implement SH
C     in its nonsymmetric form [and hence eliminate the need for
C     DIRICH], the following evolution must be conducted in
C     the MAIN PROGRAM, SUBROUTINE STATIONARYQ5, and SUBROUTINE
C     ELEVATIONQ5:
C       1. Uncomment the currently commented lines between the
C          first C+++ and the C>>> for each case of 2.)
C       2. Comment Out (or remove) lines between the C>>> and C<<<
C          comment lines.
C
C +++++++++++++++++++++++++++++++++++++++++
C + PREPARATION FOR TIME STEPPING SECTION +
C +++++++++++++++++++++++++++++++++++++++++
C
C SUBROUTINE ATMOSQ5: Initialize Arrays with values at the beginning
C                     of the simulation (i.e. time level K for the first
C                     time step).
C
C   ATMxMID(I),ATMyMID(I) - X,Y directed wind stress in pascals
C   ATEMP(I),BTEMP(I) - Terms which set the type 3 boundary condition
C                     for temperature; which is set in accordance
C                     with the expression:
C                       NdT/dz=h*(Te-T)
C                     thus:
C                       h=boundary layer coeff
C                       Te=equilibruim value of T
C                       ATEMP=h                  
C                       BTEMP=h*Te
C   PMEmid(I) - Net precipitation minus evaporation
C
C SUBROUTINE POINTSOURCEQ5:Initialize Arrays with values at the
C                          beginning of the simulation (i.e. time level
C                          K for the first time step).
C
C   SRCmid(I,JE)   - Volume source strength (m^3/sec)
C                       at horizontal node I and vertical element JE.
C   UZSRCmid(I,NNV)- Real binary flag to indicate neutrality
C                      of source being specified at horizontal node I
C                      with respect to momentum:
C                      0.0 => source is neutral wrt momentum
C                      1.0 => source is non-neutral wrt momentum
C   UZSRCmid(I,JE) - X-directed component of velocity for source
C                       at horizontal node I and vertical element JE.
C   VZSRCmid(I,JE) - Y-directed component of velocity for source
C                       at horizontal node I and vertical element JE.
C   TSRCmid(I,NNV) - Real binary flag to indicate neutrality
C                      of source being specified at horizontal node I
C                      with respect to temperature:
C                      0.0 => source is neutral wrt temperature
C                      1.0 => source is non-neutral wrt temperature
C   TSRCmid(I,JE)  - Temperature for source
C                       at horizontal node I and vertical element JE.
C   SSRCmid(I,NNV) - Real binary flag to indicate neutrality
C                      of source being specified at horizontal node I
C                      with respect to salinity:
C                      0.0 => source is neutral wrt salinity
C                      1.0 => source is non-neutral wrt salinity
C   SSRCmid(I,JE) - Salinity for source
C                       at horizontal node I and vertical element JE.
C
C MAIN PROGRAM:
C   
C   SRCSUMmid(I) - Vertical integral of SRCmid
C   SRCSUMold(I) - Set to initial conditions (i.e. SRCSUMmid(I))
C   PMEold(I) - Set to for initial conditions (i.e. PMEmid(I))
C   KDmid - Set to initial Gregorian Day
C   SECmid - Set to starting time on initial Gregorian Day
C   ITER - Set to 0
C
C +++++++++++++++++++++++++
C + TIME STEPPING SECTION +
C +++++++++++++++++++++++++
C
C MAIN PROGRAM:
C
C   SECnew - Incremented by DelT
C   KDnew - Incremented as necessary (which also affects SECnew)
C   ITER - Incremented by 1
C
C SUBROUTINES ATMOSQ5 AND POINTSOURCEQ5: Call subroutines to load
C values for arrays at the end of the current time step.  See comments
C above (or in subroutines themselves) for description of various
C array contents:
C
C   ATMxNEW(I),ATMyNEW(I),ATEMPnew(I),BTEMPnew(I),PMEnew(I)
C   SRCnew(I,J),UZSRCnew(I,J),VZSRCnew(I,J),TSRCnew(I,J),SSRCnew(I,J)
C
C MAIN PROGRAM:
C
C   SRCSUMnew(I) - Vertical integral of SRCnew(I)
C
C VARIOUS SUBROUTINES (AS NECESSARY): All at time level K:
C
C   aK(I) - Linearized bottom stress coefficients
C   BPGx(I,J),BPGy(I,J) - node based baroclinic pressure gradients
C   HRBARx(L),HRBARy(L) - element based vertically integrated 
C                         baroclinic pressure gradients
C   AGPGP(I*NFTR) - SPRSPAK STANDARD Laplacian matrix
C                  AGPGP=-<AH*grad(phi_i) . grad(phi_j)>
C                  (note: only computed if ISMAG=1)
C   AHI(I) - nodal values of horizontal eddy viscosity/diffusivity
C   VISCX(J),VISCY(J) - horizontal viscous diffusion terms
C   CONxZ(I,J),CONyZ(I,J) - node based advection plus horizontal
C                           diffusion
C   CONx(I),CONy(I) - node based vertical integrals of advection plus
C                     horizontal diffusion
C
C SUBROUTINE ELEVATIONQ5:
C
C   ZETANEW(I) - sea surface elevation at time level K+1
C   Znew(I,J) - vertical grid at time level K+1/2
C
C SUBROUTINE VERTICALQ5 (after appropriate manipulation of
C atmpspheric and point source terms to ensure forcing occurs at time 
C level K+1/2):
C
C  UZnew(I,J),VZnew(I,J) - x,y components of velocity at time level K+1
C  Q2new(I,J),Q2Lnew(I,J) - turbulence parameters at time level K+1
C  Tnew(I,J),Snew(I,J) - temperature and salinity at time level K+1
C  ENZM(I,JE),ENZH(I,JE),ENZQ(I,JE) - vertical eddy viscosities and
C                                     diffusivities at time level K+1
C
C MAIN PROGRAM:
C
C  *MID=>*OLD
C  *NEW=>*MID
C
C SUBROUTINES VERTVEL3_2 AND VERTAVG: At the end of the current time 
C step:
C
C  WZmid(I,J) - vertical velocity
C  Umid(I,J),Vmid(I,J) - vertically averaged horizontal velocity
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C ARRAY DECLARATION SECTION
C
C Objectives: Declare all arrays required by the main program.
C
C Declare all the integer program parameters via INCLUDE file
C
      INCLUDE 'DCMS.DIM'
C
C 2-D array declaration
C
      INTEGER IN(3,NEdim),INBE(4,NBEdim),INGEO(2,NBEdim),
     &LRVS(NNdim),NORM(NNdim)
      REAL ZETAOLD(NNdim),Uold(NNdim),Vold(NNdim),
     &ZETAMID(NNdim),Umid(NNdim),Vmid(NNdim),HMID(NNDIM),
     &ZETANEW(NNdim),Unew(NNdim),Vnew(NNdim),
     &ATMxMID(NNdim),ATMyMID(NNdim),ATMxNEW(NNdim),ATMyNEW(NNdim),
     &aK(NNdim),Hdown(NNdim),X(NNdim),Y(NNdim),
     &SUMx(NNdim),SUMy(NNdim),SV(NNdim),DSbdry(NNdim),DX(3,NEdim),
     &DY(3,NEdim),AR(NEdim),DELX(3),DELY(3),SVHmid(NNdim),
     &SUMx1(NNdim),SUMy1(NNdim),DSbdry1(NNdim),QP(NNdim),
     &CS(NNdim),COR(NNdim),TS(NNdim),XNOD(NNdim),YNOD(NNdim)
C
C 3-D array declaration
C
      REAL Zold(NNdim,NNVdim),Zmid(NNdim,NNVdim),Znew(NNdim,NNVdim),
     &UZnew(NNdim,NNVdim),UZmid(NNdim,NNVdim),
     &VZnew(NNdim,NNVdim),VZmid(NNdim,NNVdim)
C
C 1-D array declaration
C
      REAL ZL(NLEVdim)
C
C Array declaration for sparse matrix storage
C
      INTEGER IQ(NNdim),JQ(NFTRdim*NNdim)
      REAL PPx(NFTRdim*NNdim),PPy(NFTRdim*NNdim)
C
C Convective array declaration
C
      REAL DUZDX(NNVdim),DUZDY(NNVdim),DVZDX(NNVdim),DVZDY(NNVdim),
     &DZDX(NNVdim),DZDY(NNVdim),WZmid(NNdim,NNVdim)
C
C BCFTR = 0 ON SOLID BOUNDARIES; USED TO ZERO THE CORIOLIS EFFECT
C CONxZ, CONyZ ARE THE STORED 3D CONVECTIVE PLUS VISCOUS TERMS; 
C CONx, CONy ARE THEIR VERTICAL INTEGRALS
C VERTFLX IS THE 3D ARRAY OF VERTICAL FLUX IN THE CONVECTIVE TERMS
C AGPGP IS THE SPARSE VISCOUS MATRIX
C VISCX, VISCY ARE THE VISCOUS TERMS (ADDED INTO CONxZ, CONyZ)
C
      REAL BCFTR(NNdim),CONxZ(NNdim,NNVdim),CONyZ(NNdim,NNVdim),
     &CONx(NNdim),CONy(NNdim),VERTFLX(NNdim,NNVdim),
     &AGPGP(NFTRdim*NNdim),VISCX(NNVdim),VISCY(NNVdim)
C
C Declaration statements for BANSOLTR (symatrix) scheme
C a) The first dimension of SH(NBWE,NNdim) needs to be changed to
C    SH(NB1dim,NNdim) for the symatrix scheme
C b) DIRICH is needed only in the symatrix case
C
C+++Asymmetric storage
C+++      REAL SH(NBWE,NNdim)
C>>>Symmetric storage in symatrix
Ccgm-\/
      REAL SH(NB1dim,NNdim),DIRICH(NBWE,NRdim)
Ccgm-/\
Ccgm Symmetric and sparse storage
Ccgm+      REAL SH(NFTRdim*NNdim),DIRICH(NBWE,NRdim)
Ccgm
C<<<
C
C Variable declaration for the MY25-Galperin turbulence closure
C
      CHARACTER*72 CLOSURE
      REAL Q2new(NNdim,NNVdim),Q2Lnew(NNdim,NNVdim),
     &Q2mid(NNdim,NNVdim),Q2Lmid(NNdim,NNVdim)
C
C Variable declaration for the prognostic density calculation
C
      CHARACTER*72 PRESSURE
      REAL RHOmid(NNdim,NNVdim)
      REAL BPGx(NNdim,NNVdim),BPGy(NNdim,NNVdim)
      REAL HRBARx(NEdim),HRBARy(NEdim)
C
C Variable declaration for mass variables
C
C   Tmid,Tnew: Temperature,degrees C
C   Smid,Snew: Salinity, PSU
C
      CHARACTER*72 MASSVAR
      REAL Tmid(NNdim,NNVdim),Tnew(NNdim,NNVdim),
     &Smid(NNdim,NNVdim),Snew(NNdim,NNVdim)
C
C Array declaration for source terms
C
      REAL SRCmid(NNdim,NNVdim),SRCnew(NNdim,NNVdim),
     &     UZSRCmid(NNdim,NNVdim),VZSRCmid(NNdim,NNVdim),
     &     TSRCmid(NNdim,NNVdim),SSRCmid(NNdim,NNVdim)
      REAL SRCSUMold(NNdim),SRCSUMmid(NNdim),SRCSUMnew(NNdim)
      REAL PMEold(NNdim),PMEmid(NNdim),PMEnew(NNdim)
C
C Array declaration for the 2-D Smagorinsky type horizontal 
C viscosity/diffusivity
C
      REAL AHI(NNdim)
C
C Array declaration for the 3-D vertical diffusion coefficients
C   (element based)
C
      REAL ENZM(NNdim,NNVdim),ENZH(NNdim,NNVdim),ENZQ(NNdim,NNVdim)
C
C 2-D integer array to turn off only the convective terms at boundary
C   nodes; and everywhere if NONLIN=0
C
      INTEGER INCONV(NNdim)
C
C 2-D arrays to store the type III BC. coefficients for temperature
C
      REAL ATEMPmid(NNdim),BTEMPmid(NNdim),
     &     ATEMPnew(NNdim),BTEMPnew(NNdim)
C
C Variable declaration for I/O formatting
C
      CHARACTER*72 MESHNAME,FILicq,FILinq,COORDS
C

cJAS \/ Further declarations  
      real EKMMIN,EKHMIN,EKQMIN

cJas /\

      REAL ZERO
      DATA ZERO/0.0/
      PI=ACOS(-1.0)
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C INPUT SECTION
C
C Objectives: Read .inq file, input mesh geometry data, open and begin
C writing to echo file, and set all simulation control parameters.
C
C Set default values for standard switches and parameters
C
C  Standard physical constants (REAL)
C
      T0=10.0
      S0=32.0
      G=9.806
      OMEGA=0.00007292
      RAD=6367500.0 
C
C  Binary integer flags and constants needed due to physical assumptions
C  regarding inclusion of nonlinearities
C
      NONLIN=1
      NLBS=1
      Cd=0.005
      ISMAG=1
      AH=0.28
C
C  Physics switches (CHARACTER STRINGS) and horizontal boundary
C  condition option to indicate desired mass variable physics
C
      PRESSURE='BAROCLINIC'
      MASSVAR='TWO-PROG'
      IHHBC=2
C
C  Physics switch (CHARACTER STRING), physical constants (REAL),
C  and boundary condition options (INTEGERS) for turbulence closure
C
      CLOSURE='MY25'
      Q20=0.001
      Q2L0=0.001
      ZLOGBOT=1.0
      ZLOGTOP=1.0
      IQADVDIF=1
      IQ2TBC=2
      IQ2BBC=1
      IQ2LTBC=2
      IQ2LBBC=1
C
C  Numerical Model weighting parameters (REAL)
C
      Tau0=0.0002
      EPSN=0.50
      EPSH=1.0
      EPSQ=1.0
      THETA=0.75
C
C  Minimum parameters
C
      EKMMIN=0.00020
      EKHMIN=0.000020
      EKQMIN=0.000020
      AKMIN=0.00001
      Q2MIN= 0.000001
      Q2LMIN=0.000001
      ELLMIN=0.000001
      AHMIN=0.0
C 
C Call SUBROUTINE SETQ5 to enable the user to change standard 
C  parameters, constants, and flags.
C
      CALL SETQ5(T0,S0,G,OMEGA,NONLIN,NLBS,Cd,ISMAG,AH,AHMIN,
     &PRESSURE,MASSVAR,IHHBC,CLOSURE,Q20,Q2L0,
     &ZLOGBOT,ZLOGTOP,IQADVDIF,IQ2TBC,IQ2BBC,IQ2LTBC,IQ2LBBC,IUSTARQ2,
     &Tau0,EPSN,EPSH,EPSQ,THETA,aKMIN,Q2min,Q2Lmin,ELLmin,
     &EKMMIN,EKHMIN,EKQMIN)
C
C Enter the name of the input file (".inq") of the simulation
C
      WRITE(*,'(a)') 'Enter the name of the input file (".inq")'
      READ (*,'(a)') FILinq
      Iinq=INDEX(FILinq,' ')-1
C
C Load run information from .inq and associated input files
C
      CALL INPUTQ5(FILinq,
     &FILicq,IDD0,IMM0,IYEAR0,SEC0,IDDL,IMML,IYEARL,SECL,
     &DEGLAT,DelT,NNV,
     &MESHNAME,NN,NE,X,Y,IN,Hdown,SUMx,SUMy,DSbdry,SUMx1,SUMy1,DSbdry1,
     &INBE,NBE,INGEO,NGEO,NORM,NVN,KST,LRVS,NR,NPCN,LST,NLEV,ZL,
     &RAD,COR,CS,TS,COORDS,OMEGA)
C
C Fill nodal arrays (XNOD,YNOD) with horizontal coordinates as input 
C from the .nod file 
C
      IF(COORDS.EQ.'SPHERICAL')THEN
         FCTR=180.0/PI/RAD
      ELSEIF(COORDS.EQ.'CARTESIAN')THEN
         FCTR=1.0
      ELSE
        WRITE(*,*)'STOPPING - Unrecognized specification of coordinates'
        WRITE(*,*)'COORDS = ',COORDS
        WRITE(*,*)'COORDS  should be specified SPHERICAL or CARTESIAN'
        STOP
      ENDIF
      DO I=1,NN
         XNOD(I)=FCTR*X(I)
         YNOD(I)=FCTR*Y(I)
      ENDDO
C
C Set number of elements in the vertical
C
      NEV=NNV-1
C
C Compute half bandwidth from the element incidence list, overwrite
C  it with zero if an explicit simulation is being conducted, then
C  perform check to ensure dimensioning parameter NBWE is of sufficient
C  size
C
      CALL HALFBW(IN,3,NE,NHBW)
      IF(THETA.EQ.ZERO)NHBW=0
      NBW=2*NHBW+1
      IF(NBW.GT.NBWE)THEN
        WRITE(2,*)'EXECUTION TERMINATED DUE TO INSUFFICIENT DIMENSIONS:'
        WRITE(2,*)' NBW   =',NBW,   ',NBWE   =',NBWE
        STOP
      END IF
C
C Write plethora of settings to echo file 
C and check to ensure that settings are allowable options!
C
      CALL ECHOSETQ5(T0,S0,G,OMEGA,RAD,COR(1),COORDS,NHBW,
     &NONLIN,NLBS,Cd,ISMAG,AH,AHMIN,PRESSURE,MASSVAR,IHHBC,CLOSURE,
     &Q20,Q2L0,ZLOGBOT,ZLOGTOP,IQADVDIF,IQ2TBC,IQ2BBC,IQ2LTBC,IQ2LBBC,
     &Tau0,EPSN,EPSH,EPSQ,THETA,aKMIN,Q2min,Q2Lmin,ELLmin,NEV,
     &EKMMIN,EKHMIN,EKQMIN)
C
      WRITE(2,'(/A)')'INPUT SECTION SUCCESSFULLY COMPLETED'
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C INITIALIZATION SECTION
C
C Objectives: Initialize dependent variable arrays using either
C *COLD-START* or *HOT-START* information.  Then correct for physics
C switch settings and the desired turbulence parameterization.
C
C Initialize dependent variables
C
      CALL INITIALIZEQ5(MESHNAME,FILicq,NN,NNV,Hdown,DelT,
     &             T0,S0,Q20,Q2L0,IDD0,IMM0,IYEAR0,SEC0,
     &             ZETAOLD,Uold,Vold,Zold,ZETAMID,Umid,Vmid,
     &             Zmid,UZmid,VZmid,WZmid,Q2mid,Q2Lmid,Tmid,Smid,QP)
C
C Now apply corrections to be consistent with physics switch settings
C and desired turbulence parameterization (arrays set to reference
C values are fixed throughout the simulation).
C
C  Baroclinic Pressure Gradients
C
      IF(PRESSURE.EQ.'BAROTROPIC')THEN
         DO J=1,NNV
            DO I=1,NN
               BPGx(I,J)=0.0
               BPGy(I,J)=0.0
            ENDDO
         ENDDO
         DO L=1,NE
            HRBARx(L)=0.0
            HRBARy(L)=0.0
         ENDDO
      ENDIF
C
C  Mass Variables
C
      IF((MASSVAR.EQ.'TWO-PROG').OR.(MASSVAR.EQ.'TWO-DIAG'))THEN
         CALL EQSTATE2_2D(T0,S0,NN,NNV,TMID,SMID,RHOMID)
      ELSE IF((MASSVAR.EQ.'ONE-PROG').OR.(MASSVAR.EQ.'ONE-DIAG'))THEN
         CALL EQSTATE1_2D(T0,NN,NNV,TMID,RHOMID)
         DO J=1,NNV
            DO I=1,NN
               Smid(I,J)=S0
            ENDDO
         ENDDO
      ELSE IF((MASSVAR.EQ.'ZERO'))THEN
         DO J=1,NNV
            DO I=1,NN
               Tmid(I,J)=T0
               Smid(I,J)=S0
               RHOmid(I,J)=0.0
            ENDDO
         ENDDO
      ELSE
cJas \/  
         print *, 'MASSVAR = ', MASSVAR, ' not set correctly.'
         stop
cJas /\
      ENDIF
C
C  Constant linearized bottom stress coefficient
C
      IF(NLBS.EQ.0)THEN
         DO I=1,NN
            aK(I)=aKMIN
         ENDDO
      ENDIF
C
C  Constant vertical eddy diffusivities for momentum and heat (assuming
C  Prantdl Number=1):
C    ENZ*(I,1 to NEV) set to desired values
C    ENZ*(I,NNV) set to zero (not used during simulation)
C    Q2*(I,1 to NNV) set to minimum values (not used during simulation)
C
      IF(CLOSURE.EQ.'CONSTANT')THEN
         DO J=1,NEV
            DO I=1,NN
               ENZM(I,J)=EKMmin
               ENZH(I,J)=EKHmin
               ENZQ(I,J)=EKQmin
            ENDDO
         ENDDO
         DO I=1,NN
            ENZM(I,NNV)=0.0
            ENZH(I,NNV)=0.0
            ENZQ(I,NNV)=0.0
         ENDDO
         DO J=1,NNV
            DO I=1,NN
               Q2mid(I,J)=Q2min
               Q2Lmid(I,J)=Q2Lmin
            ENDDO
         ENDDO
      ENDIF
C                
C Load the BCFTR array to give the correct Coriolis factor of either
C   BCFTR(I) = 0 for solid boundary node with V.n = 0,
C   BCFTR(I) = 1 otherwise;for the (u + j*v) solution calculations
C
      DO I=1,NNdim
         BCFTR(I)=1.0
      ENDDO
      DO I=KST,NNdim
         NBC1=NORM(I)
         BCFTR(NBC1)=0.0
      ENDDO
C
C Load the INCONV array:
C   INCONV(I)=0 for zero convection terms at node I
C   INCONV(I)=1 convection term is included in node I
C   
C   NONLIN=0 => INCONV(I)=0 for all nodes
C   NONLIN=1 => INCONV(I)=0 for all boundary nodes
C               INCONV(I)=1 for all interior
C
      DO I=1,NNdim
         INCONV(I)=NONLIN
      ENDDO
      IF(NVN.GE.1)THEN
         DO I=1,NVN
            INCONV(NORM(I))=0
         ENDDO
      ENDIF
      IF(NNdim.GE.KST)THEN
         DO I=KST,NNdim
            INCONV(NORM(I))=0
         ENDDO
      ENDIF
      IF(NR.GE.1)THEN
         DO I=1,NR
            INCONV(LRVS(I))=0
         ENDDO
      ENDIF
      IF(NPCN.GE.1)THEN
         DO I=NR+1,NR+NPCN
            INCONV(LRVS(I))=0
         ENDDO
      ENDIF
      IF(NNdim.GE.LST)THEN
         DO I=LST, NNdim
            INCONV(LRVS(I))=0
         ENDDO
      ENDIF
      IF(NGEO.GE.1)THEN
         DO L=1,NGEO
            INCONV(INGEO(1,L))=0
            INCONV(INGEO(2,L))=0
         ENDDO
      ENDIF
C
      WRITE(2,'(/A)')'INITIALIZATION SECTION SUCCESSFULLY COMPLETED'
C
C Write timing parameters to echo file (not done until now because the
C .icq5 file provides starting time for hot starts)
C
      WRITE(2,'(/A)')'Time-stepping parameters:'
      WRITE(2,'(5X,A35,1X,I2,1X,I2,3X,I4,1X,F18.8)')
     &'Simulation start date and time     ',IDD0,IMM0,IYEAR0,SEC0
      WRITE(2,'(5X,A35,1X,I2,1X,I2,3X,I4,1X,F18.8)')
     &'Simulation end date and time       ',IDDL,IMML,IYEARL,SECL
      WRITE(2,'(5X,A35,1X,E12.4)')
     &'Time step (seconds)                ', DelT
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C START/STOP TIME SPECIFICATION SECTION
C
C Objectives: Set timing parameters for later use (these are fixed
C throughout the simulation).
C
C Starting time:
C  IDD0,IMM0,IYEAR0=starting date of the simulation (day,month,year)
C  KD0 is the day number starting date based on the
C  the Gregorian calendar (KD0=1 is January 1, 0000)
C  SEC0 starting time, in seconds after midnight, on KD0
C Stopping time:
C  (IDDL,IMML,IYEARL,KDL,SECL) are similar parameters to those above
C  86400 is the number of seconds in 1 day
C
C Get KD0/KDL from starting/stopping date and time
      CALL GDAY(IDD0,IMM0,IYEAR0,KD0)
      IF(SEC0.GE.86400.)CALL UP_DATE(KD0,SEC0)
      CALL GDAY(IDDL,IMML,IYEARL,KDL)
      IF(SECL.GE.86400.)CALL UP_DATE(KDL,SECL)
C
      WRITE(2,'(/A)')'START/STOP TIME SPECIFICATION SECTION SUCCESSFULLY
     & COMPLETED'
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C STATIONARY ARRAY CALCULATION SECTION
C
C Objectives: Compute stationary arrays (i.e. arrays which will not
C change thoughout the simulation).
C
C Depth matrix for wave equation
C     SH=(1+0.5*Tau0*DelT)*<Pi>
C       +DelT^2*0.5*THETA*G*Hdown*<grad(Pi).grad(Pj)>
C Horizontal mass matrix for horizontal diffusion, wave equation, and
C vertical structure
C     SV=<Pi>
C Horizontal derivative SPRSPAK matricies
C     PPx=<Pi*dP/dxj>
C     PPy=<Pi*dP/dyj>
C Laplacian SPRSPAK matrix for horizontal diffusion
C     AGPGP=-<AH*grad(Pi).grad(Pj)>
C     Note that AGPGP is only computed if uniform horizontal 
C     viscosity/diffusivity is used (i.e. if ISMAG=0).  Under these 
C     conditions, AHI and AGPGP will never be computed again.  Thus, AHI
C     is also filled here if ISMAG=0.
C
      CALL STATIONARYQ5(NN,NE,NHBW,X,Y,Hdown,IN,NR,LRVS,G,DelT,Tau0,
     &THETA,ISMAG,AH,AR,DX,DY,NFTR,IQ,JQ,PPx,PPy,AGPGP,SV,SH,DIRICH,CS)
      IF(ISMAG.EQ.0)THEN
         DO I=1,NN
            AHI(I)=AH
         ENDDO
      ENDIF
C
      WRITE(2,'(/A)')'STATIONARY ARRAY CALCULATION SECTION SUCCESSFULLY 
     &COMPLETED'
      WRITE(2,*)'     Sparse factor is',NFTR
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C PREPARATION FOR TIME STEPPING SECTION
C
C Objectives: Call User Subroutines to facilitate initialization and to 
C load the initial atmospheric forcing and point source information.
C Initialize timing parameters.
C
C Zero appropriate arrays and CALL ATMOSQ5 to facilitate initialization
C of local information
C
      DO I=1,NN
         ATMxMID(I)=0.0
         ATMyMID(I)=0.0
         ATEMPmid(I)=0.0
         BTEMPmid(I)=0.0
         PMEmid(I)=0.0
      ENDDO
      CALL ATMOSQ5(KD0,SEC0,0,NN,NNV,XNOD,YNOD,TMID,SMID,
     &ATMxMID,ATMyMID,ATEMPmid,BTEMPmid,PMEmid)
C
C Zero appropriate arrays, CALL POINTSOURCEQ5 to facilitate any
C initialization and to load the initial point source information, then
C compute the initial vertically integrated volumetric source rate
C
      DO J=1,NNV
         DO I=1,NN
            UZSRCmid(I,J)=0.0
            VZSRCmid(I,J)=0.0
            TSRCmid(I,J)=0.0
            SSRCmid(I,J)=0.0
            SRCmid(I,J)=0.0
         ENDDO
      ENDDO
      CALL POINTSOURCEQ5(KD0,SEC0,0,NN,NNV,Zmid,
     &SRCmid,UZSRCmid,VZSRCmid,TSRCmid,SSRCmid)
      CALL VERTSUM(SRCmid,SRCSUMmid,NN,NEV,NNdim)    
C
C CALL BCQ5, OUTPUTQ5, and RHOXYQ4 to facilitate initialization of 
C local information.  ATMxNEW and ATMyNEW are used as dummy arrays here.
C They will be zeroed again before they are used for computation herein.
C
      CALL BCQ5(KD0,SEC0,0,delt,DEGLAT,ATMxNEW,ATMyNEW,0,0)
      CALL OUTPUTQ5(KD0,SEC0,0,DelT,MESHNAME,FILinq,FILicq,
     &Cd,AHI,NN,NNV,NE,XNOD,YNOD,IN,Hdown,AKMIN,
     &ZETAMID,Umid,Vmid,Zmid,UZmid,VZmid,WZmid,Q2mid,Q2Lmid,RHOmid,
     &Tmid,Smid,ENZM,ENZH,ENZQ,ZETAOLD,Uold,Vold,Zold,QP)
      IF(PRESSURE.EQ.'BAROCLINIC')THEN
         CALL RHOXYQ4(0,G,NN,NE,X,Y,IN,PPx,PPy,IQ,JQ,NNV,
     &                 Zmid,RHOmid,NLEV,ZL,BPGx,BPGy,HRBARx,HRBARy,cs)
      ENDIF
C
C Set vertically integrated source terms and precipitation-evaporation
C for the *OLD time level (SEC0-DelT) to *MID values (SEC0) => assume
C that the terms were constant over the previous time step.
C
      DO I=1,NN
         SRCSUMold(I)=SRCSUMmid(I)
         PMEold(I)=PMEmid(I)
      ENDDO
C
C Finally, initialize timing parameters
C
      KDmid=KD0
      SECmid=SEC0
      ITER=0
C
      WRITE(2,'(/A)')'PREPARATION FOR TIME STEPPING SECTION SUCCESSFULLY
     & COMPLETED'
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C TIME STEPPING LOOP
C
  123 CONTINUE            
C
C-----------------------------------------------------------------------
C SET UP FOR THIS TIME STEP
C
C KDmid,SECmid are the timing parameters for the beginning of the
C current time step => time level K
C Set the timing parameters for the time at the end of the current
C time step         => time level K+1
C
      KDnew=KDmid
      SECnew=SECmid+DelT
      ITER=ITER+1
      IF(SECnew.GE.86400.)CALL UP_DATE(KDnew,SECnew)
C
C Zero appropriate arrays and load atmospheric forcing for the end of
C time step (i.e. at time level K+1)
C
      DO I=1,NN
         ATMxNEW(I)=0.0
         ATMyNEW(I)=0.0
         ATEMPnew(I)=0.0
         BTEMPnew(I)=0.0
         PMEnew(I)=0.0
      ENDDO
      CALL ATMOSQ5(KDnew,SECnew,ITER,NN,NNV,XNOD,YNOD,TMID,SMID,
     &ATMxNEW,ATMyNEW,ATEMPnew,BTEMPnew,PMEnew)
C
C Zero appropriate arrays, load point source information, and compute
C the vertical integral of the volumetric source rate for the end of
C the time step (i.e. at time level K+1).  The source information is
C written to the *NEW arrays for efficiency purposes.
C
      DO J=1,NNV
         DO I=1,NN
            UZnew(I,J)=0.0
            VZnew(I,J)=0.0
            Tnew(I,J)=0.0
            Snew(I,J)=0.0
            SRCnew(I,J)=0.0
         ENDDO
      ENDDO
      CALL POINTSOURCEQ5(KDnew,SECnew,ITER,NN,NNV,Zmid,
     &SRCnew,UZnew,VZnew,Tnew,Snew)
      CALL VERTSUM(SRCnew,SRCSUMnew,NN,NEV,NNdim)
C
C Evaluate the following explicitly for the current time step (i.e. at
C time level K):
C
C   Linearized Bottom Stress Coefficients
C   Baroclinic Pressure Gradients
C   Horizontal Eddy Viscosity/Diffusivity
C   Nonlinear Advection and Horizontal Diffusion of Momentum:
C
C     CONxZ(I,J),CONyZ(I,J) are the 3-D advective plus diffusion terms
C       without FL*DV/DZ => CONxZ=(+UDUDX+VDUDY-Ah*DEL^2 U,ETC).
C     The vertical part is implicit in the tri-diagonal velocity matrix.
C     CONx(I),CONy(I) are the vertical integrals of the advective
C       plus diffusion terms for the wave equation
C       => CONx=INT[CONxZ+(FL+Wmesh)*DU/DZ]dz.
C     Convective terms are turned off along all boundaries by INCONV.
C     NONLIN=0 is implemented by setting INCONV=0 for all nodes.
C
      IF(NLBS.EQ.1)THEN
         CALL QUADSTRESS(aK,Cd,aKMIN,NN,UZmid,VZmid,NNdim)
      ENDIF
      IF(PRESSURE.EQ.'BAROCLINIC')THEN
         CALL RHOXYQ4(ITER,G,NN,NE,X,Y,IN,PPx,PPy,IQ,JQ,NNV,
     &                 Zmid,RHOmid,NLEV,ZL,BPGx,BPGy,HRBARx,HRBARy,cs)
c        call printstuff(iter,nn,nnv,x,y,zmid,rhomid,
c    &                   bpgx,bpgy,hrbarx,hrbary)
      ENDIF
C
C Beginning of quoddy4_2.1_P1_main.f modification wrt quoddy4_2.1_main.f
C => new calculation of 2-D horizontal mixing.  Note that dimensioning
C    of addition real nodal array SVHmid is required.
C (DRL/JTCI/CEN 10/14/98) 
C
      IF(ISMAG.EQ.1)THEN
         DO I=1,NN
            HMID(I)=ZETAMID(I)+HDOWN(I)
         ENDDO
         CALL SMAGOR1(AGPGP,Hmid,Umid,Vmid,AR,X,Y,IN,IQ,JQ,NE,SV,NN,
     &               AHI,AH,AHMIN,NFTR,CS)
         DO I=1,NN
            SVHmid(I)=SV(I)*Hmid(I)
         ENDDO
         CALL SPRSINVMLT(AGPGP,IQ,SVHmid,NN)
      ENDIF
C
C End of quoddy4_2.1_P1_main.f modification wrt quoddy4_2.1_main.f.
C
      DO I=1,NN
         CALL SPRSMLTIN2(I,AGPGP,IQ,JQ,
     &                   UZmid,VISCX,VZmid,VISCY,NNV,NNdim)
         CALL SPRSCONV(I,PPx,PPy,IQ,JQ,UZmid,VZmid,
     &                 Zmid,DUZDX,DUZDY,DVZDX,DVZDY,DZDX,DZDY,NNV,NNdim)
         CALL CONVECTION(I,NNV,UZmid,VZmid,WZmid,Zmid,Zold,DelT,
     &           INCONV,DUZDX,DUZDY,DVZDX,DVZDY,DZDX,DZDY,VISCX,VISCY, 
     &           CONxZ,CONyZ,VERTFLX,CONx,CONy,CS)
      ENDDO
C
C-----------------------------------------------------------------------
C SOLVE THE WAVE EQUATION
C
C Objectives:Determine total depth at time level   K+1   (=>ZETANEW)
C            and vertical grid at time level       K+1/2 (=>Znew)
C            using known info at time levels       K-1   (=>*OLD)
C            and                                   K     (=>*MID)
C            with user supplied BC's at various
C            time levels
C
      CALL ELEVATIONQ5(MESHNAME,NN,NE,NHBW,NNV,DelT,G,Tau0,
     &KDnew,SECnew,ITER,DEGLAT,NR,NVN,THETA,XNOD,YNOD,IN,DX,DY,AR,Hdown,
     &ZETAOLD,Uold,Vold,ZETAMID,Umid,Vmid,Zmid,UZmid,VZmid,
     &CONx,CONy,ATMxMID,ATMyMID,aK,HRBARx,HRBARy,
     &SRCSUMold,SRCSUMmid,SRCSUMnew,PMEold,PMEmid,PMEnew,
     &LRVS,NORM,SUMx,SUMy,DSbdry,SUMx1,SUMy1,DSbdry1,QP,SV,SH,DIRICH,
Ccgm-\/
     &ZETANEW,Znew,RAD,CS,COR,TS)
Ccgm-/\
Ccgm IQ, and JQ is needed in conection with the sparse storage of SH
Ccgm+     &ZETANEW,Znew,IQ,JQ,RAD,CS,COR,TS)
Ccgm
C
C-----------------------------------------------------------------------
C SOLVE FOR THE VERTICAL STRUCTURE OF THE 3-D VARIABLES
C
C Objectives: Assemble and solve tridiagonal momentum, Q2, Q2L, T, and S
C             equations at the center of the time step, then convert to
C             values at the end of the time step.
C             => Znew, UZnew, VZnew, Q2new, Q2Lnew, Tnew, and Snew 
C                at time level K+1
C
C Overwrite *MID atmospheric and sourcerate arrays with values at K+1/2
C
      DO I=1,NN
         ATMxMID(I)=0.5*(ATMxMID(I)+ATMxNEW(I))
         ATMyMID(I)=0.5*(ATMyMID(I)+ATMyNEW(I))
         ATEMPmid(I)=0.5*(ATEMPmid(I)+ATEMPnew(I))
         BTEMPmid(I)=0.5*(BTEMPmid(I)+BTEMPnew(I))
      ENDDO
      DO J=1,NEV
         DO I=1,NN
            SRCmid(I,J)=0.5*(SRCmid(I,J)+SRCnew(I,J))
         ENDDO
      ENDDO
C
C Overwrite *NEW dependent variables with source values at K+1/2.
C Increment *SRCmid source term values such that the next time they
C are required, the *SRCmid values will contain the
C values at time level K for that time step (memory efficiency device).
C Note that *NEW(I,NNV) entries are not altered => the neutrality
C flags are not affected by this evolution.
C
      DO J=1,NEV
         DO I=1,NN
            VALMID=UZSRCmid(I,J)
            UZSRCmid(I,J)=UZnew(I,J)
            UZnew(I,J)=0.5*(VALMID+UZnew(I,J))
            VALMID=VZSRCmid(I,J)
            VZSRCmid(I,J)=VZnew(I,J)
            VZnew(I,J)=0.5*(VALMID+VZnew(I,J))
            VALMID=TSRCmid(I,J)
            TSRCmid(I,J)=Tnew(I,J)
            Tnew(I,J)=0.5*(VALMID+Tnew(I,J))
            VALMID=SSRCmid(I,J)
            SSRCmid(I,J)=Snew(I,J)
            Snew(I,J)=0.5*(VALMID+Snew(I,J))
         ENDDO
      ENDDO
C
C Compute vertical structure
C
      CALL VERTICALQ5(NN,NNV,DelT,G,Cd,CLOSURE,MASSVAR,
     &ZLOGBOT,ZLOGTOP,EPSN,KST,NR,NPCN,EPSH,IHHBC,
     &EPSQ,IQADVDIF,IQ2TBC,IQ2BBC,IQ2LTBC,IQ2LBBC,Q2min,Q2Lmin,ELLmin,
     &NORM,LRVS,Hdown,SUMx,SUMy,SUMx1,SUMy1,QP,SV,AGPGP,PPx,PPy,IQ,JQ,
     &INCONV,VERTFLX,aK,BPGx,BPGy,BCFTR,ITER,IUSTARQ2,
     &SRCmid,ATMxMID,ATMyMID,ATEMPmid,BTEMPmid,
     &ZETAMID,Zmid,UZmid,VZmid,Q2mid,Q2Lmid,RHOmid,Tmid,Smid,
     &CONxZ,CONyZ,
     &ZETANEW,Znew,UZnew,VZnew,Q2new,Q2Lnew,Tnew,Snew,
     &EKMMIN,EKHMIN,EKQMIN,ENZM,ENZH,ENZQ,RAD,CS,COR,TS)
C
C-----------------------------------------------------------------------
C INCREMENT TIMING PARAMETERS AND UPDATE/JUGGLE ARRAYS 
C
C Objectives: Increment time such that, for the next time step:
C               *OLD arrays contain information at time level K-1
C               *MID arrays contain information at time level K
C
C *OLD<=*MID:
C
      DO J=1,NNV
         DO I=1,NN
            Zold(I,J)=Zmid(I,J)
         ENDDO
      ENDDO
      DO I=1,NN
         ZETAOLD(I)=ZETAMID(I)              
         Uold(I)=Umid(I)              
         Vold(I)=Vmid(I)              
         PMEold(I)=PMEmid(I)
         SRCSUMold(I)=SRCSUMmid(I)
      ENDDO
C
C *MID<=*NEW:  Thus *MID arrays will contain information for the 
C              current time => (KDmid,SECmid), which is time level K 
C              for the next time step.
C
      KDmid=KDnew
      SECmid=SECnew
      DO J=1,NNV
         DO I=1,NN
            Zmid(I,J) =Znew(I,J)
            UZmid(I,J)=UZnew(I,J)
            VZmid(I,J)=VZnew(I,J)
            SRCmid(I,J)=SRCnew(I,J)
         ENDDO
      ENDDO
      IF(CLOSURE.EQ.'MY25')THEN
         DO J=1,NNV
            DO I=1,NN
               Q2mid(I,J)=Q2new(I,J)
               Q2Lmid(I,J)=Q2Lnew(I,J)
            ENDDO
         ENDDO
      ENDIF
      IF(MASSVAR.EQ.'TWO-PROG')THEN
         DO J=1,NNV
            DO I=1,NN
               Tmid(I,J)=Tnew(I,J)
               Smid(I,J)=Snew(I,J)
            ENDDO
         ENDDO
         CALL EQSTATE2_2D(T0,S0,NN,NNV,TMID,SMID,RHOMID)
      ELSE IF(MASSVAR.EQ.'ONE-PROG')THEN 
         DO J=1,NNV
            DO I=1,NN
               Tmid(I,J)=Tnew(I,J)
            ENDDO
         ENDDO
         CALL EQSTATE1_2D(T0,NN,NNV,TMID,RHOMID)
      ENDIF
      DO I=1,NN
         ZETAMID(I)=ZETANEW(I)       
         ATMxMID(I)=ATMxNEW(I)
         ATMyMID(I)=ATMyNEW(I)
         ATEMPmid(I)=ATEMPnew(I)
         BTEMPmid(I)=BTEMPnew(I)
         PMEmid(I)=PMEnew(I)
         SRCSUMmid(I)=SRCSUMnew(I)
      ENDDO

C Result: Arrays for the *current* time (i.e. the end of the current
C time step) are now stored in as *MID 
C
C-----------------------------------------------------------------------
C COMPUTE VERTICAL VELOCITIES AND VERTICALLY AVERAGED VELOCITIES
C AT CURRENT TIME (*MID)
C
      DO I=1,NN
         CALL SPRSCONV(I,PPx,PPy,IQ,JQ,UZmid,VZmid,
     &                 Zmid,DUZDX,DUZDY,DVZDX,DVZDY,DZDX,DZDY,NNV,NNdim)
         CALL VERTVEL3_3(I,NNV,UZmid,VZmid,WZmid,PMEmid,SRCmid,SV,
     &                Zmid,Zold,DelT,DUZDX,DVZDY,DZDX,DZDY)
      ENDDO 
      CALL VERTAVG(UZmid,Umid,Zmid,NN,NNV,NNdim)
      CALL VERTAVG(VZmid,Vmid,Zmid,NN,NNV,NNdim)
C
C-----------------------------------------------------------------------
C TIME STEP IS COMPLETE, CALL OUTPUT WITH RESULTS AT THE CURRENT TIME
C (*MID)
C
      CALL OUTPUTQ5(KDmid,SECmid,ITER,DelT,MESHNAME,FILinq,FILicq,
     &Cd,AHI,NN,NNV,NE,XNOD,YNOD,IN,Hdown,AKMIN,
     &ZETAMID,Umid,Vmid,Zmid,UZmid,VZmid,WZmid,Q2mid,Q2Lmid,RHOmid,
     &Tmid,Smid,ENZM,ENZH,ENZQ,ZETAOLD,Uold,Vold,Zold,QP)
C
C-----------------------------------------------------------------------
C  ELAPSED TIME CHECK, RETURN FOR ANOTHER TIME STEP OR STOP 
C  (AS APPROPRIATE)
C
      IF(KDmid.LT.KDL)GO TO 123
      IF((KDmid.EQ.KDL).AND.(SECmid.LT.SECL))GO TO 123
C
      WRITE(2,'(/A)')'TIME STEPPING SECTION SUCCESSFULLY COMPLETED'
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C END OF SIMULATION, COMPLETE .echo FILE AND CALL OUTPUTQ5 ONE LAST 
C TIME (WITH ITER SET TO -1)
C
      ITER=-1
      CALL OUTPUTQ5(KDmid,SECmid,ITER,DelT,MESHNAME,FILinq,FILicq,
     &Cd,AHI,NN,NNV,NE,XNOD,YNOD,IN,Hdown,AKMIN,
     &ZETAMID,Umid,Vmid,Zmid,UZmid,VZmid,WZmid,Q2mid,Q2Lmid,RHOmid,
     &Tmid,Smid,ENZM,ENZH,ENZQ,ZETAOLD,Uold,Vold,Zold,QP)
C
      CALL DMY(IDAYmid,IMNTHmid,IYEARmid,KDmid)
      WRITE(2,'(/A)')'SIMULATION SUCCESSFULLY COMPLETED'
      WRITE(2,'(/5X,A35,1X,I2,1X,I2,3X,I4,1X,F18.8)')
     &'Simulation end date and time       ',
     &IDAYmid,IMNTHmid,IYEARmid,SECmid
      CLOSE(2)
C
      STOP
      END        
C***********************************************************************
C***********************************************************************
C End of Main Program for QUODDY5 (RELEASE 1.1)
C End of Main Program for QUODDY5 (RELEASE 1.1)
C***********************************************************************
C***********************************************************************

