Skip to content

Commit

Permalink
changes to ejecta files since the SVN-GitHub switchover
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed May 10, 2017
1 parent 47e43b2 commit b8c4e80
Show file tree
Hide file tree
Showing 5 changed files with 478 additions and 167 deletions.
120 changes: 102 additions & 18 deletions src/crater/crater_subpixel_diffusion.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@
subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval,kdiffin)
use module_globals
use module_util
use module_ejecta
use module_crater, EXCEPT_THIS_ONE => crater_subpixel_diffusion
implicit none

Expand All @@ -36,11 +37,16 @@ subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval,kdiff
! Internal variables
real(DP),dimension(0:user%gridsize + 1,0:user%gridsize + 1) :: cumulative_elchange,kdiff
integer(I4B),dimension(2,0:user%gridsize + 1,0:user%gridsize + 1) :: indarray
integer(I4B) :: i,j,k,xpi,ypi,n,ntot,ioerr
integer(I4B) :: i,j,k,l,m,xpi,ypi,n,ntot,ioerr,inc,xi,xf,yi,yf
integer(I4B) :: maxhits = 1
real(DP) :: lamleft,dburial,lambda,Kbar,diam,radius,Area,avgejc,sf
real(DP),dimension(domain%pnum) :: dN,lambda_regolith,Kbar_regolith,lambda_bedrock,Kbar_bedrock

real(DP),dimension(2) :: rn
real(DP) :: superlen,rayfrac
type(cratertype) :: crater
real(DP),dimension(:,:),allocatable :: ejdistribution
integer(I4B),dimension(:,:),allocatable :: ejisray


! Create box for soften calculation (will be no bigger than the grid itself)
do j = 0,user%gridsize + 1
Expand Down Expand Up @@ -68,13 +74,6 @@ subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval,kdiff

Kbar_bedrock(i) = 0.84_DP * (0.5_DP * nflux(1,i))**4 ! Intrinsic degradation function
Kbar_regolith(i) = 0.84_DP * (0.5_DP * nflux(2,i))**4
open(unit=22,file='soften.dat',status='old',iostat=ioerr)
if (ioerr == 0) then
read(22,*) sf
close(22)
Kbar_bedrock(i) = Kbar_bedrock(i) + sf * (0.5_DP * nflux(1,i))**4 ! Extra softening
Kbar_regolith(i) = Kbar_regolith(i) + sf * (0.5_DP * nflux(2,i))**4
end if

if (user%dosoftening) then
Kbar_bedrock(i) = Kbar_bedrock(i) + user%soften_factor * (0.5_DP * nflux(1,i))**(user%soften_slope)
Expand Down Expand Up @@ -111,27 +110,112 @@ subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval,kdiff

! Now add the superdomain contribution to crater degradation
if (user%dosoftening) then
crater%ejdis = user%gridsize * user%pix * 0.5_DP
crater%continuous = crater%ejdis / DISEJB
radius = (crater%continuous / 2.348_DP)**(1.0_DP / 1.006_DP)
crater%frad = radius
crater%rad = radius
inc = nint(min(crater%ejdis / user%pix, crater%frad * user%ejecta_truncation / user%pix)) + 1
crater%xl = user%gridsize * user%pix * 0.5_DP
crater%yl = user%gridsize * user%pix * 0.5_DP
crater%xlpx = nint(crater%xl / user%pix)
crater%ylpx = nint(crater%yl / user%pix)
allocate(ejdistribution(-inc:inc,-inc:inc))
allocate(ejisray(-inc:inc,-inc:inc))
call ejecta_ray_pattern(user,surf,crater,inc,-inc,inc,-inc,inc,ejdistribution)
do j = -inc,inc
do i = -inc,inc
if (ejdistribution(i,j) > 0.0001_DP) then
ejisray(i,j) = 1
else
ejisray(i,j) = 0
end if
end do
end do
rayfrac = PI * inc**2 / (1.0_DP * count(ejisray == 1))
deallocate(ejisray,ejdistribution)

avgejc = sum(surf%ejcov) / user%gridsize**2
do i = 1,domain%pnum
if ((PI * (user%soften_size * nflux(1,i) * 0.5_DP)**2 <= (user%gridsize)**2 ).and.&
(PI * (user%soften_size * nflux(2,i) * 0.5_DP)**2 <= (user%gridsize)**2 )) cycle
do l = 1,domain%pnum
if ((PI * (user%soften_size * nflux(1,l) * 0.5_DP)**2 <= (user%gridsize)**2 ).and.&
(PI * (user%soften_size * nflux(2,l) * 0.5_DP)**2 <= (user%gridsize)**2 )) cycle

dburial = EXFAC * 0.5_DP * nflux(1,n)
if (avgejc > dburial) then
radius = nflux(2,i) * 0.5_DP
radius = nflux(2,l) * 0.5_DP
else
radius = nflux(1,i) * 0.5_DP
radius = nflux(1,l) * 0.5_DP
end if

dN(i) = nflux(3,i) * user%interval * finterval
dN(l) = nflux(3,l) * user%interval * finterval
Area = min(PI * (user%soften_size * radius)**2 , 4 * PI * user%trad**2)
Area = max(Area - (user%pix * user%gridsize)**2,0.0_DP)
if (Area == 0.0_DP) cycle
lambda = dN(i) * Area
superlen = 0.5_DP * (sqrt(Area) - user%pix * user%gridsize)

lambda = dN(l) * Area
if (lambda == 0.0_DP) cycle
k = util_poisson(lambda)

