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

@@ -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)

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)

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

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

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

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)

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