From dfbe54f3a3d9505f49d3cb5f68fda32f1e1c45c4 Mon Sep 17 00:00:00 2001 From: David Minton Date: Sun, 11 Feb 2024 10:52:50 -0500 Subject: [PATCH] Put the ejecta table values in logspace and fixed bugs due to trying to compute melt values when doregotrack was .false. --- src/crater/crater_populate.f90 | 2 -- src/ejecta/ejecta_emplace.f90 | 17 +++++++++++------ src/ejecta/ejecta_interpolate.f90 | 15 +++++++++++++-- src/ejecta/ejecta_table_define.f90 | 6 +++--- src/globals/module_globals.f90 | 2 +- src/io/io_ejecta_table.f90 | 2 +- src/regolith/regolith_melt_fraction.f90 | 6 ++++++ 7 files changed, 35 insertions(+), 15 deletions(-) diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index 59c6e20e..eb858938 100644 --- a/src/crater/crater_populate.f90 +++ b/src/crater/crater_populate.f90 @@ -315,10 +315,8 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt if (crater%ejdis > domain%smallest_ejecta) then ! Estimated size is big enough, so proceed with precise calculation if (user%doregotrack) then call ejecta_table_define(user,crater,domain,ejb,ejtble,melt) - !call ejecta_interpolate(crater,domain,crater%frad,ejb(1:ejtble),ejtble,crater%ejrim) else call ejecta_table_define(user,crater,domain,ejb,ejtble) - !call ejecta_interpolate(crater,domain,crater%frad,ejb(1:ejtble),ejtble,crater%ejrim) end if call ejecta_emplace(user,surf,crater,domain,ejb(1:ejtble),ejtble,ejbmass,& ejecta_dem,nmeltsheet,vmeltsheet) diff --git a/src/ejecta/ejecta_emplace.f90 b/src/ejecta/ejecta_emplace.f90 index 9cc16e41..b6ec3381 100644 --- a/src/ejecta/ejecta_emplace.f90 +++ b/src/ejecta/ejecta_emplace.f90 @@ -253,13 +253,16 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ maxslp = -huge(maxslp) klo = int((log(lrad) - log(crater%ejrad)) / domain%ejbres) do n = 1,MAXLOOP - call ejecta_interpolate(crater,domain,distance,ejb,ejtble,ebh,vsq=vsq,theta=ejtheta,erad=erad,melt=melt) + if (user%doregotrack) then + call ejecta_interpolate(crater,domain,distance,ejb,ejtble,ebh,vsq=vsq,theta=ejtheta,erad=erad,melt=melt) + else + call ejecta_interpolate(crater,domain,distance,ejb,ejtble,ebh,vsq=vsq,theta=ejtheta,erad=erad) + end if if (n > 1) then if ((abs(ebh0 - ebh) / ebh0) < domain%small) exit endif ebh0 = ebh - erad = exp(erad) lrange = lrad - erad baseline = ((i * crater%xslp) + (j * crater%yslp)) * user%pix @@ -378,10 +381,12 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ end do end if - if (totmelt > vmelt) then - vmeltsheet = totmelt - vmelt - else !give the craters a melt sheet of 1m - vmeltsheet = 1.0_DP * user%pix * user%pix * nmeltsheet + if (user%doregotrack) then + if (totmelt > vmelt) then + vmeltsheet = totmelt - vmelt + else !give the craters a melt sheet of 1m + vmeltsheet = 1.0_DP * user%pix * user%pix * nmeltsheet + end if end if ! Create box for soften calculation (will be no bigger than the grid itself) diff --git a/src/ejecta/ejecta_interpolate.f90 b/src/ejecta/ejecta_interpolate.f90 index c71d50a0..3463c08e 100644 --- a/src/ejecta/ejecta_interpolate.f90 +++ b/src/ejecta/ejecta_interpolate.f90 @@ -70,8 +70,19 @@ subroutine ejecta_interpolate(crater,domain,lrad,ejb,ejtble,ebh,vsq,theta,erad,m if (present(erad)) erad = ejb(k)%erad - ((ejb(k)%erad - ejb(k+1)%erad) * frac) end if ebh = exp(ebh) - if (lrad > crater%ejdis) ebh = 0._DP - + if (lrad > crater%ejdis) then + ebh = 0._DP + if (present(vsq)) vsq = 0._DP + if (present(theta)) theta = 0._DP + if (present(melt)) melt = 0._DP + if (present(erad)) erad = 0._DP + else + if (present(vsq)) vsq = exp(vsq) + if (present(theta)) theta = exp(theta) + if (present(melt)) melt = exp(melt) + if (present(erad)) erad = exp(erad) + end if + return end subroutine ejecta_interpolate diff --git a/src/ejecta/ejecta_table_define.f90 b/src/ejecta/ejecta_table_define.f90 index cf259e2d..408d9c7e 100644 --- a/src/ejecta/ejecta_table_define.f90 +++ b/src/ejecta/ejecta_table_define.f90 @@ -80,12 +80,12 @@ subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt) thick = max(crater_profile(user,crater,r),VSMALL) end if ejb(k)%thick = log(thick) - ejb(k)%vesq = vejsq - ejb(k)%angle = ejang + ejb(k)%vesq = log(vejsq) + ejb(k)%angle = log(ejang) ejb(k)%erad = log(erad) if (present(melt)) then call regolith_melt_fraction(dimp,depthb,erad,eradold,rmelt,melt) - ejb(k)%meltfrac = melt + ejb(k)%meltfrac = log(melt) end if if ((thick <= VSMALL) .or. (abs(eradold - erad) < VSMALL)) then ejtble = k diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index 157e1609..a1d38273 100644 --- a/src/globals/module_globals.f90 +++ b/src/globals/module_globals.f90 @@ -45,7 +45,7 @@ module module_globals integer(I4B), parameter :: UPPERCASE_OFFSET = iachar('A') - iachar('a') ! Miscellaneous constants: -real(DP),parameter :: VSMALL = tiny(1._DP) ! Very small number +real(DP),parameter :: VSMALL = 10*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 diff --git a/src/io/io_ejecta_table.f90 b/src/io/io_ejecta_table.f90 index 342af640..f2f178a0 100644 --- a/src/io/io_ejecta_table.f90 +++ b/src/io/io_ejecta_table.f90 @@ -40,7 +40,7 @@ subroutine io_ejecta_table(crater,domain,ejb,ejtble,filename) write(LUN,'("# ejrim = ",ES12.5, " ejdis = ",ES12.5," imp = ",ES12.5)') crater%ejrim,crater%ejdis,crater%imp write(LUN,'(A63)') '# "r (m)" "h (m)" "v (m/s)" "ang (deg)" "erad (m)"' do k=1,ejtble - write(LUN,'(5(ES13.5E3,1X))') exp(ejb(k)%lrad),exp(ejb(k)%thick),sqrt(ejb(k)%vesq),ejb(k)%angle/DEG2RAD, & + write(LUN,'(5(ES13.5E3,1X))') exp(ejb(k)%lrad),exp(ejb(k)%thick),sqrt(exp(ejb(k)%vesq)),exp(ejb(k)%angle)/DEG2RAD, & exp(ejb(k)%erad) end do close(LUN) diff --git a/src/regolith/regolith_melt_fraction.f90 b/src/regolith/regolith_melt_fraction.f90 index 4886d9fc..01dac4ac 100755 --- a/src/regolith/regolith_melt_fraction.f90 +++ b/src/regolith/regolith_melt_fraction.f90 @@ -102,5 +102,11 @@ subroutine regolith_melt_fraction(dimp,depthb,erad1,erad2,rmelt,meltfrac) write(*,*) 'regolith_melt_fraction: this is a bug!' end if + if (meltfrac < epsilon(1.0_DP)) then + meltfrac = 0.0_DP + else if (meltfrac > 1.0_DP) then + meltfrac = 1.0_DP + end if + return end subroutine regolith_melt_fraction