C     *****************************************************************
C     MAIN CODE FOR Water level modeling
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      include 'WLmodel.inc'
c      
      character*5  sNOR
      character*5  sNOL
      character*5  sVAR
      CHARACTER*2048 sTxT
      CHARACTER*32 sSerNAME(nMAXser), sSerCAP(nMAXsyn)
      real*8 rSERall(nMAX1D,2)
      real*8 rSERout(nMAXpairs, nMAXsyn), rDERout(nMAXpairs, nMAXsyn*2)
      real*8 rn(lrna)
      real*8 rRADmin 
      INTEGER nWelINT(nMAXser+1)
      INTEGER nSER, nSYN
      INTEGER ierr
      INTEGER io
      INTEGER i, j
      LOGICAL(4) bYN
c  
      io = 7
c      write(*,'(A)') ' Pathname for WLmodel files........' 
c      read(*,'(A)') sTxT
c      call TrimALL(sTxT)
c      ki = len(trim(sTxT))
c      if(ki.gt.0) bYN = ChangeDIRQQ( sTxT )
c      
   5  write(*,'(A)') '   Input file name for WLmodel........' 
      read(*,'(A)') sTxT
      call TrimALL(sTxT)
      ki = len(trim(sTxT))
c      if(ki.lt.1) goto 5
      if(ki.lt.1) sTxT = 'PEST_WLmodel.INPUT.txt  '
      open(io, file=sTxT , STATUS='OLD', err=5 )
c      
  15  write(*,'(A)') '   OUTPUT file name ........' 
      read(*,'(A)') sTxT
      call TrimALL(sTxT)
      ki = len(trim(sTxT))
c      if(ki.lt.1) goto 15
      if(ki.lt.1) sTxT = 'PEST_WLsimulated.txt '
      open(io+1, file=sTxT)
c      
  25  write(*,'(A)') '   Derivative OUTPUT file name ........' 
      read(*,'(A)') sTxT
      call TrimALL(sTxT)
      ki = len(trim(sTxT))
c      if(ki.lt.1) goto 25
      if(ki.lt.1) sTxT = 'PEST_WLjacobian.txt '
      open(io+2, file=sTxT)
      
      call DelimSET(io, sNOL, sNOR, sVAR )
c  Read time series to be processed.
      call ReadINFO(io,sNOL,sNOR,sVAR, rSERall, nWelINT, 
     +              sSerNAME, nSER, rRADmin )
c  Read each component, modify, store part, and add to synthetic water level
      call MakeSYNTHETIC(io,sNOL,sNOR,sVAR, rSERall, nWelINT, 
     +        sSerNAME, nSER, rRADmin, rSERout, rDERout, sSerCAP, nSYN )
      close(io)
c  Write results
      write(io+1,'(500(1h ,a32))')(sSerCAP(j), j=1,nSYN)  
      write(io+2,'(2I10)')(nSYN-4)*2, nWelINT(2)- 1
      do i = 1, nWelINT(2)- 1
        write(io+1,'(500(1h ,g22.16))')(rSERout(i,j), j=1,nSYN)
        write(io+2,'(500(1h ,g22.16))')(rDERout(i,j), j=1,(nSYN-4)*2)
      enddo
c
      close(io+1)
      close(io+2)
c      
      stop
      END
C ______________________________________________________________________________________________
c
      subroutine ReadINFO(io,sNOL,sNOR,sVAR, rSERall, nWelINT, 
     +                    sSerNAME, nSER, rRADmin )
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)
      include 'WLmodel.inc'
      character*5   sNOR
      character*5   sNOL
      character*5   sVAR
      CHARACTER*2048 sTxT, sTx2
      CHARACTER*32 sSerNAME(nMAXser)
      real*8 rSERall(nMAX1D,2)
      real*8 rRADmin 
      INTEGER nWelINT(nMAXser+1)
      INTEGER nSER
      INTEGER ierr
      INTEGER io
      real*8 rn(lrna)
      INTEGER i, j, iCNT, ki
c      
      rRADmin = rGET(io, sNOL, sNOR, ierr)                               !! Minimum radius for Theis transform
