Updated profile code to better connect the Pike (1977) and Fassett & Thomson (2014) models. Also improved the complex rim code to better control the shape of the scallops

src/Makefile.am

@@ -86,6 +86,7 @@ crater/crater_critical_slope.f90\
 crater/crater_degradation_function.f90\
 crater/crater_visibility.f90\
 crater/crater_get_degradation_state.f90\
+crater/crater_profile.f90\
 init/init_domain.f90\
 init/init_dist.f90\
 init/init_surf.f90\

src/crater/crater_dimensions.f90

@@ -39,11 +39,14 @@ subroutine crater_dimensions(user,crater,domain)
 
    ! Set the crater morphology
    select case(crater%morphtype)
-      case("SIMPLE","TRANSITION") ! Following Pike (1977)
-         crater%rimheight  = 0.036_DP * (crater%fcrat * 1e-3_DP)**(1.014_DP) * 1e3_DP !* crater%fcrat
-         crater%rimwidth   = 0.257_DP * (crater%fcrat * 1e-3_DP)**(1.011_DP) * 1e3_DP !* crater%fcrat
-         crater%floordepth = 0.196_DP * (crater%fcrat * 1e-3_DP)**(1.010_DP) * 1e3_DP !* crater%fcrat
-         crater%floordiam  = 0.031_DP * (crater%fcrat * 1e-3_DP)**(1.765_DP) * 1e3_DP !* crater%fcrat
+      case("SIMPLE","TRANSITION") ! A hybrid model between Pike (1977) and Fassett & Thomson (2014) 
+         !crater%rimheight  = 0.036_DP * (crater%fcrat * 1e-3_DP)**(1.014_DP) * 1e3_DP !* crater%fcrat ! Pike model
+         crater%rimheight  = 0.043_DP * (crater%fcrat * 1e-3_DP)**(1.014_DP) * 1e3_DP !* crater%fcrat  ! Closer to Fassett & Thomson
+         crater%rimwidth   = 0.257_DP * (crater%fcrat * 1e-3_DP)**(1.011_DP) * 1e3_DP !* crater%fcrat  ! Pike model
+         !crater%floordepth = 0.196_DP * (crater%fcrat * 1e-3_DP)**(1.010_DP) * 1e3_DP !* crater%fcrat ! Pike model
+         crater%floordepth = 0.224_DP * (crater%fcrat * 1e-3_DP)**(1.010_DP) * 1e3_DP !* crater%fcrat  ! Closer to Fassett & Thomson
+         !crater%floordiam  = 0.031_DP * (crater%fcrat * 1e-3_DP)**(1.765_DP) * 1e3_DP !* crater%fcrat ! Pike model
+         crater%floordiam  = 0.200_DP * (crater%fcrat * 1e-3_DP)**(1.143_DP) * 1e3_DP !* crater%fcrat  ! Fassett & Thomson for D~1km, Pike for D~20km
       case("COMPLEX","PEAKRING","MULTIRING") ! Following Pike (1977)
          crater%rimheight  = 0.236_DP * (crater%fcrat * 1e-3_DP)**(0.399_DP) * 1e3_DP !* crater%fcrat
          crater%rimwidth   = 0.467_DP * (crater%fcrat * 1e-3_DP)**(0.836_DP) * 1e3_DP !* crater%fcrat
@@ -63,6 +66,10 @@ subroutine crater_dimensions(user,crater,domain)
    ! Adjust the floor depth to measure from the pre-existing level surface, rather than the rim
    crater%floordepth  = crater%floordepth - crater%rimheight
 
+   ! Calculate the radius where the inner wall meets the original pre-existing surface
+   ! This is used to demark the location where excavation transitions to deposition
+   crater%ejrad = crater_profile_find_r_inner_wall(user,crater) * crater%frad
+
    !find rim for counting purposes
    crater%frim = RIMFAC * crater%frad
 

src/crater/crater_emplace.f90

@@ -80,8 +80,6 @@ subroutine crater_emplace(user,surf,crater,domain,deltaMtot)
    deltaMtot = 0.0_DP !ejbmass
    incsq = inc**2
 
-   ! TEMP UNTIL BETTER FORMULA
-   crater%ejrad = crater%frad
    ! This loop may not be parallelizable because of the linked list operation inside crater_form_interior
    do j=-inc,inc  ! Do the loop in pixel space
       do i=-inc,inc

src/crater/crater_form_interior.f90

