c       ************************************************************************
c       *
c       *   eva.f -
c       *
c       *   Copyright (c) 1995
c       *
c       *   ETH Zuerich
c       *   Institut fuer Molekularbiologie und Biophysik
c       *   ETH-Hoenggerberg
c       *   CH-8093 Zuerich
c       *
c       *   All Rights Reserved
c       *
c       *   Date of last modification : 95/09/15
c       *   Pathname of SCCS file     : /export/home3/cb/garant-1.0/src/SCCS/s.eva.f
c       *   SCCS identification       : 1.3
c       *
c       ************************************************************************
c     ==================================================================
      subroutine eva (expr,ires,res,cres,sres,ityp)
c
c      common /timdat/timer(100)
      character quote
      parameter (maxlin=250,maxid=32,maxvar=500,nf=32,nmax=26,
     *           maxnef=20,maxarg=20,quote='''')
      dimension inum(nmax),rnum(nmax),numtyp(nmax),idepth(maxlin),
     *          iarg(maxarg),ic(maxlin)
      complex cres,cnum(nmax)
      logical yes,take,init,opch(0:255),fnch(0:255),idch(0:255),
     *        havbra(nmax),nostr
      character expr*(*),sres*(*),snum(nmax)*(maxlin),e*(maxlin),
     *          ee*(maxlin),fnam(nf)*5,fchar*(nf+1),optra(0:255),
     *          str*(maxlin),efnam(maxnef)*(maxid),rstr*40,op*1
c
      character var*(maxid),vardef*(maxlin)
      common /vardat/nvar
     *       /vdch/var(maxvar),vardef(maxvar)
c
      save init,opch,fnch,idch,optra,ia
c
      data init/.true./
      data fnam/'abs','acos','aimag','aint','anint','asin','atan',
     *          'atan2','char','cmplx','conjg','cos','cosh','dim',
     *          'exp','ichar','index','int','len','log','log10',
     *          'max','min','mod','nint','real','sign','sin','sinh',
     *          'sqrt','tan','tanh'/
      data fchar/'ftzcdsuv"_;qxim!%a$nojkgebhpwlry'/
c
c      t0=timnow()
c     ------------------------------------------------------- initialize
      if (init) then
c        print '(A,3X,A)',(fnam(k),fchar(k:k),k=1,nf)
        ia=ichar('A')-1
        do 1 i=0,255
          op=char(i)
          fnch(i)=index(fchar,op).gt.0 .or. 
     *            (i.ge.128 .and. i.lt.128+maxnef)
          opch(i)=index('^#/|*~-+&0123456789?',op).gt.0
          idch(i)=(op.ge.'0' .and. op.le.'9') .or.
     *            (op.ge.'a' .and. op.le.'z') .or.
     *            (op.ge.'A' .and. op.le.'Z') .or. op.eq.'_'
          if (op.eq.'#') then
            optra(i)='^'
          else if (index('/|~',op).gt.0) then
            optra(i)='*'
          else if (op.eq.'-') then
            optra(i)='+'
          else if (index('12345',op).gt.0) then
            optra(i)='0'
          else if (index('89?',op).gt.0) then
            optra(i)='7'
          else
            optra(i)=op
          end if
    1     continue
        init=.false.
      end if
c     ------------------------------------------- analyze the expression
c      print '(''Expression: "'',A,''"'')',expr(1:lenstr(expr))
      lenexp=lenstr(expr)
      if (lenexp.eq.0) return
c     ----------------------------------- shortcut for positive integers
      do 5 i=1,lenexp
    5   if (expr(i:i).lt.'0' .or. expr(i:i).gt.'9') go to 7
      rstr=expr(1:lenexp)
      read (rstr,'(I20)') ires
      ityp=1
      go to 900
    7 continue
c     ----------------------------------------------- generic expression
      e=' '
      n=0
      nef=0
      l=0
      j=0
      do 100 i=1,lenexp
        if (i.gt.j .and. expr(i:i).gt.' ') then
          ityp=-i
          if (l.ge.maxlin) return
          l=l+1
c          print *,i
c         ................................................... identifier
          if ((expr(i:i).ge.'a' .and. expr(i:i).le.'z') .or.
     *        expr(i:i).eq.'_') then
            do 10 j=i,lenexp-1
   10         if (.not.idch(ichar(expr(j+1:j+1)))) go to 20
c           ................................................... variable
   20       do 25 k=nvar,1,-1
              if (expr(i:j).eq.var(k)) then
                if (n.ge.nmax) return
                n=n+1
                numtyp(n)=-k
                havbra(n)=expr(j+1:j+1).eq.'('
                e(l:l)=char(n+ia)
                go to 100
              end if
   25         continue
c           ......................................... intrinsic function
            i1=1
            i2=nf
   30       if (i1.le.i2) then
              k=(i1+i2)/2
              if (expr(i:j).eq.fnam(k)) then
                e(l:l)=fchar(k:k)
                go to 100
              else if (expr(i:j).lt.fnam(k)) then
                i2=k-1
              else
                i1=k+1
              end if
              go to 30 
            end if
c           .......................................... external function
            if (nef.ge.maxnef) return
            nef=nef+1
            efnam(nef)=expr(i:j)
            if (expr(j+1:j+1).ne.'(') then
              if (n.ge.nmax) return
              n=n+1
              numtyp(n)=-(nvar+nef)
              e(l:l)=char(n+ia)
            else
              e(l:l)=char(127+nef)
            end if
c         ................. logical constant or operation, Fortran style
          else if (expr(i:i).eq.'.' .and. 
     *             expr(i+1:i+1).ge.'a' .and. expr(i+1:i+1).le.'z') then
            if (expr(i:min(i+3,lenexp)).eq.'.eq.') then
              j=i+3
              e(l:l)='0'
            else if (expr(i:min(i+3,lenexp)).eq.'.ne.') then
              j=i+3
              e(l:l)='1'
            else if (expr(i:min(i+3,lenexp)).eq.'.lt.') then
              j=i+3
              e(l:l)='2'
            else if (expr(i:min(i+3,lenexp)).eq.'.le.') then
              j=i+3
              e(l:l)='3'
            else if (expr(i:min(i+3,lenexp)).eq.'.gt.') then
              j=i+3
              e(l:l)='4'
            else if (expr(i:min(i+3,lenexp)).eq.'.ge.') then
              j=i+3
              e(l:l)='5'
            else if (expr(i:min(i+4,lenexp)).eq.'.not.') then
              j=i+4
              e(l:l)='6'
            else if (expr(i:min(i+4,lenexp)).eq.'.and.') then
              j=i+4
              e(l:l)='7'
            else if (expr(i:min(i+3,lenexp)).eq.'.or.') then
              j=i+3
              e(l:l)='8'
            else if (expr(i:min(i+4,lenexp)).eq.'.eqv.') then
              j=i+4
              e(l:l)='9'
            else if (expr(i:min(i+5,lenexp)).eq.'.neqv.') then
              j=i+5
              e(l:l)='?'
            else if (expr(i:min(i+6,lenexp)).eq.'.false.') then
              if (n.ge.nmax) return
              n=n+1
              numtyp(n)=1
              inum(n)=0
              j=i+6
              e(l:l)=char(n+ia)
            else if (expr(i:min(i+5,lenexp)).eq.'.true.') then
              if (n.ge.nmax) return
              n=n+1
              numtyp(n)=1
              inum(n)=1
              j=i+5
              e(l:l)=char(n+ia)
            else
              return
            end if
c         ................................... logical operation, C style
          else if (expr(i:i).eq.'=' .or. expr(i:i).eq.'!' .or.
     *             expr(i:i).eq.'<' .or. expr(i:i).eq.'>' .or.
     *             expr(i:i).eq.'|' .or. expr(i:i).eq.'&') then
            if (i.eq.lenexp) return
            if (expr(i:i+1).eq.'==') then
              j=i+1
              e(l:l)='0'
            else if (expr(i:i+1).eq.'!=') then
              j=i+1
              e(l:l)='1'
            else if (expr(i:i+1).eq.'<=') then
              j=i+1
              e(l:l)='3'
            else if (expr(i:i).eq.'<') then
              j=i
              e(l:l)='2'
            else if (expr(i:i+1).eq.'>=') then
              j=i+1
              e(l:l)='5'
            else if (expr(i:i).eq.'>') then
              j=i
              e(l:l)='4' 
            else if (expr(i:i).eq.'!') then
              j=i
              e(l:l)='6'
            else if (expr(i:i+1).eq.'&&') then
              j=i+1
              e(l:l)='7'
            else if (expr(i:i+1).eq.'||') then
              j=i+1
              e(l:l)='8'
            else
              return
            end if
c         ....................................................... number
          else if (index('0123456789.eEdD',expr(i:i)).gt.0) then
            do 50 j=i,lenexp-1
              if (index(')*/,:=!<>|& ',expr(j+1:j+1)).gt.0 .or.
     *            (index('eEdD',expr(j:j)).eq.0 .and.
     *             index('+-',expr(j+1:j+1)).gt.0)) go to 60
              if (j.le.lenexp-3) then
                if (expr(j+1:j+1).eq.'.' .and. 
     *            expr(j+2:j+2).ge.'a' .and. expr(j+2:j+2).le.'z' .and. 
     *            expr(j+3:j+3).ge.'a' .and. expr(j+3:j+3).le.'z') 
     *            go to 60
              end if
   50         continue
   60       if (n.ge.nmax) return
            n=n+1
            rstr=expr(i:j)
            do 65 k=i,j
   65         if (expr(k:k).lt.'0' .or. expr(k:k).gt.'9') go to 67
            numtyp(n)=1
            read (rstr,'(BN,I20)',err=900) inum(n)
            go to 68
   67       numtyp(n)=2
            read (rstr,'(BN,F40.0)',err=900) rnum(n)
   68       e(l:l)=char(n+ia)
c         ......................................... parenthesis or comma
          else if (expr(i:i).eq.'(' .or. expr(i:i).eq.')' .or.
     *             expr(i:i).eq.',' .or. expr(i:i).eq.':') then
            j=i
            e(l:l)=expr(i:i)
c         ..................................................... operator
          else if (expr(i:i).eq.'+' .or. expr(i:i).eq.'-' .or.
     *             expr(i:i).eq.'*' .or. expr(i:i).eq.'/') then
            if (i.eq.lenexp) return
            if (expr(i+1:i+1).ne.'+' .and. expr(i+1:i+1).ne.'-' .and.
     *          expr(i+1:i+1).ne.'*' .and. expr(i+1:i+1).ne.'/') then
              j=i
              e(l:l)=expr(i:i)
            else if (expr(i:i+1).eq.'++') then
              j=i+1
              e(l:l)='+'
            else if (expr(i:i+1).eq.'+-') then
              j=i+1
              e(l:l)='-'
            else if (expr(i:i+1).eq.'-+') then
              j=i+1
              e(l:l)='-'
            else if (expr(i:i+1).eq.'--') then
              j=i+1
              e(l:l)='+'
            else if (expr(i:i+1).eq.'*+') then
              j=i+1
              e(l:l)='*'
            else if (expr(i:i+1).eq.'*-') then
              j=i+1
              e(l:l)='~'
            else if (expr(i:i+1).eq.'/+') then
              j=i+1
              e(l:l)='/'
            else if (expr(i:i+1).eq.'/-') then
              j=i+1
              e(l:l)='|'
            else if (expr(i:min(i+2,lenexp)).eq.'**+') then
              j=i+2
              e(l:l)='^'
            else if (expr(i:min(i+2,lenexp)).eq.'**-') then
              j=i+2
              e(l:l)='#'
            else if (expr(i:i+1).eq.'**') then
              j=i+1
              e(l:l)='^'
            else if (expr(i:i+1).eq.'//') then
              j=i+1
              e(l:l)='&'
            else
              return
            end if
c         .............................................. string constant
          else if (expr(i:i).eq.quote) then
            if (n.ge.nmax) return
            n=n+1
            numtyp(n)=4
            snum(n)=' '
            i2=0
            take=.false.
            do 70 j=i+1,lenexp-1
              if (expr(j:j).eq.quote) then
                take=.not.take
                if (take) then
                  if (expr(j+1:j+1).ne.quote) go to 80
                  i2=i2+1
                  snum(n)(i2:i2)=quote
                end if
              else
                i2=i2+1
                snum(n)(i2:i2)=expr(j:j)
              end if
   70         continue
            if (expr(j:j).ne.quote) return
   80       inum(n)=i2
            e(l:l)=char(n+ia)
c            print '(A)','"'//snum(n)(1:inum(n))//'"'
          else
            return
          end if
c        print '(''e, i: "'',A,''"'',I5)',expr(1:lenstr(expr)),i
        end if
  100   continue
c      print *,'n = ',n
c      timer(6)=timer(6)+timnow()-t0
c     --------- get values for variables and functions without arguments
      do 150 i=1,n
c        print *,i,numtyp(i)
        if (numtyp(i).lt.-nvar) then
          call intrin (efnam(-numtyp(i)-nvar),i,0,
     *                 inum,rnum,cnum,snum,numtyp)
        else if (numtyp(i).lt.0) then
          if (havbra(i)) then
            snum(i)=vardef(-numtyp(i))
            inum(i)=lenstr(vardef(-numtyp(i)))
            numtyp(i)=4
          else
            call varval (vardef(-numtyp(i)),
     *                   inum(i),rnum(i),cnum(i),snum(i),numtyp(i))
          end if
        end if
  150   if (numtyp(i).lt.1 .or. numtyp(i).gt.4) return
c      print '(2I5,I12,1P,3E15.6)',
c     *      (i,numtyp(i),inum(i),rnum(i),cnum(i),i=1,n)
c     ------------------------------------------ evaluate the expression
      ityp=0
      k=0
  200 if (e.ne.'A') then
c        l=lenstr(e(1:l))
c        print '(A,I5)','"'//e(1:l)//'"',l
        do 210 i=1,l
          if (e(i:i).ge.'A' .and. e(i:i).le.'Z') then
            ee(i:i)='A'
          else if (fnch(ichar(e(i:i)))) then
            ee(i:i)='a'
          else
            ee(i:i)=e(i:i)
          end if
  210     continue
c       ............................................... complex constant
        i=index(ee(1:l),'(A,A)')
        if (i.gt.1) then
          if (ee(i-1:i-1).eq.'a') i=0
        end if
        if (i.gt.0) then  
          i1=ichar(e(i+1:i+1))-ia
          i2=ichar(e(i+3:i+3))-ia
          if (numtyp(i1).ne.numtyp(i2) .or. numtyp(i1).eq.3) return
          if (numtyp(i1).eq.1) then 
            cnum(i1)=cmplx(inum(i1),inum(i2))
          else
            cnum(i1)=cmplx(rnum(i1),rnum(i2))
          end if
          numtyp(i1)=3
          e(i:i)=e(i+1:i+1)
          k=i+5
          go to 600
        end if
c       .................................................. array element
        i=index(ee(1:l),'A(A)')
        if (i.gt.0) then
          is=ichar(e(i:i))-ia
          i1=ichar(e(i+2:i+2))-ia
          iel=inum(i1)
          if (numtyp(is).ne.4 .or. numtyp(i1).ne.1 .or. iel.lt.1) return
          ls=inum(is)
          nc=0
          nostr=.true.
          do 220 j2=1,ls
            if (snum(is)(j2:j2).eq.quote) nostr=.not.nostr
            if (nostr .and. snum(is)(j2:j2).eq.',') then
              nc=nc+1
              ic(nc)=j2
            if (nc.ge.iel) go to 221
            end if
  220       continue
  221     if (iel.gt.nc+1) return
          if (iel.eq.1) then
            j1=1
          else
            j1=ic(iel-1)+1
          end if
          if (iel.le.nc) then
            j2=ic(iel)-1
          else
            j2=ls
          end if
          do 225 j=j1,j2
  225       if (snum(is)(j:j).gt.' ') go to 226
  226     if (j.gt.j2) return
          call varval (snum(is)(j:j2),
     *                 inum(is),rnum(is),cnum(is),snum(is),numtyp(is))
          if (numtyp(is).lt.1 .or. numtyp(is).gt.4) return
          k=i+4
          go to 600
        end if
c       ...................................................... substring
        i=index(ee(1:l),'A(A:A)')
        if (i.eq.0) i=index(ee(1:l),'A(A:)')
        if (i.eq.0) i=index(ee(1:l),'A(:A)')
        if (i.gt.0) then 
          is=ichar(e(i:i))-ia
          ls=inum(is)
          if (numtyp(is).ne.4 .or. ls.lt.3) return
          str=snum(is)(1:ls)
          if (str(1:1).ne.quote .or. str(ls:ls).ne.quote) return
          i2=0
          take=.false.
          do 230 j=2,ls-1
            if (str(j:j).eq.quote) then
              take=.not.take
              if (take) then
                if (str(j+1:j+1).ne.quote) go to 240
                i2=i2+1
                snum(is)(i2:i2)=quote
              end if
            else
              i2=i2+1
              snum(is)(i2:i2)=str(j:j)
            end if
  230       continue
  240     inum(is)=i2
          if (ee(i:i+5).eq.'A(A:A)') then
            i1=ichar(e(i+2:i+2))-ia
            i2=ichar(e(i+4:i+4))-ia
            if (numtyp(i1).ne.1 .or. numtyp(i2).ne.1) return
            is1=inum(i1)
            is2=inum(i2)
            k=i+6
          else if (ee(i:i+4).eq.'A(A:)') then
            i1=ichar(e(i+2:i+2))-ia
            if (numtyp(i1).ne.1) return
            is1=inum(i1)
            is2=inum(is)
            k=i+5
          else if (ee(i:i+4).eq.'A(:A)') then
            i2=ichar(e(i+3:i+3))-ia
            if (numtyp(i2).ne.1) return
            is1=1
            is2=inum(i2)
            k=i+5
          else
            return
          end if
c          if (is1.lt.1 .or. is2.gt.inum(is) .or. is1.ge.is2) return
          str=snum(is)(is1:is2)
          snum(is)=str
          inum(is)=lenstr(snum(is)(1:inum(is)))
          go to 600
        end if
c       ................................ function call with one argument
        i=index(ee(1:l),'a(A)')
        if (i.eq.0) go to 250
        narg=1
        i1=ichar(e(i+2:i+2))-ia
c       .............................................. external function
        if (ichar(e(i:i)).ge.128) then  
          iarg(1)=i1
          call intrin (efnam(ichar(e(i:i))-127),iarg,narg,
     *                 inum,rnum,cnum,snum,numtyp)
          if (numtyp(i1).lt.1 .or. numtyp(i1).gt.4) return
c       ............................................................ int
        else if (e(i:i).eq.'a') then  
          if (numtyp(i1).eq.2) then
            inum(i1)=int(rnum(i1))
          else if (numtyp(i1).eq.3) then
            inum(i1)=int(cnum(i1))
          else
            return
          end if
          numtyp(i1)=1
c       ........................................................... real
        else if (e(i:i).eq.'b') then  
          if (numtyp(i1).eq.1) then
            rnum(i1)=real(inum(i1))
          else if (numtyp(i1).eq.3) then
            rnum(i1)=real(cnum(i1))
          else
            return
          end if
          numtyp(i1)=2
c       ........................................................... aint
        else if (e(i:i).eq.'c') then  
          if (numtyp(i1).ne.2) return
          rnum(i1)=aint(rnum(i1))
c       .......................................................... anint
        else if (e(i:i).eq.'d') then  
          if (numtyp(i1).ne.2) return
          rnum(i1)=anint(rnum(i1))
c       ........................................................... nint
        else if (e(i:i).eq.'e') then  
          if (numtyp(i1).ne.2) return
          inum(i1)=nint(rnum(i1))
          numtyp(i1)=1
c       ............................................................ abs
        else if (e(i:i).eq.'f') then  
          if (numtyp(i1).eq.1) then 
            inum(i1)=abs(inum(i1))
          else if (numtyp(i1).eq.2) then
            rnum(i1)=abs(rnum(i1))
          else if (numtyp(i1).eq.3) then
            numtyp(i1)=2
            rnum(i1)=abs(cnum(i1))
          else
            return
          end if
c       ........................................................... sqrt
        else if (e(i:i).eq.'l') then  
          if (numtyp(i1).eq.2) then
            if (rnum(i1).lt.0.0) return
            rnum(i1)=sqrt(rnum(i1))
          else if (numtyp(i1).eq.3) then
            cnum(i1)=sqrt(cnum(i1))
          else
            return
          end if
c       ............................................................ exp
        else if (e(i:i).eq.'m') then  
          if (numtyp(i1).eq.1) return
          if (numtyp(i1).eq.2) then
            rnum(i1)=exp(rnum(i1))
          else if (numtyp(i1).eq.3) then
            cnum(i1)=exp(cnum(i1))
          else
            return
          end if
c       ............................................................ log
        else if (e(i:i).eq.'n') then  
          if (numtyp(i1).eq.2) then
            if (rnum(i1).le.0.0) return
            rnum(i1)=log(rnum(i1))
          else if (numtyp(i1).eq.3) then
            if (cnum(i1).eq.0.0) return
            cnum(i1)=log(cnum(i1))
          else
            return
          end if
c       .......................................................... log10
        else if (e(i:i).eq.'o') then  
          if (numtyp(i1).ne.2 .or. rnum(i1).le.0.0) return
          rnum(i1)=log10(rnum(i1))
c       ............................................................ sin
        else if (e(i:i).eq.'p') then  
          if (numtyp(i1).eq.2) then
            rnum(i1)=sin(rnum(i1))
          else if (numtyp(i1).eq.3) then
            cnum(i1)=sin(cnum(i1))
          else
            return
          end if
c       ............................................................ cos
        else if (e(i:i).eq.'q') then  
          if (numtyp(i1).eq.2) then
            rnum(i1)=cos(rnum(i1))
          else if (numtyp(i1).eq.3) then
            cnum(i1)=cos(cnum(i1))
          else
            return
          end if
c       ............................................................ tan
        else if (e(i:i).eq.'r') then  
          if (numtyp(i1).ne.2) return
          rnum(i1)=tan(rnum(i1))
c       ........................................................... asin
        else if (e(i:i).eq.'s') then  
          if (numtyp(i1).ne.2 .or. abs(rnum(i1)).gt.1.0) return
          rnum(i1)=asin(rnum(i1))
c       ........................................................... acos
        else if (e(i:i).eq.'t') then  
          if (numtyp(i1).ne.2 .or. abs(rnum(i1)).gt.1.0) return
          rnum(i1)=acos(rnum(i1))
c       ........................................................... atan
        else if (e(i:i).eq.'u') then  
          if (numtyp(i1).ne.2) return
          rnum(i1)=atan(rnum(i1))
c       ........................................................... sinh
        else if (e(i:i).eq.'w') then  
          if (numtyp(i1).ne.2) return
          rnum(i1)=sinh(rnum(i1))
c       ........................................................... cosh
        else if (e(i:i).eq.'x') then  
          if (numtyp(i1).ne.2) return
          rnum(i1)=cosh(rnum(i1))
c       ........................................................... tanh
        else if (e(i:i).eq.'y') then  
          if (numtyp(i1).ne.2) return
          rnum(i1)=tanh(rnum(i1))
c       .......................................................... aimag
        else if (e(i:i).eq.'z') then  
          if (numtyp(i1).ne.3) return
          rnum(i1)=aimag(cnum(i1))
          numtyp(i1)=2
c       .......................................................... ichar
        else if (e(i:i).eq.'!') then  
          if (numtyp(i1).ne.4 .or. inum(i1).ne.1) return
          inum(i1)=ichar(snum(i1)(1:1))
          numtyp(i1)=1
c       ........................................................... char
        else if (e(i:i).eq.'"') then  
          if (numtyp(i1).ne.1) return
          snum(i1)(1:1)=char(inum(i1))
          inum(i1)=1
          numtyp(i1)=4
c       ............................................................ len
        else if (e(i:i).eq.'$') then  
          if (numtyp(i1).ne.4) return
          numtyp(i1)=1
c       .......................................................... conjg
        else if (e(i:i).eq.';') then  
          if (numtyp(i1).ne.3) return
          cnum(i1)=conjg(cnum(i1))
        end if
        go to 300
c       ............................... function call with two arguments
  250   i=index(ee(1:l),'a(A,A)')
        if (i.eq.0) go to 260
        narg=2
        i1=ichar(e(i+2:i+2))-ia
        i2=ichar(e(i+4:i+4))-ia
c       .............................................. external function
        if (ichar(e(i:i)).ge.128) then
          iarg(1)=i1
          iarg(2)=i2
          call intrin (efnam(ichar(e(i:i))-127),iarg,narg,
     *                 inum,rnum,cnum,snum,numtyp)
          if (numtyp(i1).lt.1 .or. numtyp(i1).gt.4) return
c       ............................................................ mod
        else if (e(i:i).eq.'g') then  
          if (numtyp(i1).eq.1 .and. numtyp(i2).eq.1) then 
            if (inum(i2).eq.0) return
            inum(i1)=mod(inum(i1),inum(i2))
          else if (numtyp(i1).eq.2 .and. numtyp(i2).eq.2) then
            if (rnum(i2).eq.0) return
            rnum(i1)=mod(rnum(i1),rnum(i2))
          else
            return
          end if
c       ........................................................... sign
        else if (e(i:i).eq.'h') then  
          if (numtyp(i1).eq.1 .and. numtyp(i2).eq.1) then 
            inum(i1)=sign(inum(i1),inum(i2))
          else if (numtyp(i1).eq.2 .and. numtyp(i2).eq.2) then
            rnum(i1)=sign(rnum(i1),rnum(i2))
          else
            return
          end if
c       ............................................................ dim
        else if (e(i:i).eq.'i') then  
          if (numtyp(i1).eq.1 .and. numtyp(i2).eq.1) then 
            inum(i1)=dim(inum(i1),inum(i2))
          else if (numtyp(i1).eq.2 .and. numtyp(i2).eq.2) then
            rnum(i1)=dim(rnum(i1),rnum(i2))
          else
            return
          end if
c       .......................................................... atan2
        else if (e(i:i).eq.'v') then  
          if (numtyp(i1).ne.2 .or. numtyp(i2).ne.2 .or. 
     *        (rnum(i1).eq.0.0 .and. rnum(i2).eq.0.0)) return
          rnum(i1)=atan2(rnum(i1),rnum(i2))
c       .......................................................... cmplx
        else if (e(i:i).eq.'_') then  
          if (numtyp(i1).eq.1 .and. numtyp(i2).eq.1) then 
            cnum(i1)=cmplx(inum(i1),inum(i2))
          else if (numtyp(i1).eq.2 .and. numtyp(i2).eq.2) then
            cnum(i1)=cmplx(rnum(i1),rnum(i2))
          else
            return
          end if
          numtyp(i1)=3
c       .......................................................... index
        else if (e(i:i).eq.'%') then  
          if (numtyp(i1).ne.4 .or. numtyp(i2).ne.4) return 
          inum(i1)=index(snum(i1)(1:inum(i1)),snum(i2)(1:inum(i2)))
          numtyp(i1)=1
        else
          go to 260
        end if
        go to 300
c       ........................... function with more than two arguments
  260   i=index(ee(1:l),'a(A,A')
        if (i.eq.0) go to 310
        narg=2
        i1=ichar(e(i+2:i+2))-ia
        iarg(1)=i1
        do 270 narg=2,maxarg
          if (ee(i+narg*2-1:i+narg*2).ne.',A') go to 271
  270     iarg(narg)=ichar(e(i+narg*2:i+narg*2))-ia
        return
  271   narg=narg-1
        if (e(i+narg*2+1:i+narg*2+1).ne.')') go to 310
c       ............................................... external function
        if (ichar(e(i:i)).ge.128) then
          call intrin (efnam(ichar(e(i:i))-127),iarg,narg,
     *                 inum,rnum,cnum,snum,numtyp)
          if (numtyp(i1).lt.1 .or. numtyp(i1).gt.4) return
          go to 300
c       ........................................................ max, min
        else if (e(i:i).eq.'j' .or. e(i:i).eq.'k') then
          do 275 jarg=2,narg
            i2=iarg(jarg)
            if (numtyp(i1).eq.1 .and. numtyp(i2).eq.1) then 
              if (e(i:i).eq.'j') then
                inum(i1)=max(inum(i1),inum(i2))
              else
                inum(i1)=min(inum(i1),inum(i2))
              end if
            else if (numtyp(i1).eq.2 .and. numtyp(i2).eq.2) then
              if (e(i:i).eq.'j') then
                rnum(i1)=max(rnum(i1),rnum(i2))
              else
                rnum(i1)=min(rnum(i1),rnum(i2))
              end if
            else
              return
            end if
  275       continue
          go to 300
        end if
        go to 310
c       .................... update expression after function evaluation
  300   e(i:i)=e(i+2:i+2)
        k=i+narg*2+2
        go to 600
c       .................................................... parentheses
  310   i=index(ee(1:l),'(A)')
        if (i.gt.0) then 
          e(i:i)=e(i+1:i+1)
          k=i+3
          go to 600
        end if
c       ..................................... evaluate binary operations   
        ndepth=0
        idepth(1)=0
        do 312 i=1,l
          if (i.gt.1) idepth(i)=idepth(i-1)
          if (e(i:i).eq.'(') then
            idepth(i)=idepth(i)+1
            ndepth=max(ndepth,idepth(i))
          else if (e(i:i).eq.')') then
            idepth(i)=idepth(i)-1
          end if
  312     continue
        do 313 i=1,l
          ich=ichar(ee(i:i))
          if (opch(ich)) then
            if (idepth(i).lt.ndepth) then
              ee(i:i)='@'
            else
              ee(i:i)=optra(ich)
            end if
          end if
  313     continue
        do 320 i=l-2,1,-1
  320     if (ee(i:i+2).eq.'A^A') go to 340
        i=index(ee(1:l),'A*A')
        if (i.gt.0) go to 340 
        i=index(ee(1:l),'A+A')
        if (i.gt.0) go to 335 
        i=index(ee(1:l),'A&A')
        if (i.gt.0) go to 340 
  335   if (ee(i+1:i+1).eq.'+') then
          ib=i
          i=0
          if (ee(1:2).eq.'+A') i=1
          if (i.eq.0) i=index(ee(1:l),'(+A')
          if (i.eq.0) i=index(ee(1:l),',+A')
          if (i.gt.0) go to 400
          i=ib
        end if
  340   if (i.gt.0) then
          op=e(i+1:i+1)
          i1=ichar(e(i:i))-ia
          i2=ichar(e(i+2:i+2))-ia
c          print *,'inum1,inum2,rnum1,rnum2: ',
c     *            inum(i1),inum(i2),rnum(i1),rnum(i2)
c          print *,'before: ',e(i:i+2),ee(i:i+2),str(k:k),inum(i2)
          if (index('#|~',op).gt.0) then
            if (numtyp(i2).eq.1) then
              inum(i2)=-inum(i2)
            else if (numtyp(i2).eq.2) then
              rnum(i2)=-rnum(i2)
            else if (numtyp(i2).eq.3) then
              cnum(i2)=-cnum(i2)
            end if
          end if
c          print *,'after: ',op,inum(i2)
          if (numtyp(i1).eq.1 .and. numtyp(i2).eq.1) then
            if (op.eq.'^' .or. op.eq.'#') then
              inum(i1)=inum(i1)**inum(i2)
            else if (op.eq.'/' .or. op.eq.'|') then
              if (inum(i2).eq.0) return
              inum(i1)=inum(i1)/inum(i2)
            else if (op.eq.'*' .or. op.eq.'~') then
              inum(i1)=inum(i1)*inum(i2)
            else if (op.eq.'-') then
              inum(i1)=inum(i1)-inum(i2)
            else if (op.eq.'+') then
              inum(i1)=inum(i1)+inum(i2)
            else
              return
            end if
          else if (numtyp(i1).le.2 .and. numtyp(i2).le.2) then
            if (numtyp(i1).eq.1) rnum(i1)=inum(i1)
            if (numtyp(i2).eq.1) rnum(i2)=inum(i2)
            numtyp(i1)=2
            if (op.eq.'^' .or. op.eq.'#') then
              if (numtyp(i2).eq.1) then
                rnum(i1)=rnum(i1)**inum(i2)
              else
                rnum(i1)=rnum(i1)**rnum(i2)
              end if
            else if (op.eq.'/' .or. op.eq.'|') then
              if (rnum(i2).eq.0.0) return
              rnum(i1)=rnum(i1)/rnum(i2)
            else if (op.eq.'*' .or. op.eq.'~') then
              rnum(i1)=rnum(i1)*rnum(i2)
            else if (op.eq.'-') then
              rnum(i1)=rnum(i1)-rnum(i2)
            else if (op.eq.'+') then
              rnum(i1)=rnum(i1)+rnum(i2)
            else
              return
            end if
          else if (numtyp(i1).le.3 .and. numtyp(i2).le.3) then
            if (numtyp(i1).eq.1) then
              cnum(i1)=inum(i1)
            else if (numtyp(i1).eq.2) then
              cnum(i1)=rnum(i1)
            end if
            if (numtyp(i2).eq.1) then
              cnum(i2)=inum(i2)
            else if (numtyp(i2).eq.2) then
              cnum(i2)=rnum(i2)
            end if
            numtyp(i1)=3
            if (op.eq.'^' .or. op.eq.'#') then
              if (numtyp(i2).eq.1) then
                cnum(i1)=cnum(i1)**inum(i2)
              else
                cnum(i1)=cnum(i1)**cnum(i2)
              end if
            else if (op.eq.'/' .or. op.eq.'|') then
              if (cnum(i2).eq.0.0) return
              cnum(i1)=cnum(i1)/cnum(i2)
            else if (op.eq.'*' .or. op.eq.'~') then
              cnum(i1)=cnum(i1)*cnum(i2)
            else if (op.eq.'-') then
              cnum(i1)=cnum(i1)-cnum(i2)
            else if (op.eq.'+') then
              cnum(i1)=cnum(i1)+cnum(i2)
            else
              return
            end if
          else
            if (op.eq.'&') then
              if (numtyp(i1).ne.4 .or. numtyp(i2).ne.4) return
              snum(i1)(inum(i1)+1:inum(i1)+inum(i2))=
     *          snum(i2)(1:inum(i2))
              inum(i1)=inum(i1)+inum(i2)
            else
              return
            end if
          end if
          k=i+3
          go to 600
        end if
c       ...................................... evaluate unary operations
        i=0
        if (ee(1:2).eq.'+A') i=1
        if (i.eq.0) i=index(ee(1:l),'(+A')
        if (i.eq.0) i=index(ee(1:l),',+A')
        if (i.eq.0) i=index(ee(1:l),'0+A')
        if (i.eq.0) i=index(ee(1:l),'6+A')
        if (i.eq.0) i=index(ee(1:l),'7+A')
  400   if (i.gt.0) then
          if (ee(i:i).ne.'+') i=i+1
          if (e(i:i).eq.'-') then
            i1=ichar(e(i+1:i+1))-ia
            if (numtyp(i1).eq.1) then
              inum(i1)=-inum(i1)
            else if (numtyp(i1).eq.2) then
              rnum(i1)=-rnum(i1)
            else if (numtyp(i1).eq.3) then
              cnum(i1)=-cnum(i1)
            else
              return
            end if
          end if
          str(i+1:l)=e(i+1:l)
          e(i:l)=str(i+1:l)
          l=l-1
          go to 200
        end if
c       ................................. evaluate relational operations
        i=index(ee(1:l),'A0A')
        if (i.gt.0) then
          op=e(i+1:i+1)
          i1=ichar(e(i:i))-ia
          i2=ichar(e(i+2:i+2))-ia
c          print *,'inum1,inum2,rnum1,rnum2: ',
c     *            inum(i1),inum(i2),rnum(i1),rnum(i2)
          if (numtyp(i1).eq.1 .and.numtyp(i2).eq.1) then
            if (op.eq.'0') then
              yes=inum(i1).eq.inum(i2)
            else if (op.eq.'1') then
              yes=inum(i1).ne.inum(i2)
            else if (op.eq.'2') then
              yes=inum(i1).lt.inum(i2)
            else if (op.eq.'3') then
              yes=inum(i1).le.inum(i2)
            else if (op.eq.'4') then
              yes=inum(i1).gt.inum(i2)
            else if (op.eq.'5') then
              yes=inum(i1).ge.inum(i2)
            end if
          else if (numtyp(i1).le.2 .and.numtyp(i2).le.2) then
            if (numtyp(i1).eq.1) rnum(i1)=inum(i1)
            if (numtyp(i2).eq.1) rnum(i2)=inum(i2)
            if (op.eq.'0') then
              yes=rnum(i1).eq.rnum(i2)
            else if (op.eq.'1') then
              yes=rnum(i1).ne.rnum(i2)
            else if (op.eq.'2') then
              yes=rnum(i1).lt.rnum(i2)
            else if (op.eq.'3') then
              yes=rnum(i1).le.rnum(i2)
            else if (op.eq.'4') then
              yes=rnum(i1).gt.rnum(i2)
            else if (op.eq.'5') then
              yes=rnum(i1).ge.rnum(i2)
            end if
          else if (numtyp(i1).le.3 .and.numtyp(i2).le.3) then
            if (numtyp(i1).eq.1) then
              cnum(i1)=inum(i1)
            else if (numtyp(i1).eq.2) then
              cnum(i1)=rnum(i1)
            end if
            if (numtyp(i2).eq.1) then
              cnum(i2)=inum(i2)
            else if (numtyp(i2).eq.2) then
              cnum(i2)=rnum(i2)
            end if
            if (op.eq.'0') then
              yes=cnum(i1).eq.cnum(i2)
            else if (op.eq.'1') then
              yes=cnum(i1).ne.cnum(i2)
            else
              return
            end if
          else if (numtyp(i1).eq.4 .and.numtyp(i2).eq.4) then
            if (op.eq.'0') then
              yes=snum(i1)(1:inum(i1)).eq.snum(i2)(1:inum(i2))
            else if (op.eq.'1') then
              yes=snum(i1)(1:inum(i1)).ne.snum(i2)(1:inum(i2))
            else if (op.eq.'2') then
              yes=snum(i1)(1:inum(i1)).lt.snum(i2)(1:inum(i2))
            else if (op.eq.'3') then
              yes=snum(i1)(1:inum(i1)).le.snum(i2)(1:inum(i2))
            else if (op.eq.'4') then
              yes=snum(i1)(1:inum(i1)).gt.snum(i2)(1:inum(i2))
            else if (op.eq.'5') then
              yes=snum(i1)(1:inum(i1)).ge.snum(i2)(1:inum(i2))
            end if
          else
            return
          end if
          numtyp(i1)=1
          if (yes) then
            inum(i1)=1
          else
            inum(i1)=0
          end if
          k=i+3
          go to 600
        end if
c       ....................................... evaluate .not. operation
        i=index(ee(1:l),'6A')
        if (i.gt.0) then
          i1=ichar(e(i+1:i+1))-ia
          if (numtyp(i1).ne.1) return
          if (inum(i1).eq.0) then
            inum(i1)=1
          else 
            inum(i1)=0
          end if
          str(i+1:l)=e(i+1:l)
          e(i:l)=str(i+1:l)
          l=l-1
          go to 200
        end if
c       ............................. evaluate binary logical operations
        i=index(ee(1:l),'A7A')
        if (i.gt.0) then
          op=e(i+1:i+1)
          i1=ichar(e(i:i))-ia
          i2=ichar(e(i+2:i+2))-ia
          if (numtyp(i1).ne.1 .or. numtyp(i2).ne.1) return
          if (op.eq.'7') then
            yes=inum(i1).ne.0 .and. inum(i2).ne.0
          else if (op.eq.'8') then
            yes=inum(i1).ne.0 .or. inum(i2).ne.0
          else if (op.eq.'9') then
            yes=inum(i1).ne.0 .eqv. inum(i2).ne.0
          else if (op.eq.'?') then
            yes=inum(i1).ne.0 .neqv. inum(i2).ne.0
          end if
          if (yes) then
            inum(i1)=1
          else
            inum(i1)=0
          end if
          k=i+3
          go to 600
        end if
        return
  600   if (k.gt.l) then
          e(i+1:l)=' '
        else
          str(k:l)=e(k:l)
          e(i+1:l)=str(k:l)
        end if
        l=l-k+i+1
        go to 200
      end if
c     ---------------------------------------------------- return result
      if (numtyp(1).eq.1) then
        ires=inum(1)
      else if (numtyp(1).eq.2) then
        res=rnum(1)
      else if (numtyp(1).eq.3) then
        cres=cnum(1)
      else
        sres=snum(1)(1:inum(1))
        ires=inum(1)
      end if
      ityp=numtyp(1)
  900 continue
c     timer(4)=timer(4)+timnow()-t0
      end
c     ==================================================================
      subroutine varval (vardef,inum,rnum,cnum,snum,numtyp)
c
      character quote
      parameter (maxlin=500,quote='''')
      character*(*) vardef,snum,str*(maxlin)
      complex cnum
      logical take
c
      numtyp=0
      l=lenstr(vardef)
c      print *,i,l,vardef(1:l)
      if (l.eq.0) return
      if (vardef(1:1).eq.quote) then
        if (l.lt.3 .or. vardef(l:l).ne.quote) return
        numtyp=4
        str=' '
        i2=0
        take=.false.
        do 110 j=2,l-1
          if (vardef(j:j).eq.quote) then
            take=.not.take
            if (take) then
              if (vardef(j+1:j+1).ne.quote) go to 120
              i2=i2+1
              str(i2:i2)=quote
            end if
          else
            i2=i2+1
            str(i2:i2)=vardef(j:j)
          end if
  110     continue
  120   snum=str
        inum=i2
      else if (vardef(1:l).eq.'.true.') then
        numtyp=1
        inum=1
      else if (vardef(1:l).eq.'.false.') then
        numtyp=1
        inum=0
      else if (vardef(1:1).eq.'(') then
        if (l.lt.5 .or. vardef(l:l).ne.')') return
        j=index(vardef(1:l),',')
        if (j.lt.3 .or. j.gt.l-1) return
        str(1:40)=vardef(2:j-1)
        read (str(1:40),'(BN,F40.0)',err=900) tr
        str(1:40)=vardef(j+1:l-1)
        read (str(1:40),'(BN,F40.0)',err=900) ti
        numtyp=3
        cnum=cmplx(tr,ti)
c       print *,'"'//vardef(1:l)//'" ',tr,ti,cnum
      else         
        do 130 k=1,l
          if (vardef(k:k).eq.'.' .or. vardef(k:k).eq.'e' .or.
     *        vardef(k:k).eq.'E' .or. vardef(k:k).eq.'d' .or.
     *        vardef(k:k).eq.'D') go to 140
  130     continue
        numtyp=1
        str(1:20)=vardef(1:l)
        read (str(1:20),'(BN,I20)',err=900) inum
        return
  140   numtyp=2
        str(1:40)=vardef(1:l)
        read (str(1:40),'(BN,F40.0)',err=900) rnum
      end if
  900 return
      end
c     ==================================================================
      function ieva(expr,illexp)
      character expr*(*),s
      logical illexp
      complex c
      call eva (expr,ieva,r,c,s,ityp)
      illexp=ityp.ne.1
      end
c     ==================================================================
      function reva(expr,illexp)
      character expr*(*),s
      logical illexp
      complex c
      call eva (expr,i,reva,c,s,ityp)
      if (ityp.eq.1) reva=i
      illexp=ityp.ne.1 .and. ityp.ne.2
      end
c     ==================================================================
      complex function ceva(expr,illexp)
      character expr*(*),s
      logical illexp
      call eva (expr,i,r,ceva,s,ityp)
      if (ityp.eq.1) ceva=i
      if (ityp.eq.2) ceva=r
      illexp=ityp.lt.1 .or. ityp.gt.3
      end
c     ==================================================================
      subroutine seva (expr,res,illexp)
      character*(*) expr,res
      logical illexp
      complex c
      call eva (expr,i,r,c,res,ityp)
      illexp=ityp.ne.4
      end
c     ==================================================================
      logical function leva (expr,illexp)
      character expr*(*),s
      logical illexp
      complex c
      call eva (expr,i,r,c,s,ityp)
      illexp=ityp.ne.1
      leva=i.ne.0
      end