c
      i = 0
      iCNT = 0
      do while(i .lt. nMAXser )
        call ncread(io,sTxT,sNOL,sNOR,ierr)
        if( ierr.lt.1 ) then
          call TrimALL(sTxT)
          call UPCASE(sTxT)
          if( index(sTxT,"SERIESNAME") .gt. 0 ) then
            i = i + 1 
            ki = index(sTxT,':')+1
            sTxT = sTxT(ki:2048)
            ki = index(sTxT,' ')-1
            sSerNAME(i) = sTxT(1:ki)//'                            '
            nWelINT(i) = iCNT + 1
          else
            call qread(rn,2,sTxT,ierr)
            iCNT = iCNT + 1
            rSERall(iCNT,1) = rn(1)
            rSERall(iCNT,2) = rn(2)
            do while( iCNT.gt.0 )
              iCNT = iCNT + 1
              read(io,*,err=5) rSERall(iCNT,1), rSERall(iCNT,2)
            enddo
    5       backspace(io)  
            iCNT = iCNT - 1
        endif
        else
          nWelINT(i+1) = iCNT + 1
          nSER = i
          i = nMAXser
        endif
      end do !  end of series read loop 
c
      return
      end
C ______________________________________________________________________________________________
c
      subroutine MakeSYNTHETIC(io,sNOL,sNOR,sVAR, rSERall, nWelINT, 
     +        sSerNAME, nSER, rRADmin, rSERout,rDERout, sSerCAP, nSYN )
      include 'WLmodel.inc'
      character*5   sNOR
      character*5   sNOL
      character*5   sVAR
      CHARACTER*2048 sTxT
      character*32   sFIND, sTIDE(nMAXser), sFORM
      CHARACTER*32 sSerNAME(nMAXser), sSerCAP(nMAXsyn), sSERmult
      real*8 rSERall(nMAX1D,2), rSERout(nMAXpairs, nMAXsyn)
      real*8 rDERout(nMAXpairs, nMAXsyn*2)
      real*8 rTEMP(nMAXpairs), rDERtemp(nMAXpairs,2), rMULT(nMAXpairs) 
      real*8 rRADmin, rRAD, rTime, rDT, rDlog 
      real*8 rETime,rLat, rLong, rAlt, rAz, rPPB, rPPBdF, rPPBdB  
      INTEGER nWelINT(nMAXser+1)
      INTEGER nSER, nSYN, nMEAS
      INTEGER ierr
      INTEGER io
      real*8 rn(lrna)
      INTEGER i, j, iCNT, kb,ki, n1,n2, iNfac
      sTIDE(3) = 'FIRS                                   '
      sTIDE(4) = 'SECO                                   '
      sTIDE(6) = 'GRAV                                   ' 
      sTIDE(7) = 'TILT                                   '
      sTIDE(8) = 'DRY                                    '
      sTIDE(9) = 'SHEA                                   '
      sFORM = '(F12.3)'
      rDT = 0.01000000000000
      rDlog = exp(rDT)
