Skip to content

Commit

Permalink
Simpler crater tally
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Dec 5, 2016
1 parent 74c84bd commit 2e02ff7
Showing 1 changed file with 118 additions and 88 deletions.
206 changes: 118 additions & 88 deletions src/crater/crater_tally_observed.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,18 +55,22 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o
! tallying arrays (these contain only the number of craters on the grid)
real(SP),dimension(:,:),allocatable :: poslist
real(DP),dimension(:),allocatable :: tlist
real(DP),dimension(:),allocatable :: rimelevation,bowlelevation,bowldepth_orig
real(DP),dimension(:),allocatable :: dis,elev
integer(I4B),dimension(:),allocatable :: elevind,totpix,istart,iend
integer(I4B) :: tnum ! True number and observed number
integer(I4B) :: inc,xpi,ypi
real(DP) :: avgrim,avgbowldepth,avgbowldepth_orig,xk,Mkm1,xkdd,Mkm1dd
real(DP) :: xkdev,Mkm1dev,deviation,sigma
logical :: killable
integer(I4B) :: nrim,nbowl
integer(I4B) :: nprof,Ni
real(DP) :: totavg,avgelev,profres,avgdis,rim,bowl,rad

! Executable code


profres = 0.5_DP * user%pix

! First build up a master table of every crater for which a pixel still exists on the surface
mnum=0
mlistind = 0
Expand All @@ -84,17 +88,18 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o
end do
end do

allocate(mlist(mnum))
allocate(misrim(mnum))
allocate(melevation(mnum))
allocate(mbaseline(mnum))
allocate(mposlist(2,mnum))
allocate(mpxlist(2,mnum))
allocate(moriginal_depth_list(mnum))
allocate(mlayerlist(mnum))
allocate(ind(mnum))
allocate(istart(mnum))
allocate(iend(mnum))

allocate(mlist(max(mnum,1)))
allocate(misrim(max(mnum,1)))
allocate(melevation(max(mnum,1)))
allocate(mbaseline(max(mnum,1)))
allocate(mposlist(2,max(mnum,1)))
allocate(mpxlist(2,max(mnum,1)))
allocate(moriginal_depth_list(max(mnum,1)))
allocate(mlayerlist(max(mnum,1)))
allocate(ind(max(mnum,1)))
allocate(istart(max(mnum,1)))
allocate(iend(max(mnum,1)))

mlist = 0._DP
ind = 0
Expand Down Expand Up @@ -165,87 +170,114 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o
!$OMP SHARED(tnum,totpix,poslist,tlist,tmp_original_depth,tmp_current_depth,tmp_p,tmp_deviation_sigma,countable) &
!$OMP REDUCTION(+:nkilled) &
!$OMP REDUCTION(+:onum)
allocate(rimelevation(maxpix))
allocate(bowlelevation(maxpix))
allocate(bowldepth_orig(maxpix))
allocate(elevind(maxpix))
!allocate(dis(4*user%gridsize**2))
!allocate(elev(4*user%gridsize**2))
!allocate(elevind(4*user%gridsize**2))
! Here we go through the list of craters and take the measures of each one to
! be used by the tally subroutine in the next pass
!$OMP DO
do craternum=1,tnum
! This is the first pixel of this crater, so record the appropriate values
tlist(craternum) = mlist(ind(istart(craternum)))
poslist(:,craternum) = mposlist(:,ind(istart(craternum)))
tmp_original_depth(craternum) = 0._SP
tmp_current_depth(craternum) = 0._SP
totpix(craternum) = 1 + iend(craternum) - istart(craternum)
if (totpix(craternum) < int(0.1_DP * 0.25_DP * PI * tlist(craternum) / user%pix)) then
countable(craternum) = .false.
killable = .true.
tmp_p(craternum) = 0.0_DP
else
nrim = 0
nbowl = 0
nprof = 0
bowl = 0.0_DP
rim = 0.0_DP
inc = ceiling(0.5_DP * tlist(craternum) / user%pix * 1.25_DP)
do j = -inc, inc
do i = -inc, inc
rad = sqrt((i**2 + j**2)*1._DP) * user%pix / tlist(craternum)
if (rad < 0.1_DP) then
xpi = int(poslist(1,craternum) / user%pix) + i
ypi = int(poslist(2,craternum) / user%pix) + j
call util_periodic(xpi,ypi,user%gridsize)
bowl = bowl + surf(xpi,ypi)%dem
nbowl = nbowl + 1
else if (rad > 0.4_DP) then
xpi = int(poslist(1,craternum) / user%pix) + i
ypi = int(poslist(2,craternum) / user%pix) + j
call util_periodic(xpi,ypi,user%gridsize)
rim = rim + surf(xpi,ypi)%dem
nrim = nrim + 1
end if
end do
end do
rim = rim / nrim
bowl = bowl / nbowl
tmp_p(craternum) = (rim - bowl) / tlist(craternum)

! Get the average and std. dev. of rim and bowl elevations for this crater
rimelevation = 0._DP
bowlelevation = 0._DP
nrim = 0
nbowl = 0
do i = istart(craternum),iend(craternum)
! Find out if this pixel is a rim or a bowl, and add the appropriate tally
if (misrim(ind(i)) == 1) then
nrim = nrim + 1
rimelevation(nrim) = melevation(ind(i)) - mbaseline(ind(i))
else
nbowl = nbowl + 1
bowlelevation(nbowl) = melevation(ind(i)) - mbaseline(ind(i))
bowldepth_orig(nbowl) = moriginal_depth_list(ind(i))
end if
end do

