diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index 885f3f6b..73740b8c 100644 --- a/src/crater/crater_populate.f90 +++ b/src/crater/crater_populate.f90 @@ -200,7 +200,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt end if ! Place crater onto the surface - if (crater%fcrat > domain%smallest_crater) then + if (crater%fcrat > domain%smallest_counted_crater) then call crater_emplace(user,surf,crater,domain,ejbmass) !call crater_mass_conservation(user,surf,crater) ! mass conservation is now done in crater_emplace diff --git a/src/crater/crater_record.f90 b/src/crater/crater_record.f90 index 3fbf7c0b..d603ca79 100644 --- a/src/crater/crater_record.f90 +++ b/src/crater/crater_record.f90 @@ -54,7 +54,7 @@ subroutine crater_record(user,surf,crater) call util_periodic(xpi,ypi,user%gridsize) ! record crater in available layer - call util_add_to_layer(user,surf(xpi,ypi),crater%fcrat,crater%xl,crater%yl) + call util_add_to_layer(user,surf(xpi,ypi),crater) end if end do diff --git a/src/crater/crater_tally_observed.f90 b/src/crater/crater_tally_observed.f90 index 5b612b07..17cf1d5a 100644 --- a/src/crater/crater_tally_observed.f90 +++ b/src/crater/crater_tally_observed.f90 @@ -37,13 +37,12 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o ! Internal variables integer(I4B) :: i,j,layer,n,m,craternum,obstot + type(cratertype) :: crater ! master list arrays (these contain information on every pixel that is considered a crater) integer(I4B),dimension(user%numlayers*(user%gridsize)**2) :: mlistind real(DP),dimension(:),allocatable :: mlist,melevation real(SP),dimension(:,:),allocatable :: mposlist - real(SP),dimension(:),allocatable :: moriginal_depth_list,mbaseline - real(SP),dimension(:),allocatable :: tmp_original_depth,tmp_current_depth real(SP),dimension(:),allocatable :: tmp_depthdiam,tmp_deviation_sigma logical,dimension(:),allocatable :: countable integer(I4B),dimension(:,:),allocatable:: mpxlist @@ -61,8 +60,7 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o 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 + real(DP) :: totavg,avgelev,profres,avgdis,rim,bowl,rad,baseline ! Executable code @@ -147,17 +145,14 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o allocate(totpix(tnum)) allocate(poslist(2,tnum)) allocate(tlist(tnum)) - allocate(tmp_original_depth(tnum)) - allocate(tmp_current_depth(tnum)) allocate(tmp_depthdiam(tnum)) - allocate(tmp_deviation_sigma(tnum)) allocate(countable(tnum)) countable = .false. nkilled = 0 onum = 0 !$OMP PARALLEL DEFAULT(PRIVATE) IF(tnum > INCPAR) & !$OMP SHARED(user,surf) & - !$OMP SHARED(maxpix,istart,iend,ind,mlist,mposlist,melevation,mpxlist,mlayerlist) & + !$OMP SHARED(domain,maxpix,istart,iend,ind,mlist,mposlist,melevation,mpxlist,mlayerlist) & !$OMP SHARED(tnum,totpix,poslist,tlist,tmp_depthdiam,countable) & !$OMP REDUCTION(+:nkilled) & !$OMP REDUCTION(+:onum) @@ -167,51 +162,55 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o ! 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 + 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))) 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_depthdiam(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 + crater%fcrat = tlist(craternum) + crater%fradpx = int(crater%fcrat * 0.5_DP / user%pix) + 1 + crater%xlpx = int(poslist(1,craternum) / user%pix) + crater%ylpx = int(poslist(2,craternum) / user%pix) + call crater_averages(user,surf,crater) + nrim = 0 + nbowl = 0 + bowl = 0.0_DP + rim = 0.0_DP + inc = ceiling(0.5_DP * crater%fcrat / user%pix * 1.25_DP) + do j = -inc, inc + do i = -inc, inc + rad = sqrt((i**2 + j**2)*1._DP) * user%pix / crater%fcrat + baseline = crater%melev + ((i * crater%xslp) + (j * crater%yslp)) * user%pix + if (rad < 0.1_DP) then + xpi = crater%xlpx + i + ypi = crater%ylpx + j + call util_periodic(xpi,ypi,user%gridsize) + bowl = bowl + surf(xpi,ypi)%dem - baseline + nbowl = nbowl + 1 + else if (rad > 0.45_DP) then + xpi = crater%xlpx + i + ypi = crater%ylpx + j + call util_periodic(xpi,ypi,user%gridsize) + rim = rim + surf(xpi,ypi)%dem - baseline + nrim = nrim + 1 + end if end do - rim = rim / nrim - bowl = bowl / nbowl - tmp_depthdiam(craternum) = (rim - bowl) / tlist(craternum) + end do + rim = rim / nrim + bowl = bowl / nbowl + tmp_depthdiam(craternum) = (rim - bowl) / crater%fcrat - if (tmp_depthdiam(craternum) > 0.05_DP) then - countable(craternum) = .true. - killable = .false. - else - countable(craternum) = .false. - killable = .true. - end if + if (tmp_depthdiam(craternum) > 0.03_DP) then + countable(craternum) = .true. + killable = .false. + else + countable(craternum) = .false. + killable = .true. + end if + if ((crater%fcrat < domain%smallest_counted_crater) .or. & + (totpix(craternum) < int(0.1_DP * 0.25_DP * PI * crater%fcrat / user%pix))) then + countable(craternum) = .false. + killable = .true. end if ! TESTING if (user%testtally) then