!*****************************************************************************
!
!    Subroutine: GHIntegrate
!
!    Purpose: Calculates an integral using four-point Gauss-Hermite
!             quadrature.
!
!    Author: Ben Malin
!
!    Date: 4/09
!    Revised: 
!
!*****************************************************************************

subroutine GHIntegrate(nc, IntResult, check)

  ! inputs:
  !   nc -- Index of country.
  !   check -- Denotes whether GHIntegrate is called in policy function iteration
  !             or in simulation.
  !
  ! output:
  !   IntResult -- The value of the integral. 
  !
  ! remarks:
  !   The subroutine constructs the integral on the RHS of the Euler Equation
  !   for country nc.
  !
  
  use Params
  
  implicit none

  !Declaration of Variables
  integer,intent(in) :: nc
  real(prec),intent(out) :: IntResult 
  integer,intent(in) :: check
  integer :: i
  
  !Store values of integrand
  real(prec),dimension(GHdim) :: Integrand    
  
  !Returns Integrand on RHS of Euler Equation at specified points
  real(prec), external :: ExpectFOC
  real(prec), external :: ExpectFOC_Simul
  
  if (check == 1) then
     !Calculate values of integrand at Hermite nodes
     do i = 1, GHdim
        Integrand(i) = ExpectFOC(N, GHnodes(i,:), nc)
     enddo
  elseif (check == 3) then
     !Calculate values of integrand at Hermite nodes
     do i = 1, GHdim
        Integrand(i) = ExpectFOC_Simul(N, GHnodes(i,:), nc)
     enddo
  endif
  
  !Weight and Sum values of integrand  -- See formula (7.5.11) on pg. 275 of Judd 
  IntResult = 0.0
  
  do i = 1, GHdim
     IntResult = IntResult + GHweights(i) * Integrand(i)
  enddo

end subroutine GHIntegrate