! Count up all the rim and bowl pixels and calculate the current average rim-bowl height and shape deviation
if ((nbowl > 1) .and. (nrim > 0)) then
! Calculate mean
avgrim = sum(rimelevation(1:nrim)) / real(nrim, kind = DP)

! Sort the bowl pixels and only count the lowest (original) of the set by a fraction set by the user
call util_mrgrnk(bowldepth_orig(1:nbowl),elevind(1:nbowl))
ntot = min(ceiling(BOWLFRAC * nbowl),nbowl)
avgbowldepth = avgrim - bowlelevation(elevind(nbowl))
avgbowldepth_orig = bowldepth_orig(elevind(nbowl))
deviation = avgbowldepth - avgbowldepth_orig
sigma = 0.0_DP
do i = 2,nbowl
xk = avgrim - bowlelevation(elevind(nbowl + 1 - i))
Mkm1 = avgbowldepth
avgbowldepth = Mkm1 + (xk - Mkm1) / real(i,kind = DP)

xkdd = bowldepth_orig(elevind(nbowl + 1 - i))
Mkm1dd = avgbowldepth_orig
avgbowldepth_orig = Mkm1dd + (xkdd - Mkm1dd) / real(i,kind = DP)
! if ((i**2 + j**2) < inc**2) then
! xpi = int(poslist(1,craternum) / user%pix) + i
! ypi = int(poslist(2,craternum) / user%pix) + j
! call util_periodic(xpi,ypi,user%gridsize)
! nprof = nprof + 1
! elev(nprof) = surf(xpi,ypi)%dem
! dis(nprof) = sqrt((i**2 + j**2)*1.d0) * user%pix
! totavg = totavg + elev(nprof)
! end if
! end do
! end do
! totavg = totavg / nprof
! elev = elev - totavg
! call util_mrgrnk(dis(1:nprof),elevind(1:nprof))
! avgelev = 0.0_DP
! avgdis = 0.0_DP
! Ni = 1
! ntot = 0
! i = 0
! rim = 0.0_DP
! bowl = 0.0_DP
! nrim = 0
! nbowl = 0
! outer: do
! do
! i = i + 1
! if (i > nprof) exit outer
! if (dis(elevind(i)) > profres * Ni) then
! avgdis = avgdis / ntot / tlist(craternum)
! avgelev = avgelev / ntot / tlist(craternum)
! if (avgdis >= 0.4_DP) then
! rim = rim + avgelev
! nrim = nrim+ 1
! else if (avgdis <= 0.1_DP) then
! bowl = bowl + avgelev
! nbowl = nbowl + 1
! end if
! avgelev = elev(elevind(i))
! avgdis = dis(elevind(i))
! ntot = 1
! exit
! else
! avgelev = avgelev + elev(elevind(i))
! avgdis = avgdis + dis(elevind(i))
! ntot = ntot + 1
! end if
! end do
! Ni = Ni + 1
! end do outer

if (i == ntot) then ! Save the values for the bottom of the bowl
tmp_current_depth(craternum) = real(avgbowldepth,kind = SP)
tmp_original_depth(craternum) = real(avgbowldepth_orig,kind = SP)
end if

xkdev = xk - xkdd
Mkm1dev = deviation
deviation = Mkm1dev + (xkdev - Mkm1dev) / real(i,kind = DP)
sigma = sigma + (xkdev - Mkm1dev) * (xkdev - deviation) ! Method for computing running std. dev from
! Knuth _Art of Computer Programming_, Vol. 2 pg. 232 3rd ed
end do
sigma = sqrt(sigma / (nbowl - 1))
tmp_deviation_sigma(craternum) = real(sigma,kind = SP)
else ! This crater is uncountable
tmp_current_depth(craternum) = 1._SP
tmp_original_depth(craternum) = huge(tmp_original_depth(craternum))
tmp_deviation_sigma(craternum) = huge(tmp_deviation_sigma(craternum))
if (tmp_p(craternum) > 0.05_DP) then
countable(craternum) = .true.
killable = .false.
else
countable(craternum) = .false.
killable = .true.
end if
end if
! TESTING
if (user%testtally) then
countable = .true.
killable = .false.
end if

! Use the calibrated crater count subroutine to determine if this is a countable crater or not
if (totpix(craternum) > SMALLESTCOUNTABLE) then
call crater_tally_calibrated_count(user,tlist(craternum),tmp_current_depth(craternum),&
tmp_original_depth(craternum),tmp_deviation_sigma(craternum),&
countable(craternum),killable,tmp_p(craternum))
else
countable(craternum) = .false.
killable = .true.
end if
if (killable) then ! Obliterate the crater from the record
do i=istart(craternum),iend(craternum)
call util_remove_from_layer(surf(mpxlist(1,ind(i)),mpxlist(2,ind(i))),mlayerlist(ind(i)))
Expand All @@ -254,13 +286,11 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o
end if
if (countable(craternum)) onum = onum + 1
!deallocate(csurf)
end do
!$OMP END DO
deallocate(rimelevation)
deallocate(bowlelevation)
deallocate(bowldepth_orig)
deallocate(elevind)
!$OMP END PARALLEL
end do
!$OMP END DO
!deallocate(dis)
!deallocate(elev)
!$OMP END PARALLEL

! Bin the observed craters if required
if (present(obsdist)) then
Expand Down

0 comments on commit 2e02ff7

Please sign in to comment.