C***********************************************************************
C***********************************************************************
      PROGRAM TRUXTON6
C-----------------------------------------------------------------------
C
C Include file for parameter statements
C
      INCLUDE 'DCMS.DIM'
C
C Arrays for TRUXTON
C
C Arrays/assignments filled during initialization
C
      REAL FREQ(NFDIM)
      INTEGER INFER(NFDIM),NF1,NF2
      COMPLEX AMPRAT(NFDIM),EYE,ONE
      Real relsize(nfdim)
      complex FP(nbcdim,nbcdim)
C
      REAL XOBS(NOBSDIM),YOBS(NOBSDIM),TOBS(NOBSDIM),TOBSgmt0(NOBSDIM),
     &UBAROBS(NOBSDIM),VBAROBS(NOBSDIM),ZETAOBS(NOBSDIM),
     &DOBSMIN(NOBSDIM),DOBSMAX(NOBSDIM),DOBSTOT(NOBSDIM)
      COMPLEX UBAROBSC(NOBSDIM),VBAROBSC(NOBSDIM),ZETAOBSC(NOBSDIM)
C
      CHARACTER*72 CINVTYPE,DATATYPE
C
C Arrays filled during Least Squares Inversion
C
      COMPLEX UBARSAMP(NOBSDIM),VBARSAMP(NOBSDIM),ZETASAMP(NOBSDIM)
      Real Als(NROWDIM,NCOLDIM),Bls(NROWDIM)
      Real Zls(NCOLDIM)
      Real singvals(NCOLDIM)     
      Real Els(NROWDIM),Esamp(NROWDIM)
      REAL UBARLS(NOBSDIM),VBARLS(NOBSDIM),ZETALS(NOBSDIM),
     &UERR(NOBSDIM),VERR(NOBSDIM),ZERR(NOBSDIM)
      COMPLEX UBARLSC(NOBSDIM),VBARLSC(NOBSDIM),ZETALSC(NOBSDIM),
     &UERRC(NOBSDIM),VERRC(NOBSDIM),ZERRC(NOBSDIM)
C
C Arrays filled during Least Squares Postprocessing
C
      Real stderr(nbcdim,nfdim),CV(nbcdim,nfdim)
      complex ZCamp(nbcdim,nfdim) 
      real Cov(NCOLDIM,NCOLDIM)
      complex CCov(NBCdim*NFdim,NBCdim*NFdim)
      real se(NCOLDIM)    
      real Wls(NCOLDIM,NCOLDIM)
C
C Arrays filled during output processing
C
      real Upartial(nobsdim,nfdim),Vpartial(nobsdim,nfdim)
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C End of TRUXTON declarations; proceed with FUNDY declarations
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C Arrays filled by Subroutine INPUT
C
      INTEGER IN(3,NEDIM),INBE(4,NBEDIM),INGEO(2,NBEDIM),
     &NORM(NNDIM),LRVS(NNDIM)
      REAL X(NNDIM),Y(NNDIM),HDOWN(NNDIM),SUMX(NNDIM),SUMY(NNDIM),
     &DSBDRY(NNDIM),ZL(NLEVDIM),COR(NNDIM),CS(NNDIM)
      CHARACTER*72 FILINQ,MESHNAME,COORDS
C
C Sparse matrix derivative arrays filled by Subroutine BUILDPPXPPY
C
      INTEGER IQ(NNDIM),JQ(NFTRDIM*NNDIM)
      REAL PPX(NFTRDIM*NNDIM),PPY(NFTRDIM*NNDIM)
C
C Arrays filled by Main Program
C
      REAL XNOD(NNDIM),YNOD(NNDIM)
C
C Arrays filled by Subroutine VERTGRID
C
      REAL AK(NNDIM)
      REAL Z(NNDIM,NNVDIM),ENZ(NNDIM,NNVDIM)
C
C Arrays filled by Subroutine BAROCLINIC
C
      COMPLEX HRBARX(NEDIM),HRBARY(NEDIM)
      COMPLEX BPGX(NNDIM,NNVDIM),BPGY(NNDIM,NNVDIM)
C
C Arrays filled by Subroutine ATMOS
C
      COMPLEX ATMX(NNDIM),ATMY(NNDIM)
C
C Elevation solution arrays computed in Subroutine ELEVATION
C
      COMPLEX QG(NBWE,NNDIM)
      COMPLEX RHS(NNDIM)
