From 6373bf6b07893e99824db693241f12d9fc3654f5 Mon Sep 17 00:00:00 2001 From: daminton Date: Wed, 25 Jan 2017 15:00:03 +0000 Subject: [PATCH] Reverted back to making final craterform initially --- src/crater/crater_emplace.f90 | 10 ++++------ src/crater/crater_form_exterior.f90 | 16 ++++++++-------- src/crater/crater_form_exterior_func.f90 | 2 +- src/crater/crater_form_interior.f90 | 17 ++++++++++++----- src/crater/module_crater.f90 | 4 ++-- 5 files changed, 27 insertions(+), 22 deletions(-) diff --git a/src/crater/crater_emplace.f90 b/src/crater/crater_emplace.f90 index 2d3684f8..ca774af6 100644 --- a/src/crater/crater_emplace.f90 +++ b/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) diff --git a/src/crater/crater_form_exterior.f90 b/src/crater/crater_form_exterior.f90 index 83d196a7..76254898 100644 --- a/src/crater/crater_form_exterior.f90 +++ b/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 diff --git a/src/crater/crater_form_exterior_func.f90 b/src/crater/crater_form_exterior_func.f90 index 8ebe608d..46a94dbe 100644 --- a/src/crater/crater_form_exterior_func.f90 +++ b/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 diff --git a/src/crater/crater_form_interior.f90 b/src/crater/crater_form_interior.f90 index 222cde06..d8a1d29d 100644 --- a/src/crater/crater_form_interior.f90 +++ b/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 diff --git a/src/crater/module_crater.f90 b/src/crater/module_crater.f90 index 89e06b7f..755272f7 100644 --- a/src/crater/module_crater.f90 +++ b/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