From e55e59fb8710069035d3ac16c412b9564773b432 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 29 Apr 2019 11:55:10 -0400 Subject: [PATCH] Updated with code used in Minton et al. (2019) --- src/crater/crater_populate.f90 | 13 +- src/crater/crater_soften.f90 | 117 ++++--- src/crater/crater_soften_accumulate.f90 | 6 +- src/crater/crater_subpixel_diffusion.f90 | 284 ++++++++-------- src/ejecta/ejecta_emplace.f90 | 148 +++++--- src/ejecta/ejecta_ray_pattern.f90 | 411 ++++++++++++++++------- src/ejecta/ejecta_table_define.f90 | 2 +- src/ejecta/module_ejecta.f90 | 5 +- src/globals/module_globals.f90 | 41 ++- src/io/io_input.f90 | 30 +- src/io/io_updatePbar.f90 | 24 +- src/io/module_ejecta.f90 | 156 +++++++++ 12 files changed, 843 insertions(+), 394 deletions(-) create mode 100644 src/io/module_ejecta.f90 diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index 1eb783cb..1229cd9f 100644 --- a/src/crater/crater_populate.f90 +++ b/src/crater/crater_populate.f90 @@ -73,6 +73,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt logical :: makecrater real(DP),dimension(user%gridsize,user%gridsize) :: kdiff TARGET :: surf + integer(I4B) :: oldpbarpos ! ejecta blanket array type(ejbtype),dimension(EJBTABSIZE) :: ejb ! Ejecta blanket lookup table @@ -124,10 +125,13 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt domain%tallycoverage = 0 domain%subpixelcoverage = 0 kdiff = 0.0_DP + pbarpos = 0 + call io_updatePbar("") + oldpbarpos = 0 do while (icrater < ntotcrat) makecrater = .true. icrater = icrater + 1 - pbarpos = ceiling(real(icrater) / real(ntotcrat) * PBARRES) + pbarpos = nint(real(icrater) / real(ntotcrat) * PBARRES) ! generate random crater call crater_generate(user,crater,domain,prod,production_list,vdist,surf) if (user%testflag) write(*,*) 'Dcrat = ',crater%fcrat @@ -224,7 +228,12 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt call util_sort_layer(user,surf,crater) vistrue = vistrue + 1 nsincetally = nsincetally + 1 - if (.not.user%testflag) call io_updatePbar("") + if (.not.user%testflag) then + if (pbarpos /= oldpbarpos) then + call io_updatePbar("") + oldpbarpos = pbarpos + end if + end if !if (user%docrustal_thinning) call crust_thin(user,surf,crater,domain,mdepth) diff --git a/src/crater/crater_soften.f90 b/src/crater/crater_soften.f90 index d6943b76..e20d4f1a 100644 --- a/src/crater/crater_soften.f90 +++ b/src/crater/crater_soften.f90 @@ -31,42 +31,46 @@ subroutine crater_soften(user,surf,crater,domain) type(domaintype),intent(in) :: domain ! Internal variables - real(DP),dimension(:,:),allocatable :: cumulative_elchange,kappat + real(DP),dimension(:,:),allocatable :: cumulative_elchange,kdiff integer(I4B),dimension(:,:,:),allocatable :: indarray + real(DP),dimension(0:user%gridsize+1,0:user%gridsize+1) :: big_cel,big_kdiff + integer(I4B),dimension(2,0:user%gridsize+1,0:user%gridsize+1) :: big_indarray integer(I4B) :: maxhits - integer(I4B) :: inc,incsq,N,xpi,ypi,iradsq,i,j,ss,ioerr - real(DP) :: xp,yp,fradsq,xbar,ybar,areafrac,kappatmax,sf + integer(I4B) :: inc,incsq,N,xpi,ypi,iradsq,i,j,ss,ioerr,bigi,bigj + real(DP) :: xp,yp,fradsq,xbar,ybar,areafrac,kdiffmax,sf - !kappatmax = user%soften_factor * crater%frad**user%soften_slope - open(unit=22,file='soften.dat',status='old',iostat=ioerr) - if (ioerr /= 0) return - read(22,*) sf - close(22) - !inc = max(min(int(crater%fradpx + 2),PBCLIM * user%gridsize),2) - inc = crater%fradpx + 2 + kdiffmax = user%Kd1 * crater%frad**user%psi + inc = int(min(user%fe * crater%frad / user%pix,SQRT2 * user%gridsize)) + 3 maxhits = (1 + (inc / (user%gridsize / 2)))**2 crater%maxinc = max(crater%maxinc,inc) - incsq = inc**2 - kappatmax = (sf / PI) * crater%frad**2 + incsq = (inc - 2)**2 - allocate(kappat(-inc:inc,-inc:inc)) + allocate(kdiff(-inc:inc,-inc:inc)) allocate(indarray(2,-inc:inc,-inc:inc)) allocate(cumulative_elchange(-inc:inc,-inc:inc)) - cumulative_elchange = 0._DP + cumulative_elchange = 0.0_DP indarray = inc - kappat = 0.0_DP + kdiff = 0.0_DP ! Loop over affected matrix area do j=-inc,inc ! Do the loop in pixel space do i=-inc,inc + ! find elevation and grid point + xpi = crater%xlpx + i + ypi = crater%ylpx + j + + ! periodic boundary conditions + call util_periodic(xpi,ypi,user%gridsize) + + indarray(1,i,j) = xpi + indarray(2,i,j) = ypi + ! find distance from crater center iradsq = i*i + j*j - if (iradsq <= incsq) then - ! find elevation and grid point - xpi = crater%xlpx + i - ypi = crater%ylpx + j + ! We set the diffusion constant to be proportional to the fraction of pixel area covered by the interior of the crater + !if (iradsq <= incsq) then ! Find distance from crater center to current pixel center in real space xp = xpi * user%pix yp = ypi * user%pix @@ -74,31 +78,68 @@ subroutine crater_soften(user,surf,crater,domain) xbar = xp - crater%xl ybar = yp - crater%yl - ! periodic boundary conditions - call util_periodic(xpi,ypi,user%gridsize) - - indarray(1,i,j) = xpi - indarray(2,i,j) = ypi - - ! We set the diffusion constant to be proportional to the fraction of pixel area covered by the interior of the crater - areafrac = util_area_intersection(crater%frad*0.98_DP,xbar,ybar,user%pix) + areafrac = util_area_intersection(user%fe * crater%frad,xbar,ybar,user%pix) - kappat(i,j) = kappatmax * areafrac ! This is the extra per-crater diffusion required to match equilibrium - end if + kdiff(i,j) = kdiffmax * areafrac ! Apply user-defined diffusive degradation to degradation region + !end if end do end do !end area loopover + !write(*,*) 'crater_soften: ',crater%frad,maxval(kdiff) + if (2 * inc + 1 < user%gridsize) then + write(*,*) 'crater inc = ',inc + write(*,*) 'crater kdiff ' + write(*,*) maxval(kdiff) + call util_diffusion_solver(user,surf,2 * inc + 1,indarray,kdiff,cumulative_elchange,maxhits) + write(16,*) crater%frad + do j = -inc,inc + do i = -inc,inc + xpi = indarray(1,i,j) + ypi = indarray(2,i,j) + write(16,*) i,j,xpi,ypi,kdiff(i,j),surf(xpi,ypi)%dem,cumulative_elchange(i,j) + !surf(xpi,ypi)%dem = surf(xpi,ypi)%dem + cumulative_elchange(i,j) + !surf(xpi,ypi)%ejcov = max(surf(xpi,ypi)%ejcov + cumulative_elchange(i,j),0.0_DP) + end do + end do + else + maxhits = 1 + big_kdiff = 0.0_DP + big_cel = 0.0_DP + do bigj = 0,user%gridsize + 1 + do bigi = 0,user%gridsize + 1 + xpi = bigi + ypi = bigj + call util_periodic(xpi,ypi,user%gridsize) + big_indarray(1,bigi,bigj) = xpi + big_indarray(2,bigi,bigj) = ypi + end do + end do + + do j = -inc,inc + do i = -inc,inc + xpi = indarray(1,i,j) + ypi = indarray(2,i,j) + big_kdiff(xpi,ypi) = big_kdiff(xpi,ypi) + kdiff(i,j) + end do + end do + + big_kdiff(0,0) = big_kdiff(user%gridsize,user%gridsize) + big_kdiff(user%gridsize + 1,user%gridsize + 1) = big_kdiff(1,1) + big_kdiff(0,:) = big_kdiff(user%gridsize,:) + big_kdiff(:,0) = big_kdiff(:,user%gridsize) + big_kdiff(user%gridsize + 1,:) = big_kdiff(1,:) + big_kdiff(:,user%gridsize + 1) = big_kdiff(:,1) - call util_diffusion_solver(user,surf,2 * inc + 1,indarray,kappat,cumulative_elchange,maxhits) - do j = -inc,inc - do i = -inc,inc - xpi = indarray(1,i,j) - ypi = indarray(2,i,j) - surf(xpi,ypi)%dem = surf(xpi,ypi)%dem + cumulative_elchange(i,j) - surf(xpi,ypi)%ejcov = max(surf(xpi,ypi)%ejcov + cumulative_elchange(i,j),0.0_DP) + call util_diffusion_solver(user,surf,user%gridsize + 2,big_indarray,big_kdiff,big_cel,maxhits) + do j = 1,user%gridsize + do i = 1,user%gridsize + surf(i,j)%dem = surf(i,j)%dem + big_cel(i,j) + surf(i,j)%ejcov = max(surf(i,j)%ejcov + big_cel(i,j),0.0_DP) + end do end do - end do - deallocate(kappat,indarray,cumulative_elchange) + end if + + deallocate(kdiff,indarray,cumulative_elchange) return end subroutine crater_soften diff --git a/src/crater/crater_soften_accumulate.f90 b/src/crater/crater_soften_accumulate.f90 index 0817ff4c..0e8c24d3 100644 --- a/src/crater/crater_soften_accumulate.f90 +++ b/src/crater/crater_soften_accumulate.f90 @@ -41,9 +41,9 @@ subroutine crater_soften_accumulate(user,surf,crater,domain,kdiff) hit = .false. - kappatmax = user%soften_factor / (PI * user%soften_size**2) * crater%frad**(user%soften_slope - 2.0_DP) + kappatmax = user%Kd1 * crater%frad**(user%psi) - inc = int(min(user%soften_size * crater%frad / user%pix, user%gridsize / SQRT2)) + 2 + inc = int(min(user%fe * crater%frad / user%pix, user%gridsize / SQRT2)) + 2 crater%maxinc = max(crater%maxinc,inc) fradsq = crater%frad**2 incsq = inc**2 @@ -93,7 +93,7 @@ subroutine crater_soften_accumulate(user,surf,crater,domain,kdiff) ! periodic boundary conditions call util_periodic(xpi,ypi,user%gridsize) - areafrac = util_area_intersection(user%soften_size * crater%frad,xbar,ybar,user%pix) + areafrac = util_area_intersection(user%fe * crater%frad,xbar,ybar,user%pix) areafrac = areafrac * craterhole(xpi,ypi) if (.not.hit(xpi,ypi)) then diff --git a/src/crater/crater_subpixel_diffusion.f90 b/src/crater/crater_subpixel_diffusion.f90 index 86b53973..7f8c27a3 100644 --- a/src/crater/crater_subpixel_diffusion.f90 +++ b/src/crater/crater_subpixel_diffusion.f90 @@ -37,15 +37,18 @@ 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,l,m,xpi,ypi,n,ntot,ioerr,inc,xi,xf,yi,yf + integer(I4B) :: i,j,l,xpi,ypi,k,ktot,kk,ioerr,inc,incsq,iradsq,xi,xf,yi,yf,crat,imin,imax,jmin,jmax + integer(I8B) :: m,N 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 + real(DP) :: dburial,lambda,dKdN,diam,radius,Area,avgejc + real(DP) :: dN,diam_regolith,diam_bedrock,RR + real(DP),dimension(3) :: rn + real(DP) :: superlen,rayfrac,cutout,fe,fd type(cratertype) :: crater - real(DP),dimension(:,:),allocatable :: ejdistribution + real(DP),dimension(:,:),allocatable :: diffdistribution,ejdistribution integer(I4B),dimension(:,:),allocatable :: ejisray + real(DP) :: xbar,ybar,dD,xp,yp,areafrac,krad,supersize,lrad + integer(I8B),dimension(user%gridsize,user%gridsize) :: Ngrid ! Create box for soften calculation (will be no bigger than the grid itself) @@ -60,165 +63,158 @@ subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval,kdiff end do end do - ntot = 1 - ! calculate the subpixel diffusion probability function - Area = user%pix**2 + avgejc = sum(surf%ejcov) / domain%parea - do i = 1,domain%pnum - if ((nflux(1,i) > domain%smallest_crater).and.(nflux(2,i) > domain%smallest_crater)) exit - ntot = i - dN(i) = nflux(3,i) * user%interval * finterval + fe = FEPROX + fd = user%ejecta_truncation + if (user%dosoftening) fe = user%fe - lambda_bedrock(i) = dN(i) * Area - lambda_regolith(i) = dN(i) * Area + ! Generate both the subpixel and superdomain diffusive degradation + do k = 1,domain%pnum - 1 + dN = nflux(3,k) * user%interval * finterval - 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 - - if (user%dosoftening) then - Kbar_bedrock(i) = Kbar_bedrock(i) + user%soften_factor * (0.5_DP * nflux(1,i))**(user%soften_slope) - Kbar_regolith(i) = Kbar_regolith(i) + user%soften_factor * (0.5_DP * nflux(2,i))**(user%soften_slope) + diam_bedrock = nflux(1,k) + diam_regolith = nflux(2,k) + + dburial = 0.5_DP * EXFAC * max(diam_bedrock,diam_regolith) + if (dburial > avgejc) then + diam = diam_bedrock + else + diam = diam_regolith end if + radius = 0.5_DP * diam + + if ((fd * diam < user%pix) .or. (dN * PI * (fe * radius)**2 > 0.1_DP)) then + !Do the average degradation per pixel for the subpixel component + if (user%dosoftening) then + ! User-defined degradation function + dKdN = user%Kd1 * PI * user%fe**2 * (radius)**(2.0_DP + user%psi) / domain%parea + else + !Empirically-derived "intrinsic" degradation function from proximal ejecta redistribution + dKdN = KD1PROX * PI * FEPROX**2 * (radius)**(2.0_DP + PSIPROX) / domain%parea + end if - end do + lambda = dN * domain%parea - do j = 1,user%gridsize - do i = 1,user%gridsize - do n = 1,ntot - dburial = EXFAC * 0.5_DP * nflux(1,n) - if (surf(i,j)%ejcov > dburial) then - lambda = lambda_regolith(n) - Kbar = Kbar_regolith(n) - diam = nflux(2,n) - else - lambda = lambda_bedrock(n) - Kbar = Kbar_bedrock(n) - diam = nflux(1,n) - end if - if (diam > domain%smallest_crater) exit - k = util_poisson(lambda) - kdiff(i,j) = kdiff(i,j) + k * Kbar / Area + ! Don't parallelize the random + do j = 1,user%gridsize + do i = 1,user%gridsize + Ngrid(i,j) = util_poisson(lambda) + end do end do - end do - end do - kdiff(0,0) = kdiff(user%gridsize,user%gridsize) - kdiff(user%gridsize + 1,user%gridsize + 1) = kdiff(1,1) - kdiff(0,:) = kdiff(user%gridsize,:) - kdiff(:,0) = kdiff(:,user%gridsize) - kdiff(user%gridsize + 1,:) = kdiff(1,:) - kdiff(:,user%gridsize + 1) = kdiff(:,1) - - ! 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 + + !$OMP PARALLEL DO DEFAULT(PRIVATE) & + !$OMP SHARED(user) & + !$OMP SHARED(Ngrid,dKdN,kdiff,dN,lambda) + do j = 1,user%gridsize + do i = 1,user%gridsize + if ((Ngrid(i,j)) == 0) cycle + if ((Ngrid(i,j) > 0).and.(Ngrid(i,j) == Ngrid(i,j))) then + kdiff(i,j) = kdiff(i,j) + dKdN * Ngrid(i,j) + else ! In case of integer overflow + kdiff(i,j) = kdiff(i,j) + dKdN * lambda + end if + end do 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 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,l) * 0.5_DP - else - radius = nflux(1,l) * 0.5_DP + !$OMP END PARALLEL DO + + else if (user%dosoftening) then + ! Do the degradation as individual circles + + superlen = fd * diam + domain%side + cutout = 0.0_DP + crater%continuous = RCONT * radius**(EXPCONT) + if (diam > domain%smallest_crater) then + if (.not.user%superdomain) exit + ! Superdomain craters + cutout = domain%side + crater%continuous end if + Area = superlen**2 + if (diam < domain%biggest_crater) Area = Area - cutout**2 - 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 - superlen = 0.5_DP * (sqrt(Area) - user%pix * user%gridsize) - - lambda = dN(l) * Area - if (lambda == 0.0_DP) cycle - k = util_poisson(lambda) - 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 + lambda = dN * Area + dD = nflux(2,k+1) - nflux(2,k) + + N = util_poisson(lambda) + do m = 1, N 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 + if (diam < domain%smallest_crater) then + ! Subpixel craters with large degradation region + crater%xl = real(superlen * (rn(1) - 0.5_DP) + 0.5_DP * domain%side, kind=SP) + crater%yl = real(superlen * (rn(2) - 0.5_DP) + 0.5_DP * domain%side, kind=SP) + else if (.not.user%superdomain) then + cycle else - crater%yl = crater%yl - user%pix * user%gridsize + ! Superdomain craters + crater%xl = real((superlen - cutout) * (rn(1) - 0.5_DP), kind=SP) + if (crater%xl > 0.0_DP) crater%xl = crater%xl + cutout + crater%yl = real((superlen - cutout) * (rn(2) - 0.5_DP), kind=SP) + if (crater%yl > 0.0_DP) crater%yl = crater%yl + cutout 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) + crater%frad = 0.5_DP * (diam + dD * rn(3)) + crater%continuous = RCONT * crater%frad**(EXPCONT) + krad = fd * crater%frad + + dKdN = user%Kd1 * crater%frad**(user%psi) + inc = int(krad / user%pix) + 2 + incsq = inc**2 + + xpi = max(crater%xlpx - inc,1) + imin = xpi - crater%xlpx + xpi = min(crater%xlpx + inc,user%gridsize) + imax = xpi - crater%xlpx + ypi = max(crater%ylpx - inc,1) + jmin = ypi - crater%ylpx + ypi = min(crater%ylpx + inc,user%gridsize) + jmax = ypi - crater%ylpx + + + allocate(diffdistribution(imin:imax,jmin:jmax)) + allocate(ejdistribution(imin:imax,jmin:jmax)) + call ejecta_ray_pattern(user,surf,crater,inc,imin,imax,jmin,jmax,diffdistribution,ejdistribution) + ! Loop over affected matrix area + !$OMP PARALLEL DO DEFAULT(PRIVATE) IF(inc > INCPAR) & + !$OMP SHARED(jmin,jmax,imin,imax,kdiff,dKdN,krad,diffdistribution,ejdistribution) & + !$OMP SHARED(crater,user,surf) + do j = jmin,jmax + do i = imin,imax + xpi = crater%xlpx + i + ypi = crater%ylpx + j + kdiff(xpi,ypi) = kdiff(xpi,ypi) + dKdN * diffdistribution(i,j) + !TEMP + xp = xpi * user%pix + yp = ypi * user%pix + lrad = sqrt((xp - crater%xl)**2 + (yp - crater%yl)**2) + surf(xpi,ypi)%ejcov = surf(xpi,ypi)%ejcov + ejdistribution(i,j) & + * 0.14_DP * crater%frad**(0.74_DP) * (lrad / crater%frad)**(-3.0_DP) end do end do - deallocate(ejdistribution,ejisray) + !$OMP END PARALLEL DO + deallocate(diffdistribution,ejdistribution) end do - end do - end if - + end if + + end do + + kdiff(0,0) = kdiff(user%gridsize,user%gridsize) + kdiff(user%gridsize + 1,user%gridsize + 1) = kdiff(1,1) + kdiff(0,:) = kdiff(user%gridsize,:) + kdiff(:,0) = kdiff(:,user%gridsize) + kdiff(user%gridsize + 1,:) = kdiff(1,:) + kdiff(:,user%gridsize + 1) = kdiff(:,1) + + write(*,*) + write(*,*) 'avgkdiff = ',sum(kdiff) / (user%gridsize + 2)**2 / finterval + write(*,*) + open(unit=55,file='avgkdiff.dat',status='unknown',position='append') + write(55,*) sum(kdiff) / (user%gridsize + 2)**2 / finterval + close(55) + call util_diffusion_solver(user,surf,user%gridsize + 2,indarray,kdiff,cumulative_elchange,maxhits) do j = 1,user%gridsize diff --git a/src/ejecta/ejecta_emplace.f90 b/src/ejecta/ejecta_emplace.f90 index 837d197e..16c4a9c1 100644 --- a/src/ejecta/ejecta_emplace.f90 +++ b/src/ejecta/ejecta_emplace.f90 @@ -97,12 +97,13 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) real(DP) :: lrad,lradsq,cdepth integer(I4B),parameter :: MAXLOOP = 100 ! Maximum number of times to loop the ejecta angle correction calculation integer(I4B) :: xpi,ypi,i,j,k,n,inc,incsq,iradsq,idistorted,jdistorted - real(DP) :: xp,yp,fradsq,fradpxsq,radsq,ebh,ejdissq,ejbmass,fmasscons + real(DP) :: xp,yp,fradsq,fradpxsq,radsq,ebh,ejdissq,ejbmass,fmasscons,areafrac,xbar,ybar,krad,kdiffmax real(DP),dimension(:,:),allocatable :: cumulative_elchange,big_cumulative_elchange,kdiff,big_kdiff,cel,big_cel integer(I4B),dimension(:,:,:),allocatable :: indarray,big_indarray - real(DP),dimension(:,:),allocatable :: ejdistribution - integer(I4B) :: bigi,bigj,maxhits,nin,nnot + real(DP),dimension(:,:),allocatable :: ejdistribution,diffdistribution + integer(I4B) :: bigi,bigj,maxhits,nin,nnot,dradsq character(len=MESSAGESIZE) :: message ! message for the progress bar + ! Ray mixing model variables real(DP) :: dsc @@ -112,8 +113,10 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) ! Ejecta pattern distortion parameters real(DP) :: distance,erad,craterslope,landslope,baseline,lrange,frac,ejheight,ebh0,maxdistance + real(DP) :: maxslp real(DP) :: vsq, ejtheta - integer(I4B) :: klo,ind + integer(I4B) :: ind,klo + ! Executable code !call regolith_melt_zone(user,crater,crater%imp,crater%impvel,rm,dm) @@ -128,15 +131,26 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) inc = max(min(nint(min(PI * user%trad / user%pix, min(crater%ejdis, crater%frad * user%ejecta_truncation / user%pix))) + 1, & user%gridsize - nint(crater%continuous / user%pix)),1) maxdistance = inc * user%pix + ! Increase the box a bit to take into account possible ejecta pattern distortion due to topography - incsq = inc**2 inc = ceiling(inc * 1.5_DP) + + if (user%dosoftening) then + krad = user%ejecta_truncation * crater%frad + kdiffmax = user%Kd1 * crater%frad**(user%psi) + dradsq = int(krad / user%pix) + 3 + inc = max(inc,dradsq) + dradsq = dradsq**2 + end if + crater%maxinc = max(crater%maxinc,inc) radsq = crater%rad**2 fradsq = crater%frad**2 fradpxsq = crater%fradpx**2 ejdissq = crater%ejdis**2 + incsq = inc**2 + if (inc >= user%gridsize / 2) then if (user%testflag) then write(*,*) 'Big ejecta: fcrat =',crater%fcrat, ' Ej/S =',(crater%ejdispx*user%pix)/domain%side, ' Ejrim =', crater%ejrim @@ -147,25 +161,27 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) endif allocate(cumulative_elchange(-inc:inc,-inc:inc)) allocate(cel(-inc:inc,-inc:inc)) - allocate(kdiff(-inc-1:inc+1,-inc-1:inc+1)) + allocate(kdiff(-inc:inc,-inc:inc)) allocate(indarray(2,-inc:inc,-inc:inc)) allocate(ejdistribution(-inc:inc,-inc:inc)) + allocate(diffdistribution(-inc:inc,-inc:inc)) cumulative_elchange = 0.0_DP kdiff = 0.0_DP - indarray = inc ! initialize this array to point to a corner (this should have 0 elevation change since we're only doing work + indarray = inc - 1 ! initialize this array to point to a corner (this should have 0 elevation change since we're only doing work ! within a circle of radius irad ! *************************** Continuous Ejecta Formula *****************************! - call ejecta_ray_pattern(user,surf,crater,inc,-inc,inc,-inc,inc,ejdistribution) + call ejecta_ray_pattern(user,surf,crater,inc,-inc,inc,-inc,inc,diffdistribution,ejdistribution) ejbmass = 0.0_DP nin = 0 nnot = 0 -! !$OMP PARALLEL DO DEFAULT(PRIVATE) IF(inc > INCPAR) & -! !$OMP SHARED(user,domain,crater,surf,ejb,ejtble,mvrld,mvrldsc,nrays,n2,nef,ef,n2f,n1f,nfrays,rayf) & -! !$OMP SHARED(inc,incsq,ejdissq,fradsq,radsq,indarray,cumulative_elchange,rn) -! !$OMP REDUCTION(+:massray,massej,massrayef) + + !!$OMP PARALLEL DO DEFAULT(PRIVATE) IF(inc > INCPAR) & + !!$OMP SHARED(user,domain,crater,surf,ejb,ejtble) & + !!$OMP SHARED(inc,incsq) & + !!$OMP SHARED(cumulative_elchange,kdiff,kdiffmax,indarray,ejdistribution,diffdistribution) do j = -inc,inc do i = -inc,inc ! find distance from crater center @@ -176,9 +192,6 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) ! Find distance from crater center to current pixel center in real space xp = xpi * user%pix yp = ypi * user%pix - - lradsq = (crater%xl - xp)**2 + (crater%yl - yp)**2 - lrad = sqrt(lradsq) ! periodic boundary conditions call util_periodic(xpi,ypi,user%gridsize) @@ -186,10 +199,14 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) indarray(1,i,j) = xpi indarray(2,i,j) = ypi + lradsq = (crater%xl - xp)**2 + (crater%yl - yp)**2 + lrad = sqrt(lradsq) ! Estimate ejecta pattern distortion due to target surface angle and topography ! This must be done iteratively because the ejection distance and ejection angle vary distance = lrad + maxslp = -huge(maxslp) + klo = int((log(lrad) - log(crater%rad)) / domain%ejbres) do n = 1,MAXLOOP call ejecta_interpolate(crater,domain,distance,ejb,ejtble,ebh,vsq=vsq,theta=ejtheta,erad=erad,melt=melt) if ((n > 1).and.((abs(ebh0 - ebh) / ebh0) < domain%small)) exit @@ -200,6 +217,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) baseline = ((i * crater%xslp) + (j * crater%yslp)) * user%pix craterslope = atan(baseline / lrad) + if (craterslope > maxslp) maxslp = craterslope ejheight = erad * sin(craterslope) + crater%melev @@ -229,33 +247,66 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) iradsq = idistorted**2 + jdistorted**2 if ((iradsq > incsq).or.(distance <= crater%rad)) cycle - ebh = ejdistribution(idistorted,jdistorted) * ebh + ! we need to cut a hole out from the inside of the crater + xbar = xpi * user%pix - crater%xl + ybar = ypi * user%pix - crater%yl - if (user%doregotrack .and. ebh>1.0e-8) then - call regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,lrad,ebh,rm) - end if + areafrac = (1.0_DP - util_area_intersection(crater%rad,xbar,ybar,user%pix)) - cumulative_elchange(i,j) = cumulative_elchange(i,j) + ebh - ejbmass = ejbmass + ebh - if (ebh > 0.0_DP) then - kdiff(i,j) = user%soften_factor / (PI * user%soften_size**2) * crater%frad ** (user%soften_slope - 2.0) - nin = nin + 1 - else - nnot = nnot + 1 + ebh = areafrac * ejdistribution(idistorted,jdistorted) * ebh + cumulative_elchange(i,j) = areafrac * cumulative_elchange(i,j) + ebh + + if (user%dosoftening) then + ! Do extra diffusive degradation over ejecta region + areafrac = (1.0_DP - util_area_intersection(crater%frad,xbar,ybar,user%pix)) + kdiff(i,j) = areafrac * diffdistribution(idistorted,jdistorted) * kdiffmax end if + + + !if (user%doregotrack .and. ebh>1.0e-8_DP) then + ! call regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,lrad,ebh,rm) + !end if + + end do end do -! !$OMP END PARALLEL DO - kdiff = kdiff * (nnot + nin) / (1.0_DP * nin) + !!$OMP END PARALLEL DO + ejbmass = sum(cumulative_elchange) + + ! Create buffer to prevent infinite hole bug + kdiff(-inc,:) = 0.0_DP + kdiff(inc,:) = 0.0_DP + kdiff(:,-inc) = 0.0_DP + kdiff(:,inc) = 0.0_DP + + cumulative_elchange(-inc,:) = 0.0_DP + cumulative_elchange(inc,:) = 0.0_DP + cumulative_elchange(:,-inc) = 0.0_DP + cumulative_elchange(:,inc) = 0.0_DP ! Do mass conservation by adjusting ejecta thickness fmasscons = (-deltaMtot)/ ejbmass cumulative_elchange = cumulative_elchange * fmasscons - cel = 0.0_DP maxhits = 1 ! Create box for soften calculation (will be no bigger than the grid itself) if (2 * inc + 1 < user%gridsize) then + + ! extra soften calculation + if (user%dosoftening) then + cel = 0.0_DP + call util_diffusion_solver(user,surf,2 * inc + 1,indarray,kdiff,cel,maxhits) + do j = -inc,inc + do i = -inc,inc + xpi = indarray(1,i,j) + ypi = indarray(2,i,j) + surf(xpi,ypi)%dem = surf(xpi,ypi)%dem + cel(i,j) + surf(xpi,ypi)%ejcov = max(surf(xpi,ypi)%ejcov + cel(i,j),0.0_DP) + end do + end do + end if + call ejecta_soften(user,surf,2 * inc + 1,indarray,cumulative_elchange) + ! Add the ejecta back to the DEM do j = -inc,inc do i = -inc,inc @@ -265,22 +316,15 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) surf(xpi,ypi)%ejcov = max(surf(xpi,ypi)%ejcov + cumulative_elchange(i,j), 0.0_DP) end do end do - call util_diffusion_solver(user,surf,2 * inc + 1,indarray,kdiff,cel,maxhits) - do j = -inc,inc - do i = -inc,inc - xpi = indarray(1,i,j) - ypi = indarray(2,i,j) - surf(xpi,ypi)%dem = surf(xpi,ypi)%dem + cel(i,j) - surf(xpi,ypi)%ejcov = max(surf(xpi,ypi)%ejcov + cel(i,j),0.0_DP) - end do - end do + + else ! Ejecta wraps around the grid. ! We will therefore send in the whole grid with the total ejecta thickness added to each pixel allocate(big_cumulative_elchange(0:user%gridsize+1,0:user%gridsize+1)) allocate(big_cel(0:user%gridsize+1,0:user%gridsize+1)) allocate(big_indarray(2,0:user%gridsize+1,0:user%gridsize+1)) allocate(big_kdiff(0:user%gridsize+1,0:user%gridsize+1)) - big_kdiff = 0.0_DP + do bigj = 0,user%gridsize + 1 do bigi = 0,user%gridsize + 1 xpi = bigi @@ -292,6 +336,8 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) end do end do + big_kdiff = 0.0_DP + do j = -inc,inc do i = -inc,inc xpi = indarray(1,i,j) @@ -300,26 +346,34 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) big_kdiff(xpi,ypi) = big_kdiff(xpi,ypi) + kdiff(i,j) end do end do + + if (user%dosoftening) then + big_cel = 0.0_DP + call util_diffusion_solver(user,surf,user%gridsize + 2,big_indarray,big_kdiff,big_cel,maxhits) + + do j = 1,user%gridsize + do i = 1,user%gridsize + surf(i,j)%dem = surf(i,j)%dem + big_cel(i,j) + surf(i,j)%ejcov = max(surf(i,j)%ejcov + big_cel(i,j),0.0_DP) + end do + end do + end if + call ejecta_soften(user,surf,user%gridsize + 2,big_indarray,big_cumulative_elchange) + do i = 1,user%gridsize do j = 1,user%gridsize surf(i,j)%dem = surf(i,j)%dem + big_cumulative_elchange(i,j) surf(i,j)%ejcov = max(surf(i,j)%ejcov + big_cumulative_elchange(i,j),0.0_DP) end do end do - call util_diffusion_solver(user,surf,user%gridsize + 2,big_indarray,big_kdiff,big_cel,maxhits) - do j = 1,user%gridsize - do i = 1,user%gridsize - surf(i,j)%dem = surf(i,j)%dem + big_cel(i,j) - surf(i,j)%ejcov = max(surf(i,j)%ejcov + big_cel(i,j),0.0_DP) - end do - end do + deallocate(big_cumulative_elchange,big_indarray,big_kdiff,big_cel) end if !if (user%doregotrack) call regolith_rays(user,crater,domain,ejtble,ejb) - deallocate(cumulative_elchange,indarray,ejdistribution,kdiff,cel) + deallocate(cumulative_elchange,indarray,diffdistribution,ejdistribution,kdiff,cel) return end subroutine ejecta_emplace diff --git a/src/ejecta/ejecta_ray_pattern.f90 b/src/ejecta/ejecta_ray_pattern.f90 index b376735e..2e15dc73 100644 --- a/src/ejecta/ejecta_ray_pattern.f90 +++ b/src/ejecta/ejecta_ray_pattern.f90 @@ -31,7 +31,7 @@ !*** !********************************************************************************************************************************** -subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,ejdistribution) +subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,diffdistribution,ejdistribution) use module_globals use module_util use module_io @@ -43,14 +43,15 @@ subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,ejdistribution) ! Arguments type(usertype),intent(in) :: user type(surftype),dimension(:,:),intent(in) :: surf - type(cratertype),intent(in) :: crater + type(cratertype),intent(inout) :: crater integer(I4B),intent(in) :: inc,xi,xf,yi,yf + real(DP),dimension(xi:xf,yi:yf),intent(out) :: diffdistribution real(DP),dimension(xi:xf,yi:yf),intent(out) :: ejdistribution ! Internal variables - integer(I4B) :: nrays,i,j,k,n,nef,incsq,iradsq,xpi,ypi - real(DP) :: frac,mef,lrad,lradsq,xp,yp,binres - real(DP),dimension(1) :: rn + integer(I4B) :: nrays,i,j,k,n,nef,incsq,iradsq,xpi,ypi,ejpxsq + real(DP) :: frac,mef,lrad,lradsq,xp,yp,binres,areafrac,xbar,ybar + real(DP) :: rn real(DP) :: theta, lradp, maxdistance real(DP), parameter :: n1 = 4.0_DP real(DP) :: n2, mag @@ -58,6 +59,10 @@ subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,ejdistribution) real(DP),dimension(:),allocatable :: numinray,totnum real(DP),dimension(:),allocatable :: mefarray + + real(DP) :: C1,C2,p ! Fits to fe vs fd equation for the new ray pattern + real(DP) :: rmin,rmax,r + ! Flowery ray variables integer(I4B) :: nfrays real(DP) :: n1f, n2f, magf, lradf, rayf,rayavg @@ -65,10 +70,167 @@ subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,ejdistribution) real(DP), parameter :: b = 0.120621 !0.143 ! based on Jake's crater rays mapping studies! real(DP) :: mvrld ! median value of ray length distribution real(DP) :: mvrldsc ! median value of ray length distribution scaled by continuous ejecta ext/allent + logical :: ej + real(DP),dimension(Nraymax) :: thetari + + !TEMPORARY + + if (user%dorays) then + do i = 1,Nraymax + thetari(i) = 2 * pi * i / Nraymax + end do + call shuffle(thetari) ! randomize the ray pattern + + call random_number(rn) ! randomize the ray orientation + rmax = user%ejecta_truncation + rmin = crater%continuous / crater%frad + crater%fe = 10.0_DP ! Estimate the equivalent degradation radius + !ejdistribution = 0.0_DP + !diffdistribution = 0.0_DP + !!$OMP PARALLEL DO DEFAULT(PRIVATE) & + !!$OMP SHARED(user,crater) & + !!$OMP SHARED(xi,xf,yi,yf,rn,diffdistribution,ejdistribution,thetari,rmin) + do j = yi,yf + do i = xi,xf + xpi = crater%xlpx + i + ypi = crater%ylpx + j + + ! Find distance from crater center to current pixel center in real space + xp = xpi * user%pix + yp = ypi * user%pix + + xbar = xp - crater%xl + ybar = yp - crater%yl + + areafrac = util_area_intersection(user%ejecta_truncation * crater%frad,xbar,ybar,user%pix) + r = sqrt(xbar**2 + ybar**2) / crater%frad + theta = mod(atan2(ybar,xbar) + pi + rn * 2 * pi,2 * pi) + diffdistribution(i,j) = areafrac * pattern(theta,r,rmin,rmax,thetari,.false.) + ejdistribution(i,j) = areafrac * pattern(theta,r,rmin,rmax,thetari,.true.) + end do + end do + !!$OMP END PARALLEL DO + + + else + !Do simple circular region + incsq = inc**2 + ejdistribution = 0.0_DP + diffdistribution = 0.0_DP + !$OMP PARALLEL DO DEFAULT(PRIVATE) & + !$OMP SHARED(user,crater) & + !$OMP SHARED(inc,incsq,xi,xf,yi,yf,diffdistribution,ejdistribution) + do j = yi,yf + do i = xi,xf + iradsq = i*i + j*j + + if (iradsq < incsq) then + + xpi = crater%xlpx + i + ypi = crater%ylpx + j + + ! Find distance from crater center to current pixel center in real space + xp = xpi * user%pix + yp = ypi * user%pix + + xbar = xp - crater%xl + ybar = yp - crater%yl + areafrac = util_area_intersection(user%ejecta_truncation * crater%frad,xbar,ybar,user%pix) ! uniform circular + diffdistribution(i,j) = areafrac + ejdistribution(i,j) = areafrac + end if + end do + end do + !$OMP END PARALLEL DO + + end if + return + contains + + pure function ray(theta,thetar,r,n,w) result(ans) + implicit none + real(DP) :: ans + real(DP),intent(in) :: theta,thetar,r,w + integer(I4B),intent(in) :: n + real(DP) :: thetap,thetapp,a,b,c,dtheta + + c = w / r + b = thetar + dtheta = min(2*pi - abs(theta - b),abs(theta - b)) + a = sqrt(2 * pi) / (n * c * erf(pi / (2 *sqrt(2._DP) * c))) + ans = a * exp(-dtheta**2 / (2 * c**2)) + + return + end function ray + + pure function pattern(theta,r,rmin,rmax,thetari,ej) result(ans) + implicit none + real(DP) :: ans + real(DP),intent(in) :: r,rmin,rmax,theta + real(DP),dimension(:),intent(in) :: thetari + logical,intent(in) :: ej + real(DP) :: a,c + real(DP) :: thetar,rw,rw0,rw1 + real(DP) :: f,rtrans,length,rpeak,minray,FF + integer(I4B) :: n,i + + + minray = rmin * 3 - call random_number(rn) + if (r > rmax) then + ans = 0._DP + else if (r < 1.0_DP) then + if (ej) then + ans = 1.0_DP + else + ans = 0.0_DP + end if + else + rw0 = rmin * pi / Nraymax / 1 + rw1 = 2 * pi / Nraymax + rw = rw0 * (1._DP - (1.0_DP - rw1 / rw0) * exp(1._DP - (r / rmin)**2)) + n = max(min(floor((Nraymax**rayp - (Nraymax**rayp - 1) * log(r/minray) / log(rray/minray))**(1._DP/rayp)),Nraymax),1) ! Exponential decay of ray number with distance + ans = 0._DP + rtrans = r - 1.0_DP + c = rw / r + a = sqrt(2 * pi) / (n * c * erf(pi / (2 *sqrt(2._DP) * c))) + do i = 1,Nraymax + length = minray * exp(log(rray/minray) * ((Nraymax - i + 1)**rayp - 1_DP) / ((Nraymax**rayp - 1))) + rpeak = (length - 1_DP) * 0.5_DP + if (ej) then + FF = 1.0_DP + if (r > length) then + f = 0.0_DP + else + f = a + end if + else + FF = rayfmult * (20 / rmax)**(0.5_DP) * 0.25_DP + f = FF * fpeak * (rtrans / rpeak)**rayq * exp(1._DP / rayq * (1.0_DP - (rtrans / rpeak)**rayq)) + end if + ans = ans + ray(theta,thetari(i),r,n,rw) * f / a + end do + end if + + + end function pattern + + subroutine shuffle(a) + real(DP), intent(inout) :: a(:) + integer :: i, randpos + real(DP) :: r,temp + + do i = size(a), 2, -1 + call random_number(r) + randpos = int(r * i) + 1 + temp = a(randpos) + a(randpos) = a(i) + a(i) = temp + end do - !^^^^^^^^^^^^^^^^^^^^^^^^ + end subroutine shuffle + + !^^^^^^^^^^^^^^^^^^^^^^^^ ! *************************** Superformula Ray Model ************************************! ! *************************** Part I. Spoke and Skinny Ray Model ************************************! @@ -81,128 +243,133 @@ subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,ejdistribution) ! The average number of rays is about 10. The minimum and maximum number is 6 and 14 respectively. ! Determine parameters of Ray model based on Gielis's superformula - ! There are three parameters regarding to our desired ray shape: n1, n2, m +! ! There are three parameters regarding to our desired ray shape: n1, n2, m ! m: the repeating part of formula, and it controls the number of rays (arms). ! n1: the default value is set as 4.0 in our rays case - ! n2: n2 and n1 combining together is to control the slenderness of a ray, in general, n1/n2 is always smaller than 1 in our cases. - ! The smaller the ratio, the skinnier a ray. - - nrays = nint(8 * rn(1)) + 6 - mvrld = crater%ejdis - mvrldsc = mvrld / crater%continuous - n2 = 8.0_DP * ( log10(mvrldsc/1.5_DP) / log10(2.0_DP) ) + 2.0_DP - ! *************************** Part II. Flowery Ray Model ************************************! - nfrays = 20 - n1f = 1.0_DP - n2f = 0.5_DP - rayf = 4.0 !7.5_DP - - isray = 0.0_DP - incsq = inc**2 - maxdistance = inc * user%pix +! ! n2: n2 and n1 combining together is to control the slenderness of a ray, in general, n1/n2 is always smaller than 1 in our cases. +!! ! The smaller the ratio, the skinnier a ray. + +! nrays = nint(8 * rn(1)) + 6 +! mvrld = crater%ejdis +! mvrldsc = mvrld / crater%continuous +! n2 = 8.0_DP * ( log10(mvrldsc/1.5_DP) / log10(2.0_DP) ) + 2.0_DP +! ! *************************** Part II. Flowery Ray Model ************************************! +! nfrays = 20 +! n1f = 1.0_DP +! n2f = 0.5_DP +! rayf = 4.0 !7.5_DP +! +! isray = 0.0_DP +! incsq = inc**2 +! maxdistance = inc * user%pix ! Mass enhancement factor array - ! Don't let the bin resolution get too small for small craters. - binres = min(10.0_DP,crater%ejdis / 10) * user%pix - nef = ceiling(crater%ejdis / binres) ! bin size for enhancement factor calculation - allocate(numinray(nef),totnum(nef),mefarray(nef)) - - totnum = 0 - numinray = 0 - do j = yi,yf - do i = xi,xf - iradsq = i*i + j*j - - xpi = crater%xlpx + i - ypi = crater%ylpx + j - - ! Find distance from crater center to current pixel center in real space - xp = xpi * user%pix - yp = ypi * user%pix - - lradsq = (crater%xl - xp)**2 + (crater%yl - yp)**2 - lrad = sqrt(lradsq) - - if ((iradsq > incsq).or.(lrad <= crater%rad)) cycle - - ! Sum up the number of pixels in this bin for the mass enhancement calculation - k = ceiling(lrad / binres) - totnum(k) = totnum(k) + 1.0_DP - - ! Ray pattern setup - ! Take average of five points in the pixel to "feather" the edges of the ray - rayavg = 0.0_DP - do n = 1, 5 - select case(n) - case(1) - theta = atan2(j * 1._DP,i * 1._DP) + 2.0_DP * PI - case(2) - theta = atan2(j * 1._DP - 0.5_DP,i * 1._DP - 0.5_DP) + 2.0_DP * PI - lradsq = (crater%xl - xp - 0.5_DP * user%pix)**2 + (crater%yl - yp - 0.5_DP * user%pix)**2 - lrad = sqrt(lradsq) - case(3) - theta = atan2(j * 1._DP - 0.5_DP,i * 1._DP + 0.5_DP) + 2.0_DP * PI - lradsq = (crater%xl - xp - 0.5_DP * user%pix)**2 + (crater%yl - yp + 0.5_DP * user%pix)**2 - lrad = sqrt(lradsq) - case(4) - theta = atan2(j * 1._DP + 0.5_DP,i * 1._DP + 0.5_DP) + 2.0_DP * PI - lradsq = (crater%xl - xp + 0.5_DP * user%pix)**2 + (crater%yl - yp + 0.5_DP * user%pix)**2 - lrad = sqrt(lradsq) - case(5) - theta = atan2(j * 1._DP + 0.5_DP,i * 1._DP - 0.5_DP) + 2.0_DP * PI - lradsq = (crater%xl - xp + 0.5_DP * user%pix)**2 + (crater%yl - yp - 0.5_DP * user%pix)**2 - lrad = sqrt(lradsq) - end select - - mag = ( ( (abs(cos(nrays * theta / 4.0_DP)))**n2 + & - (abs(sin(nrays * theta / 4.0_DP)))**n2 )**(-1.0_DP/n1) ) - magf = rayf * ( ( (abs(cos(nfrays * theta / 4.0_DP)))**n2f + & - (abs(sin(nfrays * theta / 4.0_DP)))**n2f )**(-1.0_DP/n1f)) - lradp = crater%frad * mag - lradf = crater%frad * magf - lradp = lradp + lradf !max(lradp, lradf) - - if ((lrad < maxdistance) .and. (lrad < lradp)) then - rayavg = rayavg + 1.0_DP - end if - end do - rayavg = rayavg / 5 - isray(i,j) = rayavg - numinray(k) = numinray(k) + rayavg - end do - end do - if ((xi /= -inc).or.(xf /= inc).or.(yi /= -inc).or.(yf /= inc)) then - ejdistribution(xi:xf,yi:yf) = isray(xi:xf,yi:yf) - else - do k = 1,nef - if (numinray(k) == 0) then - mefarray(k) = 1.0_DP - else - mefarray(k) = totnum(k) / numinray(k) - end if - end do - ejdistribution = 1.0_DP +! ! Don't let the bin resolution get too small for small craters. +! binres = min(10.0_DP,crater%ejdis / 10) * user%pix +! nef = ceiling(crater%ejdis / binres) ! bin size for enhancement factor calculation +! allocate(numinray(nef),totnum(nef),mefarray(nef)) + +! totnum = 0 +! numinray = 0 +! do j = yi,yf +! do i = xi,xf +! iradsq = i*i + j*j + +! xpi = crater%xlpx + i +! ypi = crater%ylpx + j + +! ! Find distance from crater center to current pixel center in real space +! xp = xpi * user%pix +! yp = ypi * user%pix +! +! lradsq = (crater%xl - xp)**2 + (crater%yl - yp)**2 +! lrad = sqrt(lradsq) + +! if ((iradsq > incsq).or.(lrad <= crater%rad)) cycle +! +! ! Sum up the number of pixels in this bin for the mass enhancement calculation +! k = ceiling(lrad / binres) +! totnum(k) = totnum(k) + 1.0_DP + +! ! Ray pattern setup +! ! Take average of five points in the pixel to "feather" the edges of the ray +! rayavg = 0.0_DP +! do n = 1, 5 +! select case(n) +! case(1) +! theta = atan2(j * 1._DP,i * 1._DP) + 2.0_DP * PI +! case(2) +! theta = atan2(j * 1._DP - 0.5_DP,i * 1._DP - 0.5_DP) + 2.0_DP * PI +! lradsq = (crater%xl - xp - 0.5_DP * user%pix)**2 + (crater%yl - yp - 0.5_DP * user%pix)**2 +! lrad = sqrt(lradsq) +! case(3) +! theta = atan2(j * 1._DP - 0.5_DP,i * 1._DP + 0.5_DP) + 2.0_DP * PI +! lradsq = (crater%xl - xp - 0.5_DP * user%pix)**2 + (crater%yl - yp + 0.5_DP * user%pix)**2 +! lrad = sqrt(lradsq) +! case(4) +! theta = atan2(j * 1._DP + 0.5_DP,i * 1._DP + 0.5_DP) + 2.0_DP * PI +! lradsq = (crater%xl - xp + 0.5_DP * user%pix)**2 + (crater%yl - yp + 0.5_DP * user%pix)**2 +! lrad = sqrt(lradsq) +! case(5) +! theta = atan2(j * 1._DP + 0.5_DP,i * 1._DP - 0.5_DP) + 2.0_DP * PI +! lradsq = (crater%xl - xp + 0.5_DP * user%pix)**2 + (crater%yl - yp - 0.5_DP * user%pix)**2 +! lrad = sqrt(lradsq) +! end select +! +! mag = ( ( (abs(cos(nrays * theta / 4.0_DP)))**n2 + & +! (abs(sin(nrays * theta / 4.0_DP)))**n2 )**(-1.0_DP/n1) ) +! magf = rayf * ( ( (abs(cos(nfrays * theta / 4.0_DP)))**n2f + & +! (abs(sin(nfrays * theta / 4.0_DP)))**n2f )**(-1.0_DP/n1f)) +! lradp = crater%frad * mag +! lradf = crater%frad * magf +! lradp = lradp + lradf !max(lradp, lradf) +! +! if ((lrad < maxdistance) .and. (lrad < lradp)) then +! rayavg = rayavg + 1.0_DP +! end if +! end do +! rayavg = rayavg / 5 +! isray(i,j) = rayavg +! numinray(k) = numinray(k) + rayavg +! end do +! end do +! if ((xi /= -inc).or.(xf /= inc).or.(yi /= -inc).or.(yf /= inc)) then +! ejdistribution(xi:xf,yi:yf) = isray(xi:xf,yi:yf) +! else +! do k = 1,nef +! if (numinray(k) == 0) then +! mefarray(k) = 1.0_DP +! else +! mefarray(k) = totnum(k) / numinray(k) +! end if +! end do +! ejdistribution = 1.0_DP ! Now calculate the ejecta distribution with mass enhancement factor - do j = yi,yf - do i = xi,xf - ejdistribution(i,j) = 0.0_DP - iradsq = i*i + j*j - lrad = sqrt(lradsq) - if (lrad < crater%continuous) then - ejdistribution(i,j) = 1.0_DP - else - frac = lrad / binres - k = ceiling(frac) - ! Interpolate the mass enhancement factor between array points - mef = mefarray(k) + (frac - k * 1.0_DP) * (mefarray(k) - mefarray(k - 1)) - ejdistribution(i,j) = mef * isray(i,j) - end if - end do - end do - end if +! do j = yi,yf +! do i = xi,xf +! ejdistribution(i,j) = 0.0_DP +! iradsq = i*i + j*j +! lrad = sqrt(lradsq) +! if (lrad < crater%continuous) then +! ejdistribution(i,j) = 1.0_DP +! else +! frac = lrad / binres +! k = ceiling(frac) +! ! Interpolate the mass enhancement factor between array points +! if (k > 1) then +! mef = mefarray(k) + (frac - k * 1.0_DP) * (mefarray(k) - mefarray(k - 1)) +! else +! mef = mefarray(k) + (frac - k * 1.0_DP) * (mefarray(k + 1) - mefarray(k)) +! end if +! + ! ejdistribution(i,j) = mef * isray(i,j) + ! end if + ! end do + ! end do + !end if - deallocate(numinray,totnum,mefarray) + !deallocate(numinray,totnum,mefarray) end subroutine ejecta_ray_pattern diff --git a/src/ejecta/ejecta_table_define.f90 b/src/ejecta/ejecta_table_define.f90 index 37f48693..bbb732e3 100644 --- a/src/ejecta/ejecta_table_define.f90 +++ b/src/ejecta/ejecta_table_define.f90 @@ -42,7 +42,7 @@ subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt) ! Executable code ! Get estimate of size of ejb table - crater%continuous = 2.348_DP * crater%frad**(1.006_DP) ! Continuous ejecta distance From Moore (1974) eq. 1 + crater%continuous = RCONT * crater%frad**(EXPCONT) ! Continuous ejecta distance From Moore (1974) eq. 1 crater%ejdis = DISEJB * crater%continuous ! We go out a factor of 3 to get the discontinuous ejecta thickness domain%ejbres = (log(crater%ejdis) - log(crater%rad)) / EJBTABSIZE diff --git a/src/ejecta/module_ejecta.f90 b/src/ejecta/module_ejecta.f90 index 6cc14921..669ea8a2 100644 --- a/src/ejecta/module_ejecta.f90 +++ b/src/ejecta/module_ejecta.f90 @@ -39,13 +39,14 @@ end subroutine ejecta_emplace end interface interface - subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,ejdistribution) + subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,diffdistribution,ejdistribution) use module_globals implicit none type(usertype),intent(in) :: user type(surftype),dimension(:,:),intent(in) :: surf - type(cratertype),intent(in) :: crater + type(cratertype),intent(inout) :: crater integer(I4B),intent(in) :: inc,xi,xf,yi,yf + real(DP),dimension(xi:xf,yi:yf),intent(out) :: diffdistribution real(DP),dimension(xi:xf,yi:yf),intent(out) :: ejdistribution end subroutine ejecta_ray_pattern end interface diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index ffb8a952..a46c2fa2 100644 --- a/src/globals/module_globals.f90 +++ b/src/globals/module_globals.f90 @@ -48,15 +48,18 @@ module module_globals real(DP),parameter :: VSMALL = tiny(1._DP) ! Very small number real(DP),parameter :: LOGVSMALL = log(VSMALL) ! log of a very small number real(DP),parameter :: VBIG = huge(1._DP) ! Very big number -real(DP),parameter :: SMALLFAC = 1e-5_DP ! Smallest unit of measurement proportional to pixel size +real(DP),parameter :: SMALLFAC = 1e-5_DP ! Smallest unit of measurement proportional to pixel size integer(I4B),parameter :: MAXLAYER=20 ! Maximum number of layers (you need roughly 1-2 layers per order of magnitude of ! resolution -real(DP),parameter :: TALLYCOVERAGE = 0.01_DP ! The total area coverage to reach before a tally step is executed +real(DP),parameter :: TALLYCOVERAGE = 0.01_DP ! The total area coverage to reach before a tally step is executed real(DP),parameter :: SUBPIXELCOVERAGE = 0.025_DP ! The total area coverage to reach before a subpixel evaluate step is executed: 0.05_DP real(DP),parameter :: COOKIESIZE = 3.0_DP ! Relative size of old crater to new crater that cookie cutting is applied ! Only craters smaller than COOKIESIZE times the new crater are cookie cut real(DP),parameter :: ALPHA = 0.125_DP real(DP),parameter :: DISEJB = 100.0_DP ! The extent of discontinuous ejecta in the unit of crater radii. It is used in ejecta_table_define.f90 +real(DP),parameter :: RCONT = 2.25267_DP ! Coefficient of continuous ejecta size power law from Moore et al. (1974) - scaled from km to m +real(DP),parameter :: EXPCONT = 1.006_DP ! Exponentt of continuous ejecta size power law from Moore et al. (1974) + type regodatatype real(DP) :: thickness @@ -100,6 +103,7 @@ module module_globals real(DP) :: cxexp,cxtran ! simple to complex scaling parameters real(DP) :: kdiffterm, saccelterm ! seismic diffusion and accelleration terms real(DP) :: continuous ! Size of the continuous ejecta blanket + real(DP) :: fe ! Equivalent degradation radius ! Pixel domain properties integer(I4B) :: xlpx,ylpx ! Crater center in pixels integer(I4B) :: fcratpx,fradpx,rimdispx,ejdispx @@ -178,13 +182,15 @@ module module_globals real(DP) :: regcoh ! target surface regolith layer cohesion ! Crater diffusion input parameters - real(DP) :: soften_factor ! Kbar,0 from Minton et al. (2017) - real(DP) :: soften_slope ! power law index of diffusion model Kbar - real(DP) :: soften_size ! size of diffusion area in crater radii + real(DP) :: Kd1 ! Degradation function coefficient (from Minton et al. (2018)) + real(DP) :: psi ! Degradation function exponent (from Minton et al. (2018)) + real(DP) :: fe ! Scale factor for size of degradation region (from Minton et al. (2018)) ! Ejecta softening variables logical :: dosoftening ! Set T to use the extra crater softening model real(DP) :: ejecta_truncation ! Set the number of crater diameters to truncate the ejecta + logical :: dorays ! Set T to use ray model + logical :: superdomain ! Set T to include the superdomain ! Regolith tracking variables logical :: doregotrack ! Set T to use the regolith tracking model (EXPERIMENTAL) @@ -235,7 +241,7 @@ module module_globals ! Progress bar variables integer(I4B),parameter :: PBARRES = 100 -integer(I4B),parameter :: PBARSIZE = 40 +integer(I4B),parameter :: PBARSIZE = 50 integer(I4B),parameter :: MESSAGESIZE = 32 integer(I4B) :: pbarival integer(I4B) :: pbarpos @@ -283,11 +289,11 @@ module module_globals real(DP),parameter :: TRSIM = 1.25_DP ! ? real(DP),parameter :: EXFAC = 0.1_DP ! Excavation depth relative to transient crater diameter real(DP),parameter :: CXEXPS = 1._DP / 0.885_DP - 1.0_DP ! Complex crater scaling exponent (see Croft 1985) -real(DP),parameter :: SIMCOMKS = 16533.8_DP ! ? -real(DP),parameter :: SIMCOMPS = -1.0303_DP ! ? +real(DP),parameter :: SIMCOMKS = 16533.8_DP ! Simple-to-complex transition scaling coefficient for silicate rock +real(DP),parameter :: SIMCOMPS = -1.0303_DP ! Simple-to-complex transition scaling exponent for silicate rock real(DP),parameter :: CXEXPI = 0.155_DP ! ? -real(DP),parameter :: SIMCOMKI = 3081.39_DP ! ? -real(DP),parameter :: SIMCOMPI = -1.22486_DP ! ? +real(DP),parameter :: SIMCOMKI = 3081.39_DP ! Simple-to-complex transition scaling coefficient for ice +real(DP),parameter :: SIMCOMPI = -1.22486_DP ! Simple-to-complex transition scaling exponent for ice real(DP),parameter :: SUBPIXFAC = 0.1_DP ! Subpixel resolution (used for lookup tables and rim creation) integer(I4B),parameter :: EJBTABSIZE = 1000 ! Lookup table size real(DP),parameter :: CRITSLP = 0.7_DP ! critical slope angle @@ -295,10 +301,9 @@ module module_globals real(DP),parameter :: BOWLFRAC = 0.2_DP ! Fraction of crater interior pixels to use for the bowl-to-rim height calculation ! (calibrated for Orientale using Potter et al. 2012) real(DP),parameter :: TRNRATIO = 0.30_DP ! The ratio of the transient crater depth to the crater diameter. - -!real(DP),parameter :: SOFTEN_FACTOR = 0.60_DP ! Extra per crater diffusion constant -!#real(DP) :: SOFTEN_FACTOR = 0.60_DP ! Extra per crater diffusion constant -!real(DP),parameter :: SOFTEN_SLOPE = 1.3_DP ! Extra per crater diffusion power law slope +real(DP),parameter :: KD1PROX = 0.27_DP ! Intrinsic proximal ejecta degradation function coefficient +real(DP),parameter :: PSIPROX = 2.0_DP ! Intrinsic proximal ejecta degradation function exponent +real(DP),parameter :: FEPROX = 1.0_DP ! Intrinsic proximal ejecta degradation function size scale factor ! Seismic shaking parameters real(DP),parameter :: SEISFREQ = 20.0_DP ! seismic wave frequency @@ -310,4 +315,12 @@ module module_globals real(DP),parameter :: GFAC = 0.487_DP ! gravitational acceleration exponent real(DP),parameter :: DFAC = 0.556_DP ! impact distance exponent +! Crater ray parameters +real(DP),parameter :: rray = 16_DP +integer(I4B),parameter :: Nraymax = 16 +real(DP),parameter :: fpeak = 8000_DP ! narrow ray: rw0 propto 1/4 +real(DP),parameter :: rayp = 2.0_DP +integer(I4B),parameter :: rayq = 4 +real(DP),parameter :: rayfmult = (5)**(-4.0_DP / (1.2_DP)) + end module module_globals diff --git a/src/io/io_input.f90 b/src/io/io_input.f90 index 3a53b6d4..258a83ae 100644 --- a/src/io/io_input.f90 +++ b/src/io/io_input.f90 @@ -77,6 +77,8 @@ subroutine io_input(infile,user) user%doscour = .false. user%tallyonly = .false. user%dosoftening = .true. + user%dorays = .false. + user%superdomain = .true. user%doregotrack = .false. user%basinimp = huge(0._DP) user%maxcrat = 1.00_DP @@ -84,9 +86,9 @@ subroutine io_input(infile,user) user%testtally = .false. user%killatmaxcrater = .false. user%doporosity = .false. - user%soften_factor = 2.0_DP - user%soften_slope = 4.0_DP - user%soften_size = 30.0_DP + user%Kd1 = 0.00_DP + user%psi = 2.0_DP + user%fe = 1.0_DP user%ejecta_truncation = 10.0_DP open(unit=LUN,file=infile,status="old",iostat=ierr) @@ -274,6 +276,16 @@ subroutine io_input(infile,user) call io_get_token(line, ilength, ifirst, ilast, ierr) token = line(ifirst:ilast) read(token, *) user%dosoftening + case ("DORAYS") + ifirst = ilast + 1 + call io_get_token(line, ilength, ifirst, ilast, ierr) + token = line(ifirst:ilast) + read(token, *) user%dorays + case ("SUPERDOMAIN") + ifirst = ilast + 1 + call io_get_token(line, ilength, ifirst, ilast, ierr) + token = line(ifirst:ilast) + read(token, *) user%superdomain case ("DOREGOTRACK") ifirst = ilast + 1 call io_get_token(line, ilength, ifirst, ilast, ierr) @@ -425,21 +437,21 @@ subroutine io_input(infile,user) token = line(ifirst:ilast) read(token, *) user%shadedmaxh read(token, *) user%shadedminh - case ("SOFTEN_FACTOR") + case ("KD1") ifirst = ilast + 1 call io_get_token(line, ilength, ifirst, ilast, ierr) token = line(ifirst:ilast) - read(token, *) user%soften_factor - case ("SOFTEN_SLOPE") + read(token, *) user%Kd1 + case ("PSI") ifirst = ilast + 1 call io_get_token(line, ilength, ifirst, ilast, ierr) token = line(ifirst:ilast) - read(token, *) user%soften_slope - case ("SOFTEN_SIZE") + read(token, *) user%psi + case ("FE") ifirst = ilast + 1 call io_get_token(line, ilength, ifirst, ilast, ierr) token = line(ifirst:ilast) - read(token, *) user%soften_size + read(token, *) user%fe !************************************************************************** ! The following is for backwards compatibility with older style input files !************************************************************************** diff --git a/src/io/io_updatePbar.f90 b/src/io/io_updatePbar.f90 index f7e9f57a..7ba679d9 100644 --- a/src/io/io_updatePbar.f90 +++ b/src/io/io_updatePbar.f90 @@ -39,20 +39,20 @@ subroutine io_updatePbar(message) perc=floor(100*real(pbarpos)/real(PBARRES)) write(unit=bar(1:3),fmt="(i3)") perc endpoint = min(ceiling(real(pbarpos*PBARSIZE)/real(PBARRES)),PBARSIZE) -do k = 1,endpoint - 1 +do k = 1,endpoint !- 1 bar(6+k:6+k)=pbarchar(k:k) end do -select case(flip) -case(1) - bar(6+endpoint:6+endpoint)="/" -case(2) - bar(6+endpoint:6+endpoint)="-" -case(3) - bar(6+endpoint:6+endpoint)="\" -case(4) - bar(6+endpoint:6+endpoint)="|" -end select -flip = flip + 1 +!select case(flip) +!case(1) +! bar(6+endpoint:6+endpoint)="/" +!case(2) +! bar(6+endpoint:6+endpoint)="-" +!case(3) +! bar(6+endpoint:6+endpoint)="\" +!case(4) +! bar(6+endpoint:6+endpoint)="|" +!end select +!flip = flip + 1 if (flip > 4) flip = 1 ! print the progress bar. write(fmtlabel,'("(A1,A",I2.2,",1X,A",I2.2,",$)")') PBARSIZE+7,MESSAGESIZE diff --git a/src/io/module_ejecta.f90 b/src/io/module_ejecta.f90 new file mode 100644 index 00000000..669ea8a2 --- /dev/null +++ b/src/io/module_ejecta.f90 @@ -0,0 +1,156 @@ +!****h* ejecta/module_ejecta +! Name +! module_ejecta -- Module for ejecta +! Notes +! It includes crater ray model based on Superformula. +!*** + +!********************************************************************************************************************************** +! +! Unit Name : module_crater +! Unit Type : module +! Project : CTEM +! Language : Fortran 2003 +! +! Description : Parameters and subroutine interface blocks for crater creation, including scaling laws and Monte Carlo +! routines +! +! Notes : +! +!********************************************************************************************************************************** +module module_ejecta +use module_globals +implicit none +public +save + + interface + subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) + use module_globals + implicit none + type(usertype),intent(in) :: user + type(surftype),dimension(:,:),intent(inout) :: surf + type(cratertype),intent(inout) :: crater + type(domaintype),intent(in) :: domain + integer(I4B),intent(in) :: ejtble + type(ejbtype),dimension(ejtble),intent(in) :: ejb + real(DP),intent(in) :: deltaMtot + end subroutine ejecta_emplace + end interface + + interface + subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,diffdistribution,ejdistribution) + use module_globals + implicit none + type(usertype),intent(in) :: user + type(surftype),dimension(:,:),intent(in) :: surf + type(cratertype),intent(inout) :: crater + integer(I4B),intent(in) :: inc,xi,xf,yi,yf + real(DP),dimension(xi:xf,yi:yf),intent(out) :: diffdistribution + real(DP),dimension(xi:xf,yi:yf),intent(out) :: ejdistribution + end subroutine ejecta_ray_pattern + end interface + + interface + subroutine ejecta_rootfind(user,crater,domain,erad,lrad,vejsq,ejang,firstrun) + use module_globals + implicit none + type(usertype),intent(in) :: user + type(cratertype),intent(in) :: crater + type(domaintype),intent(in) :: domain + real(DP),intent(inout) :: erad + real(DP),intent(in) :: lrad + real(DP),intent(out) :: vejsq,ejang + logical,intent(inout) :: firstrun + end subroutine ejecta_rootfind + end interface + + interface + subroutine ejecta_blanket(user,crater,domain,erad,lrad,vejsq,ejang,firstrun) + use module_globals + implicit none + type(usertype),intent(in) :: user + type(cratertype),intent(in) :: crater + type(domaintype),intent(in) :: domain + real(DP),intent(in) :: erad + real(DP),intent(out) :: lrad,vejsq,ejang + logical,intent(inout) :: firstrun + end subroutine ejecta_blanket + end interface + + interface + subroutine ejecta_thickness(user,crater,erad1,erad2,lrad1,lrad2,thick) + use module_globals + implicit none + type(usertype),intent(in) :: user + type(cratertype),intent(in) :: crater + real(DP),intent(in) :: erad1,erad2,lrad1,lrad2 + real(DP),intent(out) :: thick + end subroutine ejecta_thickness + end interface + + interface + function ejecta_blanket_func(user,crater,domain,erad,lrad,vejsq,ejang,firstrun) result(ans) + use module_globals + implicit none + type(usertype),intent(in) :: user + type(cratertype),intent(in) :: crater + type(domaintype),intent(in) :: domain + real(DP),intent(in) :: erad,lrad + logical,intent(inout) :: firstrun + real(DP),intent(out) :: vejsq,ejang + real(DP) :: ans + end function ejecta_blanket_func + end interface + + interface + subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt) + use module_globals + implicit none + type(usertype),intent(in) :: user + type(cratertype),intent(inout) :: crater + type(domaintype),intent(inout) :: domain + type(ejbtype),dimension(EJBTABSIZE),intent(out) :: ejb + integer(I4B),intent(out) :: ejtble + real(DP),intent(out),optional :: melt + end subroutine ejecta_table_define + end interface + + interface + subroutine ejecta_interpolate(crater,domain,lrad,ejb,ejtble,ebh,vsq,theta,erad,melt) + use module_globals + implicit none + type(cratertype),intent(in) :: crater + type(domaintype),intent(in) :: domain + real(DP),intent(in) :: lrad + integer(I4B),intent(in) :: ejtble + type(ejbtype),dimension(ejtble),intent(in) :: ejb + real(DP),intent(out) :: ebh + real(DP),intent(out),optional :: vsq,theta + real(DP),intent(out),optional :: erad,melt + end subroutine ejecta_interpolate + end interface + + interface + subroutine ejecta_soften(user,surf,N,indarray,cumulative_elchange) + use module_globals + implicit none + type(usertype),intent(in) :: user + type(surftype),dimension(:,:),intent(inout) :: surf + integer(I4B),intent(in) :: N + integer(I4B),dimension(2,N,N),intent(in) :: indarray + real(DP),dimension(N,N),intent(inout) :: cumulative_elchange + end subroutine ejecta_soften + end interface + + interface + subroutine ejecta_distance_estimate(user,crater,domain,ejdis_estimate) + use module_globals + implicit none + type(usertype),intent(in) :: user + type(cratertype),intent(inout) :: crater + type(domaintype),intent(in) :: domain + real(DP),intent(out) :: ejdis_estimate + end subroutine ejecta_distance_estimate + end interface +end module