C
C Velocity solution arrays computed by Subroutine VELOCITY
C
      COMPLEX UBAR(NNDIM),VBAR(NNDIM) 
      COMPLEX UVERT(NNDIM,NNVDIM),VVERT(NNDIM,NNVDIM),
     &        WVERT(NNDIM,NNVDIM)
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 Standard physical constants (REAL)
C
      write(*,*) 'set constants'
      G=9.806
      OMEGA=0.7292E-4
      RAD=6.3675E+6
C
C Call SUBROUTINE SETF6 to enable the user to change standard 
C physical constants and mass variable character strings
C
      write(*,*) 'setf6'
      CALL SETF6(G,OMEGA,RAD)
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
C
C Load run information from .inq and associated input files
C
      CALL INPUTF6(FILINQ,OMEGA,RAD,
     &COORDS,COR,CS,DEGLAT,NNV,
     &MESHNAME,NN,NE,X,Y,IN,HDOWN,SUMX,SUMY,DSBDRY,
     &INBE,NBE,INGEO,NGEO,NORM,NVN,KST,LRVS,NR,NPCN,LST,NLEV,ZL)
C
C Compute IQ,JQ,PPX,PPY
C
      CALL BUILDPPXPPYF6(NN,NE,X,Y,IN,IQ,JQ,PPX,PPY,NFTR,CS)
      WRITE(2,'(/A)')'SPARSE DERIVATIVE MATRICIES SET'
      WRITE(2,*)'     Sparse factor is',NFTR
C
C Set input nodal coordinate arrays for communication with/to the
C user subroutines.
C
      DO I=1,NN
         XNOD(I)=X(I)
         YNOD(I)=Y(I)
      ENDDO
      IF(COORDS.EQ.'SPHERICAL')THEN
         DO I=1,NN
            XNOD(I)=X(I)*180.0/ACOS(-1.0)/RAD
            YNOD(I)=Y(I)*180.0/ACOS(-1.0)/RAD
         ENDDO
      ENDIF
C
C Call SUBROUTINE BCF6 to set the tidal frequency being solved for
C (note: UBAR and VBAR are used as dummy variables here) 
C
CEN-NOTFORT6      CALL BCF6(FREQ,UBAR,VBAR,0)
C
C Write settings to echo file 
C
      CALL ECHOSETF6(G,OMEGA,RAD,DEGLAT,COR(1),0.0)
C
      WRITE(2,'(/A)')'INPUT SECTION COMPLETED'
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C VERTICAL STRUCTURE AND ATMOSPHERIC FORCING SECTION
C
C Load vertical structure and mixing coefficients
C
      CALL VERTGRIDF6(NN,NE,XNOD,YNOD,IN,NLEV,ZL,NNV,HDOWN,
     &Z,AK,ENZ)
      WRITE(2,'(/A)')'VERTICAL GRID SPECIFICATION COMPLETED'
C
C Initialize baroclinic forcing then load nonzero values
C
      DO L=1,NE
         HRBARX(L)=CMPLX(0.0,0.0)
         HRBARY(L)=CMPLX(0.0,0.0)
      ENDDO
      DO I=1,NN
         DO J=1,NNV
            BPGX(I,J)=CMPLX(0.0,0.0)
            BPGY(I,J)=CMPLX(0.0,0.0)
         ENDDO
      ENDDO
CEN-NOTFORT6      CALL BAROCLINICF6(NN,NE,X,Y,IN,NLEV,ZL,NNV,Z,G,
CEN-NOTFORT6     &PPX,PPY,IQ,JQ,BPGX,BPGY,HRBARX,HRBARY,CS)
CEN-NOTFORT6      WRITE(2,'(/A)')'BAROCLINIC FORCING SPECIFICATION COMPLETED'
      WRITE(2,'(/A)')'NO BAROCLINIC FORCING ALLOWED'
C
C Initialize atmospheric forcing then load nonzero values
C
      DO I=1,NN
         ATMX(I)=CMPLX(0.0,0.0)
         ATMY(I)=CMPLX(0.0,0.0)
      ENDDO
