diff --git a/src/ejecta/ejecta_ray_pattern.f90 b/src/ejecta/ejecta_ray_pattern.f90 old mode 100755 new mode 100644 index 03dfaa39..32982add --- a/src/ejecta/ejecta_ray_pattern.f90 +++ b/src/ejecta/ejecta_ray_pattern.f90 @@ -1,208 +1,375 @@ -!****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 - +!****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,diffdistribution,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(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,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 + real(DP),dimension(xi:xf,yi:yf) :: isray + 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 + 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 + 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 + + 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 / 2 + 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 ************************************! + ! 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 +! 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) + +end subroutine ejecta_ray_pattern + diff --git a/src/ejecta/ejecta_rootfind.f90 b/src/ejecta/ejecta_rootfind.f90 index 9e7f952b..716da6cf 100644 --- a/src/ejecta/ejecta_rootfind.f90 +++ b/src/ejecta/ejecta_rootfind.f90 @@ -153,7 +153,7 @@ subroutine ejecta_rootfind(user,crater,domain,erad,lrad,vejsq,ejang,firstrun) end if niter = i erad = resultat - if ((erad <= crater%frad).and.(lrad >= crater%frad).and.(erad >= 0._DP)) exit + if ((erad <= crater%rad).and.(lrad >= crater%rad).and.(erad >= 0._DP)) exit factor = 0.5_DP * (factor + 1.0_DP) ! Failed. Try again with a new factor end do everything diff --git a/src/ejecta/ejecta_table_define.f90 b/src/ejecta/ejecta_table_define.f90 index 06f54905..bbb732e3 100644 --- a/src/ejecta/ejecta_table_define.f90 +++ b/src/ejecta/ejecta_table_define.f90 @@ -42,15 +42,12 @@ subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt) ! Executable code ! Get estimate of size of ejb table - if (.not.user%discontinuous) then - crater%ejdis = 2 * 2.3_DP * crater%frad**(1.006_DP) ! Continuous ejecta distance From Melosh (1989) eq. 6.3.1 - else - crater%ejdis = DISEJB * 2.3_DP * crater%frad**(1.006_DP) ! Continuous ejecta distance From Melosh (1989) eq. 6.3.1 - end if + 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 = (crater%ejdis - crater%rad) / EJBTABSIZE - lrad = crater%frad - erad = crater%rad + domain%ejbres = (log(crater%ejdis) - log(crater%rad)) / EJBTABSIZE + lrad = crater%rad !exp(log(crater%rad) !+ domain%ejbres) + erad = crater%rad ejtble = EJBTABSIZE firstrun = .true. thick = 0._DP @@ -74,21 +71,18 @@ subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt) do k = 0,EJBTABSIZE call ejecta_rootfind(user,crater,domain,erad,lrad,vejsq,ejang,firstrun) if (k >= 1) then - !if (crater%frad >= 240.0) then - ! call ejecta_thickness(user,crater,eradold,erad,lrad - domain%ejbres,lrad,thick) - !else + !call ejecta_thickness(user,crater,eradold,erad,lrad - domain%ejbres,lrad,thick) + ejb(k)%lrad = log(lrad) + ! Use McGetchin et al. 1973 for ejecta thickness if (lrad >= crater%frad) then - thick = 0.14_DP * crater%frad**(0.74_DP) * (lrad /crater%frad)**(-3.0_DP) + thick = 0.14_DP * crater%frad**(0.74_DP) * (lrad / crater%frad)**(-3.0_DP) else - thick = 0.14_DP * crater%frad**(0.74_DP) / (crater%frad - crater%rad) * (lrad - crater%rad) + thick = 0.14_DP * crater%frad**(0.74_DP) / (crater%frad - crater%rad) * (lrad - crater%rad) end if - !end if - !call ejecta_thickness(user,crater,eradold,erad,lrad - domain%ejbres,lrad,thick) - ejb(k)%lrad = log(lrad - 0.5_DP * domain%ejbres) ejb(k)%thick = log(thick) ejb(k)%vesq = vejsq ejb(k)%angle = ejang - ejb(k)%erad = erad + ejb(k)%erad = log(erad) if (present(melt)) then call regolith_melt_fraction(dimp,depthb,erad,eradold,rmelt,melt) ejb(k)%meltfrac = melt @@ -100,7 +94,7 @@ subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt) exit end if end if - lrad = lrad + domain%ejbres + lrad = exp(log(lrad) + domain%ejbres) eradold = erad end do !write(*,*) 'A MELT ZONE of ',crater%frad,' meter-sized crater: ',rmelt,'at a rim',ejb(1)%meltfrac diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index a8b1b60c..269eeb56 100644 --- a/src/globals/module_globals.f90 +++ b/src/globals/module_globals.f90 @@ -14,7 +14,7 @@ module module_globals implicit none public -character(len=*),parameter :: CTEMVER = "1.3 DEVELOPMENT" +character(len=*),parameter :: CTEMVER = "1.4 DEVELOPMENT" ! Symbolic names for kind types of 4-, 2-, and 1-byte integers: integer, parameter :: I8B = selected_int_kind(17) @@ -48,21 +48,20 @@ 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 -integer(I2B),parameter :: MAXAGEBINS=60 ! Maximum number of bins in age distribution reset by impact melting -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 :: PF = 5.0e9 ! The shock pressure exceeding the hungoit elastic limit of common geological materials in Pa. (5 GPa) -real(DP),parameter :: RAD_GP = 20.0_DP ! The maximum radial position of producing impact glass spherules within a transient crater (unit of crater radii, crater%rad) +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(SP),dimension(MAXAGEBINS) :: age ! Record impact glass spherule ages in an array (or histogram defined by MAXAGEBINS) real(DP) :: thickness real(DP) :: meltfrac real(DP) :: comp @@ -78,8 +77,9 @@ module module_globals ! Derived data type for simulated surface type surftype - real(DP),dimension(MAXLAYER) :: diam + real(DP),dimension(MAXLAYER) :: diam ! Crater diameter real(SP),dimension(MAXLAYER) :: xl,yl ! Crater center + real(SP),dimension(MAXLAYER) :: timestamp ! Crater formation time real(DP) :: ejcov ! Ejecta coverage real(DP) :: dem ! Digital elevation model real(DP) :: mantle ! Height of mantle (should be smaller than dem) @@ -92,6 +92,7 @@ module module_globals real(DP) :: imp,imprad,impvel,sinimpang,impmass ! Impactor properties ! Real domain properties real(SP) :: xl,yl ! Crater center in simulated surface size units + real(SP) :: timestamp ! The formation time of the crater real(DP) :: rad ! Transient radius real(DP) :: grad ! Strengthless material transient crater radius real(DP) :: frad ! Final crater radius @@ -104,6 +105,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 @@ -185,13 +187,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. (2019)) + real(DP) :: psi ! Degradation function exponent (from Minton et al. (2019)) + real(DP) :: fe ! Scale factor for size of degradation region (from Minton et al. (2019)) ! 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 logical :: discontinuous ! Regolith tracking variables @@ -245,7 +249,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 @@ -253,6 +257,7 @@ module module_globals ! Grid array file names character(*),parameter :: DIAMFILE = 'surface_diam.dat' +character(*),parameter :: TIMEFILE = 'surface_time.dat' character(*),parameter :: EJCOVFILE = 'surface_ejc.dat' character(*),parameter :: DEMFILE = 'surface_dem.dat' character(*),parameter :: REGOFILE = 'surface_regotop.dat' @@ -280,7 +285,7 @@ module module_globals integer(I4B),parameter :: PBCLIM = 1 ! periodic boundary condition limit integer(I4B),parameter :: SMALLESTCOUNTABLE = 5 ! Minimum pixel diameter for a crater to be considered countable real(DP),parameter :: SMALLESTEJECTA = 1.5 ! Minimum number of pixels from center of crater for an ejecta to have any surface effects -integer(I4B),parameter :: TRUECOLS = 6 ! Number of columns in the true crater count array +integer(I4B),parameter :: TRUECOLS = 7 ! Number of columns in the true crater count array integer(I4B) :: NTHREADS = 1 ! Number of OpenMP threads (reset by OpenMP if a parallel environment is detected) integer(I4B),parameter :: INCPAR = 1 ! Minimum size of inc variables before parallelization kicks in @@ -295,11 +300,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 @@ -307,10 +312,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