subroutine ggaec(rs,zet,t,ec,h)
c  called by subroutine ggaexc
c  gga91 correlation
c  input rs  : seitz radius
c  input zet : relative spin polarization
c  input t   : abs(grad d)/(d*2.*ks*g)
c  input     : correlation energy per electron (ec)
c  output h  : nonlocal part of correlation energy per electron
      implicit double precision (a-h,o-z)
      data xnu,cc0,cx,alf/15.75592d0,0.004235d0,-0.001667212d0,0.09d0/
      data c1,c2,c3,c4/0.002568d0,0.023266d0,7.389d-6,8.723d0/
      data c5,c6,a4/0.472d0,7.389d-2,100.d0/
      data thrd2/0.666666666667d0/
      pi = 4.d0*datan(1.d0)
      fk = 1.91915829d0/rs
      sk = dsqrt(4.d0*fk/pi)
      g = ((1.d0+zet)**thrd2+(1.d0-zet)**thrd2)/2.d0
      bet = xnu*cc0
      delt = 2.d0*alf/bet
      g3 = g**3
      g4 = g3*g
      pon = -delt*ec/(g3*bet)
      b = delt/(dexp(pon)-1.d0)
      b2 = b*b
      t2 = t*t
      t4 = t2*t2
      rs2 = rs*rs
      rs3 = rs2*rs
      q4 = 1.d0+b*t2
      q5 = 1.d0+b*t2+b2*t4
      q6 = c1+c2*rs+c3*rs2
      q7 = 1.d0+c4*rs+c5*rs2+c6*rs3
      cc = -cx + q6/q7
      r0 = (sk/fk)**2
      r1 = a4*r0*g4
      coeff = cc-cc0-3.d0*cx/7.d0
      r2 = xnu*coeff*g3
      r3 = dexp(-r1*t2)
      h0 = g3*(bet/delt)*dlog(1.d0+delt*q4*t2/q5)
      h1 = r3*r2*t2
      h = h0 + h1
      return
      end