CEN-NOTFORT6      CALL ATMOSF6(NN,NE,XNOD,YNOD,IN,ATMX,ATMY)
CEN-NOTFORT6      WRITE(2,'(/A)')'ATMOSPHERIC FORCING SPECIFICATION COMPLETED'
      WRITE(2,'(/A)')'NO ATMOSPHERIC FORCING ALLOWED'
      WRITE(2,'(/A)')'END OF FUNDY6 BASIC INPUT/ECHO'
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C TRUXTON SOLUTION SECTION 
C
C ++++++++++++++++
C +INITIALIZATION+
C ++++++++++++++++
C      
C Initialize variables and arrays
C
      EYE=CMPLX(0.0,1.0)
      ONE=CMPLX(1.0,0.0)
C
      NBC=NR
      do i=1,NCOLDIM
         se(i)=0.
         do j=1,NCOLDIM
            Cov(i,j)=0.
         enddo
      enddo
C
      do i=1,NROWDIM
         do j=i,NCOLDIM
            als(i,j)=0.0
         enddo
      enddo      
C
C Check for allowable BC's; if not, stop.
C
      Call BnodeCheck(NBC,NBCdim,NGEO,NVN,KST,NNDIM,NPCN,LST)
C
C Get the frequencies to be used in inversion
C
      Call GetFreqI6_3(NF,NF1,NF2,Freq,relsize,Infer,Amprat,NFDIM)
      write(2,20)
      write(2,'(/A/)')'INVERSE SPECTRUM'
      write(2,*)'Number of constituents         ', NF      
      write(2,*)'Number of primary constituents ',NF1
      write(2,*)'          i, FREQ(i), RELSIZE(i)'
      do i = 1,NF1
         write(2,*)i,FREQ(i),RELSIZE(i)
         if(infer(i).ne.i.or.Amprat(i).ne.ONE) then
            write(2,*)'error in inference data; quitting'     
            write(2,*)'I, FREQ(i), Infer(i), Amprat(i)'
            write(2,*)i,FREQ(i), Infer(i), Amprat(i)
            close(2)
            stop
         endif
      enddo      
c      
      write(2,*)'Number of inferred constituents', NF2
      if(nf2.gt.0) then
         write(2,*)'          I, FREQ(i),           Infer(i), Amprat(i)'
         do i = NF1+1, NF
            write(2,*)i,FREQ(i), Infer(i), Amprat(i)
            if(infer(i).lt.1.or.infer(i).gt.nf1) then	 
               write(2,*)'error in inference data; quitting'
               close(2)
               stop
            endif
         enddo
      endif
C
C Establish Inverse method to be used, and enter/process data 
C
C   WLS: get Rcond, Vrms, w0, w1 
C        Lspec=0 is set to avoid calculations related to fourier method
C        NCOL=2*NBC*NF1
C   SVD: get rcond
C        Lspec=0 is set to avoid calculations related to fourier method
C        NCOL=2*NBC*NF1
C   FOURIER: get Nmode, rcond, and build FP matrix.
C        Nmode is the Fourier half-width
C        Lspec=2*Nmode+1 is the length of the Fourier spectrum
C        NCOL=2*Lspec*NF1
C
      WRITE(*,*)'Enter the type of inverse desired (WLS,SVD,FOURIER)'
      READ(*,'(A)')CINVTYPE
      IF(CINVTYPE.EQ.'WLS')THEN
         write(*,*)'enter Rcond, Vrms, w0, w1'
         read(*,*) Rcond, Vrms, w0, w1
         Lspec=0
         NCOL=2*NBC*NF1
      ELSEIF(CINVTYPE.EQ.'SVD')THEN
         WRITE(*,*)'enter Rcond'
         READ(*,*)Rcond          
         Lspec=0
         NCOL=2*NBC*NF1
      ELSEIF(CINVTYPE.EQ.'FOURIER')THEN
         write(*,*)'Input Nmode, the number of Fourier modes to be ',
     &   'used around the boundary'
         read(*,*)Nmode
         WRITE(*,*)'enter Rcond'
         READ(*,*)Rcond          
         if(2*Nmode.gt.Nbc)then
            write(2,*) '2*Nmode.gt.Nbc; stopping'
            stop
         elseif(Nmode.eq.0)then
            write(2,*) 'Nmode.eq.0; stopping'
            stop
         endif
         Lspec=min(2*Nmode+1,Nbc)    
         NCOL=2*Lspec*NF1
         call BuildFP1(FP, Nbc, Nmode,Lspec)
      ELSE
         WRITE(2,*)'Invalid inverse type selected'
         STOP
      ENDIF         