c
      sSerCAP(1) = 'DATE-TIME'
      sSerCAP(2) = 'MEASURED'
      sSerCAP(3) = 'SYNTHETIC'
      sSerCAP(4) = 'DIFFERENCES'
      nMEAS = nWelINT(2) - nWelINT(1)
      nSYN = 4
      do j = 1, nMEAS
        rSERout(j,1) = rSERall(j,1)
        rSERout(j,2) = rSERall(j,2)
        rSERout(j,3) = 0.00000
        rSERout(j,4) = 0.00000
      enddo  
      i = nSYN
      do while(i .lt. nMAXsyn )
        call ncread(io,sTxT,sNOL,sNOR,ierr)
        if( ierr.lt.1 ) then
          i = i + 1 
          call TrimALL(sTxT)
          call UPCASE(sTxT)
          if( index(sTxT(1:12),'SLOPE+OFFSET') .gt. 0 ) then
            sSerCAP(i) = 'Slope+Offset'
            call qread(rn,2,sTxT,ierr)          !! rn(1) is the slope, rn(2) is the offset       
            do j = 1, nMEAS
              rDERtemp(j,1) = (rSERout(j,1)-rSERout(1,1))
              rDERtemp(j,2) = 1.00000000000000000
              rTEMP(j) = rn(1)*rDERtemp(j,1) + rn(2)
            enddo  
          elseif( index(sTxT(1:4),'STEP') .gt. 0 ) then
            sSerCAP(i) = 'STEP'
            call qread(rn,2,sTxT,ierr)          !! rn(1) is the time, rn(2) is the offset       
            call AddTEXTnNUMB(sSerCAP(i),'-',rn(2),sFORM)
            do j = 1, nMEAS
              rDERtemp(j,1) = 0.
              rDERtemp(j,2) = dIIf(rSERout(j,1).ge.rn(1), 1.D0, 0.D0)
              rTEMP(j) =  rn(2) * rDERtemp(j,2)
            enddo  
          elseif( index(sTxT(1:4),'TIDE') .gt. 0 ) then
            ki = index(sTxT,':')+1
            sTxT = sTxT(ki:2048)
            ki = index(sTxT,' ')-1
            sSerCAP(i) = sTxT(1:ki)//'TIDExxx'
            ki = index(sSerCAP(i),'xxx')
            write(sSerCAP(i)(ki:ki+2),'(i3.3)') i
            sFIND = sTxT(1:4)//'                                 '
            ki = iMATCHs(sFIND, sTIDE,9 ) 
            call qread(rn,6,sTxT,ierr)      !! rn(1)=multiplier, rn(2)=phase, rn(3)=Latitude, rn(4)=Longitude, rn(5)=Altitude, rn(6)=Azimuth
            rLat = rn(3)
            rLong = rn(4)
            rAlt = rn(5)
            rAz = dIIf(ki.gt.3, 0.D0, rn(6) ) 
            do j = 1, nMEAS
              rETime = rSERout(j,1) + rn(2)
              rPPB = AllTide(rETime, rLat, rLong, rAlt, rAz, ki)
              rTEMP(j) = rn(1)*rPPB 
              rDERtemp(j,1) = rPPB
              rPPBdB = AllTide(rETime-rDT, rLat,rLong,rAlt,rAz, ki)
              rPPBdF = AllTide(rETime+rDT, rLat,rLong,rAlt,rAz, ki)
              rDERtemp(j,2) = rn(1)*0.5*(rPPBdF - rPPBdB) / rDT
            enddo  
          elseif( index(sTxT(1:5),'GAMMA') .gt. 0 ) then
            ki = index(sTxT,' ')-1
            sSerCAP(i) = sTxT(1:ki)
            call qread(rn,4,sTxT,ierr)      !! rn(1)=Amplitude, rn(2)=k-term, rn(3)=exponent-n  , rn(4)=time-conversion
            if( ierr.gt.0 ) rn(4) = 1.0000000000000
            if( ierr.gt.1 ) rn(3) = 2.0
            iNfac = rn(3)+0.5
            ki = index(sSerCAP(i),':')+1
            sFIND = sSerCAP(i)(ki:32)//'                   '
            ki = iMATCHs(sFIND, sSerNAME, nSER ) 
            n1 = nWelINT(ki)
            n2 = nWelINT(ki+1) - 1
            ki = index(sSerCAP(i),' ')
            sSerCAP(i)(ki:ki) = '_'
            write(sSerCAP(i)(ki+1:ki+3),'(i3.3)') i
            do j = 1, nMEAS
              rTime = rSERout(j,1)
              rDERtemp(j,1) = dSUMPrecipINC(rSERall, n1,n2, 
     +                                       rn(2), iNfac, rTime)

              rTEMP(j) = rn(1)*rDERtemp(j,1)
              rPPB = rn(1)*dSUMPrecipINC(rSERall, n1,n2, 
     +                                   rn(2)*rDlog, iNfac, rTime)
              rDERtemp(j,2) = (rPPB -rTEMP(j)) / rDT / rn(2)
            enddo  
          elseif( index(sTxT(1:5),'THEIS') .gt. 0 ) then
            ki = index(sTxT,' ')-1
            sSerCAP(i) = sTxT(1:ki)
            call qread(rn,4,sTxT,ierr)      !! rn(1)=Transmissivity, rn(2)=StorageCoefficient, rn(3)=Radius , rn(4)=conversion
            if( ierr.gt.0 ) rn(4) = 1.0000000000000
            ki = index(sSerCAP(i),':')+1
            sFIND = sSerCAP(i)(ki:32)//'                   '
            ki = iMATCHs(sFIND, sSerNAME, nSER ) 
            n1 = nWelINT(ki)
            n2 = nWelINT(ki+1) - 1
            rRAD = rn(3) + rRADmin
            ki = index(sSerCAP(i),' ')
            sSerCAP(i)(ki:ki) = '_'
            write(sSerCAP(i)(ki+1:ki+3),'(i3.3)') i
            do j = 1, nMEAS
              rTime = rSERout(j,1)
              rTEMP(j) = rSUMTheis(rSERall,n1,n2, rn(1),rn(2),
     +                             rRAD, rTime, rn(4))
              rPPB = rSUMTheis(rSERall,n1,n2, rn(1)*rDlog,rn(2),
     +                             rRAD, rTime, rn(4))
              rDERtemp(j,1) = (rPPB -rTEMP(j)) / rDT / rn(1)
              rPPB = rSUMTheis(rSERall,n1,n2, rn(1),rn(2)*rDlog,
     +                             rRAD, rTime, rn(4))
              rDERtemp(j,2) = (rPPB -rTEMP(j)) / rDT / rn(2)
            enddo  
          elseif( index(sTxT(1:3),'AIR') .gt. 0 ) then
            ki = index(sTxT,' ')-1
            sSerCAP(i) = sTxT(1:ki)
            call qread(rn,3,sTxT,ierr)      !! rn(1)=Transmissivity, rn(2)=StorageCoefficient, rn(3)=Unsat length 
            ki = index(sSerCAP(i),':')+1
            sFIND = sSerCAP(i)(ki:32)//'                   '
            ki = iMATCHs(sFIND, sSerNAME, nSER ) 
            n1 = nWelINT(ki)
            n2 = nWelINT(ki+1) - 1
            do j = 1, nMEAS
              rTime = rSERout(j,1)
              rTEMP(j)= rSUMroraHD(rSERall,n1,n2,
     +                             rn(1),rn(2),rn(3),rn(3),rTime)    !!          Tran, Sy,x,a, DateNOW)
              rPPB = rSUMroraHD(rSERall,n1,n2,rn(1)*rDlog,
     +                             rn(2),rn(3),rn(3),rTime)          !!          Tran, Sy,x,a, DateNOW)
              rDERtemp(j,1) = (rPPB -rTEMP(j)) / rDT / rn(1)
              rPPB = rSUMroraHD(rSERall,n1,n2,rn(1),rn(2)*rDlog,
     +                             rn(3),rn(3),rTime)                !!          Tran, Sy,x,a, DateNOW)
              rDERtemp(j,2) = (rPPB -rTEMP(j)) / rDT / rn(2)
            enddo  
          else                                   !!  default to standard series
            kb = index(sTxT,':')+1
            ki = index(sTxT,' ')-1
            sSerCAP(i) = sTxT(kb:ki)
            call qread(rn,3,sTxT,ierr)      !! rn(1)=multiplier, rn(2)=phase, rn(3)=averaging period if present
            ki = iMATCHs(sSerCAP(i), sSerNAME, nSER ) 
            n1 = nWelINT(ki)
            n2 = nWelINT(ki+1) - 1
            call AvgSER(rSERall,n1,n2,rn(3),ierr,nWelINT(nSER+1))
            do j = 1, nMEAS
              rTime = rSERout(j,1) + rn(2)
              n1 = iMATCHrNset(rTime, rSERall, n1, n2)
              rDERtemp(j,1) = rINTERP(rTime, rSERall, n1)
              rTEMP(j) = rn(1)*rDERtemp(j,1)
              rPPB = rINTERP(rTime+rDT, rSERall, n1)
              rDERtemp(j,2) = rn(1)*(rPPB -rDERtemp(j,1)) / rDT
            enddo  
            if(ierr.eq.0) then
              call AddTEXTnNUMB(sSerCAP(i),'-Avg',rn(3),sFORM)
            else  
              sSerCAP(i) = trim(sSerCAP(i))//'-Raw'
            endif  
          endif
