Reverted back to making final craterform initially

src/crater/crater_emplace.f90

@@ -96,10 +96,12 @@ subroutine crater_emplace(user,surf,crater,domain,deltaMtot)
          call util_periodic(xpi,ypi,user%gridsize)
 
          lradsq = (crater%xl - xp)**2 + (crater%yl - yp)**2
-         if (lradsq <= crater%rad**2) then
+         if (lradsq <= crater%frad**2) then
             call crater_form_interior(user,surf(xpi,ypi),crater,lradsq,newelev,deltaMi)
-            deltaMtot = deltaMtot + deltaMi
+         else
+            call crater_form_exterior(user,surf(xpi,ypi),crater,domain,lradsq,newelev,deltaMi)
          end if
+         deltaMtot = deltaMtot + deltaMi
 
          ! do porosity computation if (user%doporosity)
          ! It is still important to consider the physical meaning of frad and rad. 
@@ -114,10 +116,6 @@ subroutine crater_emplace(user,surf,crater,domain,deltaMtot)
       end do
    end do
 
-   !call crater_form_exterior_rootfind(user,surf,crater,domain,deltaMtot)
-   lastloop = .true.
-   deltaMtot = crater_form_exterior_func(user,surf,crater,domain,RIMDROP,deltaMtot,lastloop)
-
    domain%tallycoverage = domain%tallycoverage + int(fradsq * PI / user%pix**2)
    domain%subpixelcoverage = domain%subpixelcoverage + int(fradsq * PI / user%pix**2)
 

src/crater/crater_form_exterior.f90

@@ -18,7 +18,7 @@
 !  Notes       :  
 !
 !**********************************************************************************************************************************
-subroutine crater_form_exterior(user,surfi,crater,domain,lradsq,newelev,rd,deltaMi)
+subroutine crater_form_exterior(user,surfi,crater,domain,lradsq,newelev,deltaMi)
    use module_globals
    use module_util
    use module_ejecta
@@ -30,7 +30,7 @@ subroutine crater_form_exterior(user,surfi,crater,domain,lradsq,newelev,rd,delta
    type(surftype),intent(inout) :: surfi
    type(cratertype),intent(in) :: crater
    type(domaintype),intent(in) :: domain
-   real(DP),intent(in) :: lradsq,newelev,rd
+   real(DP),intent(in) :: lradsq,newelev
    real(DP),intent(out) :: deltaMi
    integer(I4B) :: l
 
@@ -43,15 +43,15 @@ subroutine crater_form_exterior(user,surfi,crater,domain,lradsq,newelev,rd,delta
    lrad = sqrt(lradsq)
 
    ! change digital elevation map
-   cform = (1.3_DP * crater%rheight * ((crater%rad / lrad)**rd)) 
+   cform = crater%rheight * ((crater%rad / lrad)**RIMDROP) 
 
    ! elevation changes
-   if (lrad <= crater%frim) then
-      hcorr = cform * ((crater%frim - lrad) / (crater%frim - crater%rad))**20 
-      elchange = hcorr + cform 
-   else 
+   !if (lrad <= crater%frim) then
+   !   hcorr = cform * ((crater%frim - lrad) / (crater%frim - crater%rad))**20 
+   !   elchange = hcorr + cform 
+   !else 
       elchange = cform 
-   end if
+   !end if
    surfi%dem = surfi%dem + elchange 
    deltaMi = elchange
 

src/crater/crater_form_exterior_func.f90

@@ -45,7 +45,7 @@ function crater_form_exterior_func(user,surf,crater,domain,rd,deltaMtot,lastloop
    crater%maxinc = max(crater%maxinc,inc)
    incsq = inc**2
 
-   radsq = crater%rad**2
+   radsq = crater%frad**2
 
    deltaMp = 0.0_DP
 

src/crater/crater_form_interior.f90

@@ -45,11 +45,18 @@ subroutine crater_form_interior(user,surfi,crater,lradsq,newelev,deltaMi)
    !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
-
-   newdem = newelev - cform
+   !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
+   !          
+   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
+   else 
+      cform =  crater%rheight !0.188_DP - 0.187_DP * r + 0.018_DP * r**2 + 0.015_DP * r**3 
+   end if      
+   newdem = newelev + cform 
 
    pikeD = 1.044e3_DP * (crater%fcrat * 1e-3_DP)**(0.301_DP) ! Pike (1977)
    !write(*,*) melev,pikeD

src/crater/module_crater.f90

@@ -120,14 +120,14 @@ end subroutine crater_form_interior
    end interface
 
    interface
-      subroutine crater_form_exterior(user,surfi,crater,domain,lradsq,newelev,rd,deltaMi)
+      subroutine crater_form_exterior(user,surfi,crater,domain,lradsq,newelev,deltaMi)
       use module_globals
       implicit none
       type(usertype),intent(in) :: user
       type(surftype),intent(inout) :: surfi
       type(cratertype),intent(in) :: crater
       type(domaintype),intent(in) :: domain
-      real(DP),intent(in) :: newelev,lradsq,rd
+      real(DP),intent(in) :: newelev,lradsq
       real(DP),intent(out) :: deltaMi
       end subroutine crater_form_exterior
    end interface