C***********************************************************************
C***********************************************************************
C Beginning of Main Program for QUODDY4_PSEUDO (RELEASE 2.1)
C Beginning of Main Program for QUODDY4_PSEUDO (RELEASE 2.1)
C***********************************************************************
C***********************************************************************
      PROGRAM QUODDY4_PSEUDO
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 'QUODDY.DIM'
C
C 2-D array declaration
C
      INTEGER IN(3,NEDIM),INBE(4,NBEDIM),INGEO(2,NBEDIM),
     &LRVS(NNDIM),NORM(NNDIM)
      REAL HOLD(NNDIM),UOLD(NNDIM),VOLD(NNDIM),
     &HOLDS(NNDIM),UOLDS(NNDIM),VOLDS(NNDIM),
     &HMID(NNDIM),UMID(NNDIM),VMID(NNDIM),
     &HNEW(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),
     &SUMX1(NNDIM),SUMY1(NNDIM),DSBDRY1(NNDIM),QP(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)
      REAL ZOLDS(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
      REAL SH(NB1DIM,NNDIM),DIRICH(NBWE,NRDIM)
C<<<
      REAL ZERO
      DATA ZERO/0.0/
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,MESHICQ,FILINQ,CITER,INTERNAL,
     &ICQ4LIST,HEADER
      CHARACTER*72 ICQ4FILE(NNDIM)
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 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 INPUTQ4(FILINQ,
     &FILICQ,IDD0,IMM0,IYEAR0,SEC0,IDDL,IMML,IYEARL,SECL,
     &DEG,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)
C
C Standard physical constants (REAL)
C
      Cd=0.005
C
C Initialize data arrays not passed to OUTPUTQ4 in current version
C
      DO I=1,NN
         QP(I)=0.0
         AHI(I)=0.0
         DO J=1,NNV
            ENZM(I,J)=0.0
            ENZH(I,J)=0.0
            ENZQ(I,J)=0.0
         ENDDO
      ENDDO
C
C Load .icq4 file names
C
      WRITE(*,'(A)')'Enter the name of file containing the .icq4 file na
     &mes'
      READ(*,'(A)')ICQ4LIST
      WRITE(*,'(A)')'.icq4 files to be read:'
      OPEN(57,FILE=ICQ4LIST,STATUS='OLD')
      DO IFILE=1,NNDIM
         READ(57,'(A)',ERR=5)ICQ4FILE(IFILE)
         WRITE(*,*)ICQ4FILE(IFILE)
         NFILES=IFILE
      ENDDO
 5    CLOSE(57)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C TIME STEPPING LOOP
C
C Loop over icq4 files
C
      DO ITER=0,NFILES-1
C
C   Read the .icq4 file
C
         FILICQ=ICQ4FILE(ITER+1)
         CALL ICQ4READ(FILICQ,MESHICQ,NN,NNV,
     &   DELTS,IDD,IMM,IYEAR,SEC,HMID,UMID,VMID,HOLDS,UOLDS,VOLDS,
     &   ZMID,ZOLDS,UZMID,VZMID,WZMID,Q2MID,Q2LMID,TMID,SMID,NNDIM,
     &   NNVDIM)
C
C   Determine Gregorian timing parameters
C
         CALL GDAY(IDD,IMM,IYEAR,KDMID)
         SECMID=SEC
C
C   If this is the first icq4 file, open the second one and read the 
C   time in order to set the timestep and the old timing parameters
C
         IF(ITER.EQ.0)THEN
            OPEN(57,FILE=ICQ4FILE(2),STATUS='OLD')
            DO I=1,4
               READ(57,'(A)')HEADER
            ENDDO
            READ(57,*)NNOLD,NNVOLD
            READ(57,*)IDD1,IMM1,IYEAR1,SEC1
            CALL GDAY(IDD1,IMM1,IYEAR1,KD1)
            DELTOLD=(KD1-KDMID)*24.0*3600+SEC1-SECMID
            KDOLD=KDMID
            SECOLD=SECMID-DELTOLD
            CLOSE(57)
         ENDIF
C
C   Ensure the time-step between this and the most recent .icq4 file 
C   is the same as previous ones
C
         DELT=(KDMID-KDOLD)*24.0*3600+SECMID-SECOLD
         IF(NINT(DELT).NE.NINT(DELTOLD))THEN
            WRITE(*,'(2A)')'ERROR IN .icq4 FILE: ',FILICQ
            WRITE(*,'(2A)')'  TIME-STEP DOES NOT MATCH THAT OF',
     &      ' PREVIOUS FILES'
         STOP
         ENDIF
         CALL CROPSTRING(FILICQ,is1,ib1,ie1)
         WRITE(*,'(A,5I7,F7.0)')FILICQ(is1:ie1),
     &   ITER,IDD,IMM,IYEAR,KDMID,SECMID
C
C   Check that grid matches previous .icq4 file
C
         CALL CROPSTRING(MESHNAME,is1,ib1,ie1)
         CALL CROPSTRING(MESHICQ, is2,ib2,ie2)
         IFLAG=0
         IF(MESHNAME(ib1:ie1).NE.MESHICQ(ib2:ie2))THEN
            WRITE(*,'(2A)')'ERROR IN .icq4 FILE: ',FILICQ
            WRITE(*,'(2A)')'  MESH IDENTIFIED IN FILE DOES NOT MATCH',
     &      ' THAT LISTED IN .inq FILE'
            IFLAG=1
         ENDIF
         IF(NN.NE.NNOLD)THEN
            WRITE(*,'(2A)')'ERROR IN .icq4 FILE: ',FILICQ
            WRITE(*,'(2A)')'  NUMBER OF HORIZONTAL NODES DOES NOT',
     &      ' MATCH THAT LISTED IN .inq FILE'
            IFLAG=1
         ENDIF
         IF(NNV.NE.NNVOLD)THEN
            WRITE(*,'(2A)')'ERROR IN .icq4 FILE: ',FILICQ
            WRITE(*,'(2A)')'  NUMBER OF VERTICAL NODES DOES NOT',
     &      ' MATCH THAT LISTED IN .inq FILE'
            IFLAG=1
         ENDIF
         IF(IFLAG.EQ.1)STOP
C
C   Calculate necessary additions for OUTPUTQ4
C
C     All iterations
C
         DO I=1,NN
            HMID(I)=ZMID(I,NNV)-ZMID(I,1)
         ENDDO
         CALL EQSTATE2(10.0,32.0,RHOREF)
         DO I=1,NN
            DO J=1,NNV
               CALL EQSTATE2(TMID(I,J),SMID(I,J),RRMID)
               RHOMID(I,J)=(RRMID-RHOREF)/(1.E3+RHOREF)
            ENDDO
         ENDDO
C
C     Only for ITER=0
C
         IF(ITER.EQ.0)THEN
            DO I=1,NN
               HOLD(I)=HMID(I)
               UOLD(I)=UMID(I)
               VOLD(I)=VMID(I)
               DO J=1,NNV
                  ZOLD(I,J)=ZMID(I,J)
               ENDDO
            ENDDO
         ENDIF
C
C   Call OUTPUTQ4
C
         CALL OUTPUTQ4(KDMID,SECMID,ITER,DELT,MESHNAME,FILINQ,FILICQ,
     &   Cd,AHI,NN,NNV,NE,X,Y,IN,HDOWN,
     &   HMID,UMID,VMID,ZMID,UZMID,VZMID,WZMID,Q2MID,Q2LMID,RHOMID,
     &   TMID,SMID,ENZM,ENZH,ENZQ,HOLD,UOLD,VOLD,ZOLD,QP)
C
C   Index arrays
C
         KDOLD=KDMID
         SECOLD=SECMID
         NNOLD=NN
         NNVOLD=NNV
         DO I=1,NN
            HOLD(I)=HMID(I)
            UOLD(I)=UMID(I)
            VOLD(I)=VMID(I)
            DO J=1,NNV
               ZOLD(I,J)=ZMID(I,J)
            ENDDO
         ENDDO
      ENDDO
 10   STOP
      END        
C***********************************************************************
C***********************************************************************
C End of Main Program for QUODDY4_PSEUDO (RELEASE 2.1)
C End of Main Program for QUODDY4_PSEUDO (RELEASE 2.1)
C***********************************************************************
C***********************************************************************
