From 10aa7f3b5a565de4995492195b6f3cd067eba871 Mon Sep 17 00:00:00 2001 From: David Minton Date: Tue, 24 Mar 2020 22:15:17 -0400 Subject: [PATCH] Updated complex morphology parameters --- src/crater/crater_dimensions.f90 | 38 ++++++++++++----- src/crater/crater_form_interior.f90 | 48 +++++++++++++++------- src/crater/crater_generate.f90 | 3 ++ src/crater/crater_populate.f90 | 4 +- src/crater/crater_realistic_topography.f90 | 27 +++++------- src/globals/module_globals.f90 | 5 ++- 6 files changed, 80 insertions(+), 45 deletions(-) diff --git a/src/crater/crater_dimensions.f90 b/src/crater/crater_dimensions.f90 index 8c63b614..6ea7302a 100644 --- a/src/crater/crater_dimensions.f90 +++ b/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 diff --git a/src/crater/crater_form_interior.f90 b/src/crater/crater_form_interior.f90 index 889b8d01..3f5b76bc 100644 --- a/src/crater/crater_form_interior.f90 +++ b/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 diff --git a/src/crater/crater_generate.f90 b/src/crater/crater_generate.f90 index f10ccc98..cc6b9a23 100644 --- a/src/crater/crater_generate.f90 +++ b/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 diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index bf8ef50f..d1e2f03b 100644 --- a/src/crater/crater_populate.f90 +++ b/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 diff --git a/src/crater/crater_realistic_topography.f90 b/src/crater/crater_realistic_topography.f90 index b24f55b0..5c4ef4e6 100644 --- a/src/crater/crater_realistic_topography.f90 +++ b/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,7 +242,7 @@ 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 @@ -251,22 +250,18 @@ subroutine complex_wall(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi) 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 diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index 11a9d383..abb30406 100644 --- a/src/globals/module_globals.f90 +++ b/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