Updates to rename some variables and make some small adjustments to dimension parameters. Working to build a more robust initial crater morphology model.

src/crater/crater_emplace.f90

@@ -101,7 +101,7 @@ subroutine crater_emplace(user,surf,crater,domain,deltaMtot)
          lradsq = x_relative**2 + y_relative**2
 
          if (lradsq > crater%frad**2) cycle
-         call crater_form_interior(user,surf(xpi,ypi),crater,x_relative, y_relative,newelev,deltaMi)
+         call crater_form_interior(user,surf(xpi,ypi),crater,x_relative,y_relative,newelev,deltaMi)
          deltaMtot = deltaMtot + deltaMi
 
          ! do porosity computation if (user%doporosity)

src/crater/crater_form_exterior_func.f90

@@ -39,7 +39,7 @@ function crater_form_exterior_func(user,surf,crater,domain,rd,deltaMtot,lastloop
    type(surftype) :: surfi
 
    ! Now make the exterior of the crater
-   inc = int(crater%frad/user%pix*(domain%small/crater%rheight)**(-1._DP/RIMDROP)) !  Maximum distance of crater form
+   inc = int(crater%frad/user%pix*(domain%small/crater%rimheight)**(-1._DP/RIMDROP)) !  Maximum distance of crater form
    inc = max(min(max(crater%rimdispx,inc),PBCLIM*user%gridsize),1)
 
    crater%maxinc = max(crater%maxinc,inc)

src/crater/crater_populate.f90

@@ -164,7 +164,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
          ! Stamp the current time onto the crater
          crater%timestamp = real(curyear + real(icrater,kind=DP) / real(ntotcrat,kind=DP) * user%interval,kind=SP)
          ! Find the visible crater parameters
-         call crater_find_visible(user,crater,domain)
+         call crater_dimensions(user,crater,domain)
 
          ! Crater is big enough to keep, so record it into the true distribution 
          ntrue = ntrue + 1

src/crater/crater_realistic_topography.f90

@@ -39,60 +39,58 @@ subroutine crater_realistic_topography(user,surf,crater,domain,deltaMtot)
    type(surftype),dimension(:,:),intent(inout) :: surf
    type(cratertype),intent(inout) :: crater
    type(domaintype),intent(in) :: domain
-   real(DP),intent(out) :: deltaMtot
+   real(DP),intent(inout) :: deltaMtot
 
    ! Internal variables
-   real(DP) :: lradsq,newelev, x_relative, y_relative 
+   real(DP) :: lradsq,x_relative, y_relative 
    integer(I4B) :: xpi,ypi,i,j,inc,incsq,iradsq
-   real(DP) :: xp,yp,fradsq,deltaMi,rimheight
+   real(DP) :: r,xp,yp,fradsq,deltaMi,rimheight
    logical :: lastloop
    real(DP), dimension(2) :: rn
 
+   interface
+      subroutine complex_floor(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi)
+      use module_globals
+      implicit none
+      type(usertype),intent(in) :: user
+      type(surftype),intent(inout) :: surfi
+      type(cratertype),intent(in) :: crater
+      real(DP),intent(in) :: r,x_relative, y_relative 
+      real(DP),dimension(:),intent(in) :: rn
+      real(DP),intent(out) :: deltaMi
+      end subroutine complex_floor
 
-   ! Topographic noise parameters
-   real(DP) :: xynoise, znoise
-   integer(I4B) :: octave
-   real(DP) :: noise
-
-
-   ! Complex crater floor
-   integer(I4B), parameter :: complex_floor_num_octaves  = 8   ! Number of Perlin noise octaves
-   integer(I4B), parameter :: complex_floor_offset = 10000 ! Scales the random xy-offset so that each crater's random noise is unique 
-   real(DP), parameter :: complex_floor_xy_noise_fac = 2.0_DP  ! Spatial "size" of noise features at the first octave
-   real(DP), parameter :: complex_floor_noise_height = 1e-3_DP ! Vertical height of noise features as a function of crater radius at the first octave
-   real(DP), parameter :: complex_floor_freq = 2.0_DP     ! Spatial size scale factor multiplier at each octave level
-   real(DP), parameter :: complex_floor_pers = 0.80_DP  ! The relative size scaling at each octave level
-
-   ! Complex crater peak
-   integer(I4B), parameter :: complex_peak_num_octaves  = 8   ! Number of Perlin noise octaves
-   integer(I4B), parameter :: complex_peak_offset = 1000 ! Scales the random xy-offset so that each crater's random noise is unique 
-   real(DP), parameter :: complex_peak_rad = 0.20_DP      ! Radial size of central peak relative to frad
-   real(DP), parameter :: complex_peak_xy_noise_fac = 16.0_DP  ! Spatial "size" of noise features at the first octave
-   real(DP), parameter :: complex_peak_freq = 2.0_DP     ! Spatial size scale factor multiplier at each octave level
-   real(DP), parameter :: complex_peak_pers = 0.80_DP  ! The relative size scaling at each octave level
-
-   ! Complex crater wall
-   integer(I4B), parameter :: complex_wall_num_octaves  = 3   ! Number of Perlin noise octaves
-   integer(I4B), parameter :: complex_wall_offset = 2000 ! Scales the random xy-offset so that each crater's random noise is unique 
-   real(DP), parameter :: complex_wall_xy_noise_fac = 2.0_DP  ! Spatial "size" of noise features at the first octave
-   real(DP), parameter :: complex_wall_freq = 2.0_DP     ! Spatial size scale factor multiplier at each octave level
-   real(DP), parameter :: complex_wall_pers = 0.80_DP  ! The relative size scaling at each octave level
-   real(DP), parameter :: complex_wall_noise_height = 0.10_DP ! Vertical height of noise features as a function of crater radius at the first octave
-   real(DP), parameter :: complex_wall_slope = 25._DP * DEG2RAD
-   integer(I4B), parameter :: complex_wall_terracefac = 16 ! Spatial size factor for terraces relative to crater radius
-   integer(I4B) :: complex_wall_terrace_num
+      subroutine complex_peak(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi)
+      use module_globals
+      implicit none
+      type(usertype),intent(in) :: user
+      type(surftype),intent(inout) :: surfi
+      type(cratertype),intent(in) :: crater
+      real(DP),intent(in) :: r,x_relative, y_relative 
+      real(DP),dimension(:),intent(in) :: rn
+      real(DP),intent(out) :: deltaMi
+      end subroutine complex_peak
 
+      subroutine complex_wall(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi)
+      use module_globals
+      implicit none
+      type(usertype),intent(in) :: user
+      type(surftype),intent(inout) :: surfi
+      type(cratertype),intent(in) :: crater
+      real(DP),intent(in) :: r,x_relative, y_relative 
+      real(DP),dimension(:),intent(in) :: rn
+      real(DP),intent(out) :: deltaMi
+      end subroutine complex_wall
+   end interface
 
    ! Executable code
-
    call random_number(rn)
 
    ! determine area to effect
    ! First make the interior of the crater
    inc = max(min(crater%fradpx,PBCLIM*user%gridsize),1) + 1
    crater%maxinc = max(crater%maxinc,inc)
    fradsq = crater%frad**2
-   deltaMtot = 0.0_DP 
    incsq = inc**2
 
    do j=-inc,inc  ! Do the loop in pixel space
@@ -101,7 +99,6 @@ subroutine crater_realistic_topography(user,surf,crater,domain,deltaMtot)
          ! find distance from crater center
 
          ! find elevation and grid point
-         newelev = crater%melev + ((i * crater%xslp) + (j * crater%yslp)) * user%pix
          xpi = crater%xlpx + i
          ypi = crater%ylpx + j
 
@@ -116,51 +113,178 @@ subroutine crater_realistic_topography(user,surf,crater,domain,deltaMtot)
          lradsq = x_relative**2 + y_relative**2
 
          if (lradsq > crater%frad**2) cycle
+
+         r = sqrt(x_relative**2+y_relative**2) / crater%frad
 
-
-
+         if (crater%fcrat > crater%cxtran * 2) then
+            call complex_floor(user,surf(xpi,ypi),crater,r,x_relative,y_relative,rn,deltaMi)
+            call complex_peak (user,surf(xpi,ypi),crater,r,x_relative,y_relative,rn,deltaMi)
+            call complex_wall (user,surf(xpi,ypi),crater,r,x_relative,y_relative,rn,deltaMi)
+         end if
 
          deltaMtot = deltaMtot + deltaMi
 
       end do
    end do
 
 
-!            ! Make the floor topographically "noisy"
-!            noise = 0.0_DP
-!            do octave = 1, complex_floor_num_octaves
-!               xynoise = complex_floor_xy_noise_fac * complex_floor_freq ** octave / crater%fcrat
-!               znoise = complex_floor_noise_height * complex_floor_pers ** (octave - 1) * crater%fcrat
-!               noise = noise + util_perlin_noise(xynoise * x_relative + complex_floor_offset * rn(1), &
-!                                                 xynoise * y_relative + complex_floor_offset * rn(2)) * znoise
-!            end do
-!            newdem = newdem + max(noise,0.0_DP)
-!
-!            ! Add in "noisy" central peak
-!            if (r < complex_peak_rad) then
-!               Hpeak = 0.032e3_DP * (crater%fcrat * 1e-3_DP)**(0.900_DP) ! Pike (1977)
-!               noise = 0.0_DP
-!               do octave = 1, complex_peak_num_octaves 
-!                  xynoise = complex_peak_xy_noise_fac * complex_peak_freq ** octave / crater%fcrat
-!                  znoise = Hpeak * (1.0_DP - r / complex_peak_rad) * complex_peak_pers ** (octave - 1)
-!                  noise = noise + util_perlin_noise(xynoise * x_relative + complex_peak_offset * rn(1), &
-!                                                    xynoise * y_relative + complex_peak_offset * rn(2)) * znoise
-!               end do 
-!               newdem = newdem + max(noise,0.0_DP) 
-!            end if
-!         else ! Terraced walls 
-!            complex_wall_terrace_num = int(complex_wall_terracefac * r) + 1
-!            noise = 0.0_DP
-!            do octave = 1, complex_wall_num_octaves
-!               xynoise = complex_wall_xy_noise_fac * complex_wall_freq ** octave / crater%fcrat
-!               znoise = complex_wall_noise_height * complex_wall_pers ** (octave - 1) * crater%fcrat
-!               noise = noise + util_perlin_noise(xynoise * x_relative + complex_wall_terrace_num * complex_wall_offset * rn(1), &
-!                                                 xynoise * y_relative + complex_wall_terrace_num * complex_wall_offset * rn(2)) &
-!                                                 * znoise
-!            end do
-!            newdem = max(newdem + noise,crater%melev - pikeD) !/ cos(complex_wall_slope)
-
-
    return
 end subroutine crater_realistic_topography
 
+subroutine complex_peak(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi)
+   use module_globals
+   use module_util
+   use module_crater, EXCEPT_THIS_ONE => crater_realistic_topography
+   implicit none
+
+   ! Arguments
+   type(usertype),intent(in) :: user
+   type(surftype),intent(inout) :: surfi
+   type(cratertype),intent(in) :: crater
+   real(DP),intent(in) :: r,x_relative, y_relative 
+   real(DP),dimension(:),intent(in) :: rn
+   real(DP),intent(out) :: deltaMi
+
+   ! Complex crater peak
+   real(DP), parameter :: complex_peak_rad = 0.20_DP      ! Radial size of central peak relative to frad
+   integer(I4B), parameter :: complex_peak_num_octaves  = 8   ! Number of Perlin noise octaves
+   integer(I4B), parameter :: complex_peak_offset = 1000 ! Scales the random xy-offset so that each crater's random noise is unique 
+   real(DP), parameter :: complex_peak_xy_noise_fac = 16.0_DP  ! Spatial "size" of noise features at the first octave
+   real(DP), parameter :: complex_peak_freq = 2.0_DP     ! Spatial size scale factor multiplier at each octave level
+   real(DP), parameter :: complex_peak_pers = 0.80_DP  ! The relative size scaling at each octave level
+
+   ! Internal variables
+   real(DP) :: newdem,elchange,pikeD,Hpeak
+
+   ! Topographic noise parameters
+   real(DP) :: xynoise, znoise
+   integer(I4B) :: octave
+   real(DP) :: noise
+
+   if (r > complex_peak_rad) return
+   ! Add in "noisy" central peak
+   newdem = surfi%dem
+   Hpeak = 0.032e3_DP * (crater%fcrat * 1e-3_DP)**(0.900_DP) ! Pike (1977)
+   noise = 0.0_DP
+   do octave = 1, complex_peak_num_octaves 
+      xynoise = complex_peak_xy_noise_fac * complex_peak_freq ** octave / crater%fcrat
+      znoise = Hpeak * (1.0_DP - r / complex_peak_rad) * complex_peak_pers ** (octave - 1)
+      noise = noise + util_perlin_noise(xynoise * x_relative + complex_peak_offset * rn(1), &
+                                        xynoise * y_relative + complex_peak_offset * rn(2)) * znoise
+   end do 
+   newdem = newdem + max(noise,0.0_DP) 
+
+   elchange  = newdem - surfi%dem
+   deltaMi = elchange
+   surfi%dem = newdem
+
+end subroutine complex_peak
+
+subroutine complex_floor(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi)
+   use module_globals
+   use module_util
+   use module_crater, EXCEPT_THIS_ONE => crater_realistic_topography
+   implicit none
+
+   ! Arguments
+   type(usertype),intent(in) :: user
+   type(surftype),intent(inout) :: surfi
+   type(cratertype),intent(in) :: crater
+   real(DP),intent(in) :: r,x_relative, y_relative 
+   real(DP),dimension(:),intent(in) :: rn
+   real(DP),intent(out) :: deltaMi
+
+   ! Complex crater floor parameters
+   integer(I4B), parameter :: complex_floor_num_octaves  = 8   ! Number of Perlin noise octaves
+   integer(I4B), parameter :: complex_floor_offset = 10000 ! Scales the random xy-offset so that each crater's random noise is unique 
+   real(DP), parameter :: complex_floor_xy_noise_fac = 2.0_DP  ! Spatial "size" of noise features at the first octave
+   real(DP), parameter :: complex_floor_noise_height = 1e-3_DP ! Vertical height of noise features as a function of crater radius at the first octave
+   real(DP), parameter :: complex_floor_freq = 2.0_DP     ! Spatial size scale factor multiplier at each octave level
+   real(DP), parameter :: complex_floor_pers = 0.80_DP  ! The relative size scaling at each octave level
+
+   ! Internal variables
+   real(DP) :: newdem,elchange,pikeD,Hpeak
+
+   ! Topographic noise parameters
+   real(DP) :: xynoise, znoise
+   integer(I4B) :: octave
+   real(DP) :: noise
+
+   newdem = surfi%dem
+   ! Make the floor topographically "noisy"
+   noise = 0.0_DP
+   do octave = 1, complex_floor_num_octaves
+      xynoise = complex_floor_xy_noise_fac * complex_floor_freq ** octave / crater%fcrat
+      znoise = complex_floor_noise_height * complex_floor_pers ** (octave - 1) * crater%fcrat
+      noise = noise + util_perlin_noise(xynoise * x_relative + complex_floor_offset * rn(1), &
+                                        xynoise * y_relative + complex_floor_offset * rn(2)) * znoise
+   end do
+   newdem = newdem + max(noise,0.0_DP)
+
+   elchange  = newdem - surfi%dem
+   deltaMi = elchange
+   surfi%dem = newdem
+
+end subroutine complex_floor
+
+
+subroutine complex_wall(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi)
+   use module_globals
+   use module_util
+   use module_crater, EXCEPT_THIS_ONE => crater_realistic_topography
+   implicit none
+
+   ! Arguments
+   type(usertype),intent(in) :: user
+   type(surftype),intent(inout) :: surfi
+   type(cratertype),intent(in) :: crater
+   real(DP),intent(in) :: r,x_relative, y_relative 
+   real(DP),dimension(:),intent(in) :: rn
+   real(DP),intent(out) :: deltaMi
+
+   ! Internal variables
+   real(DP) :: newdem,elchange,pikeD,pikeFloor,Hpeak
+
+   ! Topographic noise parameters
+   real(DP) :: xynoise, znoise
+   integer(I4B) :: octave
+   real(DP) :: noise
+
+   ! Complex crater wall
+   integer(I4B), parameter :: complex_wall_num_octaves  = 3   ! Number of Perlin noise octaves
+   integer(I4B), parameter :: complex_wall_offset = 2000 ! Scales the random xy-offset so that each crater's random noise is unique 
+   real(DP), parameter :: complex_wall_xy_noise_fac = 2.0_DP  ! Spatial "size" of noise features at the first octave
+   real(DP), parameter :: complex_wall_freq = 2.0_DP     ! Spatial size scale factor multiplier at each octave level
+   real(DP), parameter :: complex_wall_pers = 0.80_DP  ! The relative size scaling at each octave level
+   real(DP), parameter :: complex_wall_noise_height = 0.10_DP ! Vertical height of noise features as a function of crater radius at the first octave
+   real(DP), parameter :: complex_wall_slope = 25._DP * DEG2RAD
+   integer(I4B), parameter :: complex_wall_terracefac = 32 ! Spatial size factor for terraces relative to crater radius
+   integer(I4B) :: complex_wall_terrace_num
+
+
+   newdem = surfi%dem
+   pikeD = min(1.044e3_DP * (crater%fcrat * 1e-3_DP)**(0.301_DP), user%deplimit) ! Pike (1977) depth/diameter ratio of complex craters
+   !pikeFloor = 
+   pikeFloor = 0.7_DP
+   if (r < pikeFloor) return
+   complex_wall_terrace_num = int(complex_wall_terracefac * r) + 1
+   noise = 0.0_DP
+   do octave = 1, complex_wall_num_octaves
+      xynoise = complex_wall_xy_noise_fac * complex_wall_freq ** octave / crater%fcrat
+      znoise = complex_wall_noise_height * complex_wall_pers ** (octave - 1) * crater%fcrat
+      noise = noise + util_perlin_noise(xynoise * x_relative + complex_wall_terrace_num * complex_wall_offset * rn(1), &
+                                        xynoise * y_relative + complex_wall_terrace_num * complex_wall_offset * rn(2)) &
+                                        * znoise
+   end do
+   newdem = max(newdem + noise,crater%melev - pikeD) !/ cos(complex_wall_slope)
+
+
+   elchange  = newdem - surfi%dem
+   deltaMi = elchange
+   surfi%dem = newdem
+
+end subroutine complex_wall
+
+
+
+

src/ejecta/ejecta_emplace.f90

@@ -94,7 +94,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot)
    real(DP),intent(in) :: deltaMtot
 
    ! Internal variables
-   real(DP) :: lrad,lradsq,cdepth
+   real(DP) :: lrad,lradsq
    integer(I4B),parameter :: MAXLOOP = 100 ! Maximum number of times to loop the ejecta angle correction calculation
    integer(I4B) :: xpi,ypi,i,j,k,n,inc,incsq,iradsq,idistorted,jdistorted
    real(DP) :: xp,yp,fradsq,fradpxsq,radsq,ebh,ejdissq,ejbmass,fmasscons,areafrac,xbar,ybar,krad,kdiffmax
@@ -121,9 +121,8 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot)
 
    !call regolith_melt_zone(user,crater,crater%imp,crater%impvel,rm,dm)
 
-   cdepth = DDRATIO * crater%fcrat
-   crater%vdepth = crater%ejrim + cdepth
-   crater%vrim   = crater%ejrim + crater%rheight
+   crater%vdepth = crater%ejrim + crater%floordepth
+   crater%vrim   = crater%ejrim + crater%rimheight
 
    if (crater%ejdis <= crater%rad) return