Skip to content

Commit

Permalink
change double intrinsics to generics
Browse files Browse the repository at this point in the history
  • Loading branch information
@belericant
belericant committed Mar 10, 2023
1 parent 9820936 commit 6d73bf8
Show file tree
Hide file tree
Showing 11 changed files with 102 additions and 102 deletions.
2 changes: 1 addition & 1 deletion src/common/m_helper.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,7 +116,7 @@ subroutine s_comp_n_from_cons(vftmp, nRtmp, ntmp, weight)
real(wp), dimension(nb) :: weight

nR3 = dot_product(weight, nRtmp**3._wp)
ntmp = DSQRT((4._wp*pi/3._wp)*nR3/vftmp)
ntmp = SQRT((4._wp*pi/3._wp)*nR3/vftmp)

end subroutine s_comp_n_from_cons

Expand Down
28 changes: 14 additions & 14 deletions src/post_process/m_global_parameters.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -679,7 +679,7 @@ subroutine s_initialize_nonpoly

temp = 293.15_wp
D_m = (0.242_wp * (10._wp ** -(4)))
uu = DSQRT(pl0/rhol0)
uu = SQRT(pl0/rhol0)

omega_ref = 3._wp*k_poly*Ca + 2._wp*(3._wp*k_poly - 1._wp)/Web

Expand All @@ -688,10 +688,10 @@ subroutine s_initialize_nonpoly
R_n = Ru/M_n
R_v = Ru/M_v
! phi_vn & phi_nv (phi_nn = phi_vv = 1)
phi_vn = (1._wp + DSQRT(mu_v/mu_n)*(M_n/M_v)**(0.25_wp))**2 &
/(DSQRT(8._wp)*DSQRT(1._wp + M_v/M_n))
phi_nv = (1._wp + DSQRT(mu_n/mu_v)*(M_v/M_n)**(0.25_wp))**2 &
/(DSQRT(8._wp)*DSQRT(1._wp + M_n/M_v))
phi_vn = (1._wp + SQRT(mu_v/mu_n)*(M_n/M_v)**(0.25_wp))**2 &
/(SQRT(8._wp)*SQRT(1._wp + M_v/M_n))
phi_nv = (1._wp + SQRT(mu_n/mu_v)*(M_v/M_n)**(0.25_wp))**2 &
/(SQRT(8._wp)*SQRT(1._wp + M_n/M_v))
! internal bubble pressure
pb0 = pl0 + 2._wp*ss/(R0ref*R0)

Expand Down Expand Up @@ -726,7 +726,7 @@ subroutine s_initialize_nonpoly
! keeps a constant (cold liquid assumption)
Tw = 1._wp
! natural frequencies
omegaN = DSQRT(3._wp*k_poly*Ca + 2._wp*(3._wp*k_poly - 1._wp)/(Web*R0))/R0
omegaN = SQRT(3._wp*k_poly*Ca + 2._wp*(3._wp*k_poly - 1._wp)/(Web*R0))/R0

pl0 = 1._wp
do ir = 1, Nb
Expand Down Expand Up @@ -861,31 +861,31 @@ subroutine s_simpson(Npt)
!R0mx = 1.3D0

sd = poly_sigma
R0mn = 0.8_wp*DEXP(-2.8_wp*sd)
R0mx = 0.2_wp*DEXP(9.5_wp*sd) + 1._wp
R0mn = 0.8_wp*EXP(-2.8_wp*sd)
R0mx = 0.2_wp*EXP(9.5_wp*sd) + 1._wp

! phi = ln( R0 ) & return R0
do ir = 1, Npt
phi(ir) = DLOG(R0mn) &
+ dble(ir - 1)*DLOG(R0mx/R0mn)/dble(Npt - 1)
R0(ir) = DEXP(phi(ir))
phi(ir) = LOG(R0mn) &
+ dble(ir - 1)*LOG(R0mx/R0mn)/dble(Npt - 1)
R0(ir) = EXP(phi(ir))
end do
dphi = phi(2) - phi(1)

! weights for quadrature using Simpson's rule
do ir = 2, Npt - 1
! Gaussian
tmp = DEXP(-0.5_wp*(phi(ir)/sd)**2)/DSQRT(2._wp*pi)/sd
tmp = EXP(-0.5_wp*(phi(ir)/sd)**2)/SQRT(2._wp*pi)/sd
if (mod(ir, 2) == 0) then
weight(ir) = tmp*4._wp*dphi/3._wp
else
weight(ir) = tmp*2._wp*dphi/3._wp
end if
end do