C
C ++++++++++++++++++++++++++++
C +BUILD LEAST SQUARES SYSTEM+
C ++++++++++++++++++++++++++++
C
      write(2,20)
      write(2,'(/A/)')'BUILDING THE LEAST SQUARES SYSTEM' 
C
C Read the observations to be fit during the inversion
C    NROW=number of rows for the least squares system
C    BLS=the RHS for the WLS system: ALS*ZLS=BLS+ELS
C
      WRITE(*,*)'Enter data input data type: (i.e. UXV,UXZ,UXVC,UXZC)'
      read(*,'(A)')DATATYPE
      IF(DATATYPE.EQ.'UXZC')THEN
         CALL GETOBS_UXZC(NOBSDIM,NOBS,XOBS,YOBS,ZETAOBSC,
     &   IYEAR,NNDIM,NEDIM,IN,XNOD,YNOD,NN,NE)
         NROW=2*NOBS
         DO I=1,NOBS
            BLS(2*I-1)=REAL(ZETAOBSC(I))
            BLS(2*I)  =AIMAG(ZETAOBSC(I))
            UBAROBS(I)=CABS(ZETAOBSC(I))
            VBAROBS(I)=0.0
            TOBS(I)=0.0
            TOBSGMT0(I)=0.0
            DOBSMIN(I)=0.0
            DOBSMAX(I)=0.0
            DOBSTOT(I)=0.0
         ENDDO
         IF(NF.NE.1)THEN
            WRITE(2,*)'NF>1 is not appropriate for complex data'
            STOP
         ENDIF
      ELSEIF(DATATYPE.EQ.'UXZ')THEN
         CALL GETOBS_UXZ(NOBSDIM,NOBS,TOBSGMT0,XOBS,YOBS,ZETAOBS,
     &   TOBS,IYEAR,NNDIM,NEDIM,IN,XNOD,YNOD,NN,NE)
         NROW=NOBS
         DO I=1,NOBS
            BLS(I)=ZETAOBS(I)
            UBAROBS(I)=ZETAOBS(I)
            VBAROBS(I)=0.0
            DOBSMIN(I)=0.0
            DOBSMAX(I)=0.0
            DOBSTOT(I)=0.0
         ENDDO
      ELSEIF(DATATYPE.EQ.'UXVC')THEN
         CALL GETOBS_UXVC(NOBSDIM,NOBS,XOBS,YOBS,IYEAR,
     &   UBAROBSC,VBAROBSC,DOBSMIN,DOBSMAX,DOBSTOT,
     &   NNDIM,NEDIM,IN,XNOD,YNOD,NN,NE)
         NROW=4*NOBS
         DO I=1,NOBS
            BLS(4*I-3)=REAL(UBAROBSC(I))
            BLS(4*I-2)=AIMAG(UBAROBSC(I))
            BLS(4*I-1)=REAL(VBAROBSC(I))
            BLS(4*I  )=AIMAG(VBAROBSC(I))
            UBAROBS(I)=CABS(UBAROBSC(I))
            VBAROBS(I)=CABS(VBAROBSC(I))
            TOBS(I)=0.0
            TOBSGMT0(I)=0.0
         ENDDO
         IF(NF.NE.1)THEN
            WRITE(2,*)'NF>1 is not appropriate for complex data'
            STOP
         ENDIF
      ELSEIF(DATATYPE.EQ.'UXV')THEN
         CALL GETOBS_UXV(NOBSDIM,NOBS,Xobs,Yobs,Tobs,Tobsgmt0,IYEAR,
     &   UBAROBS,VBAROBS,DOBSMIN,DOBSMAX,DOBSTOT,
     &   nndim,nedim,in,XNOD,YNOD,nn,ne)
         NROW=2*NOBS
         DO I=1,NOBS
            BLS(2*I-1)=UBAROBS(I)
            BLS(2*I)  =VBAROBS(I)
         ENDDO
      ELSE
         WRITE(2,*)'Invalid option for data type: ',DATATYPE
         STOP
      ENDIF
C 
C Loop over frequencies to build the least squares matrix: ALS
C
      do ifreq = 1, NF    
         write(2,*)'Frequency ',ifreq,' = ',freq(ifreq)
         FREQI = freq(ifreq)
