crater rim outline and elevation

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

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
 

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

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_distance<rim_distance .and. radial_distance>crater%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
+
+
+
+