SUBROUTINE sample_distribution(params,g,eta,go,etao)
TYPE(KORC_PARAMS), INTENT(IN) :: params
REAL(rp), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: g
REAL(rp), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: eta
REAL(rp), INTENT(OUT) :: go
REAL(rp), INTENT(OUT) :: etao
REAL(rp) :: go_root
REAL(rp) :: etao_root
REAL(rp), DIMENSION(:), ALLOCATABLE :: p
REAL(rp) :: p_buffer
REAL(rp) :: p_test
REAL(rp) :: eta_buffer
REAL(rp) :: eta_test
REAL(rp) :: ratio
REAL(rp) :: rand_unif
REAL(rp), DIMENSION(:), ALLOCATABLE :: p_samples
REAL(rp), DIMENSION(:), ALLOCATABLE :: eta_samples
REAL(rp), DIMENSION(:), ALLOCATABLE :: p_tmp
REAL(rp), DIMENSION(:), ALLOCATABLE :: eta_tmp
REAL(rp) :: minmax
REAL(rp) :: min_p
REAL(rp) :: max_p
REAL(rp) :: min_pitch_angle
REAL(rp) :: max_pitch_angle
REAL(rp) :: deta
REAL(rp) :: dp
LOGICAL :: lp
LOGICAL :: leta
INTEGER :: num_accepted
INTEGER :: ii
INTEGER :: jj
INTEGER :: ppp
INTEGER :: nsamples
INTEGER :: mpierr
ppp = SIZE(g)
nsamples = ppp*params%mpi_params%nmpi
ALLOCATE(p(ppp))
deta = (pdf_params%max_pitch_angle - &
pdf_params%min_pitch_angle)/100.0_rp
dp = (pdf_params%max_p - pdf_params%min_p)/100.0_rp
do jj=1_idef,INT(minmax_buffer_size,idef)
minmax = pdf_params%min_p - REAL(jj,rp)*dp
if (minmax.GT.0.0_rp) then
min_p = minmax
end if
end do
max_p = pdf_params%max_p + minmax_buffer_size*dp
if (pdf_params%min_pitch_angle.GE.korc_zero) then
do jj=1_idef,INT(minmax_buffer_size,idef)
minmax = pdf_params%min_pitch_angle - REAL(jj,rp)*deta
if (minmax.GT.0.0_rp) then
min_pitch_angle = minmax
end if
end do
else
min_pitch_angle = 0.0_rp
end if
do jj=1_idef,INT(minmax_buffer_size,idef)
minmax = pdf_params%max_pitch_angle + REAL(jj,rp)*deta
if (minmax.LE.90.0_rp) then
max_pitch_angle = minmax
else
max_pitch_angle = pdf_params%max_pitch_angle
EXIT
end if
end do
if (params%mpi_params%rank.EQ.0_idef) then
ALLOCATE(p_samples(nsamples))
! Number of samples to distribute among all MPI processes
ALLOCATE(eta_samples(nsamples))
! Number of samples to distribute among all MPI processes
ALLOCATE(p_tmp(nsamples))
! Number of samples to distribute among all MPI processes
ALLOCATE(eta_tmp(nsamples))
! Number of samples to distribute among all MPI processes
!* * * Transient * * *!
call RANDOM_SEED()
call RANDOM_NUMBER(rand_unif)
eta_buffer = pdf_params%min_pitch_angle + &
(pdf_params%max_pitch_angle - pdf_params%min_pitch_angle)*rand_unif
call RANDOM_NUMBER(rand_unif)
p_buffer = pdf_params%min_p + (pdf_params%max_p - &
pdf_params%min_p)*rand_unif
ii=2_idef
do while (ii .LE. 1000_idef)
eta_test = eta_buffer + random_norm(0.0_rp,deta)
do while ((ABS(eta_test) .GT. pdf_params%max_pitch_angle).OR. &
(ABS(eta_test) .LT. pdf_params%min_pitch_angle))
eta_test = eta_buffer + random_norm(0.0_rp,deta)
end do
p_test = p_buffer + random_norm(0.0_rp,dp)
do while ((p_test.LT.pdf_params%min_p).OR.(p_test.GT. &
pdf_params%max_p))
p_test = p_buffer + random_norm(0.0_rp,dp)
end do
ratio = fRE(eta_test,p_test)/fRE(eta_buffer,p_buffer)
if (ratio .GE. 1.0_rp) then
p_buffer = p_test
eta_buffer = eta_test
ii = ii + 1_idef
else
call RANDOM_NUMBER(rand_unif)
if (rand_unif .LT. ratio) then
p_buffer = p_test
eta_buffer = eta_test
ii = ii + 1_idef
end if
end if
end do
!* * * Transient * * *!
eta_tmp(1) = eta_buffer
p_tmp(1) = p_buffer
num_accepted = 0_idef
do while(num_accepted.LT.nsamples)
ii=2_idef
do while (ii .LE. nsamples)
eta_test = eta_tmp(ii-1) + random_norm(0.0_rp,deta)
do while ((ABS(eta_test) .GT. max_pitch_angle).OR. &
(ABS(eta_test) .LT. min_pitch_angle))
eta_test = eta_tmp(ii-1) + random_norm(0.0_rp,deta)
end do
p_test = p_tmp(ii-1) + random_norm(0.0_rp,dp)
do while ((p_test.LT.min_p).OR.(p_test.GT.max_p))
p_test = p_tmp(ii-1) + random_norm(0.0_rp,dp)
end do
ratio = fRE(eta_test,p_test)/fRE(eta_tmp(ii-1),p_tmp(ii-1))
if (ratio .GE. 1.0_rp) then
p_tmp(ii) = p_test
eta_tmp(ii) = eta_test
ii = ii + 1_idef
else
call RANDOM_NUMBER(rand_unif)
if (rand_unif .LT. ratio) then
p_tmp(ii) = p_test
eta_tmp(ii) = eta_test
ii = ii + 1_idef
end if
end if
end do
eta_tmp = ABS(eta_tmp)
ii = 1_idef
do while ( (ii.LT.nsamples).AND.(num_accepted.LT.nsamples) )
lp = (p_tmp(ii).LE.pdf_params%max_p).AND.(p_tmp(ii).GE. &
pdf_params%min_p)
leta = (eta_tmp(ii).LE.pdf_params%max_pitch_angle).AND. &
(eta_tmp(ii).GE.pdf_params%min_pitch_angle)
if (lp.AND.leta) then
num_accepted = num_accepted + 1_idef
p_samples(num_accepted) = p_tmp(ii)
eta_samples(num_accepted) = eta_tmp(ii)
end if
ii = ii + 1_idef
end do
end do
go = SUM(SQRT(1.0_rp + p_samples**2))/nsamples
etao = SUM(eta_samples)/nsamples
end if
CALL MPI_SCATTER(p_samples,ppp,MPI_REAL8,p,ppp,MPI_REAL8,0,MPI_COMM_WORLD,mpierr)
CALL MPI_SCATTER(eta_samples,ppp,MPI_REAL8,eta,ppp,MPI_REAL8,0,MPI_COMM_WORLD,mpierr)
CALL MPI_BCAST(go,1,MPI_REAL8,0,MPI_COMM_WORLD,mpierr)
CALL MPI_BCAST(etao,1,MPI_REAL8,0,MPI_COMM_WORLD,mpierr)
call MPI_BARRIER(MPI_COMM_WORLD,mpierr)
g = SQRT(1.0_rp + p**2)
DEALLOCATE(p)
if (params%mpi_params%rank.EQ.0_idef) then
DEALLOCATE(p_samples)
DEALLOCATE(eta_samples)
DEALLOCATE(p_tmp)
DEALLOCATE(eta_tmp)
end if
END SUBROUTINE sample_distribution