      subroutine h2oex(pot,rij)
C
C     H2O - FITPOTENTIAL - EXCITED STATE (1B1)
C     CALCULATED BY STAEMMLER + PALMA;  FIT BY SCHINKE
C
C     ABSTAENDE IN ANGSTROEM,  ENERGIE IN EV
C
      implicit real*8 (a-h,o-z)
C
      dimension r(3),re(3),alpha(3),c(100),phi(100)
      dimension rij(3)
      real*8 angtoau, autoev
      data angtoau/0.529177249d0/, autoev/27.2113961d0/
C
C
      data doh,betoh,roh/4.621d0 , 2.294d0 , 0.971d0/
C
      rmorse(x)=doh*(exp(-betoh*(x-roh))-1.0d0)**2 - doh
C
      xsi(x,y)=(1.0d0-exp(-gamma*x**2)) * (1.0d0-exp(-gamma*y**2))
C
C     EINLESEN DER DATEN
C
      data norder/ 50/
      data re/      1.0000d0,      1.0000d0,      1.0000d0/
      data alpha/      1.2500d0,      0.7500d0,      1.2500d0/
      data gamma/      0.2500d0/
      data c(  1)/ -0.93562626d+00/
      data c(  2)/ -0.33992392d+01/
      data c(  3)/ -0.37923375d+01/
      data c(  4)/  0.35229926d+01/
      data c(  5)/  0.20508566d+01/
      data c(  6)/  0.33290770d+00/
      data c(  7)/  0.64682997d+00/
      data c(  8)/ -0.51447712d+02/
      data c(  9)/ -0.49046646d+01/
      data c( 10)/  0.24435891d+02/
      data c( 11)/  0.76337272d+00/
      data c( 12)/  0.10582765d+02/
      data c( 13)/  0.73669551d+01/
      data c( 14)/  0.64329848d+02/
      data c( 15)/  0.89923222d+01/
      data c( 16)/ -0.89690335d+02/
      data c( 17)/  0.64529244d+02/
      data c( 18)/ -0.53894246d+02/
      data c( 19)/ -0.59093483d+00/
      data c( 20)/ -0.12138697d+02/
      data c( 21)/  0.65697698d+01/
      data c( 22)/  0.59787268d+02/
      data c( 23)/ -0.29842680d+02/
      data c( 24)/  0.30342123d+02/
      data c( 25)/ -0.29963985d+02/
      data c( 26)/  0.24344415d+02/
      data c( 27)/ -0.62823352d+02/
      data c( 28)/  0.18190455d+03/
      data c( 29)/  0.10358976d+02/
      data c( 30)/ -0.39393237d+02/
      data c( 31)/  0.39698652d+01/
      data c( 32)/  0.12165069d+02/
      data c( 33)/  0.51422343d+01/
      data c( 34)/ -0.47491777d+01/
      data c( 35)/ -0.15323205d+01/
      data c( 36)/  0.92191236d+01/
      data c( 37)/ -0.74335345d+01/
      data c( 38)/  0.15738989d+03/
      data c( 39)/  0.33299372d+01/
      data c( 40)/ -0.35652609d+01/
      data c( 41)/ -0.15416913d+03/
      data c( 42)/ -0.17472317d+01/
      data c( 43)/  0.52101774d+02/
      data c( 44)/  0.16650884d+03/
      data c( 45)/ -0.23987576d+02/
      data c( 46)/ -0.69990914d+02/
      data c( 47)/  0.19884456d+02/
      data c( 48)/  0.40429481d+02/
      data c( 49)/ -0.97198548d+01/
      data c( 50)/  0.18539893d+01/
C
c  The following conventions exist for the internuclear distances:
c     r(1) = r(OH),  r(2) = r(HH),  r(3) = r(OH)
c  Convert from aph to staemmler-schinke coordinate convention
c  Convert distances from au to angstroms
      r(1) = rij(3)*angtoau
      r(2) = rij(1)*angtoau
      r(3) = rij(2)*angtoau
c
      rmax=dmax1(r(1),r(3))/0.529177249d0
      if(rmax.ge.5.5d0)goto 10
C
C
      call func(r,phi,norder,re,alpha)
C
      pot=0.0d0
      do 1 i=1,norder
 1    pot=pot+c(i)*phi(i)
      pot=pot+rmorse(r(1))*xsi(r(2),r(3))
     *       +rmorse(r(3))*xsi(r(2),r(1))
C
      if(rmax.le.3.5d0)goto 20
C
      vint=rint(r)
      coscos=cos(1.5707963d0*(rmax-3.5d0))
      gw=0.5d0+0.5d0*coscos
      fw=0.5d0-0.5d0*coscos
      pot=pot*gw+vint*fw
      goto 20
C
 10   continue
      pot=rint(r)
C
C
c  Shift the zero of energy by De(OH)
 20   pot=(pot + 4.621d0)/autoev
      return
      end
