c     File contains utility functions and subroutines for
c     string manipulation and various other operations
      integer function strlen(str)
      implicit none
c
c      finds length to the last non-space character
c
      character*(*)    str
      character*1 c1,c2
      integer length
c
      length = len(str)
      do while (length .gt. 0 .and. (str(length:length).eq.' '
     ;           .or. str(length:length).eq.""))
        length = length - 1
      end do
      strlen = length
      return
      end
c
c     Return the index indicating the start of str1
c     within str2
      integer function indexStr(str1,str2)
      implicit none
      character*(*) str1,str2
c
      integer len1,len2,i,j,k
      character c1,c2
      integer strlen
      external strlen
c
      indexStr=0
      len2=strlen(str2)
      len1=strlen(str1)
      if (len1.gt.len2) then
        return
      elseif (len1.eq.len2) then
        if (LLE(str1,str2).and.LGE(str1,str2)) then
c         two strings are identical return 1
          indexStr=1
          return
        else
          return
        endif
      else
        do i=1,len2
          c2=str2(i:i)
          c1=str1(1:1)
          if (c1.eq.c2) then
            do j=2,len1
              k=i+j
              if (k.gt.len2) then
                goto 12  !stop: the entire str1 did not occur in str2
              endif
              c1=str1(j:j)
              c2=str2(k:k)
              if (c1.ne.c2) then
                goto 10  !break
              endif
            enddo
            indexStr=i
            return
          endif
10      enddo
      endif
12    return
      end
c
c     Compare two strings if they are the same
c     return TRUE
      logical function cmpstr(target,testStr)
      implicit none
      character*(*) testStr,target
c
      character	c1,c2
      logical   done,irange
      integer   strlen,lentar,lentest,i
      external  irange
      external  strlen
c
      lentar = strlen(target)
      lentest= strlen(testStr)
      if (lentar .ne. lentest) then
        cmpstr=.FALSE.
        return
      end if
      i = 1
      done = .FALSE.
      c1 = testStr(i:i)
      c2 = target(i:i)
      cmpstr=.TRUE.
      do while ( .not. done .and. 
     ;           irange(c1) .and. irange(c2) 
     ;           .and. (i .le. lentar) )
       if ( c1.ne.c2 ) then
         cmpstr = .FALSE. 
         done = .TRUE.
        else
         i = i + 1
         c1 = testStr(i:i)
         c2 = target(i:i)
        end if
      end do
      return
      end

      logical function irange(c)
      implicit none
      character    c
c
      integer        intc,izero,icharz
c
      intc = ICHAR(c)
      izero = ICHAR('0')
      icharz = ICHAR('z')
      if (intc .ge. izero .and. intc.le.icharz) then
       irange = .TRUE.
      else
       irange = .FALSE.
      end if
      return
      end

c     copy vector
      subroutine veccpr(istart,iend,vecfrm,vecto,ierr)
      implicit none
      integer           istart,iend,ierr
      double precision  vecfrm(*),vecto(*)
c
      integer i
c
      if (istart .gt. iend) then
        ierr = 1
        return
      end if
      do i=istart,iend
        vecto(i) = vecfrm(i)
      end do
      return
      end
c
c     copy integer vector
      subroutine veccpi(istart,iend,vecfrm,vecto,ierr)
      implicit none
      integer        istart,iend,ierr
      integer        vecfrm(*),vecto(*)
c
      integer i
c
      if ((istart .gt. iend) .or. ierr .gt. 0) then
       ierr = 1
       return
      end if
      do i=istart,iend
       vecto(i) = vecfrm(i)
      end do
      return
      end
c
c     initialization vector to a single value
      subroutine vecstd(istart,iend,vec,value)
      implicit none
      include 'defs.fi'
      integer        istart,iend
      double precision    vec(*),value
c
      integer i
c
      if (iend .lt. istart) then 
        write(nErrLog,*)
     ;  'Attempt to initialize a vector - incorrect bounds'
        stop
      end if
      do i=istart,iend
       vec(i) = value
      end do
      return
      end
c
c     get string tokens from a line of data input in str
      integer function getTok(str,Tok,ldt)
      implicit none
      character*(*)       str
      character*20        Tok(12)
      integer             ldt
c
      integer            strlen,i,j,n
      integer            actLen

c     set array of token strings to all blanks

      do j=1,ldt
        Tok(j)(1:20)='                    '
      end do
c
      i=1
      actLen=strlen(str)
      do while (str(i:i).eq.' ' .or. str(i:i).eq.'\t')
       i=i+1
      end do
      n=0
      do while (i.le.actLen)
       j=1
       n=n+1
