Added code to test a variable fe hypothesis

MintonGroup · Sep 10, 2019 · 0687dbc · 0687dbc
1 parent 1986bd8
commit 0687dbc
Show file tree

Hide file tree

Showing 4 changed files with 58 additions and 8 deletions.
diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90
@@ -190,12 +190,6 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
          ! Do seismic shaking
          if (user%doseismic) call seismic_shake(user,surf,crater,domain)
 
-         ! Generate interior anomalous diffusion
-         !if (user%dosoftening) call crater_soften(user,surf,crater,domain)
-
-         ! Generate distal anomalous diffusion
-         !if (user%dosoftening) call crater_soften_accumulate(user,surf,crater,domain,kdiff)
-
          ! find the average height and slope at crater location
          call crater_averages(user,surf,crater)
 

diff --git a/src/crater/crater_soften.f90 b/src/crater/crater_soften.f90
@@ -40,8 +40,18 @@ subroutine crater_soften(user,surf,crater,domain)
    integer(I4B) :: inc,incsq,N,xpi,ypi,iradsq,i,j,ss,ioerr,bigi,bigj
    real(DP) :: xp,yp,fradsq,xbar,ybar,areafrac,kdiffmax,sf
 
+   real(DP) :: mfe,bfe
+
+   ! TEST VARIABLE FE MODEL
+   mfe = (user%femax - user%femin) / (log10(user%rmaxfe) - log10(user%rminfe))
+   bfe = 0.5_DP * ((user%femax + user%femin) - mfe * (log10(user%rmaxfe) + log10(user%rminfe)))
+
+   crater%fe = mfe * log10(crater%frad) + bfe
+   crater%fe = min(user%femax,crater%fe) 
+   !^^^^^^^^^^^^^^^^^^^^^^^^^^
+
    kdiffmax = user%Kd1 * crater%frad**user%psi
-   inc = int(min(user%fe * crater%frad / user%pix,SQRT2 * user%gridsize)) + 3 
+   inc = int(min(crater%fe * crater%frad / user%pix,SQRT2 * user%gridsize)) + 3 
    maxhits = (1 + (inc / (user%gridsize / 2)))**2
    crater%maxinc = max(crater%maxinc,inc)
    incsq = (inc - 2)**2
@@ -78,7 +88,7 @@ subroutine crater_soften(user,surf,crater,domain)
             xbar = xp - crater%xl 
             ybar = yp - crater%yl
 
-            areafrac = util_area_intersection(user%fe * crater%frad,xbar,ybar,user%pix)
+            areafrac = util_area_intersection(crater%fe * crater%frad,xbar,ybar,user%pix)
 
             kdiff(i,j) = kdiffmax * areafrac  ! Apply user-defined diffusive degradation to degradation region
          !end if

diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90
@@ -187,6 +187,13 @@ module module_globals
    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))
+
+   ! Testing variable fe model
+   logical  :: dovariablefe ! Testing a variable fe model
+   real(DP) :: femin  ! fe at largest crater in linear interpolation
+   real(DP) :: femax  ! fe at smallest crate in linear interpolation
+   real(DP) :: rminfe   ! crater size with lowest fe 
+   real(DP) :: rmaxfe   ! crater size with highest fe 
 
    ! Ejecta softening variables
    logical           :: dosoftening  ! Set T to use the extra crater softening model

diff --git a/src/io/io_input.f90 b/src/io/io_input.f90
@@ -89,7 +89,12 @@ subroutine io_input(infile,user)
    user%Kd1 = 0.00_DP
    user%psi = 2.0_DP
    user%fe  = 1.0_DP
+   user%dovariablefe = .false.
    user%ejecta_truncation = 10.0_DP
+   user%femax = 10.0_DP
+   user%femin = 2.0_DP
+   user%rmaxfe = 100_DP
+   user%rminfe = 400e3_DP
 
    open(unit=LUN,file=infile,status="old",iostat=ierr)
    if (ierr /= 0) then
@@ -452,6 +457,40 @@ subroutine io_input(infile,user)
             call io_get_token(line, ilength, ifirst, ilast, ierr)
             token = line(ifirst:ilast)
             read(token, *) user%fe
+
+!Used to test variable fe model
+         case ("DOVARIABLEFE")
+            ifirst = ilast + 1
+            call io_get_token(line, ilength, ifirst, ilast, ierr)
+            token = line(ifirst:ilast)
+            read(token, *) user%dovariablefe
+         case ("FEMIN")
+            ifirst = ilast + 1
+            call io_get_token(line, ilength, ifirst, ilast, ierr)
+            token = line(ifirst:ilast)
+            read(token, *) user%femin
+         case ("FEMAX")
+            ifirst = ilast + 1
+            call io_get_token(line, ilength, ifirst, ilast, ierr)
+            token = line(ifirst:ilast)
+            read(token, *) user%femax
+         case ("RMAXFE")
+            ifirst = ilast + 1
+            call io_get_token(line, ilength, ifirst, ilast, ierr)
+            token = line(ifirst:ilast)
+            read(token, *) user%rmaxfe
+         case ("RMINFE")
+            ifirst = ilast + 1
+            call io_get_token(line, ilength, ifirst, ilast, ierr)
+            token = line(ifirst:ilast)
+            read(token, *) user%rminfe
+
+
+
+
+
+
+
          !**************************************************************************
          ! The following is for backwards compatibility with older style input files
          !**************************************************************************