diff --git a/examples/morphology_test_cases/.DS_Store b/examples/morphology_test_cases/.DS_Store deleted file mode 100644 index 3a3743e9..00000000 Binary files a/examples/morphology_test_cases/.DS_Store and /dev/null differ diff --git a/python/.idea/.gitignore b/python/.idea/.gitignore new file mode 100644 index 00000000..26d33521 --- /dev/null +++ b/python/.idea/.gitignore @@ -0,0 +1,3 @@ +# Default ignored files +/shelf/ +/workspace.xml diff --git a/python/ctem/ctem/NPF.py b/python/ctem/ctem/NPF.py new file mode 100644 index 00000000..cb7ee363 --- /dev/null +++ b/python/ctem/ctem/NPF.py @@ -0,0 +1,104 @@ +import numpy as np +# The Neukum production function +# Ivanov, Neukum, and Hartmann (2001) SSR v. 96 pp. 55-86 +def N1lunar(T): + return 5.44e-14 * (np.exp(6.93*T) - 1) + 8.38e-4 * T + + +def Nlunar(D): + # Lunar crater SFD + aL00 = -3.0876 + aL01 = -3.557528 + aL02 = 0.781027 + aL03 = 1.021521 + aL04 = -0.156012 + aL05 = -0.444058 + aL06 = 0.019977 + aL07 = 0.086850 + aL08 = -0.005874 + aL09 = -0.006809 + aL10 = 8.25e-4 + aL11 = 5.54e-5 + + return np.where((D > 0.01) & (D < 1000.0), + 10**(aL00 + + aL01 * np.log10(D) ** 1 + + aL02 * np.log10(D) ** 2 + + aL03 * np.log10(D) ** 3 + + aL04 * np.log10(D) ** 4 + + aL05 * np.log10(D) ** 5 + + aL06 * np.log10(D) ** 6 + + aL07 * np.log10(D) ** 7 + + aL08 * np.log10(D) ** 8 + + aL09 * np.log10(D) ** 9 + + aL10 * np.log10(D) ** 10 + + aL11 * np.log10(D) ** 11),float('nan')) + +def Rproj(D): + #Projectile SFD + aP00 = 0.0 + aP01 = -1.375458 + aP02 = 1.272521e-1 + aP03 = -1.282166 + aP04 = -3.074558e-1 + aP05 = 4.149280e-1 + aP06 = 1.910668e-1 + aP07 = -4.260980e-2 + aP08 = -3.976305e-2 + aP09 = -3.180179e-3 + aP10 = 2.799369e-3 + aP11 = 6.892223e-4 + aP12 = 2.614385e-6 + aP13 = -1.416178e-5 + aP14 = -1.191124e-6 + + return np.where((D > 1e-4) & (D < 300), + 10**(aP00 + + aP01 * np.log10(D)**1 + + aP02 * np.log10(D)**2 + + aP03 * np.log10(D)**3 + + aP04 * np.log10(D)**4 + + aP05 * np.log10(D)**5 + + aP06 * np.log10(D)**6 + + aP07 * np.log10(D)**7 + + aP08 * np.log10(D)**8 + + aP09 * np.log10(D)**9 + + aP10 * np.log10(D)**10 + + aP11 * np.log10(D)**11 + + aP12 * np.log10(D)**12 + + aP13 * np.log10(D)**13 + + aP14 * np.log10(D)**14),float('nan')) + +def N1mars(T): + # Mars crater SFD + # Ivanov (2001) SSR v. 96 pp. 87-104 + return 2.68e-14 * (np.exp(6.93 * T) - 1) + 4.13e-4 * T + +def Nmars(D): + aM00 = -3.384 + aM01 = -3.197 + aM02 = 1.257 + aM03 = 0.7915 + aM04 = -0.4861 + aM05 = -0.3630 + aM06 = 0.1016 + aM07 = 6.756e-2 + aM08 = -1.181e-2 + aM09 = -4.753e-3 + aM10 = 6.233e-4 + aM11 = 5.805e-5 + + return np.where((D > 0.01) & (D < 1000.), + 10**(aM00 + + aM01 * np.log10(D)**1 + + aM02 * np.log10(D)**2 + + aM03 * np.log10(D)**3 + + aM04 * np.log10(D)**4 + + aM05 * np.log10(D)**5 + + aM06 * np.log10(D)**6 + + aM07 * np.log10(D)**7 + + aM08 * np.log10(D)**8 + + aM09 * np.log10(D)**9 + + aM10 * np.log10(D)**10 + + aM11 * np.log10(D)**11), + np.full_like(D, np.nan)) diff --git a/python/readme.txt b/python/ctem/readme.txt similarity index 96% rename from python/readme.txt rename to python/ctem/readme.txt index e89c42cd..4a127168 100644 --- a/python/readme.txt +++ b/python/ctem/readme.txt @@ -1,22 +1,22 @@ -1) Install a distribution of Python 2.x to your machine. This source code - has been tested on Anaconda Python 4.1.1. - -2) Insure required modules are installed: - numpy - os - scipy - shutil - subprocess - -3) Place the .py files in this directory into the directory you wish - to execute from, with the required initialization files. For this, - the IDL directory may be copied, and ctem.in replaced by the one found - in this directory. - -4) Inside Anaconda Python, the program can be run through the Spyder2 GUI. - Outside of Anaconda, the progam can be begun by configuring the ctem.in - file, then typing - - python < ctem_driver.py - - at the command line. +1) Install a distribution of Python 2.x to your machine. This source code + has been tested on Anaconda Python 4.1.1. + +2) Insure required modules are installed: + numpy + os + scipy + shutil + subprocess + +3) Place the .py files in this directory into the directory you wish + to execute from, with the required initialization files. For this, + the IDL directory may be copied, and ctem.in replaced by the one found + in this directory. + +4) Inside Anaconda Python, the program can be run through the Spyder2 GUI. + Outside of Anaconda, the progam can be begun by configuring the ctem.in + file, then typing + + python < ctem_driver.py + + at the command line. diff --git a/src/crater/crater_dimensions.f90 b/src/crater/crater_dimensions.f90 index f790cb3a..78fd064b 100644 --- a/src/crater/crater_dimensions.f90 +++ b/src/crater/crater_dimensions.f90 @@ -68,8 +68,8 @@ subroutine crater_dimensions(user,crater,domain) ! Calculate the radius where the inner wall meets the original pre-existing surface ! This is used to demark the location where excavation transitions to deposition - crater%ejrad = max(crater_profile_find_r_inner_wall(user,crater) * crater%frad, crater%rad) crater%ejrim = 0.14_DP * (crater%fcrat * 0.5_DP)**(0.74_DP) ! McGetchin et al. (1973) Thickness of ejecta at rim + crater%ejrad = max(crater_profile_find_r_inner_wall(user,crater) * crater%frad, crater%rad) !find rim for counting purposes crater%frim = RIMFAC * crater%frad diff --git a/src/crater/crater_realistic_topography.f90 b/src/crater/crater_realistic_topography.f90 index ed1b29c4..115601ff 100644 --- a/src/crater/crater_realistic_topography.f90 +++ b/src/crater/crater_realistic_topography.f90 @@ -314,10 +314,10 @@ subroutine complex_wall_texture(user,surf,crater,domain,deltaMtot) integer(I4B), parameter :: num_octaves = 10 ! Number of Perlin noise octaves integer(I4B), parameter :: offset = 4000 ! Scales the random xy-offset so that each crater's random noise is unique real(DP), parameter :: xy_noise_fac = 3.0_DP ! Spatial "size" of noise features at the first octave - real(DP), parameter :: noise_height = 7.0e-3_DP ! Spatial "size" of noise features at the first octave + real(DP), parameter :: noise_height = 3.0e-3_DP ! Spatial "size" of noise features at the first octave real(DP), parameter :: freq = 2.0_DP ! Spatial size scale factor multiplier at each octave level real(DP), parameter :: pers = 1.20_DP ! The relative size scaling at each octave level - real(DP), parameter :: outer_wall_size = 2.1_DP + real(DP), parameter :: outer_wall_size = 1.2_DP !Executable code call random_number(rn) @@ -344,7 +344,7 @@ subroutine complex_wall_texture(user,surf,crater,domain,deltaMtot) areafrac = 1.0 - util_area_intersection(0.3_DP * crater%floordiam,xbar,ybar,user%pix) areafrac = areafrac * util_area_intersection(outer_wall_size * crater%frad,xbar,ybar,user%pix) areafrac = areafrac * min((r / flr)**12,1.0_DP) ! Smooth out interface between wall and floor - areafrac = areafrac * max(min(2._DP - r,1.0_DP),0.0_DP) ! Smooth out region outside of the rim + areafrac = areafrac * max(min(outer_wall_size- r,1.0_DP),0.0_DP) ! Smooth out region outside of the rim ! Add some roughness to the walls noise = 0.0_DP @@ -356,7 +356,6 @@ subroutine complex_wall_texture(user,surf,crater,domain,deltaMtot) end do newdem = newdem + noise * areafrac if (r < flr) newdem = max(newdem,crater%melev - crater%floordepth) - if (r > 1.1_DP) newdem = max(newdem,newdem + areafrac * crater%ejrim * r**(-3)) elchange = newdem - surf(xpi,ypi)%dem deltaMtot = deltaMtot + elchange @@ -405,6 +404,7 @@ subroutine complex_terrace(user,surf,crater,deltaMtot) real(DP) :: router ! The radius of the terrace outer edge real(DP) :: rinner ! The radius of the terrace outer edge real(DP) :: upshift,dfloor + real(DP) :: h_scallop_profile ! Profile of the scallop wall integer(I4B) :: num_oct_tfloor real(DP) :: noise_height_tfloor real(DP) :: freq_tfloor @@ -424,7 +424,9 @@ subroutine complex_terrace(user,surf,crater,deltaMtot) nscallops = 16 scallop_p = 0.6_DP scallop_width = 0.20_DP - rimfloor = crater%melev + h_scallop_profile = 2.0_DP + rimfloor = crater%melev + 0.5_DP * crater%rimheight + upshift = 0.0_DP num_oct_tfloor = 4 noise_height_tfloor = crater%floordepth / nterraces @@ -441,13 +443,12 @@ subroutine complex_terrace(user,surf,crater,deltaMtot) !^^^^^^^^^^^^^^^ rad = 2.0_DP * crater%frad - upshift = 0._DP inc = max(min(nint(rad / user%pix),PBCLIM*user%gridsize),1) + 1 crater%maxinc = max(crater%maxinc,inc) flr = crater%floordiam / crater%fcrat - do terrace = 1, nterraces + do terrace = 1, nterraces router = flr + (1._DP - flr) * (terrace / real(nterraces,kind=DP))**(terracefac) ! The radius of the outer edge of the terrace rinner = flr + (1._DP - flr) * ((terrace - 1) / real(nterraces,kind=DP))**(terracefac) ! The radius of the inner edge of the terrace @@ -484,7 +485,11 @@ subroutine complex_terrace(user,surf,crater,deltaMtot) isterrace = 1.0_DP - noise < hprof - tprof = crater_profile(user,crater,rinner) * (r/rinner)**2 + tnoise ! This is the floor profile that replaces the old one at each terrace + if (r < 1.0_DP) then + tprof = crater_profile(user,crater,rinner) * (r/rinner)**h_scallop_profile + tnoise ! This is the floor profile that replaces the old one at each terrace + else + tprof = (crater_profile(user,crater,rinner) + crater%ejrim * (rinner)**(-EJPROFILE)) + tnoise ! This is the floor profile that replaces the old one at each terrace + end if if (isterrace) then newdem = min(tprof,newdem) @@ -517,7 +522,7 @@ subroutine complex_terrace(user,surf,crater,deltaMtot) r = sqrt(xbar**2 + ybar**2) / crater%frad ! Make scalloped rim - znoise = scallop_width**(1._DP / (2 * scallop_p)) + znoise = (scallop_width)**(1._DP / (2 * scallop_p)) xynoise = (nscallops / PI) / crater%fcrat dnoise = util_perlin_noise(xynoise * xbar + offset * rn(1), & xynoise * ybar + offset * rn(2)) * znoise @@ -563,7 +568,7 @@ subroutine complex_terrace(user,surf,crater,deltaMtot) r = sqrt(xbar**2 + ybar**2) / crater%frad if (r > flr) then - newdem = newdem + upshift * max(min(3.0_DP - 2 * r,1.0_DP),0.0_DP) + newdem = newdem + upshift * max(min(3.0_DP - 2 * r**2,1.0_DP),0.0_DP) elchange = newdem - surf(xpi,ypi)%dem deltaMtot = deltaMtot + elchange surf(xpi,ypi)%dem = newdem diff --git a/src/util/util_diffusion_solver.f90 b/src/util/util_diffusion_solver.f90 index b9c23f45..7189dbd7 100644 --- a/src/util/util_diffusion_solver.f90 +++ b/src/util/util_diffusion_solver.f90 @@ -22,6 +22,7 @@ subroutine util_diffusion_solver(user,surf,N,indarray,kdiff,cumulative_elchange,maxhits) use module_globals use module_util, EXCEPT_THIS_ONE => util_diffusion_solver + use, intrinsic :: ieee_exceptions implicit none ! Arguments @@ -40,7 +41,12 @@ subroutine util_diffusion_solver(user,surf,N,indarray,kdiff,cumulative_elchange, integer(I4B) :: mi,mip1,mim1,mj,mjp1,mjm1,nloops real(DP),dimension(N,N) :: elchange integer(I4B),dimension(6,N,N) :: indarray_extended + logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes, fpe_quiet_modes + logical, dimension(size(IEEE_USUAL)) :: fpe_flag + call ieee_get_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily + fpe_quiet_modes(:) = .false. + call ieee_set_halting_mode(IEEE_ALL,fpe_quiet_modes) cumulative_elchange = 0.0_DP kdiffmax = maxval(kdiff) @@ -122,6 +128,8 @@ subroutine util_diffusion_solver(user,surf,N,indarray,kdiff,cumulative_elchange, cumulative_elchange = cumulative_elchange + elchange end do + call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) ! Restore the original FPE halting modes + return end subroutine util_diffusion_solver