if (k > 0) kdiff = kdiff + k * user%soften_factor / (PI * user%soften_size**2) * radius**(user%soften_slope - 2.0_DP)
do m = 1, k

! generate random crater on the superdomain
! Set up size of crater and ejecta blanket
crater%frad = radius
crater%rad = radius
crater%continuous = 2.348_DP * crater%frad**(1.006_DP)
crater%ejdis = DISEJB * crater%continuous
inc = nint(min(crater%ejdis / user%pix, crater%frad * user%ejecta_truncation / user%pix)) + 1

! find the x and y position of the crater
call random_number(rn)
crater%xl = real(superlen * (2 * rn(1) - 1.0_DP), kind=SP)
crater%yl = real(superlen * (2 * rn(2) - 1.0_DP), kind=SP)

if (crater%xl > 0.0_SP) then
crater%xl = crater%xl + user%pix * user%gridsize
else
crater%xl = crater%xl - user%pix * user%gridsize
end if
if (crater%yl > 0.0_SP) then
crater%yl = crater%yl + user%pix * user%gridsize
else
crater%yl = crater%yl - user%pix * user%gridsize
end if

crater%xlpx = nint(crater%xl / user%pix)
crater%ylpx = nint(crater%yl / user%pix)

xi = 1 - crater%xlpx
xf = user%gridsize - crater%xlpx
yi = 1 - crater%ylpx
yf = user%gridsize - crater%ylpx

allocate(ejdistribution(xi:xf,yi:yf))
allocate(ejisray(xi:xf,yi:yf))
! Now generate ray pattern
call ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,ejdistribution)
do j = yi,yf
do i = xi,xf
if (ejdistribution(i,j) > 0.0001_DP) then
ejisray(i,j) = 1
else
ejisray(i,j) = 0
end if
end do
end do

do j = 1, user%gridsize
do i = 1, user%gridsize
xpi = i - crater%xlpx
ypi = j - crater%ylpx
if ((abs(xpi) > inc) .or. (abs(ypi) > inc)) cycle
if (ejisray(xpi,ypi) == 0) cycle
kdiff(i,j) = kdiff(i,j) + &
rayfrac * user%soften_factor / (PI * user%soften_size**2) * crater%frad**(user%soften_slope - 2.0_DP)
end do
end do
deallocate(ejdistribution,ejisray)
end do
end do
end if

Expand Down
100 changes: 100 additions & 0 deletions src/crater/crater_tally_calibrated_count.f90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,100 @@
!**********************************************************************************************************************************
!
! Unit Name : crater_tally_calibrated_count
! Unit Type : subroutine
! Project : CTEM
! Language : Fortran 2003
!
! Description : Uses a human-calibrated counting model to determine the countability of a given crater
!
!
! Input
! Arguments :
!
! Output
! Arguments :
!
!
! Notes :
!
!**********************************************************************************************************************************
subroutine crater_tally_calibrated_count(user,diameter,current_depth,original_depth,deviation_sigma,&
countable,killable,p)
use module_globals
use module_util
use module_crater, EXCEPT_THIS_ONE => crater_tally_calibrated_count
implicit none

! Arguments
type(usertype),intent(in) :: user
real(DP),intent(in) :: diameter
real(SP),intent(in) :: current_depth,original_depth,deviation_sigma
logical,intent(out) :: countable,killable
real(SP),intent(out) :: p

! Internal variables
real(DP) :: Rd,Rsig,pd,psig,a,b,c,d,Dtran,complex_correction

Rd = current_depth/original_depth
Rsig = log10(deviation_sigma/current_depth)

! Executable code
select case(user%countingmodel)
case("FASSETT")
a = -0.0237520259849034_DP
b = 0.208739267563691_DP
c = 0.992603771065306_DP
if (Rd < 0.081895370584072_DP) then
pd = 0._DP
else if (Rd > 0.9158503765541_DP) then
pd = 1._DP
else
pd = a + b * Rd + c * Rd**2
end if
pd = min(max(pd, 0.0_DP), 1._DP)

a = -0.0856409196596715_DP
b = -1.16657321813299_DP
psig = a + b * Rsig
psig = min(max(psig, 0.0_DP), 1._DP)

p = real(min(pd,psig),kind = SP)
p = min(max(p, 0.0_SP), 1._SP)
case("HOWL")
a = -0.48_DP
b = -0.55_DP
c = 8.0_DP
d = -4.5_DP
!p = (Rsig - a * log(Rd)) / b - c * (diameter / user%pix)**d - 0.5_DP
p = 1.0_DP
if (Rd < 0.2_DP) p = 0.0_DP
end select

select case (user%mat)
case ("ROCK")
Dtran = 2 * SIMCOMKS * user%gaccel**SIMCOMPS
case ("ICE")
Dtran = 2 * SIMCOMKI * user%gaccel**SIMCOMPI
case default
Dtran = 2 * SIMCOMKS * user%gaccel**SIMCOMPS ! Use rock as default just in case it dien't get set
end select

complex_correction = min(Dtran / diameter, 1.0_DP)

if (p > 0.50_SP * complex_correction) then ! Kill the crater if its p-score is too low
countable = .true.
killable = .false.
else
countable = .false.
killable = .true.
end if

! TESTING
if (user%testtally) then
countable = .true.
killable = .false.
end if

return
end subroutine crater_tally_calibrated_count

Loading

0 comments on commit b8c4e80

Please sign in to comment.