subroutine model(lpar,par,  gamma,ss,psil,psiy,psix,psic,info)
!
! Model   Compute equilibrium for the BCA prototype model with 
!         unit root technology.
!       

!         Ellen McGrattan, 6-2-06
!         Revised, ERM, 10-2-09


  implicit none

  integer, parameter               :: lwork=5000
  integer, intent(in)              :: lpar
  integer                          :: i,j,l,info
  integer, dimension(lwork)        :: ipiv
  real, dimension(lpar),intent(in) :: par
  real, dimension(lpar)            :: param
  real, dimension(lwork)           :: work
  real, dimension(4,1)             :: sbar,p0,gam
  real, dimension(4,4)             :: imat,p,q,atem,ptem
  real                             :: irate,beta,delta,gn,psi,sigma,theta,zs, &
                                      tauls,tauxs,gs,kls,ls,a,b,ks,cs,ys,xs,  &
                                      rp,rm,a0,a1,a2,lam1,lam2,gamk,philh,    &
                                      philk,philz,phill,philx,philg,philkp,   &
                                      phiyk,phiyz,phiyl,phiyx,phiyg,phiykp,   &
                                      phixk,phixz,phixl,phixx,phixg,phixkp,   &
                                      phick,phicz,phicl,phicx,phicg,phickp
  real, dimension(1,4)             :: b0,b1,btem
  real, dimension(11)              :: z,del,zp,zm,dr
  real, dimension(6),intent(out)   :: gamma,ss,psil,psiy,psix,psic

  !___________________________________________________________________________
  ! PARAMETERS
  !
  info     = 0
  irate    = par(1)
  delta    = par(2)
  gn       = par(3)
  psi      = par(4)
  sigma    = par(5)
  theta    = par(6)
  l        = 11
  imat     = 0.
  do i=1,4
    sbar(i,1) = par(i+6)
    imat(i,i) = 1. 
    do j=1,4
      p(i,j)  = par(l)
      l       = l+1
    enddo
  enddo
  do i=1,4
    do j=1,4
      q(i,j)  = par(l)
      l       = l+1
    enddo
  enddo
  ptem        = imat-p
  p0          = matmul(ptem,sbar)

  !____________________________________________________________________________
  ! COMPUTE EQUILIBRIUM
  ! a. Steady state
  !
  zs       = exp(sbar(1,1))
  beta     = zs**sigma/(1.+irate)
  tauls    = sbar(2,1)
  tauxs    = sbar(3,1)
  gs       = exp(sbar(4,1))
  kls      = ((1+tauxs)*(1-beta*zs**(-sigma)*(1-delta))/                      &
              (beta*zs**(-sigma)*theta))**(1/(theta-1))*zs
  a        = kls**(theta-1)*zs**(-theta)-(1+gn)+(1-delta)/zs
  b        = (1-tauls)*(1-theta)*kls**theta*zs**(-theta)/psi
  ks       = (b+gs)/(a+b/kls)
  cs       = a*ks-gs
  ls       = ks/kls
  ys       = (ks/zs)**theta*ls**(1-theta)
  xs       = ys-cs-gs

  !
  ! b. Call subroutine with residuals:
  !

  z        = (/ log(ks), log(ks), log(ks), log(zs), log(zs), tauls, tauls,    &
                tauxs, tauxs, log(gs), log(gs) /)
  do i=1,11
    del(i) = max(abs(z(i)*.00001),.00000001)
  enddo
                                                                                
  param    = par
  param(1) = beta
  do i=1,11
    zp     = z
    zm     = z
    zp(i)  = z(i)+del(i)
    zm(i)  = z(i)-del(i)
    call reswedgeur(param,zp, rp)
    call reswedgeur(param,zm, rm)
    
    dr(i)  = (rp-rm)/(2.*del(i))

  enddo

  !
  ! c. Solution:  log k[t+1] = gamma0 + gamk* log k[t] + gamma* S[t]
  !

  a0       = dr(1)
  a1       = dr(2)
  a2       = dr(3)

  l        = 1
  do i=4,10,2
    b0(1,l)= dr(i)
    b1(1,l)= dr(i+1)
    l      = l+1
  enddo

  lam1     = .5*(-a1+sqrt(a1*a1-4.*a0*a2))/a0
  lam2     = .5*(-a1-sqrt(a1*a1-4.*a0*a2))/a0
  if (abs(lam1)<1) then
    gamk   = lam1
  else
    gamk   = lam2
  endif
  atem  =  (gamk*a0+a1)*imat+a0*transpose(p)
  call dgetrf(4,4,atem,4,ipiv,info)
  if (info.ne.0) return
  call dgetri(4,atem,4,ipiv,work,lwork,info)
  if (info.ne.0) return
  btem       = matmul(b0,p)+b1
  gam        = -matmul(atem,transpose(btem))
  gamma(1)   = gamk
  gamma(2:5) = gam(1:4,1)
  gamma(6)   = (1.-gamk)*log(ks)-gam(1,1)*log(zs)-gam(2,1)*tauls-            & 
               gam(3,1)*tauxs-gam(4,1)*log(gs)

  philh      =-ys*(1.-theta)*(psi+theta*(1.-tauls)*(1.-ls)/ls+1.-tauls)
  philk      = (psi*ys*theta+psi*(1.-delta)*ks/zs-                           &
                theta*(1.-theta)*(1.-tauls)*ys*(1-ls)/ls)/philh
  philz      =-philk
  phill      = ((1.-theta)*ys*(1-ls)/ls)/philh
  philx      = 0.
  philg      = -psi*gs/philh
  philkp     = -psi*(1.+gn)*ks/philh
  phiyk      = theta+(1-theta)*philk
  phiyz      = -phiyk
  phiyl      = (1.-theta)*phill
  phiyx      = 0.
  phiyg      = (1.-theta)*philg
  phiykp     = (1.-theta)*philkp
  phixk      = -ks/(zs*xs)*(1.-delta)
  phixkp     = ks/xs*(1+gn)
  phixz      = -phixk
  phixl      = 0.
  phixx      = 0.
  phixg      = 0.
  phick      = theta*ys/cs+(1.-theta)*philk*ys/cs+(1.-delta)*ks/(cs*zs)
  phicz      = -phick
  phicl      = (1.-theta)*phill*ys/cs
  phicx      = 0.
  phicg      = (1.-theta)*philg*ys/cs-gs/cs
  phickp     = (1.-theta)*philkp*ys/cs-(1.+gn)*ks/cs

                                                                                
  psil(1)     = philk+philkp*gamk
  psil(2)     = philz+philkp*gam(1,1)
  psil(3)     = phill+philkp*gam(2,1)
  psil(4)     = philx+philkp*gam(3,1)
  psil(5)     = philg+philkp*gam(4,1)
  psil(6)     = log(ls)-psil(1)*log(ks)-psil(2)*log(zs)-psil(3)*tauls        &
                       -psil(4)*tauxs-psil(5)*log(gs)

  psiy(1)     = phiyk+phiykp*gamk
  psiy(2)     = phiyz+phiykp*gam(1,1)
  psiy(3)     = phiyl+phiykp*gam(2,1)
  psiy(4)     = phiyx+phiykp*gam(3,1)
  psiy(5)     = phiyg+phiykp*gam(4,1)
  psiy(6)     = log(ys)-psiy(1)*log(ks)-psiy(2)*log(zs)-psiy(3)*tauls        &
                       -psiy(4)*tauxs-psiy(5)*log(gs)

  psix(1)     = phixk+phixkp*gamk
  psix(2)     = phixz+phixkp*gam(1,1)
  psix(3)     = phixl+phixkp*gam(2,1)
  psix(4)     = phixx+phixkp*gam(3,1)
  psix(5)     = phixg+phixkp*gam(4,1)
  psix(6)     = log(xs)-psix(1)*log(ks)-psix(2)*log(zs)-psix(3)*tauls        &
                       -psix(4)*tauxs-psix(5)*log(gs)

  psic(1)     = phick+phickp*gamk
  psic(2)     = phicz+phickp*gam(1,1)
  psic(3)     = phicl+phickp*gam(2,1)
  psic(4)     = phicx+phickp*gam(3,1)
  psic(5)     = phicg+phickp*gam(4,1)
  psic(6)     = log(cs)-psic(1)*log(ks)-psic(2)*log(zs)-psic(3)*tauls        &
                       -psic(4)*tauxs-psic(5)*log(gs)

  ss          = (/ log(ks), log(zs), tauls, tauxs, log(gs), 1. /)

end subroutine model