c          
c   Create paired multiplier array, fill with 1s if a multiplier array was not specified           
          kb = index(sTxT,'MULTIPLYSERIES:') 
          if(kb .gt. 0) then
            sTxT = sTxT(kb:2048)
            kb = index(sTxT,':')+1
            ki = index(sTxT,' ')-1
            sSERmult = sTxT(kb:ki)
            ki = iMATCHs(sSERmult, sSerNAME, nSER ) 
            n1 = nWelINT(ki)
            n2 = nWelINT(ki+1) - 1
c            call AvgSER(rSERall,n1,n2,rn(3),ierr,nWelINT(nSER+1))  !!  Maybe at a future date
            do j = 1, nMEAS
              rTime = rSERout(j,1)  
              n1 = iMATCHrNset(rTime, rSERall, n1, n2)
              rMULT(j) = rINTERP(rTime, rSERall, n1)
            enddo  
          else
            do j = 1, nMEAS
              rMULT(j) = 1.000000000000000000
            enddo  
          endif  
c   Save component and add to synthetic water level          
          do j = 1, nMEAS
            rSERout(j,i) = rTEMP(j)*rMULT(j)                                     !!  Component of water level model   
            rSERout(j,3) = rSERout(j,3) + rSERout(j,i)                           !!  Add component to water level model   
            rDERout(j,(i-nSYN)*2-1) = rDERtemp(j,1)*rMULT(j)                     !!  Derivative of component with respect to amplitude
            rDERout(j,(i-nSYN)*2)   = rDERtemp(j,2)*rMULT(j)                     !!  Derivative of component with respect to Phase
          enddo  
        else
          nSYN = i
          i = nMAXsyn
        endif
      end do !  end of series read loop 
