From 0d6094796034ebf7d322db474d8ec85b7fd7ffc4 Mon Sep 17 00:00:00 2001 From: David Minton Date: Thu, 26 Sep 2019 12:49:32 -0400 Subject: [PATCH] Re-wrote tally to explicitly look at degradation state to conform to the equations in Minton et al. (2019) --- src/Makefile.am | 6 +- src/crater/crater_degradation_function.f90 | 2 +- src/crater/crater_get_degradation_state.f90 | 89 +++++++++++++++++++++ src/crater/crater_tally_observed.f90 | 82 ++++--------------- src/crater/crater_visibility.f90 | 50 ++++++++++++ src/crater/module_crater.f90 | 26 ++++++ 6 files changed, 186 insertions(+), 69 deletions(-) create mode 100644 src/crater/crater_get_degradation_state.f90 create mode 100644 src/crater/crater_visibility.f90 diff --git a/src/Makefile.am b/src/Makefile.am index 1d53fddc..b3b8df93 100755 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -6,9 +6,9 @@ bin_PROGRAMS = CTEM #AM_FCFLAGS = -O3 -p -g -openmp -debug all -traceback -CB -assume byterecl -m64 -heap-arrays -FR #gfortran optimized flags -AM_FCFLAGS = -O3 -fopenmp -ffree-form -g -fbounds-check -fbacktrace +#AM_FCFLAGS = -O3 -fopenmp -ffree-form -g -fbounds-check -fbacktrace #gfortran debug flags -#AM_FCFLAGS = -O0 -g -fopenmp -fbounds-check -Wall -Warray-bounds -Warray-temporaries -Wimplicit-interface -ffree-form +AM_FCFLAGS = -O0 -g -fopenmp -fbounds-check -Wall -Warray-bounds -Warray-temporaries -Wimplicit-interface -ffree-form CTEM_SOURCES = globals/module_globals.f90\ util/module_util.f90\ @@ -82,6 +82,8 @@ crater/crater_subpixel_diffusion.f90\ crater/crater_make_list.f90\ crater/crater_critical_slope.f90\ crater/crater_degradation_function.f90\ +crater/crater_visibility.f90\ +crater/crater_get_degradation_state.f90\ init/init_domain.f90\ init/init_dist.f90\ init/init_surf.f90\ diff --git a/src/crater/crater_degradation_function.f90 b/src/crater/crater_degradation_function.f90 index 37f6e511..f4fd9d57 100644 --- a/src/crater/crater_degradation_function.f90 +++ b/src/crater/crater_degradation_function.f90 @@ -36,8 +36,8 @@ function crater_degradation_function(user,r) result(Kd) delta = 1.0_DP A = user%Kd1 * r_break**(alpha_1) / (1_DP + (1.d0 / delta) * (alpha_1 - alpha_2) / alpha_1)**2 - Kd = A * (r / r_break)**(alpha_1) * (0.5_DP * (1.0_DP + (r / r_break)**(1.0_DP / delta)))**((alpha_2 - alpha_1 ) / delta) + return end function crater_degradation_function diff --git a/src/crater/crater_get_degradation_state.f90 b/src/crater/crater_get_degradation_state.f90 new file mode 100644 index 00000000..37356104 --- /dev/null +++ b/src/crater/crater_get_degradation_state.f90 @@ -0,0 +1,89 @@ +! Project : CTEM +! Language : Fortran 2003 +! +! Description : Determines the current degradation state of the crater +! +! +! Input +! Arguments : +! +! Output +! Arguments : +! +! +! Notes : +! +!********************************************************************************************************************************** + + +function crater_get_degradation_state(user,surf,crater,dd) result(Kval) + use module_globals + use module_util + use module_crater, EXCEPT_THIS_ONE => crater_get_degradation_state + implicit none + + ! Arguments + type(usertype),intent(in) :: user + type(surftype),dimension(:,:),intent(in) :: surf + type(cratertype),intent(in) :: crater + real(SP),intent(out) :: dd + ! Return value + real(DP) :: Kval + + ! Internal variables + integer(I4B) :: inc,xpi,ypi,i,j + integer(I4B) :: nrim,nbowl,nouter + real(DP) :: rim,bowl,outer,rad,baseline,ddinit + + ! Definitions of crater dimensions + real(DP),parameter :: OCUTOFF = 5.5e-2_DP + real(DP),parameter :: RIMDI = 1.0_DP + real(DP),parameter :: RIMDO = 1.2_DP + real(DP),parameter :: BOWLD = 0.2_DP + real(DP),parameter :: OUTERD = 2.0_DP + + ! Fit for complex crater diffusion model + real(DP),parameter :: Bconst = 8.0_DP + real(DP),parameter :: dDpikeC = 1.725_DP ! CTEM version of the Pike (1977) depth-to-diam constant + real(DP),parameter :: ddsimple = 0.21362307923564000_DP ! CTEM initial depth-to-diameter of simple craterws + + ddinit = min(ddsimple, dDpikeC * (crater%fcrat * 1e-3_DP)**(0.301_DP - 1._DP)) + + inc = int(OUTERD * crater%frad / user%pix) + 1 + + + nrim = 0 + nbowl = 0 + nouter = 0 + bowl = 0.0_DP + rim = 0.0_DP + outer = 0.0_DP + + do j = -inc, inc + do i = -inc, inc + rad = sqrt((i**2 + j**2)*1._DP) * user%pix / crater%frad + baseline = crater%melev + ((i * crater%xslp) + (j * crater%yslp)) * user%pix + xpi = crater%xlpx + i + ypi = crater%ylpx + j + call util_periodic(xpi,ypi,user%gridsize) + if (rad <= BOWLD) then + bowl = bowl + surf(xpi,ypi)%dem - baseline + nbowl = nbowl + 1 + else if ((rad >= RIMDI).and.(rad < RIMDO)) then + rim = rim + surf(xpi,ypi)%dem - baseline + nrim = nrim + 1 + else if (rad > RIMDO) then + outer = outer + surf(xpi,ypi)%dem - baseline + nouter = nouter + 1 + end if + end do + end do + rim = rim / nrim + bowl = bowl / nbowl + outer = outer / nouter + dd = (rim - bowl) / crater%fcrat + Kval = - 1._SP / Bconst * log(dd / ddinit) + return + +end function crater_get_degradation_state + diff --git a/src/crater/crater_tally_observed.f90 b/src/crater/crater_tally_observed.f90 index af8211eb..e163bc73 100644 --- a/src/crater/crater_tally_observed.f90 +++ b/src/crater/crater_tally_observed.f90 @@ -45,6 +45,7 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o real(SP),dimension(:,:),allocatable :: mposlist real(SP),dimension(:),allocatable :: tmp_depthdiam logical,dimension(:),allocatable :: countable + real(DP),dimension(:),allocatable :: Kval integer(I4B),dimension(:,:),allocatable:: mpxlist integer(I4B),dimension(:),allocatable :: ind,mlayerlist integer(I4B) :: q,mnum,imnum,maxpix,npix,ntot @@ -55,17 +56,7 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o real(DP),dimension(:),allocatable :: dis,elev integer(I4B),dimension(:),allocatable :: totpix,istart,iend integer(I4B) :: tnum ! True number and observed number - integer(I4B) :: inc,xpi,ypi - logical :: killable - integer(I4B) :: nrim,nbowl,nouter - real(DP) :: rim,bowl,outer,rad,baseline,dd - ! Counting parameters from Howl study - real(DP),parameter :: DDCUTOFF = 5.0e-2_DP - real(DP),parameter :: OCUTOFF = 5.5e-2_DP - real(DP),parameter :: RIMDI = 1.0_DP - real(DP),parameter :: RIMDO = 1.2_DP - real(DP),parameter :: BOWLD = 0.2_DP - real(DP),parameter :: OUTERD = 2.0_DP + ! Executable code @@ -148,17 +139,18 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o allocate(tlist(tnum)) allocate(tmp_depthdiam(tnum)) allocate(countable(tnum)) + allocate(Kval(tnum)) countable = .false. nkilled = 0 onum = 0 + ! 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 PARALLEL DEFAULT(PRIVATE) IF(tnum > INCPAR) & !$OMP SHARED(user,surf) & !$OMP SHARED(domain,maxpix,istart,iend,ind,mlist,mposlist,mpxlist,mlayerlist) & - !$OMP SHARED(tnum,totpix,poslist,tlist,tmp_depthdiam,countable) & + !$OMP SHARED(tnum,totpix,poslist,tlist,tmp_depthdiam,countable,Kval) & !$OMP REDUCTION(+:nkilled) & !$OMP REDUCTION(+:onum) - ! 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 @@ -170,71 +162,28 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o crater%fradpx = int(crater%frad / 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 - nouter = 0 - bowl = 0.0_DP - rim = 0.0_DP - outer = 0.0_DP - inc = int(OUTERD * crater%frad / user%pix) + 1 - do j = -inc, inc - do i = -inc, inc - rad = sqrt((i**2 + j**2)*1._DP) * user%pix / crater%frad - baseline = crater%melev + ((i * crater%xslp) + (j * crater%yslp)) * user%pix - xpi = crater%xlpx + i - ypi = crater%ylpx + j - call util_periodic(xpi,ypi,user%gridsize) - if (rad <= BOWLD) then - bowl = bowl + surf(xpi,ypi)%dem - baseline - nbowl = nbowl + 1 - else if ((rad >= RIMDI).and.(rad < RIMDO)) then - rim = rim + surf(xpi,ypi)%dem - baseline - nrim = nrim + 1 - else if (rad > RIMDO) then - outer = outer + surf(xpi,ypi)%dem - baseline - nouter = nouter + 1 - end if - end do - end do - rim = rim / nrim - bowl = bowl / nbowl - outer = outer / nouter - tmp_depthdiam(craternum) = (rim - bowl) / crater%fcrat - !if (crater%fcrat < 20e3_DP) then - ! dd = DDCUTOFF - !else - ! dd = DDCUTOFF - (crater%fcrat - 20e3_DP) * 2e-7 - !end if - dd = min(DDCUTOFF, 24.897 * crater%fcrat ** (-0.632545)) + call crater_averages(user,surf,crater) - if (((tmp_depthdiam(craternum) > dd).and.((outer - rim) / crater%fcrat < OCUTOFF)).and.& - (nrim /= 0).and.(nbowl /= 0).and.(nouter /= 0)) 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. + else + Kval(craternum) = crater_get_degradation_state(user,surf,crater,tmp_depthdiam(craternum)) + countable(craternum) = crater_visibility(user,crater,Kval(craternum)) end if - ! TESTING - if (user%testtally) then - countable = .true. - killable = .false. + + if (user%testtally) then ! used for testing the tally system + countable(craternum) = .true. end if - if (killable) then ! Obliterate the crater from the record + + if (.not.countable(craternum)) 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))) end do nkilled = nkilled + 1 end if if (countable(craternum)) onum = onum + 1 - !deallocate(csurf) end do !$OMP END DO !$OMP END PARALLEL @@ -293,6 +242,7 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o deallocate(tmp_depthdiam) deallocate(ind) deallocate(countable) + deallocate(Kval) return end subroutine crater_tally_observed diff --git a/src/crater/crater_visibility.f90 b/src/crater/crater_visibility.f90 new file mode 100644 index 00000000..234b7bcc --- /dev/null +++ b/src/crater/crater_visibility.f90 @@ -0,0 +1,50 @@ +! Project : CTEM +! Language : Fortran 2003 +! +! Description : Calculates the visibility function from Minton et al. (2019) +! +! +! Input +! Arguments : +! +! Output +! Arguments : +! +! +! Notes : +! +!********************************************************************************************************************************** + + +function crater_visibility(user,crater,Kval) result(iscountable) + use module_globals + use module_util + use module_crater, EXCEPT_THIS_ONE => crater_visibility + implicit none + + ! Arguments + type(usertype),intent(in) :: user + type(cratertype),intent(in) :: crater + real(DP),intent(in) :: Kval + + ! Result variable + logical :: iscountable + + ! Internal variables + real(DP) :: Kv + real(DP),parameter :: Kv1 = 0.17 + real(DP),parameter :: gam = 2.0 + + ! Counting parameter for simple craters from Bryan Howl's study in Minton et al. (2019) + + Kv = Kv1 * (crater%frad)**gam + if (Kval >= kv) then + iscountable = .false. + else + iscountable = .true. + end if + + return + +end function crater_visibility + diff --git a/src/crater/module_crater.f90 b/src/crater/module_crater.f90 index b96fb0dd..719b96b8 100644 --- a/src/crater/module_crater.f90 +++ b/src/crater/module_crater.f90 @@ -271,6 +271,32 @@ function crater_degradation_function(user,r) result(Kd) end function crater_degradation_function end interface + interface + function crater_get_degradation_state(user,surf,crater,dd) result(Kval) + use module_globals + implicit none + ! Arguments + type(usertype),intent(in) :: user + type(surftype),dimension(:,:),intent(in) :: surf + type(cratertype),intent(in) :: crater + real(SP),intent(out) :: dd + ! Return value + real(DP) :: Kval + end function crater_get_degradation_state + end interface + + interface + function crater_visibility(user,crater,Kval) result(iscountable) + use module_globals + implicit none + ! Arguments + type(usertype),intent(in) :: user + type(cratertype),intent(in) :: crater + real(DP),intent(in) :: Kval + ! Result variable + logical :: iscountable + end function crater_visibility + end interface end module