From 1739f9396d4619b418a64b0d96fa04fbf78e3aa9 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 3 Nov 2021 14:05:13 -0400 Subject: [PATCH] Fixed problem dealing with periodic boundary for the realistic morphology functions. --- src/crater/crater_realistic_topography.f90 | 52 +++++++++------------- 1 file changed, 21 insertions(+), 31 deletions(-) diff --git a/src/crater/crater_realistic_topography.f90 b/src/crater/crater_realistic_topography.f90 index 7031867a..a34c0f98 100644 --- a/src/crater/crater_realistic_topography.f90 +++ b/src/crater/crater_realistic_topography.f90 @@ -184,13 +184,13 @@ subroutine complex_peak(user,surf,crater,deltaMtot) xpi = crater%xlpx + i ypi = crater%ylpx + j + xbar = xpi * user%pix - crater%xl + ybar = ypi * user%pix - crater%yl + ! periodic boundary conditions call util_periodic(xpi,ypi,user%gridsize) newdem = surf(xpi,ypi)%dem - xbar = xpi * user%pix - crater%xl - ybar = ypi * user%pix - crater%yl - areafrac = util_area_intersection(0.5_DP * crater%floordiam,xbar,ybar,user%pix) r = sqrt(xbar**2 + ybar**2) / crater%frad @@ -255,13 +255,13 @@ subroutine complex_floor(user,surf,crater,deltaMtot) xpi = crater%xlpx + i ypi = crater%ylpx + j + xbar = xpi * user%pix - crater%xl + ybar = ypi * user%pix - crater%yl + ! periodic boundary conditions call util_periodic(xpi,ypi,user%gridsize) newdem = surf(xpi,ypi)%dem - xbar = xpi * user%pix - crater%xl - ybar = ypi * user%pix - crater%yl - areafrac = util_area_intersection(0.5_DP * crater%floordiam,xbar,ybar,user%pix) ! Make the floor topographically "noisy" @@ -330,13 +330,13 @@ subroutine complex_wall_texture(user,surf,crater,domain,deltaMtot) xpi = crater%xlpx + i ypi = crater%ylpx + j + xbar = xpi * user%pix - crater%xl + ybar = ypi * user%pix - crater%yl + ! periodic boundary conditions call util_periodic(xpi,ypi,user%gridsize) newdem = surf(xpi,ypi)%dem - xbar = xpi * user%pix - crater%xl - ybar = ypi * user%pix - crater%yl - r = sqrt(xbar**2 + ybar**2) / crater%frad areafrac = 1.0 - util_area_intersection(0.3_DP * crater%floordiam,xbar,ybar,user%pix) @@ -429,7 +429,6 @@ subroutine complex_terrace(user,surf,crater,deltaMtot) pers_tfloor = 0.5_DP xy_size_tfloor = 5.0_DP / crater%fcrat - ! Lansberg values !terracefac = 1.0_DP !nterraces = 16 @@ -438,10 +437,6 @@ subroutine complex_terrace(user,surf,crater,deltaMtot) !scallop_width = 0.10_DP !^^^^^^^^^^^^^^^ - - - - rad = 2.0_DP * crater%frad upshift = 0._DP @@ -458,13 +453,13 @@ subroutine complex_terrace(user,surf,crater,deltaMtot) xpi = crater%xlpx + i ypi = crater%ylpx + j + xbar = xpi * user%pix - crater%xl + ybar = ypi * user%pix - crater%yl + ! periodic boundary conditions call util_periodic(xpi,ypi,user%gridsize) newdem = surf(xpi,ypi)%dem - xbar = xpi * user%pix - crater%xl - ybar = ypi * user%pix - crater%yl - r = sqrt(xbar**2 + ybar**2) / crater%frad ! Make scalloped terraces @@ -484,10 +479,8 @@ subroutine complex_terrace(user,surf,crater,deltaMtot) xynoise * ybar + offset * rn(2))* znoise end do - isterrace = 1.0_DP - noise < hprof - tprof = crater_profile(user,crater,rinner) * (r/rinner)**2 + tnoise ! This is the floor profile that replaces the old one at each terrace if (isterrace) then @@ -511,13 +504,13 @@ subroutine complex_terrace(user,surf,crater,deltaMtot) xpi = crater%xlpx + i ypi = crater%ylpx + j + xbar = xpi * user%pix - crater%xl + ybar = ypi * user%pix - crater%yl + ! periodic boundary conditions call util_periodic(xpi,ypi,user%gridsize) newdem = surf(xpi,ypi)%dem - xbar = xpi * user%pix - crater%xl - ybar = ypi * user%pix - crater%yl - r = sqrt(xbar**2 + ybar**2) / crater%frad ! Make scalloped rim @@ -557,13 +550,13 @@ subroutine complex_terrace(user,surf,crater,deltaMtot) xpi = crater%xlpx + i ypi = crater%ylpx + j + xbar = xpi * user%pix - crater%xl + ybar = ypi * user%pix - crater%yl + ! periodic boundary conditions call util_periodic(xpi,ypi,user%gridsize) newdem = surf(xpi,ypi)%dem - xbar = xpi * user%pix - crater%xl - ybar = ypi * user%pix - crater%yl - r = sqrt(xbar**2 + ybar**2) / crater%frad if (r > flr) then @@ -675,15 +668,14 @@ subroutine ejecta_texture(user,surf,crater,deltaMtot,inc,ejecta_dem) ! Add the base texture to the ejecta proportional to the thickness - do j = -inc,inc do i = -inc,inc xpi = indarray(1,i,j) ypi = indarray(2,i,j) - xbar = xpi * user%pix - crater%xl - ybar = ypi * user%pix - crater%yl + xbar = (crater%xlpx + i) * user%pix - crater%xl + ybar = (crater%ylpx + j) * user%pix - crater%yl r = sqrt(xbar**2 + ybar**2) / crater%frad phi = atan2(ybar,xbar) @@ -697,9 +689,7 @@ subroutine ejecta_texture(user,surf,crater,deltaMtot,inc,ejecta_dem) areafrac = areafrac * (1.0_DP - max(min(2._DP - r,1.0_DP),0.0_DP)) ! Blend in with the wall texture - ! Make the splat pattern - splatnoise = 0.0_DP do octave = 1,nsplat_octaves xysplat = (nsplats / PI) * freq ** (octave -1) / crater%fcrat @@ -795,7 +785,7 @@ subroutine crater_realistic_slope_texture(user,critical,inc,critarray) noise = noise + util_perlin_noise(xynoise * xbar + offset * rn(1), & xynoise * ybar + offset * rn(2)) * znoise end do - !write(*,*) i,j,noise + critarray(i,j) = max(critical * (1.0_DP + noise),0.0_DP) end do end do