c
c       Get next token
c
       do while (str(i:i).ne.' ' .and. str(i:i).ne.'\t' 
     ;   .and. i.le.actLen)
         Tok(n)(j:j)=str(i:i)
         i=i+1
         j=j+1
       end do
c
c       Skip over blanks and tabs
c
       do while ((str(i:i).eq.' ' .or. str(i:i).eq.'\t')
     ;   .and. i.le.actLen)
         i=i+1
        end do
      end do
      getTok=n
      return
      end
c
c     find the layer number of node NodeNum
      integer function layOfNod(NumNP,layerDep,dz,NodeNum,NLayers,
     ;                          ierr)
      implicit none
      double precision layerDep(NLayers),dz
      integer NumNP,nodeNum,NLayers,ierr
c
      double precision depthNod
      integer i
c
      layOfNod=0
      if (NLayers.eq.1) then
        layOfNod=1
        return
      endif
c     special case for bottom node, the layer number must
c     be 1
      if (NodeNum.eq.1) then
        layOfNod=1
        return
      endif
      depthNod=(NumNP-NodeNum)*dz
      do i=NLayers,1,-1
        if (layerDep(i).ge.depthNod) then
          layOfNod=i
          return
        endif
      enddo
c If the nodal depth wasn't found return error code of zero
      return
      end

      double precision function layThick(NLayers,layNum,top,layerDep)
      implicit none
      integer layNum,NLayers
      double precision top,layerDep
      dimension layerDep(NLayers)
      layThick=top - layerDep(layNum)
      return
      end
c
      integer function nodOfZ(NumNP,nodeZ,depth,ierr)
      implicit none
      integer NumNP,ierr
      double precision nodeZ,depth
      dimension nodeZ(NumNP)
c
      integer i
      if (depth.gt.nodeZ(1).or.depth.lt.nodeZ(NumNP)) then
        ierr=1
        nodOfZ=0
      else
        do i=NumNP,1,-1
c         go through nodeZ from z=0 until nodeZ is ge depth
c         return the node number corresponding to depth
          if (nodeZ(i).ge.depth) then
            nodOfZ=i
            goto 100
          endif
        enddo
      endif
100   return
      end
c
c Extract the Par(NTPar) vector from ParD for this material
      subroutine layerPar(ParD,Par,NTPar,iMat)
      implicit none
      integer NTPar,iMat
      double precision ParD(NTPar,*),Par(NTPar)
c
      integer i
c
      do i=1,NTPar
        Par(i) = ParD(i,iMat)
      enddo
      return
      end

c     DateCv is a date conversion utility that accepts a cumulative 
c     number of days starting from Jan. 1 of some year and returns
c     the date in yyyy,mm,dd. 
c     Input: day       - the cumulative number of days
c            styear    - starting year
c            cal       - an array of days in each month of the year
c                        cal(1:12) normal year
c                        cal(13:24) leap year
c
c     Return: yyyy     - current year
c             mm       - month        
c             dd       - day of month
c
      subroutine datecv(day,styear,yyyy,mm,dd)
      implicit none
      include 'defs.fi'
      integer    day,styear,yyyy,mm,dd
c
      logical        isleap
c
      integer        accum,year,calIndx,dInYear,targday,days,
     ;               dayPerM
      integer*4    cal(24)
      data (cal(i),i=1,12) /31,28,31,30,31,30,31,31,30,31,30,31/
      data (cal(i),i=13,24)/31,29,31,30,31,30,31,31,30,31,30,31/

      if (isleap(styear)) then
        dInYear = 366
      else
        dInYear = 365
      end if
      year = styear
      accum = 0
10    continue
      year = year + 1
      accum = accum + dInYear
      if (isleap(year)) then
        dInYear = 366
      else
        dInYear = 365
      end if
      if (accum.le.day) goto 10
      year = year - 1
      yyyy = year
      if (isleap(year)) then
        calIndx = 12
        accum = accum - 366 
      else
        calIndx = 0
        accum = accum - 365
      endif
      mm = 1
      days = cal(1)
      dayPerM = days
      targday = day - accum
      do while (days .lt. targday)
        if (days .gt. 366) then
          write (nErrLog,*) 'Error in calendar specification'
          stop
        end if
        mm = mm + 1
        dayPerM = cal(calIndx+mm)
        days = days + dayPerM 
      enddo
      dayPerM=cal(calIndx+mm)
      days = days - dayPerM
      dd = targday - days
      return
      end

      logical function isleap(year)
      implicit none
      integer  year
