!*********************************************************************
! program    : mpi_samp6.f90
! programmer : makoto nakajima
! description: run a simple monte carlo simulation.
! date       : april 30, 2006
!*********************************************************************
program main

	!******* prohibit implicit variable declaration *******
	implicit none     !prohibit implicit declaration of variables

	!******* activate mpi library *******
	include 'mpif.h'  !include mpi library

	!******* external procedure *******
	real(8),external:: optdecfun

	!******* variables related mpi *******
	integer:: ierror !return error message from mpi subroutines
	integer:: id     !identification number of each processor
	integer:: nproc  !total number of processors

	!******* model parameters *******
	integer,parameter:: ns=2 !number of shocks
	real(8),dimension(ns):: sval  !value of shocks
	real(8),dimension(ns,ns):: stran !transition matrix

	!******* other variable declaration *******
	integer,parameter:: master_id=0   !node which works as the frontend.
	integer,parameter:: n_all=1000000 !total number of agents used.
	integer,parameter:: nt=5000 !length of periods for simulation.
	integer:: t0     !index for period
	integer:: n_each !number of agents assigned to each node
	integer:: i0     !index for individual agent
	integer:: s0,s1  !index for individual's s
	real(8):: k_each !sum of capital at each node
	real(8):: k_all  !sum of capital across all nodes
	real(8):: l_each !sum of labor at each node
	real(8):: l_all  !sum of labor across all nodes
	real(8):: rnum   !draw of random number
	real(8),dimension(:),allocatable:: aind0,aind1 !asset of each agent
	integer,dimension(:),allocatable:: sind0,sind1 !s of each agent

	!******* initialization of mpi environment *******
	call mpi_init(ierror)

	!******* obtain id for each node *******
	call mpi_comm_rank(mpi_comm_world, id, ierror)

	!******* obtains the number of nodes *******
	call mpi_comm_size(mpi_comm_world, nproc, ierror)

	!******* set shock *******
	sval(1)=-0.20_8
	sval(2)= 0.20_8
	stran(1,1)=0.980_8
	stran(1,2)=1.0_8-stran(1,1)
	stran(2,2)=0.980_8
	stran(2,1)=1.0_8-stran(2,2)

	!******* compute agents assigned to each node *******
	do i0=1,n_all
		if (i0*nproc>=n_all) exit
	end do
	if (id/=nproc) then
		n_each=i0
	else
		n_each=n_all-(nproc-1)*i0
	end if

	!******* allocate arrays *******
	allocate(aind0(n_each),aind1(n_each),sind0(n_each),sind1(n_each))

	!******* set initial condition for each agent *******
	aind0(:)=0.0_8
	sind0(:)=1

	!******* loop for simulation *******
	do t0=1,nt

		!******* gather information every 100 periods
		if (mod(t0,100)==0) then
			!******* compute aggregate data at each node *******
			k_each=0.0_8
			l_each=0.0_8
			do i0=1,n_each
				k_each=k_each+aind0(i0)
				l_each=l_each+exp(sval(sind0(i0)))
			end do

			!******* sum up all capital to master node *******
			call mpi_reduce(k_each,k_all,1,mpi_double_precision,&
			mpi_sum,master_id,mpi_comm_world,ierror)

			!******* sum up all labor to master node *******
			call mpi_reduce(l_each,l_all,1,mpi_double_precision,&
			mpi_sum,master_id,mpi_comm_world,ierror)

			!******* take average *******
			k_all=k_all/real(n_all,8)
			l_all=l_all/real(n_all,8)

			!******* master prints out *******
			if (id==master_id) then
				write(*,"('period=',i5,' k=',f10.5,' l=',f10.5,&
				' k/l=',f10.5)") t0,k_all,l_all,k_all/l_all
			end if
		end if

		!******* loop for each agent *******
		do i0=1,n_each
			!******* update s *******
			s0=sind0(i0)
			call random_number(rnum)
			do s1=1,ns
				if (sum(stran(s0,1:s1))>=rnum) exit
			end do
			sind1(i0)=s1
			!******* update a *******
			aind1(i0)=optdecfun(sind0(i0),aind0(i0))
		end do

		!******* update individual states *******
		aind0=aind1
		sind0=sind1
	end do

	!******* finalize mpi environment *******
	call mpi_finalize(ierror)

end program main

!*********************************************************************
! function   : optdecfun
! description: (exogenously provided) optimal decision rule.
!              in a typical model economy, this is obtained by solving
!              agent's optimization problem.
!*********************************************************************
function optdecfun(s0,a0r) result (a1r)
	implicit none
	real(8):: a1r !output (asset in the next period)
	real(8),intent(in):: a0r !current asset
	integer,intent(in):: s0  !current shock

	!******* compute asset in the next period *******
	if (s0==1) then
		a1r=0.80_8*a0r
	else if (s0==2) then
		a1r=1.0_8+0.90_8*a0r
	end if

	!******* end of function *******
	return

end function optdecfun