tmp = DEXP(-0.5_wp*(phi(1)/sd)**2)/DSQRT(2._wp*pi)/sd
tmp = EXP(-0.5_wp*(phi(1)/sd)**2)/SQRT(2._wp*pi)/sd
weight(1) = tmp*dphi/3._wp
tmp = DEXP(-0.5_wp*(phi(Npt)/sd)**2)/DSQRT(2._wp*pi)/sd
tmp = EXP(-0.5_wp*(phi(Npt)/sd)**2)/SQRT(2._wp*pi)/sd
weight(Npt) = tmp*dphi/3._wp

end subroutine s_simpson
Expand Down
18 changes: 9 additions & 9 deletions src/pre_process/m_assign_variables.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -321,10 +321,10 @@ subroutine s_assign_patch_species_primitive_variables_bubbles(patch_id, j, k, l,
q_prim_vf(bub_idx%fullmom(i, 0, 2))%sf(j, k, l) = muV**2 + sigV**2
else if (dist_type == 2) then
q_prim_vf(bub_idx%fullmom(i, 0, 0))%sf(j, k, l) = 1._wp
q_prim_vf(bub_idx%fullmom(i, 1, 0))%sf(j, k, l) = dexp((sigR**2)/2._wp)*muR
q_prim_vf(bub_idx%fullmom(i, 1, 0))%sf(j, k, l) = exp((sigR**2)/2._wp)*muR
q_prim_vf(bub_idx%fullmom(i, 0, 1))%sf(j, k, l) = muV
q_prim_vf(bub_idx%fullmom(i, 2, 0))%sf(j, k, l) = dexp((sigR**2)*2)*(muR**2)
q_prim_vf(bub_idx%fullmom(i, 1, 1))%sf(j, k, l) = dexp((sigR**2)/2)*muR*muV
q_prim_vf(bub_idx%fullmom(i, 2, 0))%sf(j, k, l) = exp((sigR**2)*2)*(muR**2)
q_prim_vf(bub_idx%fullmom(i, 1, 1))%sf(j, k, l) = exp((sigR**2)/2)*muR*muV
q_prim_vf(bub_idx%fullmom(i, 0, 2))%sf(j, k, l) = muV**2 + sigV**2
end if

Expand Down Expand Up @@ -400,10 +400,10 @@ subroutine s_assign_patch_species_primitive_variables_bubbles(patch_id, j, k, l,
q_prim_vf(bub_idx%fullmom(i, 0, 2))%sf(j, k, l) = muV**2 + sigV**2
else if (dist_type == 2) then
q_prim_vf(bub_idx%fullmom(i, 0, 0))%sf(j, k, l) = 1._wp
q_prim_vf(bub_idx%fullmom(i, 1, 0))%sf(j, k, l) = dexp((sigR**2)/2._wp)*muR
q_prim_vf(bub_idx%fullmom(i, 1, 0))%sf(j, k, l) = exp((sigR**2)/2._wp)*muR
q_prim_vf(bub_idx%fullmom(i, 0, 1))%sf(j, k, l) = muV
q_prim_vf(bub_idx%fullmom(i, 2, 0))%sf(j, k, l) = dexp((sigR**2)*2._wp)*(muR**2)
q_prim_vf(bub_idx%fullmom(i, 1, 1))%sf(j, k, l) = dexp((sigR**2)/2._wp)*muR*muV
q_prim_vf(bub_idx%fullmom(i, 2, 0))%sf(j, k, l) = exp((sigR**2)*2._wp)*(muR**2)
q_prim_vf(bub_idx%fullmom(i, 1, 1))%sf(j, k, l) = exp((sigR**2)/2._wp)*muR*muV
q_prim_vf(bub_idx%fullmom(i, 0, 2))%sf(j, k, l) = muV**2 + sigV**2
end if
else
Expand Down Expand Up @@ -507,10 +507,10 @@ subroutine s_assign_patch_species_primitive_variables_bubbles(patch_id, j, k, l,
q_prim_vf(bub_idx%fullmom(i, 0, 2))%sf(j, k, l) = muV**2 + sigV**2
else if (dist_type == 2) then
q_prim_vf(bub_idx%fullmom(i, 0, 0))%sf(j, k, l) = 1._wp
q_prim_vf(bub_idx%fullmom(i, 1, 0))%sf(j, k, l) = dexp((sigR**2)/2._wp)*muR
q_prim_vf(bub_idx%fullmom(i, 1, 0))%sf(j, k, l) = exp((sigR**2)/2._wp)*muR
q_prim_vf(bub_idx%fullmom(i, 0, 1))%sf(j, k, l) = muV
q_prim_vf(bub_idx%fullmom(i, 2, 0))%sf(j, k, l) = dexp((sigR**2)*2._wp)*(muR**2)
q_prim_vf(bub_idx%fullmom(i, 1, 1))%sf(j, k, l) = dexp((sigR**2)/2._wp)*muR*muV
q_prim_vf(bub_idx%fullmom(i, 2, 0))%sf(j, k, l) = exp((sigR**2)*2._wp)*(muR**2)
q_prim_vf(bub_idx%fullmom(i, 1, 1))%sf(j, k, l) = exp((sigR**2)/2._wp)*muR*muV
q_prim_vf(bub_idx%fullmom(i, 0, 2))%sf(j, k, l) = muV**2 + sigV**2
end if
else
Expand Down
28 changes: 14 additions & 14 deletions src/pre_process/m_global_parameters.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -673,7 +673,7 @@ contains

temp = 293.15_wp
D_m = (0.242_wp * (10._wp ** -(4)))
uu = DSQRT(pl0/rhol0)
uu = SQRT(pl0/rhol0)

omega_ref = 3._wp*k_poly*Ca + 2._wp*(3._wp*k_poly - 1._wp)/Web

Expand All @@ -682,10 +682,10 @@ contains
R_n = Ru/M_n
R_v = Ru/M_v
! phi_vn & phi_nv (phi_nn = phi_vv = 1)
phi_vn = (1._wp + DSQRT(mu_v/mu_n)*(M_n/M_v)**(0.25_wp))**2 &
/(DSQRT(8._wp)*DSQRT(1._wp + M_v/M_n))
phi_nv = (1._wp + DSQRT(mu_n/mu_v)*(M_v/M_n)**(0.25_wp))**2 &
/(DSQRT(8._wp)*DSQRT(1._wp + M_n/M_v))
phi_vn = (1._wp + SQRT(mu_v/mu_n)*(M_n/M_v)**(0.25_wp))**2 &
/(SQRT(8._wp)*SQRT(1._wp + M_v/M_n))
phi_nv = (1._wp + SQRT(mu_n/mu_v)*(M_v/M_n)**(0.25_wp))**2 &
/(SQRT(8._wp)*SQRT(1._wp + M_n/M_v))
! internal bubble pressure
pb0 = pl0 + 2._wp*ss/(R0ref*R0)

Expand Down Expand Up @@ -722,7 +722,7 @@ contains
! keeps a constant (cold liquid assumption)
Tw = 1._wp
! natural frequencies
omegaN = DSQRT(3._wp*k_poly*Ca + 2._wp*(3._wp*k_poly - 1._wp)/(Web*R0))/R0
omegaN = SQRT(3._wp*k_poly*Ca + 2._wp*(3._wp*k_poly - 1._wp)/(Web*R0))/R0

pl0 = 1._wp
do ir = 1, Nb
Expand Down Expand Up @@ -847,30 +847,30 @@ contains
!R0mx = 150.D0

sd = poly_sigma
R0mn = 0.8_wp*DEXP(-2.8_wp*sd)
R0mx = 0.2_wp*DEXP(9.5_wp*sd) + 1._wp
R0mn = 0.8_wp*EXP(-2.8_wp*sd)
R0mx = 0.2_wp*EXP(9.5_wp*sd) + 1._wp

! phi = ln( R0 ) & return R0
do ir = 1, nb
phi(ir) = DLOG(R0mn) &
+ dble(ir - 1)*DLOG(R0mx/R0mn)/dble(nb - 1)
R0(ir) = DEXP(phi(ir))
phi(ir) = LOG(R0mn) &
+ dble(ir - 1)*LOG(R0mx/R0mn)/dble(nb - 1)
R0(ir) = EXP(phi(ir))
end do
dphi = phi(2) - phi(1)

! weights for quadrature using Simpson's rule
do ir = 2, nb - 1
! Gaussian
tmp = DEXP(-0.5_wp*(phi(ir)/sd)**2)/DSQRT(2._wp*pi)/sd
tmp = EXP(-0.5_wp*(phi(ir)/sd)**2)/SQRT(2._wp*pi)/sd
if (mod(ir, 2) == 0) then
weight(ir) = tmp*4._wp*dphi/3._wp
else
weight(ir) = tmp*2._wp*dphi/3._wp
end if
end do
tmp = DEXP(-0.5_wp*(phi(1)/sd)**2)/DSQRT(2._wp*pi)/sd
tmp = EXP(-0.5_wp*(phi(1)/sd)**2)/SQRT(2._wp*pi)/sd
weight(1) = tmp*dphi/3._wp
tmp = DEXP(-0.5_wp*(phi(nb)/sd)**2)/DSQRT(2._wp*pi)/sd
tmp = EXP(-0.5_wp*(phi(nb)/sd)**2)/SQRT(2._wp*pi)/sd
weight(nb) = tmp*dphi/3._wp
end subroutine s_simpson
Expand Down
14 changes: 7 additions & 7 deletions src/pre_process/m_patches.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -284,7 +284,7 @@ subroutine s_varcircle(patch_id, patch_id_fp, q_prim_vf) ! ---------------------
! the current patch are assigned to this cell.
do j = 0, n
do i = 0, m
myr = dsqrt((x_cc(i) - x_centroid)**2 &
myr = sqrt((x_cc(i) - x_centroid)**2 &
+ (y_cc(j) - y_centroid)**2)

if (myr <= radius + thickness/2._wp .and. &
Expand All @@ -295,7 +295,7 @@ subroutine s_varcircle(patch_id, patch_id_fp, q_prim_vf) ! ---------------------
eta, q_prim_vf, patch_id_fp)

q_prim_vf(alf_idx)%sf(i, j, 0) = patch_icpp(patch_id)%alpha(1)* &
dexp(-0.5_wp*((myr - radius)**2._wp)/(thickness/3._wp)**2._wp)
exp(-0.5_wp*((myr - radius)**2._wp)/(thickness/3._wp)**2._wp)
end if

end do
Expand Down Expand Up @@ -340,7 +340,7 @@ subroutine s_3dvarcircle(patch_id, patch_id_fp, q_prim_vf) ! -------------------
do k = 0, p
do j = 0, n
do i = 0, m
myr = dsqrt((x_cc(i) - x_centroid)**2 &
myr = sqrt((x_cc(i) - x_centroid)**2 &
+ (y_cc(j) - y_centroid)**2)

if (myr <= radius + thickness/2._wp .and. &
Expand All @@ -351,7 +351,7 @@ subroutine s_3dvarcircle(patch_id, patch_id_fp, q_prim_vf) ! -------------------
eta, q_prim_vf, patch_id_fp)

q_prim_vf(alf_idx)%sf(i, j, k) = patch_icpp(patch_id)%alpha(1)* &
dexp(-0.5_wp*((myr - radius)**2._wp)/(thickness/3._wp)**2._wp)
exp(-0.5_wp*((myr - radius)**2._wp)/(thickness/3._wp)**2._wp)
end if

end do
Expand Down Expand Up @@ -724,7 +724,7 @@ subroutine s_1D_analytical(patch_id, patch_id_fp, q_prim_vf) ! -----------------
!what variables to alter
!bump in pressure
q_prim_vf(E_idx)%sf(i, 0, 0) = q_prim_vf(E_idx)%sf(i, 0, 0)* &
(1._wp + 0.2_wp*dexp(-1._wp*((x_cb(i) - x_centroid)**2._wp)/(2._wp*0.005_wp)))
(1._wp + 0.2_wp*exp(-1._wp*((x_cb(i) - x_centroid)**2._wp)/(2._wp*0.005_wp)))

!bump in void fraction
!q_prim_vf(adv_idx%beg)%sf(i,0,0) = q_prim_vf(adv_idx%beg)%sf(i,0,0) * &
Expand Down Expand Up @@ -883,11 +883,11 @@ subroutine s_2D_analytical(patch_id, patch_id_fp, q_prim_vf) ! -----------------
!what variables to alter
!x-y bump in pressure
q_prim_vf(E_idx)%sf(i, j, 0) = q_prim_vf(E_idx)%sf(i, j, 0)* &
(1._wp + 0.2_wp*dexp(-1._wp*((x_cb(i) - x_centroid)**2._wp + (y_cb(j) - y_centroid)**2._wp)/(2._wp*0.005_wp)))
(1._wp + 0.2_wp*exp(-1._wp*((x_cb(i) - x_centroid)**2._wp + (y_cb(j) - y_centroid)**2._wp)/(2._wp*0.005_wp)))

!x-bump
!q_prim_vf(E_idx)%sf(i, j, 0) = q_prim_vf(E_idx)%sf(i, j, 0)* &
!(1._wp + 0.2_wp*dexp(-1._wp*((x_cb(i) - x_centroid)**2._wp)/(2._wp*0.005_wp)))
!(1._wp + 0.2_wp*exp(-1._wp*((x_cb(i) - x_centroid)**2._wp)/(2._wp*0.005_wp)))

!bump in void fraction
!q_prim_vf(adv_idx%beg)%sf(i,j,0) = q_prim_vf(adv_idx%beg)%sf(i,j,0) * &
Expand Down
6 changes: 3 additions & 3 deletions src/simulation/m_bubbles.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -231,8 +231,8 @@ contains
! Keller-Miksis bubbles
Cpinf = myP
Cpbw = f_cpbw_KM(R0(q), myR, myV, pb)
! c_gas = dsqrt( n_tait*(Cpbw+B_tait) / myRho)
c_liquid = DSQRT(n_tait*(myP + B_tait)/(myRho*(1._wp - alf)))
! c_gas = sqrt( n_tait*(Cpbw+B_tait) / myRho)
c_liquid = SQRT(n_tait*(myP + B_tait)/(myRho*(1._wp - alf)))
rddot = f_rddot_KM(pbdot, Cpinf, Cpbw, myRho, myR, myV, R0(q), c_liquid)
else if (bubble_model == 3) then
! Rayleigh-Plesset bubbles
Expand Down Expand Up @@ -333,7 +333,7 @@ contains
tmp = (fCpinf/(1._wp + fBtait) + 1._wp)**((fntait - 1._wp)/fntait)
tmp = fntait*(1._wp + fBtait)*tmp

f_cgas = DSQRT(tmp + (fntait - 1._wp)*fH)
f_cgas = SQRT(tmp + (fntait - 1._wp)*fH)

end function f_cgas

Expand Down
14 changes: 7 additions & 7 deletions src/simulation/m_data_output.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -902,7 +902,7 @@ contains
nR3 = nR3 + weight(s)*(nR(s)**3._wp)
end do

nbub = DSQRT((4._wp*pi/3._wp)*nR3/alf)
nbub = SQRT((4._wp*pi/3._wp)*nR3/alf)

#ifdef DEBUG
print *, 'In probe, nbub: ', nbub
Expand Down Expand Up @@ -993,7 +993,7 @@ contains
nR3 = nR3 + weight(s)*(nR(s)**3._wp)
end do

nbub = DSQRT((4._wp*pi/3._wp)*nR3/alf)
nbub = SQRT((4._wp*pi/3._wp)*nR3/alf)

R(:) = nR(:)/nbub
Rdot(:) = nRdot(:)/nbub
Expand Down Expand Up @@ -1239,7 +1239,7 @@ contains
int_pres = int_pres + (pres - 1._wp)**2._wp
end if
end do
int_pres = dsqrt(int_pres/(1._wp*npts))
int_pres = sqrt(int_pres/(1._wp*npts))

if (num_procs > 1) then
tmp = int_pres
Expand Down Expand Up @@ -1272,16 +1272,16 @@ contains
trigger = .false.
if (i == 1) then
!inner portion
if (dsqrt(x_cb(j)**2._wp + y_cb(k)**2._wp) < (rad - 0.5_wp*thickness)) &
if (sqrt(x_cb(j)**2._wp + y_cb(k)**2._wp) < (rad - 0.5_wp*thickness)) &
trigger = .true.
elseif (i == 2) then
!net region
if (dsqrt(x_cb(j)**2._wp + y_cb(k)**2._wp) > (rad - 0.5_wp*thickness) .and. &
dsqrt(x_cb(j)**2._wp + y_cb(k)**2._wp) < (rad + 0.5_wp*thickness)) &
if (sqrt(x_cb(j)**2._wp + y_cb(k)**2._wp) > (rad - 0.5_wp*thickness) .and. &
sqrt(x_cb(j)**2._wp + y_cb(k)**2._wp) < (rad + 0.5_wp*thickness)) &
trigger = .true.
elseif (i == 3) then
!everything else
if (dsqrt(x_cb(j)**2._wp + y_cb(k)**2._wp) > (rad + 0.5_wp*thickness)) &
if (sqrt(x_cb(j)**2._wp + y_cb(k)**2._wp) > (rad + 0.5_wp*thickness)) &
trigger = .true.
end if

Expand Down
28 changes: 14 additions & 14 deletions src/simulation/m_global_parameters.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -874,7 +874,7 @@ contains

temp = 293.15_wp
D_m = (0.242_wp * (10._wp ** -(4)))
uu = DSQRT(pl0/rhol0)
uu = SQRT(pl0/rhol0)

omega_ref = 3._wp*k_poly*Ca + 2._wp*(3._wp*k_poly - 1._wp)/Web

Expand All @@ -883,10 +883,10 @@ contains
R_n = Ru/M_n
R_v = Ru/M_v
! phi_vn & phi_nv (phi_nn = phi_vv = 1)
phi_vn = (1._wp + DSQRT(mu_v/mu_n)*(M_n/M_v)**(0.25_wp))**2 &
/(DSQRT(8._wp)*DSQRT(1._wp + M_v/M_n))
phi_nv = (1._wp + DSQRT(mu_n/mu_v)*(M_v/M_n)**(0.25_wp))**2 &
/(DSQRT(8._wp)*DSQRT(1._wp + M_n/M_v))
phi_vn = (1._wp + SQRT(mu_v/mu_n)*(M_n/M_v)**(0.25_wp))**2 &
/(SQRT(8._wp)*SQRT(1._wp + M_v/M_n))
phi_nv = (1._wp + SQRT(mu_n/mu_v)*(M_v/M_n)**(0.25_wp))**2 &
/(SQRT(8._wp)*SQRT(1._wp + M_n/M_v))
! internal bubble pressure
pb0 = pl0 + 2._wp*ss/(R0ref*R0)

Expand Down Expand Up @@ -922,7 +922,7 @@ contains
! keeps a constant (cold liquid assumption)
Tw = 1._wp
! natural frequencies
omegaN = DSQRT(3._wp*k_poly*Ca + 2._wp*(3._wp*k_poly - 1._wp)/(Web*R0))/R0
omegaN = SQRT(3._wp*k_poly*Ca + 2._wp*(3._wp*k_poly - 1._wp)/(Web*R0))/R0

pl0 = 1._wp
do ir = 1, Nb
Expand Down Expand Up @@ -1053,30 +1053,30 @@ contains
!R0mx = 150.D0

sd = poly_sigma
R0mn = 0.8_wp*DEXP(-2.8_wp*sd)
R0mx = 0.2_wp*DEXP(9.5_wp*sd) + 1._wp
R0mn = 0.8_wp*EXP(-2.8_wp*sd)
R0mx = 0.2_wp*EXP(9.5_wp*sd) + 1._wp

! phi = ln( R0 ) & return R0
do ir = 1, nb
phi(ir) = DLOG(R0mn) &
+ dble(ir - 1)*DLOG(R0mx/R0mn)/dble(nb - 1)
R0(ir) = DEXP(phi(ir))
phi(ir) = LOG(R0mn) &
+ dble(ir - 1)*LOG(R0mx/R0mn)/dble(nb - 1)
R0(ir) = EXP(phi(ir))
end do
dphi = phi(2) - phi(1)

! weights for quadrature using Simpson's rule
do ir = 2, nb - 1
! Gaussian
tmp = DEXP(-0.5_wp*(phi(ir)/sd)**2)/DSQRT(2._wp*pi)/sd
tmp = EXP(-0.5_wp*(phi(ir)/sd)**2)/SQRT(2._wp*pi)/sd
if (mod(ir, 2) == 0) then
weight(ir) = tmp*4._wp*dphi/3._wp
else
weight(ir) = tmp*2._wp*dphi/3._wp
end if
end do
tmp = DEXP(-0.5_wp*(phi(1)/sd)**2)/DSQRT(2._wp*pi)/sd
tmp = EXP(-0.5_wp*(phi(1)/sd)**2)/SQRT(2._wp*pi)/sd
weight(1) = tmp*dphi/3._wp
tmp = DEXP(-0.5_wp*(phi(nb)/sd)**2)/DSQRT(2._wp*pi)/sd
tmp = EXP(-0.5_wp*(phi(nb)/sd)**2)/SQRT(2._wp*pi)/sd
weight(nb) = tmp*dphi/3._wp
end subroutine s_simpson
Expand Down
Loading

0 comments on commit 6d73bf8

Please sign in to comment.