C      
C   Elevaton matrix for this frequency:
C   Elevation5x brings back LU decompositon of QG; 
C     no RHS info for TYPE I BCs
C     zeroes out QG and RHS before assembling
C
         IFRHS=0
         CALL ELEVATIONF6(NN,NE,X,Y,IN,HDOWN,SUMX,SUMY,DSBDRY,
     &   INGEO,NGEO,NORM,NVN,LRVS,NR,LST,G,COR,FREQI,CS,
     &   NNV,Z,AK,ENZ,BPGX,BPGY,HRBARX,HRBARY,ATMX,ATMY,IFRHS,
     &   QG,NHBW,RHS)
C
C   BEGIN basic loop over Type I BC's for unit repsonses
C   If CINVTYPE.EQ.'FOURIER' then loop is over fourier modes of the bc's
C
         limit=nbc
         if(CINVTYPE.EQ.'FOURIER')limit=Lspec
         do IBC=1,limit
C      
C     Set the RHS's and solve for elevation
C  
            Do I=1,NN
               RHS(I) = cmplx(0.0,0.0)
            Enddo
            if(CINVTYPE.EQ.'FOURIER') then 
               do i=1, nbc
                  rhs(lrvs(i))=fp(i,ibc)
               enddo
            else     
               RHS(LRVS(IBC))= cmplx(1.0,0.0)
            endif
            CALL CSOLVET(2,QG,RHS,NN,NHBW,NNDIM,NBWE)
C
C     Sample solution at observational locations and
C     insert appropriate rows into the Least Squares LHS matrix: ALS
C
            IF(DATATYPE.EQ.'UXZC')THEN
               CALL SAMPLER_Z(IBC,NN,NE,XNOD,YNOD,IN,RHS,
     &         NOBS,XOBS,YOBS,ZETASAMP)
               CALL BUILDALS_UXZC(NOBS,ZETASAMP,
     &         FREQI,INFER,IFREQ,AMPRAT,LSPEC,NBC,IBC,ALS)
            ELSEIF(DATATYPE.EQ.'UXZ')THEN
               CALL SAMPLER_Z(IBC,NN,NE,XNOD,YNOD,IN,RHS,
     &         NOBS,XOBS,YOBS,ZETASAMP)
               CALL BUILDALS_UXZ(NOBS,ZETASAMP,TOBS,
     &         FREQI,INFER,IFREQ,AMPRAT,LSPEC,NBC,IBC,ALS)
            ELSEIF(DATATYPE.EQ.'UXVC')THEN
               IFWVERT=0
               CALL VELOCITYF6(NN,HDOWN,SUMX,SUMY,
     &         INGEO,NGEO,NORM,NVN,KST,LRVS,NR,NPCN,LST,G,COR,FREQI,
     &         NNV,Z,AK,ENZ,BPGX,BPGY,ATMX,ATMY,RHS,IQ,JQ,PPX,PPY,
     &         IFWVERT,UVERT,VVERT,WVERT,UBAR,VBAR)
               CALL SAMPLER_V(IBC,NN,NE,XNOD,YNOD,IN,NNV,Z,UVERT,VVERT,
     &         XOBS,YOBS,NOBS,DOBSMIN,DOBSMAX,DOBSTOT,UBARSAMP,VBARSAMP)
               CALL BUILDALS_UXVC(NOBS,UBARSAMP,VBARSAMP,
     &         FREQI,INFER,IFREQ,AMPRAT,LSPEC,NBC,IBC,ALS)
            ELSEIF(DATATYPE.EQ.'UXV')THEN
               IFWVERT=0
               CALL VELOCITYF6(NN,HDOWN,SUMX,SUMY,
     &         INGEO,NGEO,NORM,NVN,KST,LRVS,NR,NPCN,LST,G,COR,FREQI,
     &         NNV,Z,AK,ENZ,BPGX,BPGY,ATMX,ATMY,RHS,IQ,JQ,PPX,PPY,
     &         IFWVERT,UVERT,VVERT,WVERT,UBAR,VBAR)
               CALL SAMPLER_V(IBC,NN,NE,XNOD,YNOD,IN,NNV,Z,UVERT,VVERT,
     &         XOBS,YOBS,NOBS,DOBSMIN,DOBSMAX,DOBSTOT,UBARSAMP,VBARSAMP)
               CALL BUILDALS_UXV(NOBS,UBARSAMP,VBARSAMP,TOBS,
     &         FREQI,INFER,IFREQ,AMPRAT,LSPEC,NBC,IBC,ALS)
            ELSE
               WRITE(2,*)'Invalid option for data type: ',DATATYPE
               STOP
            ENDIF
         ENDDO
      ENDDO
