!***********************************************************************************
!
!    FUNCTION: ExpectFOC_Simul
!
!    PURPOSE: Specifies integrand on the rhs of the euler-equation.
!              Will be integrated by MonIntegrate subroutine, when called in
!              Simul subroutine.
!
!    AUTHOR:  Ben Malin 
!
!    DATE:  6/2004
!    REVISED: 4/2009
!***********************************************************************************

function ExpectFOC_Simul(dim1, techshock, nc)

  ! inputs:
  !    dim1 -- number of countries: N
  !    techshock -- holds the technological innovation for each of N countries and an aggregate shock
  !    nc -- index denotes which country's euler equation is being evaluated
  !
  ! output:
  ! 
  ! remarks:
  !
  
  use Params
  use StaticFOC
  implicit none
  
  !Variable Declaration
  integer, intent(in) :: dim1
  real(prec), dimension(dim1+1), intent(in) :: techshock
  integer, intent(in) :: nc
  real(prec) :: ExpectFOC_Simul
  
  real(prec) :: constom1, labtom1, captomtom, constomnc, labtomnc
  real(prec), dimension(N) :: newstate, techstate
  real(prec),dimension(2**(q-d)+1,d) :: Ttemp
  integer :: i,j
  real(prec) :: MPKtom, MPCtom, density
  
  !Returns policy choice (either consumption or capital) using Chebyshev polynomial for given state
  real(prec), external::ChebValue
  
  !Variables used in nonlinear equation solver to solve for country 1's labor supply
  integer, parameter :: itmaxA = 300
  integer, parameter :: countA = 1
  real(prec) :: RelerrA, fnormA
  real(prec), dimension(countA) :: guessA, finalA
  
  !Variables used in nonlinear equation solver to solve for cons and labor of country nc
  integer, parameter :: itmaxB = 300
  integer, parameter :: countB = 2
  real(prec) :: RelerrB, fnormB
  real(prec), dimension(countB) :: guessB, finalB
  
  !IMSL routine used for pausing program in nle solver has an error
  integer, external :: IERCD
  
  !Initialize inputs to nle solver, DNEQNF
  RelerrA = 0.0001
  fnormA = 0
  guessA = 1.0
  finalA = 1.0
  
  !Initialize inputs to nle solver, DNEQNF
  RelerrB = 0.0001
  fnormB = 0
  guessB(1) = A_tfp
  guessB(2) = 0.0
  finalB = 1.0
  
  !Construct z shock 
  do i = 1,dim1
     newstate(i) = rho * Simtech(i) + sigma * sqrt(2.0) * (techshock(i) + techshock(dim1+1))  !Constucts tech shock on [zmin, zmax] grid
  enddo
  
  do i = 1, dim1
     techstate(i) = (newstate(i) - zmin)*2.0/(zmax-zmin) - 1.0	!Converts tech shock to [-1,1] grid
  enddo

  !Construct Chebyshev polynomials from capital stocks
  Ttemp(:,1:N) = unitTtemp(:,:)
  
  !Construct Chebyshev polynomials from technological shocks
  do i = 1,dim1
     Ttemp(1,i+N) = 1.0
     Ttemp(2,i+N) = techstate(i)
     do j = 3,5		!Since q = d+2, 2**(q-d)+1 = 5
        Ttemp(j,i+N) = 2.0 * techstate(i) * Ttemp(j-1,i+N) - Ttemp(j-2,i+N)
     end do
  enddo
  
  !Compute tomorrow's consumption, labor and capital stock in t+2 using Chebyshev polynomial
  constom1 = ChebValue(Ttemp,CoeffsCons1(:,:,1),CoeffsCons2(:,:,:,:,1),CoeffsCons3(:,:,1),YCons0(1))
  captomtom = ChebValue(Ttemp,CoeffsCap1(:,:,nc),CoeffsCap2(:,:,:,:,nc),CoeffsCap3(:,:,nc),YCap0(nc))
  
  !Compute labor supply of country 1
  ExCons1 = constom1
  ExZ1 = newstate(1)
  
  call DNEQNF(FConsLabOne, RelerrA, countA, itmaxA, guessA, finalA, fnormA)
  
  if (IERCD() == 3 .or. IERCD() == 1) then
     print*, 'FConsLabOne fails to converge in ExpectFOC'
     pause
  end if
  
  labtom1 = finalA(1)
  
  !Error check to make sure solution is sensible
  if (labtom1 <= 10.0**(-5)) then
     print*, 'negative labor in ExpectFOC for country', 1, labtom1
     pause
  elseif (labtom1 > Le) then
     print*, 'labor supply above Le in ExpectFOC for country', 1, labtom1
     pause
  endif
  
  !Compute consumption and labor supply of country nc
  FOCi = nc
  ExLab1 = labtom1
  ExZnc = newstate(nc)
  guessB(1) = log(ExCons1)
  call DNEQNF(FConsLabOther, RelerrB, countB, itmaxB, guessB, finalB, fnormB)
  
  if (IERCD() == 3 .or. IERCD() == 1) then
     print*, 'FConsLabOther fails to converge in ExpectFOC'
     pause
  end if
  
  constomnc = exp(finalB(1))
  labtomnc = exp(finalB(2))
  
  !Error check to make sure solution is sensible
  if (constomnc <= 10.0**(-5)) then
     print*, 'consumption below epsilon in ExpectFOC for country', nc, constomnc
     pause
  endif
  if (labtomnc <= 10.0**(-5)) then
     print*, 'labor below epsilon in ExpectFOC for country', nc, labtomnc
     pause
  elseif (labtomnc > Le) then
     print*, 'labor supply above Le in ExpectFOC for country', nc, labtomnc
     pause
  endif
  
  !Compute MPKtom
  MPKtom = 1.0 + Fk(ktemp(nc),labtomnc,newstate(nc),nc) - AdjCost_k(ktemp(nc),captomtom)
  
  !Compute MPCtom
  MPCtom = Pareto(1) * Uc(constom1,labtom1,1)
  
  !Compute Density
  density = pie ** (-(dim1+1) * 0.5)
  
  !Specify Integrand
  ExpectFOC_Simul = MPKtom * MPCtom * density

end function ExpectFOC_Simul
