function tanh (x)
c april 1977 edition.  w. fullerton, c3, los alamos scientific lab.
      dimension tanhcs(17)
      external alog, csevl, exp, inits, r1mach, sqrt
c
c series for tanh       on the interval  0.          to  1.00000d+00
c                                        with weighted error   9.88e-18
c                                         log weighted error  17.01
c                               significant figures required  16.25
c                                    decimal places required  17.62
c
      data tanhcs( 1) /   -.2582875664 3634710e0 /
      data tanhcs( 2) /   -.1183610633 0053497e0 /
      data tanhcs( 3) /    .0098694426 48006398e0 /
      data tanhcs( 4) /   -.0008357986 62344582e0 /
      data tanhcs( 5) /    .0000709043 21198943e0 /
      data tanhcs( 6) /   -.0000060164 24318120e0 /
      data tanhcs( 7) /    .0000005105 24190800e0 /
      data tanhcs( 8) /   -.0000000433 20729077e0 /
      data tanhcs( 9) /    .0000000036 75999055e0 /
      data tanhcs(10) /   -.0000000003 11928496e0 /
      data tanhcs(11) /    .0000000000 26468828e0 /
      data tanhcs(12) /   -.0000000000 02246023e0 /
      data tanhcs(13) /    .0000000000 00190587e0 /
      data tanhcs(14) /   -.0000000000 00016172e0 /
      data tanhcs(15) /    .0000000000 00001372e0 /
      data tanhcs(16) /   -.0000000000 00000116e0 /
      data tanhcs(17) /    .0000000000 00000009e0 /
c
      data nterms, sqeps, xmax /0, 0.0, 0.0/
c
      if (nterms.ne.0) go to 10
      nterms = inits (tanhcs, 17, 0.1*r1mach(3))
      sqeps = sqrt (3.0*r1mach(3))
      xmax = -0.5*alog(r1mach(3))
c
 10   y = abs(x)
      if (y.gt.1.) go to 20
c
      tanh = x
      if (y.gt.sqeps) tanh = x*(1. + csevl (2.*x*x-1., tanhcs, nterms))
      return
c
 20   if (y.gt.xmax) go to 30
      y = exp(y)
      tanh = sign ((y-1./y)/(y+1./y), x)
      return
c
 30   tanh = sign (1.0, x)
c
      return
      end