C
C +++++++++++++++++++++++++++++
C +INVERT LEAST SQUARES SYSTEM+
C +++++++++++++++++++++++++++++
C
      write(2,20)
      write(2,'(/A/)')'INVERTING THE SYSTEM'
C      
C ZLS is the solution vector
C Els is the least square residual error vector
C    Els=Als*Zls-Bls i.e. (model - observation)
C se is the normalized standard error for Zls.  
C    se is SE/(RMS Velocity error)
C Cov is the covariance matrix, normalized by the mean square velocity error.
C    se(i)=sqrt(cov(i,i))
c RMSE is the RMS velocity error.
C
C All arrays are real i.e. they treat x- and y- errors and Re and -Im parts of
C   Zls as separate real variables.
C
      if(CINVTYPE.EQ.'WLS')then
         write(2,*) 'Solution by Weighted Least Squares: '
         write(2,*) 'Rcond, Vrms, w0, w1:', Rcond, Vrms, w0, w1 
         call BuildW6_3(Wls,nobs,nbc,nf1,Vrms,w0,w1,relsize,
     &       LRVS,XNOD,YNOD,inbe,nbe)
         call InvertWLS3_2(nrow,ncol,als,zls,bls,els,
     &                wls,se,Cov,rcond,irank,RMSE,singvals)            
      else 
         write(2,*) 'Solution by SVD with Rcond=',Rcond     
         call Invert3_2(nrow,ncol,als,zls,bls,els,se,Cov,
     &                           rcond,irank,RMSE,singvals)       
      endif
C 
C Recover the model solution at the observational locations
C    (UBARLS,VBARLS) = the modeled solution at the observation points; 
C    (UERR,VERR)     = the modeled error at the observation points;
C
      IF(DATATYPE.EQ.'UXZC')THEN
         DO I=1,NOBS
            ZETALSC(I)=CMPLX(BLS(I*2-1)+ELS(I*2-1),BLS(I*2)+ELS(I*2))
            ZERRC(I)=CMPLX(-ELS(I*2-1),-ELS(I*2))
            UBARLS(I)=CABS(ZETALSC(I))
            VBARLS(I)=0.0
            UERR(I)=CABS(ZERRC(I))
            VERR(I)=0.0
         ENDDO
      ELSEIF(DATATYPE.EQ.'UXZ')THEN
         DO I=1,NOBS
            ZETALS(I)=BLS(I)+ELS(I)
            ZERR(I)=-ELS(I)
            UBARLS(I)=ZETALS(I)
            VBARLS(I)=0.0
            UERR(I)=ZERR(I)
            VERR(I)=0.0
         ENDDO
      ELSEIF(DATATYPE.EQ.'UXVC')THEN
         DO I=1,NOBS
            I4=4*I
            UBARLSC(I)=CMPLX(BLS(I4-3)+ELS(I4-3),BLS(I4-2)+ELS(I4-2))
            VBARLSC(I)=CMPLX(BLS(I4-1)+ELS(I4-1),BLS(I4  )+ELS(I4  ))
            UERRC(I)=CMPLX(-ELS(I4-3),-ELS(I4-2))
            VERRC(I)=CMPLX(-ELS(I4-1),-ELS(I4  ))
            UBARLS(I)=CABS(UBARLSC(I))
            VBARLS(I)=CABS(VBARLSC(I))
            UERR(I)=CABS(UERRC(I))
            VERR(I)=CABS(VERRC(I))
         ENDDO
      ELSEIF(DATATYPE.EQ.'UXV')THEN
         DO I=1,NOBS
            UBARLS(I)=BLS(2*I-1)+ELS(2*I-1)
            VBARLS(I)=BLS(2*I)  +ELS(2*I)
            UERR(I)=-ELS(2*I-1)
            VERR(I)=-ELS(2*I)
         ENDDO
      ELSE
         WRITE(2,*)'Invalid option for data type: ',DATATYPE
         STOP
      ENDIF
