Realistic topography much improved for complex craters similar to Copernicus. Simple craters are currently broken.

src/crater/crater_emplace.f90

@@ -79,6 +79,9 @@ subroutine crater_emplace(user,surf,crater,domain,deltaMtot)
    fradsq = crater%frad**2
    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

@@ -27,7 +27,7 @@ subroutine crater_form_interior(user,surfi,crater,x_relative, y_relative ,newele
    ! Arguments
    type(usertype),intent(in) :: user
    type(surftype),intent(inout) :: surfi
-   type(cratertype),intent(in) :: crater
+   type(cratertype),intent(inout) :: crater
    real(DP),intent(in) :: x_relative, y_relative 
    real(DP),intent(in) :: newelev
    real(DP),intent(out) :: deltaMi
@@ -66,7 +66,6 @@ subroutine crater_form_interior(user,surfi,crater,x_relative, y_relative ,newele
 
 
    ! 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
@@ -77,6 +76,8 @@ subroutine crater_form_interior(user,surfi,crater,x_relative, y_relative ,newele
    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
 
    !if (r < r_floor) then
    !   cform = -simple_depth_diam  * crater%fcrat

src/crater/crater_populate.f90

@@ -196,11 +196,6 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
          ! Place crater onto the surface
          call crater_emplace(user,surf,crater,domain,ejbmass)
 
-
-         if (user%dorealistic) call crater_realistic_topography(user,surf,crater,domain,ejbmass)
-
-
-
          call ejecta_distance_estimate(user,crater,domain,crater%ejdis) ! Fast but imprecise estimate of the total ejecta distance
                                                                         ! For very steep size distributions, only a fraction of the
                                                                         ! craters are retained. The full ejecta_table_define function
@@ -220,6 +215,8 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
          else
             ejtble = 0
          end if
+
+         if (user%dorealistic) call crater_realistic_topography(user,surf,crater,domain,ejbmass)
 
          ! Collapse any remaining unstable slopes
          if (user%docollapse) call crater_slope_collapse(user,surf,crater,domain,ejbmass)