c
      if (mod(year,400).eq.0) then
        isleap = .TRUE.
      else if (mod(year,100).eq.0) then
        isleap = .FALSE.
      else if (mod(year,4).eq.0) then
        isleap = .TRUE.
      else
        isleap = .FALSE.
      end if
      return
      end

      integer function elapsDays(yyyy1,mm1,dd1,yyyy2,mm2,dd2,ierr)
      implicit none
      integer yyyy1,mm1,dd1,yyyy2,mm2,dd2
c Two dates supplied as integer values of year, month, day of month
c Returns the number of days between them: date2 - date1
c date2 must be later than date1
      integer    cal(24),cumulDay,mm,year,i,ierr
      logical isleap
      data (cal(i),i=1,12) /31,28,31,30,31,30,31,31,30,31,30,31/
      data (cal(i),i=13,24)/31,29,31,30,31,30,31,31,30,31,30,31/
      ierr=0
      if (yyyy1.eq.yyyy2) then !simple case
        cumulDay=dd2
        do mm=mm2-1,mm1,-1
          if (isleap(yyyy1)) then
            cumulDay=cumulDay + cal(mm+12)
          else
            cumulDay=cumulDay + cal(mm)
          endif
        enddo
        elapsDays=cumulDay-dd1
        return
      else if (yyyy2.gt.yyyy1) then
        cumulDay=dd2
        do mm=mm2-1,1,-1
          if (isleap(yyyy2)) then
            cumulDay=cumulDay + cal(mm+12)
          else
            cumulDay=cumulDay + cal(mm)
          endif
        enddo
        do year=yyyy2-1,yyyy1,-1
          if (isleap(year)) then
            cumulDay=cumulDay + 366
          else
            cumulDay=cumulDay + 365
          endif
        enddo
        do mm=mm1-1,1,-1
          if (isleap(yyyy1)) then
            cumulDay=cumulDay-cal(mm+12)
          else
            cumulDay=cumulDay-cal(mm)
          endif
        enddo
        cumulDay=cumulDay-dd1
      else
        ierr=1
        return
      endif
      elapsDays=cumulDay
      return
      end

      double precision function NTReduc(t,dateInit,TVarCoef)
      implicit none
      include "defs.fi"
      double precision t,TVarCoef(nMonths)
      character*10 dateInit
c
      integer stYear,yyyy,mm,dd,ierr
      double precision zero,one
      data zero,one /0.0D00,1.0D00/
c
c     Does t represent less than one year?
      call dateParts(dateInit,mm,dd,stYear,ierr)
      call datecv(int(t),stYear,yyyy,mm,dd)
      NTReduc = TVarCoef(mm)
      return
      end

      subroutine closer(nunit)
      implicit none
      integer nunit
      close(nunit)
      return
      end

      subroutine advance(nunit,line,feof)
      implicit none
      character*(*) line
      integer nunit,ierr
      logical found,feof
c
      feof=.FALSE.
      found=.FALSE.
      do while (.not.found)
        read(nUnit,'(A)',END=800,ERR=901) line
        if (line(1:1).ne.'*') then
          found=.TRUE.
        endif
      enddo
      return
800   feof=.TRUE.
      return
901   ierr=1
      return
      end

      integer function getDoy(dateVec)
      implicit none
      integer dateVec(3)
c
      integer monthDoy(12),i,year,res
      logical isleap
      external isleap
      data monthDoy /0,31,59,90,120,151,181,212,243,273,304,334/
      res = monthDoy(dateVec(2))+dateVec(3)
      year=dateVec(1)
      if (isleap(year) .and. dateVec(2).gt.2) then
        res = res + 1
      endif
      getDoy=res
      return
      end

      double precision function getMonth(monthStr)
      implicit none
      character*20 monthStr
      logical cmpstr
      external cmpstr
      if (cmpstr('January',monthStr)) then
        getMonth=1.D00
      elseif (cmpstr('February',monthStr)) then
        getMonth=2.D00
      elseif (cmpstr('March',monthStr)) then
        getMonth=3.D00
      elseif (cmpstr('April',monthStr)) then
        getMonth=4.D00
      elseif (cmpstr('May',monthStr)) then
        getMonth=5.D00
      elseif (cmpstr('June',monthStr)) then
        getMonth=6.D00
      elseif (cmpstr('July',monthStr)) then
        getMonth=7.D00
      elseif (cmpstr('August',monthStr)) then
        getMonth=8.D00
      elseif (cmpstr('September',monthStr)) then
        getMonth=9.D00
      elseif (cmpstr('October',monthStr)) then
        getMonth=1.0D01
      elseif (cmpstr('November',monthStr)) then
        getMonth=1.1D01
      elseif (cmpstr('December',monthStr)) then
        getMonth=1.2D01
      else
        getMonth=0.0D00
      endif
      return
      end