C
C ++++++++++++++++++++++++++++++++++
C +CONSTRUCT ARRAY OF ELEVATION BCS+ 
C ++++++++++++++++++++++++++++++++++
C      
C Get the complex covariance matrix CCov from its real parts Cov
C    
      call CmplxCov(CCov, Cov, ncol)
C
C Recover the inferred constitutents: augment Zls, CCov and se;  
C 
      if(nf1.lt.nf)
     &call Inference(NBC,lspec,nf,nf1,amprat,infer,Zls, se,Ccov)
C
C Build ZCamp(i,k)=complex amplitude of BC at node i, freq k.
C    (If Fourier inverse, this involves a transform to physical space.)
C
      call CmplxAmp(Zls,ZCamp,NBC,nf,freq,lspec,FP)
C
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +PROCESS STANDARD ERROR AND COEFFICIENT OF VARIATION RESULTS+ 
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C           
C Interpret and write the std errors and covariances:
C   CCov is the complex covariance of BC's in physical space
c   Cov is the real covariance in method space (i.e. phys. or fourier)
C Convert CCov from Fourier to physical space if needed 
C
      write(2,20)
      if(CINVTYPE.EQ.'FOURIER')then
         call fourier_se1(Ccov,nbc,nf,zls,lspec,fp,
     &                    freq,rcond,irank,ncol,RMSE)
         call CCovK2J(CCov, nbc, nf, lspec, fp)
      else 
         call svd_se1(Ccov,nbc,nf,zls,
     &                    freq,rcond,irank,ncol,RMSE)
      endif
C
C Compute srandard errors and coefficients of variation
C   stderr(i,k)=standard error of ZCamp(i,k) in phys. space
C   CV=se*RMSE/Cabs(Z)=coefficient of variation 
C 
      call stderrCV(stderr, CV, Ccov, ZCamp, RMSE, nbc, nf)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C TRUXTON OUTPUT SECTION 
C
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +RECONSTRUCT AND OUTPUT FUNDY DOMAIN SOLUTION AT EACH FREQUENCY+
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      write(2,20)
      write(2,'(/A/)')'RECONSTRUCTING DOMAIN SOLUTION WITH INVERSE BCS'
C
C Solve the composite system using Zls; call a standard Fundy5 Output
C subroutine for each frequency.
C      
      Do IFREQ = 1, NF
	 write(2,*)'   frequency ',ifreq,' = ',freq(ifreq)
	 FREQI = freq(ifreq)
C
C   Reconstruct QG at this freq.
C
         IFRHS=0
         CALL ELEVATIONF6(NN,NE,X,Y,IN,HDOWN,SUMX,SUMY,DSBDRY,
     &   INGEO,NGEO,NORM,NVN,LRVS,NR,LST,G,COR,FREQI,CS,
     &   NNV,Z,AK,ENZ,BPGX,BPGY,HRBARX,HRBARY,ATMX,ATMY,IFRHS,
     &   QG,NHBW,RHS)
C
C   Assign RHS's at this freq.
C
         Do I=1,NN
            RHS(I) = Cmplx(0.0,0.0)
         Enddo
         Do IBC=1,NBC
            RHS(LRVS(IBC))= ZCamp(ibc,ifreq)
         enddo
C
C   Solve for elevation and velocity at this freq.
C
         CALL CSOLVET(2,QG,RHS,NN,NHBW,NNDIM,NBWE)      
         IFWVERT=1
         CALL VELOCITYF6(NN,HDOWN,SUMX,SUMY,
     &   INGEO,NGEO,NORM,NVN,KST,LRVS,NR,NPCN,LST,G,COR,FREQI,
     &   NNV,Z,AK,ENZ,BPGX,BPGY,ATMX,ATMY,RHS,IQ,JQ,PPX,PPY,IFWVERT,
     &   UVERT,VVERT,WVERT,UBAR,VBAR)
