Updated complex morphology parameters

src/crater/crater_dimensions.f90

@@ -30,21 +30,38 @@ subroutine crater_dimensions(user,crater,domain)
    ! Internal variables
    real(DP) :: lrad,cform
 
-   real(DP),parameter :: DDRATIO = 0.19_DP        ! Depth-diameter ratio
-   real(DP),parameter :: RDRATIO = 0.030_DP       ! Rim height to diameter ratio
+!   real(DP),parameter :: DDRATIO = 0.19_DP        ! Depth-diameter ratio
+!   real(DP),parameter :: RDRATIO = 0.030_DP       ! Rim height to diameter ratio
    real(DP),parameter :: RIMFAC = 1.5_DP          ! Ratio of radius used for counting craters to the final rim radius
 
+   !TODO: Add in transitional and multi-ringed basin crater morphology
+   !TODO: Implement Monte Carlo model for standard errors for all dimension parameters
+
+   ! 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("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
+         crater%floordepth = 1.044_DP * (crater%fcrat * 1e-3_DP)**(0.301_DP) * 1e3_DP !* crater%fcrat
+         crater%floordiam  = 0.187_DP * (crater%fcrat * 1e-3_DP)**(1.249_DP) * 1e3_DP !* crater%fcrat
+         crater%peakheight = 0.032_DP * (crater%fcrat * 1e-3_DP)**(0.900_DP) * 1e3_DP !* crater%fcrat
+   end select
+
 
    ! Executable code
-   crater%floordepth = DDRATIO * crater%fcrat
-   if (crater%fcrat <= crater%cxtran) then
-      crater%rimheight = RDRATIO * crater%fcrat
-   else
-      crater%rimheight = (RDRATIO*crater%cxtran)+(RDRATIO*((crater%fcrat-crater%cxtran)**(0.399_DP)))
-   endif
+   !crater%floordepth = DDRATIO * crater%fcrat
+   !if (crater%fcrat <= crater%cxtran) then
+   !else
+   !   crater%rimheight = (RDRATIO*crater%cxtran)+(RDRATIO*((crater%fcrat-crater%cxtran)**(0.399_DP)))
+   !endif
 
-   crater%vcorr  = crater%floordepth - crater%rimheight
-   crater%parab  = crater%floordepth / ((crater%frad)**2)
+   ! Adjust the floor depth to measure from the pre-existing level surface, rather than the rim
+   crater%floordepth  = crater%floordepth - crater%rimheight
 
    !find rim for counting purposes
    crater%frim = RIMFAC * crater%frad
@@ -64,6 +81,7 @@ subroutine crater_dimensions(user,crater,domain)
    crater%fradpx = int(crater%frad/user%pix) + 1
    crater%rimdispx = int(crater%rimdis/user%pix)  + 1
 
+
    return
 end subroutine crater_dimensions
 

src/crater/crater_form_interior.f90

@@ -33,7 +33,7 @@ subroutine crater_form_interior(user,surfi,crater,x_relative, y_relative ,newele
    real(DP),intent(out) :: deltaMi
 
    ! Internal variables
-   real(DP) :: cform,newdem,elchange,pikeD,r, Hpeak
+   real(DP) :: cform,newdem,elchange,r
    integer(I4B) :: layer
 
    ! A list for popped data 
@@ -54,31 +54,49 @@ 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
 
-   ! Executable code
 
-   !change digital elevation map
+   ! Executable code
    r = sqrt(x_relative**2+y_relative**2) / crater%frad
-   ! Use empirical crater form from Fassett et al. 2014
-   if (r < r_floor) then
-      cform = -simple_depth_diam  * crater%fcrat
-   else if (r < r_rim) then
-      cform =  (inner_c0 + inner_c1 * r + inner_c2 * r**2 + inner_c3 * r**3) * crater%fcrat 
-   else 
-      cform =  (outer_c0 + outer_c1 * r + outer_c2 * r**2 + outer_c3 * r**3) * crater%fcrat
-   end if      
+
+   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
+
+   !if (r < r_floor) then
+   !   cform = -simple_depth_diam  * crater%fcrat
+   !else if (r < r_rim) then
+   !   cform =  (inner_c0 + inner_c1 * r + inner_c2 * r**2 + inner_c3 * r**3) * crater%fcrat 
+   !else 
+   !   cform =  (outer_c0 + outer_c1 * r + outer_c2 * r**2 + outer_c3 * r**3) * crater%fcrat
+   !end if      
    newdem = newelev + cform 
 
-   pikeD = min(1.044e3_DP * (crater%fcrat * 1e-3_DP)**(0.301_DP), user%deplimit) ! Pike (1977) depth/diameter ratio of complex craters
-   if (crater%fcrat > crater%cxtran * 2) then 
+   !if (crater%fcrat > crater%cxtran * 2) then 
+   if (crater%morphtype == "COMPLEX") then
       ! Make this crater complex
 
       !if (r < r_rim) then
-         if (newdem < (crater%melev - pikeD)) then ! Flatten out the bottom of the crater
+         if (newdem < (crater%melev - crater%floordepth)) then ! Flatten out the bottom of the crater
             do layer = 1,user%numlayers ! Remove all pre-existing craters from this current pixel
                call util_remove_from_layer(surfi,layer)
             end do
-            newdem = crater%melev - pikeD 
+            newdem = crater%melev - crater%floordepth
 
 
          end if

src/crater/crater_generate.f90

@@ -175,11 +175,14 @@ subroutine crater_generate(user,crater,domain,prod,production_list,vdist,surf)
 
    trfin = 2 * TRSIM * crater%rad
    if (trfin <= crater%cxtran) then
+      crater%morphtype = "SIMPLE"
       crater%fcrat = trfin   ! Simple Crater
    else  ! Complex crater
+      crater%morphtype = "COMPLEX"
       crater%fcrat = trfin * ((trfin/crater%cxtran)**crater%cxexp) ! Complex Crater
    end if
    if (crater%imp >= user%basinimp) then
+      crater%morphtype = "MULTIRING" 
       ! This section is temporary until a better basin scaling law can be implemented. Potter et al. (2012) GRL v39. p 18203
       if ((crater%xl < (0.5_DP * domain%side) ).or.(user%testflag).or.(domain%initialize)) then ! pick TP1 for left-hand hemisphere
          crater%fcrat = 0.1354e3_DP * (0.5_DP * trfin * 1e-3_DP)**(1.389_DP)  ! TP2

src/crater/crater_populate.f90

@@ -243,7 +243,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
             imax = crater%imp
             cmax = crater%fcrat
             rhmax = crater%vrim
-            if (crater%vcorr <= user%deplimit) then
+            if (crater%floordepth <= user%deplimit) then
                rmax = crater%vdepth
             else
                rmax = user%deplimit + crater%vrim 
@@ -254,7 +254,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
             imin = crater%imp
             cmin = crater%fcrat
             rhmin = crater%vrim
-            if (crater%vcorr <= user%deplimit) then
+            if (crater%floordepth <= user%deplimit) then
                rmin = crater%vdepth
             else
                rmin = user%deplimit + crater%vrim

src/crater/crater_realistic_topography.f90

@@ -154,7 +154,7 @@ subroutine complex_peak(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi)
    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
+   real(DP) :: newdem,elchange
 
    ! Topographic noise parameters
    real(DP) :: xynoise, znoise
@@ -164,11 +164,10 @@ subroutine complex_peak(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi)
    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)
+      znoise = crater%peakheight * (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 
@@ -203,7 +202,7 @@ subroutine complex_floor(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi)
    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
+   real(DP) :: newdem,elchange
 
    ! Topographic noise parameters
    real(DP) :: xynoise, znoise
@@ -243,30 +242,26 @@ subroutine complex_wall(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi)
    real(DP),intent(out) :: deltaMi
 
    ! Internal variables
-   real(DP) :: newdem,elchange,pikeD,pikeFloor,Hpeak
+   real(DP) :: newdem,elchange
 
    ! 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_num_octaves  = 4   ! 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_xy_noise_fac = 1.5_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
+   real(DP), parameter :: complex_wall_pers = 0.70_DP  ! The relative size scaling at each octave level
+   real(DP), parameter :: complex_wall_noise_height = 0.02_DP ! Vertical height of noise features as a function of crater radius at the first octave
+   integer(I4B), parameter :: complex_wall_terracefac = 8 ! 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
+   if (r < crater%floordiam / crater%fcrat) return
    complex_wall_terrace_num = int(complex_wall_terracefac * r) + 1
    noise = 0.0_DP
    do octave = 1, complex_wall_num_octaves
@@ -276,7 +271,7 @@ subroutine complex_wall(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi)
                                         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)
+   newdem = max(newdem + noise,crater%melev - crater%floordepth) 
 
 
    elchange  = newdem - surfi%dem

src/globals/module_globals.f90

@@ -97,8 +97,8 @@ module module_globals
    real(DP) :: grad                ! Strengthless material transient crater radius
    real(DP) :: frad                ! Final crater radius
    real(DP) :: fcrat               ! Final crater diameter
-   real(DP) :: vdepth,vrim,vcorr   ! parameter for parabolic crater form
-   real(DP) :: frim,parab          ! parameter for parabolic crater form
+   real(DP) :: vdepth,vrim         ! parameter for parabolic crater form
+   real(DP) :: frim                ! parameter for parabolic crater form
    real(DP) :: rimdis              ! crater form radius (bowl + upturned rim)
    real(DP) :: ejdis               ! ejecta max distance
    real(DP) :: ejrim               ! ejecta height at crater rim
@@ -114,6 +114,7 @@ module module_globals
    real(DP) :: melev,xslp,yslp     ! Mean elevation and slopes at pre-existing impact site
 
    !Crater dimension information  - See Pike (1977) Impact and Explosion Cratering, Fig. 1
+   character(STRMAX) :: morphtype ! Type of crater. One of: "SIMPLE", "TRANSITION", "COMPLEX", "PEAKRING", "MULTIRING"
    real(DP) :: rimheight
    real(DP) :: rimwidth 
    real(DP) :: floordepth