Overhauled crater formation routine. All raised rims are now part of …

…the ejecta blanket
MintonGroup · Feb 9, 2017 · 98dd642 · 98dd642
1 parent 797aa69
commit 98dd642
Showing 1 changed file with 5 additions and 9 deletions.
diff --git a/src/crater/crater_form_interior.f90 b/src/crater/crater_form_interior.f90
@@ -42,24 +42,18 @@ subroutine crater_form_interior(user,surfi,crater,lradsq,newelev,deltaMi)
    ! Executable code
 
    !change digital elevation map
-   !cform = crater%vcorr - (crater%parab * lradsq)
-   !if (lradsq < crater%rad**2) then
    r = sqrt(lradsq) / crater%frad
-   !h = 2.5_DP*crater%rad !0.5_DP*crater%fcrat
-   !circrad = 0.5_DP * h + (2*crater%rad)**2 / (8 * h)
-   !cform = sqrt(circrad**2 - lradsq) + (circrad - h) + 0.15_DP*crater%frad
-   !          
+   ! Use empirical crater form from Fassett et al. 2014
    if (r < 0.2_DP) then
       cform = -0.181_DP * crater%fcrat
    else if (r < 0.98_DP) then
-      cform =  (-0.229_DP + 0.228_DP * r + 0.083_DP * r**2 - 0.039_DP * r**3) * crater%fcrat
+      cform =  (-0.229_DP + 0.228_DP * r + 0.083_DP * r**2 - 0.039_DP * r**3) * crater%fcrat 
    else 
-      cform =  crater%rheight !0.188_DP - 0.187_DP * r + 0.018_DP * r**2 + 0.015_DP * r**3 
+      cform =  (0.188_DP - 0.187_DP * r + 0.018_DP * r**2 + 0.015_DP * r**3) * crater%fcrat
    end if      
    newdem = newelev + cform 
 
    pikeD = 1.044e3_DP * (crater%fcrat * 1e-3_DP)**(0.301_DP) ! Pike (1977)
-   !write(*,*) melev,pikeD
    if ((crater%fcrat > crater%cxtran * 2) .and. newdem < (crater%melev - pikeD)) then
       newdem = crater%melev - pikeD ! Flatten out the bottom of the crater
    end if
@@ -69,6 +63,8 @@ subroutine crater_form_interior(user,surfi,crater,lradsq,newelev,deltaMi)
          call util_remove_from_layer(surfi,layer)
       end do
    end if
+
+   newdem = min(newdem,surfi%dem) ! Only allow excavation, no deposition
    elchange  = newdem - surfi%dem
    deltaMi = elchange
    surfi%dem = newdem