module dgp
	use basicmodule
	implicit none

	contains
	

	subroutine dgp_U
		integer	:: ir,i
		real	:: s2
		
		do ir=1,n
			i=clubid(ir)	
			s2=(1-rhos(ir)**2)*kappa2(ir)*om2
			call rnnoa(Us(:,ir))
			Us(:,ir)=rhos(ir)*CFs(:,i)+sqrt(s2)*matmul(cholSu(:,:,cindu(ir)),Us(:,ir))
			Us(0,ir)=Us(0,ir)+muU
			Xs(:,ir)=f+Us(:,ir)
		enddo
	end subroutine

	subroutine dgp_muU
		real	:: v(1)
		call rnnoa(v)
		muU=v(1)/sqrt(muprec_prior)
	endsubroutine
	
	subroutine dgp_CF
		integer	:: i,k
		real	:: s2
		do i=1,nclubs
			k=CoCid(i)
			call rnnoa(CFs(:,i))
			s2=(1-CoCrhos(i)**2)*om2*clubs2(i)
			CFs(:,i)=CoCrhos(i)*CoCFs(:,k)+sqrt(s2)*matmul(cholSu(:,:,cindCF(i)),CFs(:,i))
		enddo
	end subroutine
	
	subroutine dgp_CoCF
		integer	:: k
		do k=1,nCoCs
			call rnnoa(CoCFs(:,k))
			CoCFs(:,k)=sqrt(om2*CoCs2(k))*matmul(cholSu(:,:,cindCoC(k)),CoCFs(:,k))
		enddo
	end subroutine

	subroutine dgp_kappa2grid
		integer	:: ir,j
		real	:: unif(n),cdf(ns2grid)
		cdf=s2dist
		do j=2,ns2grid
			cdf(j)=cdf(j-1)+cdf(j)
		enddo
		call rnun(unif)
		do ir=1,n
			kappa2(ir)=s2grid(count(unif(ir)>cdf)+1)
		enddo
	end subroutine

	subroutine dgp_clubs2grid
		integer	:: i,j
		real	:: unif(nclubs),cdf(ns2grid)
		cdf=CFs2dist
		do j=2,ns2grid
			cdf(j)=cdf(j-1)+cdf(j)
		enddo
		call rnun(unif)
		do i=1,nclubs
			clubs2(i)=s2grid(count(unif(i)>cdf)+1)
		enddo
	end subroutine

	subroutine dgp_CoCs2grid
		integer	:: k,j
		real	:: unif(nCoCs),cdf(ns2grid)
		cdf=CoCs2dist
		do j=2,ns2grid
			cdf(j)=cdf(j-1)+cdf(j)
		enddo
		call rnun(unif)
		do k=1,nCoCs
			CoCs2(k)=s2grid(count(unif(k)>cdf)+1)
		enddo
	end subroutine

	subroutine dgp_s2dist
		use rngam_int
		call rngam(alphatot_prior/ns2grid,s2dist)
		s2dist=s2dist/sum(s2dist)
		call rngam(alphatot_prior/ns2grid,CFs2dist)
		CFs2dist=CFs2dist/sum(CFs2dist)
		call rngam(alphatot_prior/ns2grid,CoCs2dist)
		CoCs2dist=CoCs2dist/sum(CoCs2dist)
	end subroutine

	subroutine dgp_om2
		real	:: ssum,snu,v(1)
		ssum=om2prior
		snu=om2nuprior
		call rnchi(snu,v)
		om2=ssum/v(1)
	end subroutine
		
	subroutine dgp_rho
		integer	:: ir,j
		real	:: unif(n),cdf(nrhogrid)
		cdf=rhodist
		do j=2,nrhogrid
			cdf(j)=cdf(j-1)+cdf(j)
		enddo
		call rnun(unif)
		do ir=1,n
			rhos(ir)=rhogrid(count(unif(ir)>cdf)+1)
		enddo
	end subroutine

	subroutine dgp_CoCrho
		integer	:: i,j
		real	:: unif(nclubs),cdf(nrhogrid)
		cdf=CoCrhodist
		do j=2,nrhogrid
			cdf(j)=cdf(j-1)+cdf(j)
		enddo
		call rnun(unif)
		do i=1,nclubs
			CoCrhos(i)=rhogrid(count(unif(i)>cdf)+1)
		enddo
	end subroutine
	
	subroutine dgp_cindu
		integer	:: j,ir
		real	:: unif(n),cdf(nth)
		cdf=thdist
		do j=2,nth
			cdf(j)=cdf(j-1)+cdf(j)
		enddo
		call rnun(unif)
		do ir=1,n
			cindu(ir)=count(unif(ir)>cdf)+1
		enddo
	end subroutine

	subroutine dgp_cindCF
		integer	:: i,j
		real	:: unif(nclubs),cdf(nth)
		cdf=CFthdist
		do j=2,nth
			cdf(j)=cdf(j-1)+cdf(j)
		enddo
		call rnun(unif)
		do i=1,nclubs
			cindCF(i)=count(unif(i)>cdf)+1
		enddo
	end subroutine

	subroutine dgp_cindCoC
		integer	:: k,j
		real	:: unif(nCoCs),cdf(nth)
		cdf=CoCthdist
		do j=2,nth
			cdf(j)=cdf(j-1)+cdf(j)
		enddo
		call rnun(unif)
		do k=1,nCoCs
			cindCoC(k)=count(unif(k)>cdf)+1
		enddo
	end subroutine

	subroutine dgp_clubid
		call rnund(nclubs,clubid)
	end subroutine
	
	subroutine dgp_CoCid
		call rnund(nCoCs,CoCid)
	end subroutine

	subroutine dgp_thdist
		use rngam_int
		call rngam(alphatot_prior/size(thdist),thdist)
		thdist=thdist/sum(thdist)
	end subroutine
	
	subroutine dgp_CFthdist
		use rngam_int
		call rngam(alphatot_prior/size(CFthdist),CFthdist)
		CFthdist=CFthdist/sum(CFthdist)
	end subroutine

	subroutine dgp_CoCthdist
		use rngam_int
		call rngam(alphatot_prior/size(CoCthdist),CoCthdist)
		CoCthdist=CoCthdist/sum(CoCthdist)
	end subroutine
	
	subroutine dgp_rhodist
		use rngam_int
		call rngam(alphatot_prior/size(rhodist),rhodist)
		rhodist=rhodist/sum(rhodist)
	end subroutine

	subroutine dgp_CoCrhodist
		use rngam_int
		call rngam(alphatot_prior/size(Cocrhodist),CoCrhodist)
		CoCrhodist=CoCrhodist/sum(CoCrhodist)
	end subroutine

	subroutine dgp_F
		real	:: fa(0:q)
		call rnnoa(fm)
		fm=sqrt(sigma2_F(1))*matmul(cholfm(:,:,cindf),fm)
		call rnnoa(fa)
		f=fm+sqrt(sigma2_F(2))*matmul(cholfa(:,:),fa)
		f(0:1)=f(0:1)+mt_f
	end subroutine

	subroutine dgp_mt_f
		call rnnoa(mt_f)
		mt_f=mt_f/sqrt(muprec_prior)
	end subroutine

	subroutine dgp_cindF
		real	:: unif(1)
		call rnun(unif)
		cindf=unif(1)*nf+1
	end subroutine
		
	subroutine dgp_sigma2_F
		real	:: ssum
		real	:: snu
		real	:: v(1),unif(1)
		real	:: cdf(ns2fm)
		integer	:: l,j
		
		do l=1,ns2fm
			cdf(l)=merge(l,ns2fm+1-l,2*l<=ns2fm)
		enddo
		do j=2,ns2fm
			cdf(j)=cdf(j-1)+cdf(j)
		enddo
		cdf=cdf/cdf(ns2fm)
		call rnun(unif)
		sigma2_F(1)=s2fmtab(count(unif(1)>cdf)+1)
		
		ssum=sigma2_Faprior
		snu=sigma2_Fanuprior
		call rnchi(snu,v)
		sigma2_F(2)=ssum/v(1)
	end subroutine	

end module