diff --git a/src/crater/crater_subpixel_diffusion.f90 b/src/crater/crater_subpixel_diffusion.f90 index 2a97e6f3..86b53973 100644 --- a/src/crater/crater_subpixel_diffusion.f90 +++ b/src/crater/crater_subpixel_diffusion.f90 @@ -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 @@ -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 @@ -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) @@ -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 diff --git a/src/crater/crater_tally_calibrated_count.f90 b/src/crater/crater_tally_calibrated_count.f90 new file mode 100644 index 00000000..8ba3e744 --- /dev/null +++ b/src/crater/crater_tally_calibrated_count.f90 @@ -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 + diff --git a/src/ejecta/ejecta_emplace.f90 b/src/ejecta/ejecta_emplace.f90 index 4f9380b9..837d197e 100644 --- a/src/ejecta/ejecta_emplace.f90 +++ b/src/ejecta/ejecta_emplace.f90 @@ -96,34 +96,14 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) ! Internal variables 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 - real(DP) :: xp,yp,fradsq,fradpxsq,radsq,ebh,ejdissq,continuous,ejbmass,fmasscons - real(DP),dimension(:,:),allocatable :: cumulative_elchange,big_cumulative_elchange + 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),dimension(:,:),allocatable :: cumulative_elchange,big_cumulative_elchange,kdiff,big_kdiff,cel,big_cel integer(I4B),dimension(:,:,:),allocatable :: indarray,big_indarray - integer(I4B) :: bigi,bigj + real(DP),dimension(:,:),allocatable :: ejdistribution + integer(I4B) :: bigi,bigj,maxhits,nin,nnot character(len=MESSAGESIZE) :: message ! message for the progress bar - !Ray model parameters and variables by Gielis superformula - integer(I4B) :: nrays - real(DP),dimension(15) :: rn - real(DP) :: theta, lradp - real(DP), parameter :: n1 = 4.0_DP - real(DP) :: n2, mag - ! Ray Mass conservation - real(DP), parameter :: a = 16.8799 !a = 11.8126 ! Fitting parameters for a relation between ray length and radius of crater - 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 extent - - ! Enhanced factor test - real(DP) :: xef, yef, thetamax - integer(I4B) :: xefpi, yefpi, nef, ray_pix - integer(I4B), dimension(:), allocatable :: sf, tot - real(DP), dimension(:), allocatable :: ef - - ! Flowery ray variables - integer(I4B) :: nfrays - real(DP) :: n1f, n2f, magf, lradf, rayf - real(DP) :: vsq, ejtheta + ! Ray mixing model variables real(DP) :: dsc @@ -132,6 +112,7 @@ 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) :: vsq, ejtheta integer(I4B) :: klo,ind ! Executable code @@ -144,16 +125,16 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) if (crater%ejdis <= crater%rad) return ! determine area to effect - continuous = 2.348_DP * crater%frad**(1.006_DP) - inc = max(min(nint(crater%frad * user%ejecta_truncation / user%pix) + 1,PBCLIM*user%gridsize),1) + 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) crater%maxinc = max(crater%maxinc,inc) radsq = crater%rad**2 fradsq = crater%frad**2 fradpxsq = crater%fradpx**2 - incsq = inc**2 ejdissq = crater%ejdis**2 if (inc >= user%gridsize / 2) then @@ -165,101 +146,25 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) end if 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(indarray(2,-inc:inc,-inc:inc)) + allocate(ejdistribution(-inc:inc,-inc:inc)) - call random_number(rn) - cumulative_elchange = 0._DP + 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 ! within a circle of radius irad ! *************************** Continuous Ejecta Formula *****************************! - !^^^^^^^^^^^^^^^^^^^^^^^^ - - ! *************************** Superformula Ray Model ************************************! - ! *************************** Part I. Spoke and Skinny Ray Model ************************************! - ! From fitting Jake Elliot's ray mapping data, it is a linear function that describes the relationship between the median value of ray length - ! distribution and the radius of rayed craters (in unit of kilometers) - ! Also, we need to scale it with the continuous ejecta's extent for ray model -!# mvrld = a * (crater%frad/1000.0)**(b) -!# mvrldsc = mvrld / (continuous/crater%frad) - ! It appears that no strong correlation between the number of rays and size of craters. - ! 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 - ! 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(15)) + 6 - mvrld = crater%ejdis - mvrldsc = mvrld / continuous - n2 = 8.0_DP * ( log10(mvrldsc) / log10(2.0_DP) ) + 2.0_DP - ! *************************** Part II. Flowery Ray Model ************************************! - nfrays = 20 - n1f = 1.0_DP - n2f = 0.5_DP - rayf = 5.0 !7.5_DP - - !write(*,*) mvrld, mvrldsc, n2 - - ! Enhanced factor test: Build a enhanced factor quick lookup table - ! Use a circular sector with an angle containing a half ray, which is pi/nrays. First, we determine the looping area for - ! this sector, then the look-up table will be built while looping over each pixel and then making it to go to the right bin. - ! The size of bin is one pixel. - xef = mvrld - thetamax = PI/dble(nrays) - yef = mvrld * sin(thetamax) - xefpi = nint(xef / user%pix) - yefpi = nint(yef / user%pix) - nef = ceiling( (mvrld) / user%pix) !ceiling( (mvrld - crater%rad) / user%pix) - allocate(sf(nef)) - allocate(tot(nef)) - allocate(ef(nef)) - tot = 0 - sf = 0 - - do j = yefpi, 0, -1 - do i = xefpi, 0, -1 - xp = (i + crater%xlpx) * user%pix - yp = (j + crater%ylpx) * user%pix - lradsq = (xp - crater%xl)**2 + (yp - crater%yl)**2 - lrad = sqrt(lradsq) - ! inside or outside a ray? - theta = atan2(j * 1._DP,i * 1._DP) - 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 = max(lradp, lradf) - !if ( lrad > continuous .and. lrad < mvrld .and. theta < thetamax .and. theta>0._DP) then - if (lrad >= crater%rad .and. lrad < mvrld .and. theta < thetamax .and. theta>0._DP) then - ray_pix = ceiling( (lrad - crater%rad) / user%pix) - tot(ray_pix) = tot(ray_pix) + 1 - if (lrad < lradp) sf(ray_pix) = sf(ray_pix) + 1 - end if - end do - end do + call ejecta_ray_pattern(user,surf,crater,inc,-inc,inc,-inc,inc,ejdistribution) - ! Smooth out the enhanced factor lookup table - do i=1,nef - ef(i) = max( dble(tot(i)) / dble(sf(i)), 1.0_DP) - end do - - do i=1,nef - if (ef(i) /= ef(i) .or. ef(i) > VBIG) then - ef(i) = dble(tot(i)) / 1.0 - end if - end do - - deallocate(sf) - deallocate(tot) - 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,continuous) +! !$OMP SHARED(inc,incsq,ejdissq,fradsq,radsq,indarray,cumulative_elchange,rn) ! !$OMP REDUCTION(+:massray,massej,massrayef) do j = -inc,inc do i = -inc,inc @@ -280,17 +185,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) indarray(1,i,j) = xpi indarray(2,i,j) = ypi - if ((iradsq > incsq).or.(lrad <= crater%rad)) cycle - - ! Ray pattern setup - theta = atan2(j * 1._DP,i * 1._DP) + 2.0_DP * PI - 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 = max(lradp, lradf) + ! Estimate ejecta pattern distortion due to target surface angle and topography ! This must be done iteratively because the ejection distance and ejection angle vary @@ -325,37 +220,42 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) end do if (vsq < 0.0_DP) cycle - if ((distance < maxdistance) .and. (distance < lradp)) then ! Inside a ray or continuous ejecta - - - ! Ray mass conservation - ray_pix = max(min(ceiling((lrad - crater%rad) / user%pix),1),nef) - ebh = ebh * ef(ray_pix) - - 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 - - else ! Outside a ray - ebh = 0._DP + if (distance /= distance) cycle + idistorted = int(i * distance / lrad) + if (abs(idistorted) > inc) cycle + jdistorted = int(j * distance / lrad) + if (abs(jdistorted) > inc) cycle + + iradsq = idistorted**2 + jdistorted**2 + if ((iradsq > incsq).or.(distance <= crater%rad)) cycle + + ebh = ejdistribution(idistorted,jdistorted) * ebh + + 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 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 + end if end do end do ! !$OMP END PARALLEL DO - - deallocate(ef) + kdiff = kdiff * (nnot + nin) / (1.0_DP * nin) ! Do mass conservation by adjusting ejecta thickness fmasscons = (-deltaMtot)/ ejbmass - !write(*,*) 'fmasscons = ',fmasscons 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 + 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 @@ -365,11 +265,22 @@ 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 ejecta_soften(user,surf,2 * inc + 1,indarray,cumulative_elchange) + 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 @@ -386,21 +297,29 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) xpi = indarray(1,i,j) ypi = indarray(2,i,j) big_cumulative_elchange(xpi,ypi) = big_cumulative_elchange(xpi,ypi) + cumulative_elchange(i,j) + big_kdiff(xpi,ypi) = big_kdiff(xpi,ypi) + kdiff(i,j) end do end do + 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 ejecta_soften(user,surf,user%gridsize + 2,big_indarray,big_cumulative_elchange) - deallocate(big_cumulative_elchange,big_indarray) + 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) + deallocate(cumulative_elchange,indarray,ejdistribution,kdiff,cel) return end subroutine ejecta_emplace diff --git a/src/ejecta/ejecta_ray_pattern.f90 b/src/ejecta/ejecta_ray_pattern.f90 new file mode 100644 index 00000000..b376735e --- /dev/null +++ b/src/ejecta/ejecta_ray_pattern.f90 @@ -0,0 +1,208 @@ +!****f* ejecta/ejecta_ray_pattern +! Name +! ejecta_ray_pattern -- Calculate ejecta ray pattern +! SYNOPSIS +! This uses +! * module_globals +! * module_util +! * module_io +! * module_crater +! * module_regolith +! * module_ejecta +! +! call ejecta_ray_pattern(user,surf,crater,inc,ejdistribution) +! +! DESCRIPTION +! +! Models the discontinuous ray pattern based on the Superformula. +! Citation: Gielis, J. "A Generic Geometric Transformation that Unifies a Wide Range of Natural +! and Abstract Shapes." Amer. J. Botany 90, 333-338, 2003. +! +! ARGUMENTS +! Input +! * user -- User input parameters +! * surf -- Surface ggrid +! * crater -- Crater dimension container +! * domain -- Simulation domain variable container +! +! Output +! * ejdist -- logical array containing the pixels that contain ejecta +! +!*** + +!********************************************************************************************************************************** +subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,ejdistribution) + use module_globals + use module_util + use module_io + use module_crater + use module_regolith + use module_ejecta, EXCEPT_THIS_ONE => ejecta_ray_pattern + implicit none + + ! Arguments + type(usertype),intent(in) :: user + type(surftype),dimension(:,:),intent(in) :: surf + type(cratertype),intent(in) :: crater + integer(I4B),intent(in) :: inc,xi,xf,yi,yf + 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 + real(DP) :: theta, lradp, maxdistance + real(DP), parameter :: n1 = 4.0_DP + real(DP) :: n2, mag + real(DP),dimension(xi:xf,yi:yf) :: isray + real(DP),dimension(:),allocatable :: numinray,totnum + real(DP),dimension(:),allocatable :: mefarray + + ! Flowery ray variables + integer(I4B) :: nfrays + real(DP) :: n1f, n2f, magf, lradf, rayf,rayavg + real(DP), parameter :: a = 16.8799 !a = 11.8126 ! Fitting parameters for a relation between ray length and radius of crater + 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 + + call random_number(rn) + + !^^^^^^^^^^^^^^^^^^^^^^^^ + + ! *************************** Superformula Ray Model ************************************! + ! *************************** Part I. Spoke and Skinny Ray Model ************************************! + ! From fitting Jake Elliot's ray mapping data, it is a linear function that describes the relationship between the median value of ray length + ! distribution and the radius of rayed craters (in unit of kilometers) + ! Also, we need to scale it with the continuous ejecta's extent for ray model +!# mvrld = a * (crater%frad/1000.0)**(b) +!# mvrldsc = mvrld / (continuous/crater%frad) + ! It appears that no strong correlation between the number of rays and size of craters. + ! 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 + ! 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 + + ! 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 + + ! 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 + + 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 e64878de..37f48693 100644 --- a/src/ejecta/ejecta_table_define.f90 +++ b/src/ejecta/ejecta_table_define.f90 @@ -42,7 +42,8 @@ subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt) ! Executable code ! Get estimate of size of ejb table - crater%ejdis = DISEJB * 2.348_DP * crater%frad**(1.006_DP) ! Continuous ejecta distance From Moore (1974) eq. 1 + crater%continuous = 2.348_DP * crater%frad**(1.006_DP) ! 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 lrad = crater%rad !exp(log(crater%rad) !+ domain%ejbres) @@ -98,7 +99,6 @@ subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt) end do !write(*,*) 'A MELT ZONE of ',crater%frad,' meter-sized crater: ',rmelt,'at a rim',ejb(1)%meltfrac ! Get pixel space distance - crater%ejdis = crater%ejdis / 3.0_DP crater%ejdispx = nint(crater%ejdis / user%pix) return