c
      do j = 1, nMEAS
        rSERout(j,4) = rSERout(j,3) - rSERout(j,2)                      !! Compute residuals  
      enddo  
c
      return
      end
C ______________________________________________________________________________________________
	integer function iMATCHrNset(rFIND, rARRlook, nLOOK1,nLOOK2)
      include 'WLmodel.inc'
      integer nLOOK1, nLOOK2
      integer i
      integer iDIR, iDIRin
      real*8 rFIND
      real*8 rARRlook(nMAX1D,2)
c
      if( rFIND .lt. rARRlook(nLOOK1,1)) then
        iMATCHrNset = nLOOK1
      elseif( rFIND .ge. rARRlook(nLOOK2,1)) then
        iMATCHrNset = nLOOK2
      else  
        i = nLOOK1-1
        do while (i < nLOOK2) 
          i = i+1
          if( rFIND .lt. rARRlook(i,1)) then
            iMATCHrNset = i-1
            i = nLOOK2 
          endif
        end do
      endif
c
      return
	end
C ______________________________________________________________________________________________
c
      subroutine AvgSER(rSERall,n1,n2,rAVG,ierr,n1off)
      include 'WLmodel.inc'
      real*8 rSERall(nMAX1D,2)
      real*8 rTEMP(nMAXpairs), rAVG, rSUM1, rSUM2
      real*8 TimeNOW, rHalf 
      INTEGER nWelINT(nMAXser+1)
      INTEGER nSER
      INTEGER ierr, n, iCNT, n1off, nCNT
      INTEGER i, j,  n1,n2, n1A, n2A
c
      if( ierr.eq.0 ) then
        rHalf = rAVG*0.5 
        nCNT = 0
        n1A = n1
        n2A = n1
        i = n1
        do while( i.lt.n2 )
          rSUM1 = 0.D0
          rSUM2 = 0.D0
          iCNT = 0
          TimeNOW = rSERall(i,1) - rHalf
          n1A = iMATCHrNset(TimeNOW, rSERall, n1A, n2)
          TimeNOW = rSERall(i,1) + rHalf
          n2A = iMATCHrNset(TimeNOW, rSERall, n2A, n2)
          do j = n1A,n2A
            iCNT = iCNT + 1 
            rSUM1 = rSUM1 + rSERall(j,1)
            rSUM2 = rSUM2 + rSERall(j,2)
          enddo
          if( iCNT.gt.0 ) then
            rSERall(n1off+nCNT,1) = rSUM1 / iCNT
            rSERall(n1off+nCNT,2) = rSUM2 / iCNT
          endif
          i = i + 1
          if( n1A.gt.n1 ) nCNT = nCNT + 1  
          if( n2A.ge.n2 ) i = n2  
        end do
        n1 = n1off
        n2 = n1off + nCNT - 1
      endif
c       
      return
      end
C ______________________________________________________________________________________________
c
      real*8 function rINTERP(rTime, rSERall, n1)
      include 'WLmodel.inc'
      real*8 rSERall(nMAX1D,2)
      real*8 rTime, rDX, rDY, rDT 
      INTEGER n1
c
      rINTERP = rSERall(n1,2) 
      rDX = rSERall(n1+1,1) - rSERall(n1,1)
      if( rDX.gt. 1.0E-16 ) then
        rDY = rSERall(n1+1,2) - rSERall(n1,2)
        rDT = rTIME - rSERall(n1,1)
        rINTERP = rINTERP + rDT * rDY / rDX 
      endif       
c       
      return
      end
      