Removed code that makes the ejecta go infinite distance for velocitie…

…s above escape, as it wreaks havoc on the root finder
MintonGroup · Oct 28, 2021 · 961b8f4 · 961b8f4
1 parent 1ae66e8
commit 961b8f4
Show file tree

Hide file tree

Showing 3 changed files with 5 additions and 14 deletions.
diff --git a/src/Makefile.am b/src/Makefile.am
@@ -1,12 +1,12 @@
 bin_PROGRAMS = CTEM 
 #ifort optimized flags
-#AM_FCFLAGS =  -O3 -openmp -parallel -xHost -ipo -assume byterecl -m64 -heap-arrays -FR 
-#AM_FCFLAGS =  -O3 -openmp -parallel -assume byterecl -m64 -heap-arrays -FR 
+#AM_FCFLAGS =  -O3 -qopenmp -parallel -xHost -ipo -assume byterecl -m64 -heap-arrays -FR 
+AM_FCFLAGS =  -O3 -qopenmp -parallel -assume byterecl -m64 -heap-arrays -FR 
 #ifort debug flags
-AM_FCFLAGS = -O3 -p -g -qopenmp -debug all -traceback -CB -assume byterecl -m64 -heap-arrays -FR
+#AM_FCFLAGS = -O0 -p -g -debug all -traceback -CB -assume byterecl -m64 -heap-arrays -FR
 
 #gfortran optimized flags
-AM_FCFLAGS =  -O3 -fopenmp -ffree-form  -g -fbounds-check -fbacktrace
+#AM_FCFLAGS =  -O3 -fopenmp -ffree-form  -g -fbounds-check -fbacktrace
 #gfortran debug flags
 #AM_FCFLAGS =  -O0 -g -fopenmp -fbounds-check -Wall -Warray-bounds -Warray-temporaries -Wimplicit-interface -ffree-form
 

diff --git a/src/ejecta/ejecta_blanket.f90 b/src/ejecta/ejecta_blanket.f90
@@ -76,10 +76,6 @@ subroutine ejecta_blanket(user,crater,domain,erad,lrad,vejsq,ejang,firstrun)
    vej1sq = max(ce1sq * abs(user%gaccel * crater%grad) * (crater%grad / erad)**p,0._DP)
    gterm = ce1sq * (user%gaccel * erad)
    vejsq = max(vej1sq - gterm - yterm,0._DP)
-   if (vejsq >= domain%vescsq) then 
-      lrad = VBIG
-      return
-   end if
 
    ! Ejection angle
    ejang=DEG2RAD * (55._DP - 20 * (erad / crater%grad))

diff --git a/src/ejecta/ejecta_table_define.f90 b/src/ejecta/ejecta_table_define.f90
@@ -46,7 +46,7 @@ subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt)
    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)
+   lrad = crater%rad 
    erad = crater%rad 
    ejtble = EJBTABSIZE
    firstrun = .true.
@@ -66,12 +66,9 @@ subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt)
       call regolith_melt_zone(user,crater,dimp,vimp,rmelt,depthb)
    end if
 
-   !write(*,*) '     lrad/Df                   vej                       ebh                       melt & 
-   !           fraction              melt thickness'
    do k = 0,EJBTABSIZE
       call ejecta_rootfind(user,crater,domain,erad,lrad,vejsq,ejang,firstrun)
       if (k >= 1) then
-         !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
@@ -87,7 +84,6 @@ subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt)
             call regolith_melt_fraction(dimp,depthb,erad,eradold,rmelt,melt)
             ejb(k)%meltfrac = melt
          end if
-         !write(*,*) lrad/crater%rad,erad/crater%rad,sqrt(vejsq),thick !,melt,thick*melt
          if ((thick <= VSMALL) .or. (abs(eradold - erad) < VSMALL)) then
             ejtble = k
             crater%ejdis = lrad
@@ -97,7 +93,6 @@ subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt)
       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
    ! Get pixel space distance
    crater%ejdispx = nint(crater%ejdis / user%pix)