@@ -54,30 +54,29 @@ subroutine crater_form_interior(user,surfi,crater,x_relative, y_relative ,newele
    real(DP),parameter :: outer_c2 =  0.018_DP 
    real(DP),parameter :: outer_c3 =  0.015_DP
 
-   real(DP) :: c0,c1,c2,c3,flrad,rh,fld
+   !real(DP) :: c0,c1,c2,c3,flrad,rh,fld
 
 
    ! Executable code
    r = sqrt(x_relative**2+y_relative**2) / crater%frad
 
-   rh = crater%rimheight 
-   fld = -crater%floordepth 
-   flrad = 0.5_DP * crater%floordiam / crater%frad 
+   !rh = crater%rimheight 
+   !fld = -crater%floordepth 
+   !flrad = 0.5_DP * crater%floordiam / crater%frad 
 
 
    ! Use polynomial crater profile similar to that of Fassett et al. (2014), but the parameters are set by the crater dimensions
-   c1 = (fld - rh) / (flrad + flrad**2 / 3._DP - flrad**3 / 6._DP - 7._DP / 6._DP)
-   c0 = rh - (7._DP / 6._DP) * c1
-   c2 = c1 / 3._DP
-   c3 = -c2 / 2._DP
-
-   if (r < flrad) then
-      cform = fld 
-   else
-      cform = c0 + c1 * r + c2 * r**2 + c3 * r**3 
-   end if
-   ! TEMP UNTIL I WRITE THE SOLUTION PROPERLY
-   if (cform > 0.0_DP .and. r * crater%frad < crater%ejrad) crater%ejrad = r * crater%frad
+   !c1 = (fld - rh) / (flrad + flrad**2 / 3._DP - flrad**3 / 6._DP - 7._DP / 6._DP)
+   !c0 = rh - (7._DP / 6._DP) * c1
+   !c2 = c1 / 3._DP
+   !c3 = -c2 / 2._DP
+!
+!   if (r < flrad) then
+!      cform = fld 
+!   else
+!      cform = c0 + c1 * r + c2 * r**2 + c3 * r**3 
+!   end if
+
 
    !if (r < r_floor) then
    !   cform = -simple_depth_diam  * crater%fcrat
@@ -86,6 +85,8 @@ subroutine crater_form_interior(user,surfi,crater,x_relative, y_relative ,newele
    !else 
    !   cform =  (outer_c0 + outer_c1 * r + outer_c2 * r**2 + outer_c3 * r**3) * crater%fcrat
    !end if      
+
+   cform = crater_profile(user,crater,r)
    newdem = newelev + cform 
 
    !if (crater%fcrat > crater%cxtran * 2) then 

src/crater/crater_profile.f90

@@ -27,29 +27,134 @@ function crater_profile(user,crater,r) result(h)
    type(usertype),intent(in) :: user
    type(cratertype),intent(in) :: crater
    real(DP),intent(in) :: r
+
+   ! Result variable
    real(DP) :: h
 
    ! Internal variables
    real(DP) :: flrad,c0,c1,c2,c3
+   real(DP),parameter :: A = 4._DP / 11._DP
+   real(DP),parameter :: B = -32._DP / 187._DP
 
    ! Executable code
-   flrad = 0.5_DP * crater%floordiam / crater%frad 
+   flrad = crater%floordiam / crater%fcrat
 
    ! Use polynomial crater profile similar to that of Fassett et al. (2014), but the parameters are set by the crater dimensions
-   c1 = (-crater%floordepth - crater%rimheight) / (flrad + flrad**2 / 3._DP - flrad**3 / 6._DP - 7._DP / 6._DP)
-   c0 = crater%rimheight - (7._DP / 6._DP) * c1
-   c2 = c1 / 3._DP
-   c3 = -c2 / 2._DP
+   c1 = (-crater%floordepth - crater%rimheight) / (flrad - 1._DP + A * (flrad**2 - 1._DP) + B * (flrad**3 - 1._DP))
+   c0 = crater%rimheight - c1 * (1._DP + A + B)
+   c2 = A * c1
+   c3 = B * c1
 
    if (r < flrad) then
       h = -crater%floordepth 
-   if (r > crater%frad)
-      h = crater%rimheight * ((crater%frad / lrad)**RIMDROP)
+   else if (r > 1.0_DP) then
+      h = crater%rimheight * (r**(-RIMDROP))
    else
       h = c0 + c1 * r + c2 * r**2 + c3 * r**3 
    end if
 
+   return
+end function crater_profile
+
+
+function crater_profile_find_r_inner_wall(user,crater) result(r_inner_wall)
+   !Uses Brent's method to solve for the point in the crater profile where the inner wall of the crater profile meets the original surface
+   use module_globals
+   use module_util
+   use module_crater, EXCEPT_THIS_ONE => crater_profile_find_r_inner_wall
+   implicit none
+
+   ! Arguments
+   type(usertype),intent(in) :: user
+   type(cratertype),intent(in) :: crater
+
+   ! Result variable
+   real(DP) :: r_inner_wall
+
+   integer(I4B),parameter :: maxIterations=100
+   real(DP) :: valueAtRoot
+
+   real(DP),parameter :: Tolerance = 1e-8_DP
+   real(DP),parameter :: FPP = 1.e-11_DP
+   real(DP),parameter :: nearzero = 1.e-20_DP
+
+   integer :: niter
 
+   real(DP) :: resultat,AA,BB,CC,DD,EE,FA,FB,FC,Tol1,PP,QQ,RR,SS,xm,x1,x2
+   integer :: i, done
+
+   i = 0
+   done = 0
+
+   x1 = 0.0_DP
+   x2 = 1.0_DP
+
+   AA = x1
+   BB = x2
+   FA = crater_profile(user,crater,AA) 
+   FB = crater_profile(user,crater,BB) 
+
+   FC = FB;
+   do while (done.eq.0.and.i < maxIterations)
+      if (.not.util_rootbracketed(FC,FB)) then
+         CC = AA; FC = FA; DD = BB - AA; EE = DD
+      endif
+      if (abs(FC) < abs(FB)) then
+         AA = BB; BB = CC; CC = AA
+         FA = FB; FB = FC; FC = FA
+      endif
+      Tol1 = 2 * FPP * abs(BB) + 0.5_DP * Tolerance
+      xm = 0.5_DP * (CC-BB)
+      if ((abs(xm) <= Tol1).or.(abs(FA) < nearzero)) then
+        ! A root has been found
+        resultat = BB;
+        done = 1
+        valueAtRoot = crater_profile(user,crater,resultat) 
+      else 
+        if ((abs(EE) >= Tol1).and.(abs(FA) > abs(FB))) then
+          SS = FB/ FA;
+          if (abs(AA - CC) < nearzero) then
+            PP = 2 * xm * SS;
+            QQ = 1._DP - SS;
+          else 
+            QQ = FA/FC;
+            RR = FB /FC;
+            PP = SS * (2 * xm * QQ * (QQ - RR) - (BB-AA) * (RR - 1._DP));
+            QQ = (QQ - 1._DP) * (RR - 1._DP) * (SS - 1._DP);
+          endif
+          if (PP > nearzero) QQ = -QQ;
+          PP = abs(PP);
+          if ((2 * PP) < min(3 * xm * QQ - abs(Tol1 * QQ),abs( EE * QQ))) then
+            EE = DD;  DD = PP/QQ;
+          else 
+            DD = xm;   EE = DD;
+          endif
+        else 
+          DD = xm;
+          EE = DD;
+        endif
+        AA = BB;
+        FA = FB;
+        if (abs(DD) > Tol1) then 
+          BB = BB + DD;
+        else 
+          if (xm > 0._DP) then 
+            BB = BB + abs(Tol1)
+          else 
+            BB = BB - abs(Tol1)
+          endif
+        endif
+        FB = crater_profile(user,crater,BB)
+        i=i+1
+      endif
+   end do
+   if (i >= maxIterations) then
+      write(*,*) 'Too many iterations in profile solver.'
+      r_inner_wall = 1._DP
+      return
+   end if
+   niter = i
+   r_inner_wall = resultat
    return
-end subroutine crater_profile
+end function crater_profile_find_r_inner_wall
 

src/crater/crater_realistic_topography.f90

@@ -106,7 +106,7 @@ end subroutine complex_rim
 
    call complex_rim(user,surf,crater,rn,deltaMtot)
    !call complex_terrace(user,surf,crater,rn,deltaMtot)
-   call crater_slope_collapse(user,surf,crater,domain,(0.35_DP * user%pix)**2,deltaMtot)
+   !call crater_slope_collapse(user,surf,crater,domain,(0.35_DP * user%pix)**2,deltaMtot)
    !call complex_wall_texture(user,surf,crater,rn,deltaMtot)
    call complex_floor(user,surf,crater,rn,deltaMtot)
    call complex_peak(user,surf,crater,rn,deltaMtot)
@@ -343,22 +343,22 @@ subroutine complex_rim(user,surf,crater,rn,deltaMtot)
    real(DP) :: xynoise, znoise
    real(DP) :: noise,dnoise
 
-   ! Complex crater wall
-   integer(I4B)        :: offset       ! Scales the random xy-offset so that each crater's random noise is unique 
-   real(DP)            :: xy_noise_fac ! Spatial "size" of noise features at the first octave
-   real(DP)            :: freq         ! Spatial size scale factor multiplier at each octave level
-   real(DP)            :: pers         ! The relative size scaling at each octave level
-   real(DP)            :: noise_height ! Vertical height of noise features as a function of crater radius at the first octave
+   ! Complex crater rim
 
+   integer(I4B)            :: nscallops             ! Approximate number of scallop features on edge of crater
+   integer(I4B),parameter  :: offset = 3000         ! Scales the random xy-offset so that each crater's random noise is unique 
+   real(DP),parameter      :: noise_height = 0.15_DP ! Vertical height of noise features as a function of crater radius at the first octave
+                                                    ! Higher values make the scallop features broader
+   real(DP),parameter      :: noise_slope = 0.6_DP  ! The slope of the power law function that defines the noise shape for scallop features
+                                                    ! Higher values makes the inner sides of the scallop features more rounded, lower values 
+                                                    ! make them have sharper points at the inflections
+
+   nscallops = 16
    ! determine area to effect
    ! First make the interior of the crater
    rad = 1.25_DP * crater%frad
 
    ! Create diffusion noise for scalloped rim
-   offset = 3000
-   xy_noise_fac = 5.00_DP
-   noise_height = 1.00_DP
-
 
    inc = max(min(nint(rad / user%pix),PBCLIM*user%gridsize),1) + 1
    crater%maxinc = max(crater%maxinc,inc)
@@ -378,28 +378,22 @@ subroutine complex_rim(user,surf,crater,rn,deltaMtot)
          r = sqrt(xbar**2 + ybar**2) / crater%frad
 
          ! Make scalloped rim
-         znoise = noise_height * crater%fcrat
-         xynoise = xy_noise_fac / crater%fcrat
+         znoise = noise_height**(1._DP / (2 * noise_slope))
+         xynoise = (nscallops / PI) / crater%fcrat
          dnoise = util_perlin_noise(xynoise * xbar + offset * rn(1), &
                                     xynoise * ybar + offset * rn(2)) * znoise
-         noise = sqrt(sqrt(dnoise**2))
+         noise = (dnoise**2)**noise_slope
 
+         hprof = r**(-1)
+         tprof = 0.5_DP * crater%ejrim + crater%melev 
 
-         if (r >= 1.0) then
-            hprof = crater%ejrim * r**(-3) + crater%melev
-            tprof = 0.5_DP * crater%ejrim + crater%melev !* ((2.0_DP - r)**(-4) - 1.0_DP) + crater%melev
-         else
-            hprof = crater%ejrim * (2.0_DP - r)**(-2) + crater%melev
-            tprof = 0.5_DP * crater%ejrim + crater%melev ! (r**(-4) - 1.0_DP) + crater%melev
-         end if
-
-         if (crater%melev + crater%ejrim - noise < hprof) then
-            newdem = tprof - noise
+         if (1.0_DP - noise < hprof) then
+            newdem = min(tprof,newdem) 
+            elchange  = newdem - surf(xpi,ypi)%dem
+            deltaMtot = deltaMtot + elchange
+            surf(xpi,ypi)%dem = newdem
          end if
 
-         elchange  = newdem - surf(xpi,ypi)%dem
-         deltaMtot = deltaMtot + elchange
-         surf(xpi,ypi)%dem = newdem
       end do
    end do
 

src/crater/module_crater.f90

@@ -17,15 +17,15 @@ module module_crater
 public 
 save
 
-interface
-   subroutine crater_mass_conservation(user,surf,crater)
-   use module_globals
-   implicit none
-   type(usertype),intent(in) :: user
-   type(surftype),dimension(:,:),intent(inout) :: surf
-   type(cratertype),intent(in)  :: crater
-   end subroutine crater_mass_conservation
-end interface
+   interface
+      subroutine crater_mass_conservation(user,surf,crater)
+      use module_globals
+      implicit none
+      type(usertype),intent(in) :: user
+      type(surftype),dimension(:,:),intent(inout) :: surf
+      type(cratertype),intent(in)  :: crater
+      end subroutine crater_mass_conservation
+   end interface
 
 
    interface
@@ -313,4 +313,26 @@ end function crater_visibility
    end interface
 
 
+   interface
+      function crater_profile(user,crater,r) result(h)
+      use module_globals
+      implicit none
+      type(usertype),intent(in) :: user
+      type(cratertype),intent(in) :: crater
+      real(DP),intent(in) :: r
+      real(DP) :: h
+      end function crater_profile
+   end interface
+
+
+   interface
+      function crater_profile_find_r_inner_wall(user,crater) result(r_inner_wall)
+      use module_globals
+      implicit none
+      type(usertype),intent(in) :: user
+      type(cratertype),intent(in) :: crater
+      real(DP) :: r_inner_wall
+      end function crater_profile_find_r_inner_wall
+   end interface
+
 end module