diff --git a/src/Makefile.am b/src/Makefile.am index 90ac9873..69c0278e 100755 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -41,6 +41,7 @@ util/util_traverse_pop.f90\ util/util_destroy_list.f90\ util/util_init_list.f90\ util/util_perlin_noise.f90\ +util/util_random_number_normal.f90\ io/io_read_const.f90\ io/io_get_token.f90\ io/io_input.f90\ @@ -126,5 +127,6 @@ regolith/regolith_shock_damage_zone.f90\ regolith/regolith_shock_damage.f90\ porosity/porosity_form_interior.f90\ realistic/realistic_perlin_noise.f90\ +realistic/realistic_crater_topography.f90\ main/CTEM.f90 CLEANFILES = *.mod diff --git a/src/crater/crater_dimensions.f90 b/src/crater/crater_dimensions.f90 index 78fd064b..2fe68690 100644 --- a/src/crater/crater_dimensions.f90 +++ b/src/crater/crater_dimensions.f90 @@ -71,6 +71,7 @@ subroutine crater_dimensions(user,crater,domain) crater%ejrim = 0.14_DP * (crater%fcrat * 0.5_DP)**(0.74_DP) ! McGetchin et al. (1973) Thickness of ejecta at rim crater%ejrad = max(crater_profile_find_r_inner_wall(user,crater) * crater%frad, crater%rad) + ! print *,'in _dimension.f90',crater%rimheight,crater%ejrim !find rim for counting purposes crater%frim = RIMFAC * crater%frad diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index 7ed11620..7729c868 100644 --- a/src/crater/crater_populate.f90 +++ b/src/crater/crater_populate.f90 @@ -24,6 +24,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt use module_io use module_ejecta use module_util + use module_realistic !use module_crust use module_regolith ! simulate regolith mixing use module_crater, EXCEPT_THIS_ONE => crater_populate @@ -261,7 +262,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt ejtble = 0 end if - if (user%dorealistic) call crater_realistic_topography(user,surf,crater,domain,ejecta_dem) + if (user%dorealistic) call realistic_crater_topography(user,surf,crater,domain,ejecta_dem) deallocate(ejecta_dem) ! Collapse any remaining unstable slopes diff --git a/src/crater/crater_realistic_topography.f90 b/src/crater/crater_realistic_topography.f90 index 07a1cf5e..9ce22cc0 100644 --- a/src/crater/crater_realistic_topography.f90 +++ b/src/crater/crater_realistic_topography.f90 @@ -85,6 +85,18 @@ subroutine complex_terrace(user,surf,crater,deltaMtot) real(DP),intent(inout) :: deltaMtot end subroutine complex_terrace + + ! subroutine realistic_rim(user,surf,crater,deltaMtot) + ! use module_globals + ! use module_util + ! use module_crater + ! implicit none + ! type(usertype),intent(in) :: user + ! type(surftype),dimension(:,:),intent(inout) :: surf + ! type(cratertype),intent(inout) :: crater + ! real(DP),intent(inout) :: deltaMtot + ! end subroutine realistic_rim + subroutine ejecta_texture(user,surf,crater,deltaMtot,inc,ejecta_dem) use module_globals implicit none @@ -111,7 +123,9 @@ end subroutine crater_realistic_slope_texture deltaMtot = 0.0_DP select case(crater%morphtype) case("COMPLEX","PEAKRING","MULTIRING") - call complex_terrace(user,surf,crater,deltaMtot) + + !call realistic_rim(user,surf,crater,deltaMtot) + !call complex_terrace(user,surf,crater,deltaMtot) !call complex_wall_texture(user,surf,crater,domain,deltaMtot) !call complex_floor(user,surf,crater,deltaMtot) !call complex_peak(user,surf,crater,deltaMtot) @@ -801,3 +815,328 @@ subroutine crater_realistic_slope_texture(user,critical_value,inc,critarray) return end subroutine crater_realistic_slope_texture + + + + + + + + + + + + +! new subroutines added by jundu on 9/28/2022 + + + + +! subroutine realistic_rim(user,surf,crater,deltaMtot) +! !Makes terraced walls by applying noisy topographic diffusion in discrete radial zones along the wall +! use module_globals +! use module_util +! use module_crater +! implicit none + +! ! in and out +! type(usertype),intent(in) :: user +! type(surftype),dimension(:,:),intent(inout) :: surf +! type(cratertype),intent(inout) :: crater +! real(DP),intent(inout) :: deltaMtot + +! ! internal +! type(psdtype) :: psd_distance,psd_elevation ! PSDs of rim distance and elevation +! real(DP) :: newdem,elchange,rad_roi,xbar,ybar +! integer(I4B) :: xpi,ypi,i,j,inc +! real(DP) :: arc_length,radial_distance +! real(DP) :: rim_distance,rim_distance_delta,rim_elevation,rim_elevation_delta ! rim distance and elevation and their variations. PSD only characterizes the variation. +! real(DP) :: diameter_in_m,gridsize +! real(DP) :: rim_elevation_old +! integer(I4B),parameter :: infile = 14 +! integer(I4B),parameter :: outfile = 15 +! real(DP),dimension(:), allocatable :: amplitude_distance,wavelength_distance,phase_distance +! real(DP),dimension(:), allocatable :: amplitude_elevation,wavelength_elevation,phase_elevation + +! ! executable +! gridsize=user%pix +! diameter_in_m=crater%fcrat +! rad_roi = 2.0_DP * crater%frad +! inc = max(min(nint(rad_roi / user%pix),PBCLIM*user%gridsize),1) + 1 +! crater%maxinc = max(crater%maxinc,inc) + +! ! PSDs of rim distance +! psd_distance%diameter_in_km=diameter_in_m/1000 +! psd_distance%num_vertices=5000 +! psd_distance%num_sine=psd_distance%num_vertices/2 +! psd_distance%input_x_max=PI*psd_distance%diameter_in_km +! psd_distance%diameter_in_km_trans=20.0_DP +! psd_distance%bp2_x_k_b_sigma=[0.06294416243654823, -0.23388324873096447, 0.0056111111111111145, 0.9127777777777776, 0.06082223144489516] +! psd_distance%bp2_y_k_b_sigma=[0.17451428571428573, 0.5097142857142856, 0.0015652173913043492, 3.968695652173913, 0.39645422747641196] +! psd_distance%bp3_y_k_b_sigma=[0.13316666666666666, 1.836666666666667, 0.0004552129221732701, 4.4908957415565345, 0.27156105098740746] +! psd_distance%bp4_y_k_b_sigma=[0.200561797752809, -2.4412359550561797, 0.005084745762711865, 1.4683050847457628, 0.349619629586608] +! psd_distance%slope1_k_b_sigma=[0.004916013437849945, 2.9333146696528556, -9999.9, -9999.9,0.36717269233927413] + + +! ! PSDs of rim elevation +! psd_elevation%diameter_in_km=diameter_in_m/1000 +! psd_elevation%num_vertices=5000 +! psd_elevation%num_sine=psd_elevation%num_vertices/2 +! psd_elevation%input_x_max=PI*psd_elevation%diameter_in_km +! psd_elevation%diameter_in_km_trans=20.0_DP +! psd_elevation%bp2_x_k_b_sigma=[0.05533333333333333, -0.10666666666666666, 0.0040595399188092015, 0.918809201623816,0.11430090458471147] +! psd_elevation%bp2_y_k_b_sigma=[0.12594736842105264, -0.2859473684210526, 0.003978349120433018, 2.1534330175913396, 0.40875902323519153] +! psd_elevation%bp3_y_k_b_sigma=[0.0, 4.0, 0.0, 4.0,0.7435034294351548] +! psd_elevation%bp4_y_k_b_sigma(1)=-9999.9 +! psd_elevation%slope1_k_b_sigma=[0.09919786096256683, 1.7070427807486632, -0.008135593220338983, 3.8537118644067796,0.42722290048928596] + +! call Calculate_am_wl_phase_from_diameter(psd_distance,amplitude_distance,wavelength_distance,phase_distance) +! call Calculate_am_wl_phase_from_diameter(psd_elevation,amplitude_elevation,wavelength_elevation,phase_elevation) + + +! rim_elevation_old=(crater%rimheight - crater%ejrim) + +! do j = -inc,inc +! do i = -inc,inc +! xpi = crater%xlpx + i +! ypi = crater%ylpx + j +! xbar = xpi * user%pix - crater%xl +! ybar = ypi * user%pix - crater%yl + +! arc_length=(atan2(real(j,8),real(i,8))+PI)*diameter_in_m/2.0_DP/1000.0_DP +! radial_distance= sqrt(xbar**2 + ybar**2) + +! call Create_rim(arc_length,psd_distance,amplitude_distance,wavelength_distance,phase_distance,rim_distance_delta) +! call Create_rim(arc_length,psd_elevation,amplitude_elevation,wavelength_elevation,phase_elevation,rim_elevation_delta) + +! rim_distance=rim_distance_delta+diameter_in_m/2.0_DP +! rim_elevation=rim_elevation_delta+rim_elevation_old +! crater%rimheight=rim_elevation + +! if (radial_distance>rim_distance) then +! newdem=crater_profile(user,crater,radial_distance/rim_distance) +! elchange = newdem-surf(xpi,ypi)%dem +! deltaMtot = deltaMtot + elchange +! surf(xpi,ypi)%dem = newdem +! endif + +! if (radial_distancecrater%floordiam/2.0_DP) then +! newdem=crater_profile(user,crater,radial_distance/rim_distance) +! elchange = newdem-surf(xpi,ypi)%dem +! deltaMtot = deltaMtot + elchange +! surf(xpi,ypi)%dem = newdem +! endif + +! end do +! end do + + +! end subroutine realistic_rim + + ! Universal shapde model from PSD + +! subroutine Calculate_am_wl_phase_from_diameter(psd_1D,amplitude,wavelength,phase) +! use module_globals +! use module_crater +! implicit none +! ! in and out +! type(psdtype),intent(inout) :: psd_1D +! real(DP),dimension(:), allocatable,intent(out) :: amplitude,wavelength,phase +! ! internal +! integer(I4B) :: i +! real(DP) :: phase_random,random_number_normal +! real(DP),dimension(:), allocatable :: psd +! integer(I4B),parameter :: infile = 14 +! integer(I4B),parameter :: outfile = 15 + + +! ! excutable + + +! call Calculate_breakpoint_slope_from_diameter(psd_1D) +! call Calculate_targetPSD_from_breakpoint_slope(psd_1D,wavelength,psd) +! call Calculate_am_wl_phase_from_targetPSD(psd_1D,wavelength,psd,amplitude,phase) + + +! open(outfile,file="output_psd.txt",status='replace') +! do i = 1, size(wavelength) +! write(outfile,'(ES18.10,1(",",ES18.10))') wavelength(i) ,psd(i) +! end do +! close(outfile) + +! end subroutine Calculate_am_wl_phase_from_diameter + +! subroutine Calculate_breakpoint_slope_from_diameter(psd_1D) +! use module_globals +! use module_crater +! implicit none +! ! in and out +! type(psdtype),intent(inout) :: psd_1D + +! ! excutable + +! if ( psd_1D%diameter_in_km < psd_1D%diameter_in_km_trans ) then + +! psd_1D%bp2_x=psd_1D%bp2_x_k_b_sigma(1)*psd_1D%diameter_in_km+psd_1D%bp2_x_k_b_sigma(2) +! psd_1D%bp2_y=psd_1D%bp2_y_k_b_sigma(1)*psd_1D%diameter_in_km+psd_1D%bp2_y_k_b_sigma(2) +! psd_1D%bp3_y=psd_1D%bp3_y_k_b_sigma(1)*psd_1D%diameter_in_km+psd_1D%bp3_y_k_b_sigma(2) +! psd_1D%bp4_y=psd_1D%bp4_y_k_b_sigma(1)*psd_1D%diameter_in_km+psd_1D%bp4_y_k_b_sigma(2) + +! else + +! psd_1D%bp2_x=psd_1D%bp2_x_k_b_sigma(3)*psd_1D%diameter_in_km+psd_1D%bp2_x_k_b_sigma(4) +! psd_1D%bp2_y=psd_1D%bp2_y_k_b_sigma(3)*psd_1D%diameter_in_km+psd_1D%bp2_y_k_b_sigma(4) +! psd_1D%bp3_y=psd_1D%bp3_y_k_b_sigma(3)*psd_1D%diameter_in_km+psd_1D%bp3_y_k_b_sigma(4) +! psd_1D%bp4_y=psd_1D%bp4_y_k_b_sigma(3)*psd_1D%diameter_in_km+psd_1D%bp4_y_k_b_sigma(4) + +! endif + + +! if (psd_1D%slope1_k_b_sigma(3)>-100.0_DP) then +! if ( psd_1D%diameter_in_km < psd_1D%diameter_in_km_trans ) then +! psd_1D%slope1=psd_1D%slope1_k_b_sigma(1)*psd_1D%diameter_in_km+psd_1D%slope1_k_b_sigma(2) +! else +! psd_1D%slope1=psd_1D%slope1_k_b_sigma(3)*psd_1D%diameter_in_km+psd_1D%slope1_k_b_sigma(4) +! endif +! else +! psd_1D%slope1=psd_1D%slope1_k_b_sigma(1)*psd_1D%diameter_in_km+psd_1D%slope1_k_b_sigma(2) +! endif + + +! end subroutine Calculate_breakpoint_slope_from_diameter + +! subroutine Calculate_targetPSD_from_breakpoint_slope(psd_1D,wavelength,psd) +! use module_globals +! use module_crater +! implicit none +! !in and out +! type(psdtype), intent(in) :: psd_1D +! real(DP) ,dimension(:),allocatable,intent(out) :: wavelength,psd +! ! internal +! integer(I4B) :: i,bp2_x_index,bp3_x_index +! real(DP) :: bp3_x,bp4_x +! real(DP) :: slope_12,intercept_12,slope_23,intercept_23,slope_34,intercept_34 +! real(DP) :: random_number_normal + + +! ! excutable + +! allocate(wavelength(psd_1D%num_sine)) +! allocate(psd(psd_1D%num_sine)) +! !---------------------------------------------------------------------------------------------------------------------- +! do i = 1, psd_1D%num_sine +! wavelength(i) = psd_1D%input_x_max / (psd_1D%num_sine-i+1) +! end do +! !---------------------------------------------------------------------------------------------------------------------- + +! if (psd_1D%bp4_y_k_b_sigma(1)>-100.0_DP ) then + +! ! 3 segments, 4 breakpoints +! bp3_x=log10(psd_1D%input_x_max/2.0_DP) +! bp4_x=log10(psd_1D%input_x_max) + +! do i = 1, psd_1D%num_sine +! bp2_x_index = i +! if (wavelength(i) > 10**psd_1D%bp2_x) exit +! end do +! do i = 1, psd_1D%num_sine +! bp3_x_index = i +! if (wavelength(i) > 10**bp3_x) exit +! end do +! !---------------------------------------------------------------------------------------------------------------------- +! slope_12=psd_1D%slope1 +! intercept_12 = psd_1D%bp2_y - slope_12 * psd_1D%bp2_x + +! slope_23=(psd_1D%bp3_y-psd_1D%bp2_y)/(bp3_x-psd_1D%bp2_x) +! intercept_23 = psd_1D%bp3_y - slope_23 * bp3_x + +! slope_34=(psd_1D%bp4_y-psd_1D%bp3_y)/(bp4_x-bp3_x) +! intercept_34 = psd_1D%bp4_y - slope_34 * bp4_x +! !---------------------------------------------------------------------------------------------------------------------- +! do i = 1, bp2_x_index-1 +! psd(i) = 10** ( slope_12*log10(wavelength(i)) + intercept_12 ) +! end do +! do i = bp2_x_index, bp3_x_index-1 +! psd(i) = 10** ( slope_23*log10(wavelength(i)) + intercept_23 ) +! end do +! do i = bp3_x_index, psd_1D%num_sine +! psd(i) = 10** ( slope_34*log10(wavelength(i)) + intercept_34 ) +! end do + +! else + +! ! 2 segments, 3 breakpoints +! bp3_x=log10(psd_1D%input_x_max) + +! do i = 1, psd_1D%num_sine +! bp2_x_index = i +! if (wavelength(i) > 10**psd_1D%bp2_x) exit +! end do +! !---------------------------------------------------------------------------------------------------------------------- +! slope_12=psd_1D%slope1 +! intercept_12 = psd_1D%bp2_y - slope_12 * psd_1D%bp2_x + +! slope_23=(psd_1D%bp3_y-psd_1D%bp2_y)/(bp3_x-psd_1D%bp2_x) +! intercept_23 = psd_1D%bp3_y - slope_23 * bp3_x +! !---------------------------------------------------------------------------------------------------------------------- +! do i = 1, bp2_x_index-1 +! psd(i) = 10** ( slope_12*log10(wavelength(i)) + intercept_12 ) +! end do +! do i = bp2_x_index, psd_1D%num_sine +! psd(i) = 10** ( slope_23*log10(wavelength(i)) + intercept_23 ) +! end do +! endif +! end subroutine Calculate_targetPSD_from_breakpoint_slope + +! subroutine Calculate_am_wl_phase_from_targetPSD(psd_1D,wavelength,psd,amplitude,phase) +! use module_globals +! use module_crater +! implicit none +! !in and out +! type(psdtype), intent(in) :: psd_1D +! real(DP) ,dimension(:),intent(in) :: wavelength,psd +! real(DP) ,dimension(:),allocatable,intent(out) :: amplitude,phase +! ! internal +! integer(I4B) :: i +! real(DP) :: phase_random +! ! excutable +! allocate(amplitude(psd_1D%num_sine)) +! allocate(phase(psd_1D%num_sine)) + +! !---------------------------------------------------------------------------------------------------------------------- +! do i = 1, psd_1D%num_sine +! amplitude(i)= sqrt(psd(i) * psd_1D%input_x_max / (psd_1D%num_vertices ** 2)) +! call RANDOM_NUMBER(phase_random) +! phase(i)=wavelength(i) *phase_random +! end do + +! end subroutine Calculate_am_wl_phase_from_targetPSD + +! subroutine Create_rim(arc_length,psd_1D,amplitude,wavelength,phase,rim_parameter) + +! use module_globals +! implicit none +! ! in and out +! real(DP),intent(in) :: arc_length +! type(psdtype),intent(in) :: psd_1D +! real(DP),dimension(:),intent(in) :: amplitude +! real(DP),dimension(:),intent(in) :: wavelength +! real(DP),dimension(:),intent(in) :: phase +! real(DP),intent(out) :: rim_parameter +! ! internal +! integer(I4B) :: i +! real(DP) :: rim_parameter_ind + +! ! excutable +! rim_parameter = 0.0_DP +! do i = 1, psd_1D%num_sine +! rim_parameter_ind= amplitude(i) * sin( 2 * PI* 1/wavelength(i) * ( arc_length- phase(i) ) ) +! rim_parameter=rim_parameter+rim_parameter_ind +! end do +! rim_parameter=rim_parameter*1000.0_DP + +! end subroutine Create_rim + + + + \ No newline at end of file diff --git a/src/crater/module_crater.f90 b/src/crater/module_crater.f90 index b946c25b..a6be21fc 100644 --- a/src/crater/module_crater.f90 +++ b/src/crater/module_crater.f90 @@ -339,4 +339,62 @@ subroutine crater_superdomain(user,surf,age,age_resolution,prod,nflux,domain,fin end subroutine crater_superdomain end interface + + + +! new subroutines added by jundu on 10/15/2022 + + +! Rim_crest: + + ! interface + + ! subroutine Calculate_am_wl_phase_from_diameter(psd_1D,amplitude,wavelength,phase) + ! use module_globals + ! implicit none + ! ! in and out + ! type(psdtype),intent(inout) :: psd_1D + ! real(DP),dimension(:), allocatable,intent(out) :: amplitude,wavelength,phase + ! end subroutine Calculate_am_wl_phase_from_diameter + + ! subroutine Calculate_breakpoint_slope_from_diameter(psd_1D) + ! use module_globals + ! implicit none + ! ! in and out + ! type(psdtype),intent(inout) :: psd_1D + ! end subroutine Calculate_breakpoint_slope_from_diameter + + ! subroutine Calculate_targetPSD_from_breakpoint_slope(psd_1D,wavelength,psd) + ! use module_globals + ! implicit none + ! !in and out + ! type(psdtype), intent(in) :: psd_1D + ! real(DP) ,dimension(:),allocatable,intent(out) :: wavelength,psd + ! end subroutine Calculate_targetPSD_from_breakpoint_slope + + ! subroutine Calculate_am_wl_phase_from_targetPSD(psd_1D,wavelength,psd,amplitude,phase) + ! use module_globals + ! implicit none + ! !in and out + ! type(psdtype), intent(in) :: psd_1D + ! real(DP) ,dimension(:),intent(in) :: wavelength,psd + ! real(DP) ,dimension(:),allocatable,intent(out) :: amplitude,phase + ! end subroutine Calculate_am_wl_phase_from_targetPSD + + + ! subroutine Create_rim(arc_length,psd_1D,amplitude,wavelength,phase,rim_parameter) + ! use module_globals + ! implicit none + ! ! in and out + ! real(DP),intent(in) :: arc_length + ! type(psdtype), intent(in) :: psd_1D + ! real(DP),dimension(:),intent(in) :: amplitude + ! real(DP),dimension(:),intent(in) :: wavelength + ! real(DP),dimension(:),intent(in) :: phase + ! real(DP),intent(out) :: rim_parameter + ! end subroutine Create_rim + + ! end interface + + end module diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index 8013b051..c2006765 100644 --- a/src/globals/module_globals.f90 +++ b/src/globals/module_globals.f90 @@ -263,6 +263,32 @@ module module_globals ! real(SP) :: bedrock ! Fraction of bedrock contained in mixture end type ejbtype + + +! added by jundu on 10/17/2022 +!---------------------------------------------------------- +! power spectral density +type psdtype + real(DP) :: diameter_in_km ! diameter in km + integer(I4B) :: num_vertices ! how many vertices along the linear feature? + integer(I4B) :: num_sine ! how many sine waves needed at a given vertex? + real(DP) :: input_x_max ! period + real(DP) :: diameter_in_km_trans + real(DP) :: bp2_x + real(DP) :: bp2_y + real(DP) :: bp3_y + real(DP) :: bp4_y + real(DP) :: slope1 + real(DP),dimension(5) :: bp2_x_k_b_sigma + real(DP),dimension(5) :: bp2_y_k_b_sigma + real(DP),dimension(5) :: bp3_y_k_b_sigma + real(DP),dimension(5) :: bp4_y_k_b_sigma + real(DP),dimension(5) :: slope1_k_b_sigma + real(DP),dimension(5) :: misfit_k_b_sigma +end type psdtype +!---------------------------------------------------------- + + ! Progress bar variables integer(I4B),parameter :: PBARRES = 100 integer(I4B),parameter :: PBARSIZE = 50 diff --git a/src/realistic/module_realistic.f90 b/src/realistic/module_realistic.f90 index a5c9ace1..2af1c464 100644 --- a/src/realistic/module_realistic.f90 +++ b/src/realistic/module_realistic.f90 @@ -1,4 +1,5 @@ module module_realistic +use module_globals implicit none public save @@ -106,5 +107,100 @@ subroutine realistic_ejecta_texture(user,surf,crater,deltaMtot,inc,ejecta_dem) end subroutine realistic_ejecta_texture end interface + + +! new interfaces added by jundu on 10/22/2022 + +interface + subroutine realistic_crater_topography(user,surf,crater,domain,ejecta_dem) + use module_globals + use module_util + use module_crater + implicit none + ! in and out + type(usertype),intent(in) :: user + type(surftype),dimension(:,:),intent(inout) :: surf + type(cratertype),intent(inout) :: crater + type(domaintype),intent(in) :: domain + real(DP),dimension(:,:),intent(inout) :: ejecta_dem + end subroutine realistic_crater_topography +end interface + + +interface + subroutine realistic_rim(user,surf,crater,deltaMtot) + use module_globals + use module_util + use module_crater + implicit none + ! in and out + type(usertype),intent(in) :: user + type(surftype),dimension(:,:),intent(inout) :: surf + type(cratertype),intent(inout) :: crater + real(DP),intent(inout) :: deltaMtot + end subroutine realistic_rim +end interface + + + +! Rim_crest: +interface + + subroutine Calculate_am_wl_phase_from_diameter(psd_1D,amplitude,wavelength,phase) + use module_globals + implicit none + ! in and out + type(psdtype),intent(inout) :: psd_1D + real(DP),dimension(:), allocatable,intent(out) :: amplitude,wavelength,phase + end subroutine Calculate_am_wl_phase_from_diameter + + subroutine Calculate_breakpoint_slope_from_diameter(psd_1D) + use module_globals + use module_util + implicit none + ! in and out + type(psdtype),intent(inout) :: psd_1D + end subroutine Calculate_breakpoint_slope_from_diameter + + subroutine Calculate_targetPSD_from_breakpoint_slope(psd_1D,wavelength,psd) + use module_globals + use module_util + implicit none + !in and out + type(psdtype), intent(in) :: psd_1D + real(DP) ,dimension(:),allocatable,intent(out) :: wavelength,psd + end subroutine Calculate_targetPSD_from_breakpoint_slope + + subroutine Calculate_am_wl_phase_from_targetPSD(psd_1D,wavelength,psd,amplitude,phase) + use module_globals + implicit none + !in and out + type(psdtype), intent(in) :: psd_1D + real(DP) ,dimension(:),intent(in) :: wavelength,psd + real(DP) ,dimension(:),allocatable,intent(out) :: amplitude,phase + end subroutine Calculate_am_wl_phase_from_targetPSD + + subroutine Create_rim(arc_length,psd_1D,amplitude,wavelength,phase,rim_parameter) + use module_globals + implicit none + ! in and out + real(DP),intent(in) :: arc_length + type(psdtype), intent(in) :: psd_1D + real(DP),dimension(:),intent(in) :: amplitude + real(DP),dimension(:),intent(in) :: wavelength + real(DP),dimension(:),intent(in) :: phase + real(DP),intent(out) :: rim_parameter + end subroutine Create_rim + +end interface + + + + + + + + + end module module_realistic diff --git a/src/util/module_util.f90 b/src/util/module_util.f90 index eb0236e8..46e47e21 100644 --- a/src/util/module_util.f90 +++ b/src/util/module_util.f90 @@ -235,5 +235,33 @@ function util_perlin_noise(x,y,z) result (noise) end function util_perlin_noise end interface + +! added by jundu on 10/25/2022 +! generate random number with a normal distribution + +interface + subroutine util_random_number_uniform(u) + use module_globals + implicit none + real(DP),intent(out) :: u + end subroutine util_random_number_uniform + subroutine util_random_number_normal(x) + use module_globals + implicit none + real(DP),intent(out) :: x + end subroutine util_random_number_normal +end interface + + + + + + + + + + + + end module