C
C   Sample solution at observation locs.
C
         J=IFREQ
         IF(DATATYPE.EQ.'UXZC')THEN
            CALL SAMPLER_Z(2,NN,NE,XNOD,YNOD,IN,RHS,
     &      NOBS,XOBS,YOBS,ZETASAMP)
            DO I=1,NOBS
               UPARTIAL(I,J)=CABS(ZETASAMP(I))
               VPARTIAL(I,J)=0.0
            ENDDO
         ELSEIF(DATATYPE.EQ.'UXZ')THEN
            CALL SAMPLER_Z(2,NN,NE,XNOD,YNOD,IN,RHS,
     &      NOBS,XOBS,YOBS,ZETASAMP)
            DO I=1,NOBS
               UPARTIAL(I,J)=REAL(ZETASAMP(I)*CEXP(EYE*FREQI*TOBS(I)))
               VPARTIAL(I,J)=0.0
            ENDDO
         ELSEIF(DATATYPE.EQ.'UXVC')THEN
            CALL SAMPLER_V(2,NN,NE,XNOD,YNOD,IN,NNV,Z,UVERT,VVERT,
     &      XOBS,YOBS,NOBS,DOBSMIN,DOBSMAX,DOBSTOT,UBARSAMP,VBARSAMP)
            DO I=1,NOBS
               UPARTIAL(I,J)=CABS(UBARSAMP(I))
               VPARTIAL(I,J)=CABS(VBARSAMP(I))
            ENDDO
         ELSEIF(DATATYPE.EQ.'UXV')THEN
            CALL SAMPLER_V(2,NN,NE,XNOD,YNOD,IN,NNV,Z,UVERT,VVERT,
     &      XOBS,YOBS,NOBS,DOBSMIN,DOBSMAX,DOBSTOT,UBARSAMP,VBARSAMP)
            DO I=1,NOBS
               UPARTIAL(I,J)=REAL(UBARSAMP(I)*CEXP(EYE*FREQI*TOBS(I)))
               VPARTIAL(I,J)=REAL(VBARSAMP(I)*CEXP(EYE*FREQI*TOBS(I)))
            ENDDO
         ELSE
            WRITE(2,*)'Invalid option for data type: ',DATATYPE
            STOP
         ENDIF
C
C   Write the results at this freq.  
C
         CALL OUTPUTF6(MESHNAME,NN,NE,XNOD,YNOD,IN,NNV,Z,
     &   RHS,UBAR,VBAR,UVERT,VVERT,WVERT,FREQI)
      ENDDO
C
C +++++++++++++++++++++++
C +OUTPUT TRUXTON ARRAYS+
C +++++++++++++++++++++++
C
C Call OUTPUTT6
C      
      CALL OUTPUTT6(MESHNAME,NOBS,XOBS,YOBS,TOBS,TOBSGMT0,IYEAR,
     &UBAROBS,VBAROBS,DOBSMIN,DOBSMAX,DOBSTOT,CINVTYPE,DATATYPE,
     &NF,NF1,NF2,FREQ,INFER,AMPRAT, 
     &NBC,LRVS,ZCAMP,STDERR,CV,CCOV,RMSE,RCOND,IRANK,NCOL,SINGVALS,
     &UBARLS,VBARLS,UERR,VERR,UPARTIAL,VPARTIAL)

CALF 02/04/02 Call TruxOut6 to output uxv output 
      write(*,*)"RHERHERHERHERHE ====== TruxOut6==START==="

      CALL TruxOut6(NOBS,Xobs,Yobs,Tobs,Tobsgmt0,IYEAR,
     & Ubarobs,Vbarobs,dobsmax,dobsmin,dobstot,
     & Nf, Nf1, Nf2, freq, infer, amprat,
     & NBC,ZCamp,LRVS,CCov,
     & Ubarls, Vbarls, Upartial, Vpartial,
     & Uerr, Verr, RMSE,rcond,irank,NCOL,MESHNAME)
      write(*,*)"RHERHERHERHERHE ====== TruxOut6==STOP==="
CALF 02/04/02 end

C
C Write successful completion message to ECHO file
C
      write(2,20)
      write(2,'(/A)')'TRUXTON RUN SUCCESSFULLY COMPLETED'
      write(2,20)
      CLOSE(2)
C
C Format statements
C   
   20 format(/72('*')/72('*'))   
 110  FORMAT('Set in Subroutine VERTGRID:',
     &/5X,'The number of nodes in the vertical = ',I6)
C
C End of Main Program
C
      STOP
      END
C
C***********************************************************************
C***********************************************************************
C END OF MAIN ROUTINE
C END OF MAIN ROUTINE
C***********************************************************************
C***********************************************************************
