removed variable fe code

MintonGroup · Sep 14, 2019 · 1019372 · 1019372
1 parent 0687dbc
commit 1019372
Show file tree

Hide file tree

Showing 9 changed files with 141 additions and 44 deletions.
diff --git a/src/Makefile.am b/src/Makefile.am
@@ -81,6 +81,7 @@ crater/crater_soften_accumulate.f90\
 crater/crater_subpixel_diffusion.f90\
 crater/crater_make_list.f90\
 crater/crater_critical_slope.f90\
+crater/crater_degradation_function.f90\
 init/init_domain.f90\
 init/init_dist.f90\
 init/init_surf.f90\

diff --git a/src/Makefile.in b/src/Makefile.in
@@ -149,6 +149,7 @@ am_CTEM_OBJECTS = globals/module_globals.$(OBJEXT) \
 	crater/crater_subpixel_diffusion.$(OBJEXT) \
 	crater/crater_make_list.$(OBJEXT) \
 	crater/crater_critical_slope.$(OBJEXT) \
+	crater/crater_degradation_function.$(OBJEXT) \
 	init/init_domain.$(OBJEXT) init/init_dist.$(OBJEXT) \
 	init/init_surf.$(OBJEXT) init/init_regolith_stack.$(OBJEXT) \
 	init/init_porosity_stack.$(OBJEXT) \
@@ -387,6 +388,7 @@ crater/crater_soften_accumulate.f90\
 crater/crater_subpixel_diffusion.f90\
 crater/crater_make_list.f90\
 crater/crater_critical_slope.f90\
+crater/crater_degradation_function.f90\
 init/init_domain.f90\
 init/init_dist.f90\
 init/init_surf.f90\
@@ -685,6 +687,8 @@ crater/crater_make_list.$(OBJEXT): crater/$(am__dirstamp) \
 	crater/$(DEPDIR)/$(am__dirstamp)
 crater/crater_critical_slope.$(OBJEXT): crater/$(am__dirstamp) \
 	crater/$(DEPDIR)/$(am__dirstamp)
+crater/crater_degradation_function.$(OBJEXT): crater/$(am__dirstamp) \
+	crater/$(DEPDIR)/$(am__dirstamp)
 init/init_domain.$(OBJEXT): init/$(am__dirstamp) \
 	init/$(DEPDIR)/$(am__dirstamp)
 init/init_dist.$(OBJEXT): init/$(am__dirstamp) \

diff --git a/src/crater/crater_degradation_function.f90 b/src/crater/crater_degradation_function.f90
@@ -0,0 +1,46 @@
+!  Project     : CTEM
+!  Language    : Fortran 2003
+!
+!  Description : Calculates the degradation function from Minton et al. (2019)
+!  
+!
+!  Input
+!    Arguments : 
+!
+!  Output
+!    Arguments :
+!           
+! 
+!  Notes       :  
+!
+!**********************************************************************************************************************************
+
+
+function crater_degradation_function(user,crater) result(Kd)
+   use module_globals
+   use module_util
+   use module_crater, EXCEPT_THIS_ONE => crater_degradation_function
+   implicit none
+
+  ! Arguments
+   type(usertype),intent(in) :: user
+   type(cratertype),intent(in) :: crater
+   real(DP) :: Kd
+   real(DP) :: A,r_break,alpha_1,alpha_2,delta,r
+
+   r = crater%frad
+
+   ! Testing values that match both mare scale and highlands scale
+   alpha_1 = user%psi
+   alpha_2 = 1.1_DP
+
+   r_break = 0.10e3_DP
+   delta = 1.0_DP
+
+   A = user%Kd1 * r_break**(alpha_1) / (1_DP + (alpha_1 - alpha_2) / alpha_1)**2 
+
+   Kd =  A * (r / r_break)**(alpha_1) * (0.5_DP * (1.0_DP + (r / r_break)**(1.0_DP / delta)))**((alpha_2 - alpha_1 ) / delta)
+   return
+
+end function crater_degradation_function
+
diff --git a/src/crater/crater_generate.f90 b/src/crater/crater_generate.f90
@@ -39,6 +39,10 @@ subroutine crater_generate(user,crater,domain,prod,production_list,vdist,surf)
    integer(I4B)             :: k,khi,klo,Nk
    integer(I8B)             :: random_index,numremaining,nabove,nabovep1
 
+   !TESTING FOR VARIABLE FE MODERL
+   real(DP) :: mfe,bfe
+
+
    ! Get all six random numbers we need in one call
    if (.not.domain%initialize) call random_number(rn)
 
@@ -190,6 +194,20 @@ subroutine crater_generate(user,crater,domain,prod,production_list,vdist,surf)
 
    ! Get pixel space values
    crater%fcratpx = nint(crater%fcrat / user%pix)
+
+   ! TEST VARIABLE FE MODEL
+   if (user%dovariablefe) then
+      !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 = max(min(user%femax,crater%fe),0.d0)
+      crater%fe = max(user%femin * (crater%frad*1e-3_DP)**(0.27_DP),user%femin)
+      crater%fe = min(crater%fe * crater%frad, 3.08e6_DP) / crater%frad
+   else
+      crater%fe = user%fe
+   end if
+   !^^^^^^^^^^^^^^^^^^^^^^^^^^
    return
 end subroutine crater_generate
 
diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90
@@ -274,7 +274,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
             call io_updatePbar(message)
             craters_since_subpixel = icrater - icrater_last_subpixel
             finterval = craters_since_subpixel / real(ntotcrat,kind=DP)
-            call crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval,kdiff)
+            call crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiff)
             icrater_last_subpixel = icrater
          end if
          ! Intermediate tally step 

diff --git a/src/crater/crater_soften.f90 b/src/crater/crater_soften.f90
@@ -42,15 +42,10 @@ subroutine crater_soften(user,surf,crater,domain)
 
    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
+   !kdiffmax = user%Kd1 * crater%frad**user%psi
+   kdiffmax = crater_degradation_function(user,crater)
    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)
@@ -101,12 +96,12 @@ subroutine crater_soften(user,surf,crater,domain)
       write(*,*) 'crater kdiff '
       write(*,*) maxval(kdiff)
       call util_diffusion_solver(user,surf,2 * inc + 1,indarray,kdiff,cumulative_elchange,maxhits)
-      write(16,*) crater%frad
+      !write(16,*) crater%frad
       do j = -inc,inc
          do i = -inc,inc
             xpi = indarray(1,i,j)
             ypi = indarray(2,i,j)
-            write(16,*) i,j,xpi,ypi,kdiff(i,j),surf(xpi,ypi)%dem,cumulative_elchange(i,j)
+            !write(16,*) i,j,xpi,ypi,kdiff(i,j),surf(xpi,ypi)%dem,cumulative_elchange(i,j)
             !surf(xpi,ypi)%dem = surf(xpi,ypi)%dem + cumulative_elchange(i,j)
             !surf(xpi,ypi)%ejcov = max(surf(xpi,ypi)%ejcov + cumulative_elchange(i,j),0.0_DP)
          end do

diff --git a/src/crater/crater_subpixel_diffusion.f90 b/src/crater/crater_subpixel_diffusion.f90
@@ -19,7 +19,7 @@
 !  Notes       :  
 !
 !**********************************************************************************************************************************
-subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval,kdiffin)
+subroutine crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiffin)
    use module_globals
    use module_util
    use module_ejecta
@@ -29,26 +29,27 @@ subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval,kdiff
    ! Arguments
    type(usertype),intent(in) :: user
    type(surftype),dimension(:,:),intent(inout) :: surf
-   real(DP),dimension(:,:),intent(in) :: prod,nflux 
+   real(DP),dimension(:,:),intent(in) :: nflux 
    type(domaintype),intent(in) :: domain
    real(DP),intent(in) :: finterval
    real(DP),dimension(:,:),intent(inout) :: kdiffin
 
    ! Internal variables
    real(DP),dimension(0:user%gridsize + 1,0:user%gridsize + 1) :: cumulative_elchange,kdiff
    integer(I4B),dimension(2,0:user%gridsize + 1,0:user%gridsize + 1) :: indarray
-   integer(I4B) :: i,j,l,xpi,ypi,k,ktot,kk,ioerr,inc,incsq,iradsq,xi,xf,yi,yf,crat,imin,imax,jmin,jmax
+   integer(I4B) :: i,j,xpi,ypi,k,inc,incsq,imin,imax,jmin,jmax
    integer(I8B) :: m,N
    integer(I4B) :: maxhits = 1
    real(DP) :: dburial,lambda,dKdN,diam,radius,Area,avgejc
-   real(DP) :: dN,diam_regolith,diam_bedrock,RR
+   real(DP) :: dN,diam_regolith,diam_bedrock
    real(DP),dimension(3)    :: rn  
-   real(DP) :: superlen,rayfrac,cutout,fe,fd
+   real(DP) :: superlen,fe,fd
+   real(SP) :: cutout
    type(cratertype) :: crater
    real(DP),dimension(:,:),allocatable :: diffdistribution,ejdistribution
-   integer(I4B),dimension(:,:),allocatable :: ejisray
-   real(DP) :: xbar,ybar,dD,xp,yp,areafrac,krad,supersize,lrad
+   real(DP) :: xbar,ybar,dD,xp,yp,areafrac,krad,lrad,ebh
    integer(I8B),dimension(user%gridsize,user%gridsize) :: Ngrid
+   real(DP) ::  mfe,bfe
 
 
    ! Create box for soften calculation (will be no bigger than the grid itself)
@@ -67,7 +68,7 @@ subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval,kdiff
 
    fe = FEPROX
    fd = user%ejecta_truncation
-   if (user%dosoftening) fe = user%fe
+   if (user%dosoftening) fe = crater%fe
 
    ! Generate both the subpixel and superdomain diffusive degradation
    do k = 1,domain%pnum - 1
@@ -86,13 +87,15 @@ subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval,kdiff
 
       if ((fd * diam < user%pix) .or. (dN * PI * (fe * radius)**2 > 0.1_DP)) then 
       !Do the average degradation per pixel for the subpixel component
+
+         dKdN = 0.0_DP 
+         if (diam < user%pix)  dKdN = dKdN + KD1PROX * PI * FEPROX**2 * (radius)**(2.0_DP + PSIPROX) / domain%parea
          if (user%dosoftening) then 
          ! User-defined degradation function
-            dKdN = user%Kd1 * PI * user%fe**2 * (radius)**(2.0_DP + user%psi) / domain%parea
-         else 
-         !Empirically-derived "intrinsic" degradation function from proximal ejecta redistribution
-            dKdN = KD1PROX * PI * FEPROX**2 * (radius)**(2.0_DP + PSIPROX) / domain%parea
+            !dKdN = dKdN + user%Kd1 * PI * fe**2 * (radius)**(2.0_DP + user%psi) / domain%parea
+            dKdN = dKdN + PI * fe**2 * radius**2 * crater_degradation_function(user,crater) / domain%parea
          end if
+         !Empirically-derived "intrinsic" degradation function from proximal ejecta redistribution
 
          lambda = dN * domain%parea
 
@@ -122,15 +125,15 @@ subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval,kdiff
       ! Do the degradation as individual circles
 
          superlen = fd * diam + domain%side
-         cutout = 0.0_DP
+         cutout = 0.0_SP
          crater%continuous = RCONT * radius**(EXPCONT) 
          if (diam > domain%smallest_crater) then
             if (.not.user%superdomain) exit
             ! Superdomain craters 
-            cutout = domain%side + crater%continuous
+            cutout = real(domain%side + crater%continuous, kind=SP)
          end if
          Area = superlen**2
-         if (diam < domain%biggest_crater) Area = Area - cutout**2
+         if (diam < domain%biggest_crater) Area = Area - (1._DP * cutout)**2
 
          lambda = dN * Area
          dD = nflux(2,k+1) - nflux(2,k)
@@ -148,9 +151,9 @@ subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval,kdiff
             else
                ! Superdomain craters 
                crater%xl = real((superlen - cutout) * (rn(1) - 0.5_DP), kind=SP) 
-               if (crater%xl > 0.0_DP) crater%xl = crater%xl + cutout
+               if (crater%xl > 0.0_SP) crater%xl = crater%xl + cutout
                crater%yl = real((superlen - cutout) * (rn(2) - 0.5_DP), kind=SP) 
-               if (crater%yl > 0.0_DP) crater%yl = crater%yl + cutout
+               if (crater%yl > 0.0_SP) crater%yl = crater%yl + cutout
             end if
 
             crater%xlpx = nint(crater%xl / user%pix)
@@ -159,8 +162,21 @@ subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval,kdiff
             crater%frad = 0.5_DP * (diam + dD * rn(3))
             crater%continuous = RCONT * crater%frad**(EXPCONT) 
             krad = fd * crater%frad
+
+            if (user%dovariablefe) then
+               !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 = max(min(user%femax,crater%fe),0.0_DP) 
+               crater%fe = max(user%femin * (crater%frad*1e-3_DP)**(0.27_DP),user%femin)
+               crater%fe = min(crater%fe * crater%frad, 3.08e6_DP) / crater%frad
+            else
+               crater%fe = user%fe
+            end if
+
 
-            dKdN = user%Kd1 * crater%frad**(user%psi)
+            !dKdN = user%Kd1 * crater%frad**(user%psi)
+            dKdN = crater_degradation_function(user,crater)
             inc  = int(krad / user%pix) + 2
             incsq = inc**2
 
@@ -185,13 +201,17 @@ subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval,kdiff
                do i = imin,imax
                   xpi = crater%xlpx + i
                   ypi = crater%ylpx + j
-                  kdiff(xpi,ypi) = kdiff(xpi,ypi) + dKdN * diffdistribution(i,j)
+                  xbar = xpi * user%pix - crater%xl 
+                  ybar = ypi * user%pix - crater%yl
+                  areafrac = util_area_intersection(crater%fe * crater%frad,xbar,ybar,user%pix)
+                  kdiff(xpi,ypi) = kdiff(xpi,ypi) + dKdN * diffdistribution(i,j) * areafrac
                   !TEMP
                   xp = xpi * user%pix
                   yp = ypi * user%pix
                   lrad = sqrt((xp - crater%xl)**2 + (yp - crater%yl)**2)
-                  surf(xpi,ypi)%ejcov = surf(xpi,ypi)%ejcov + ejdistribution(i,j) &
-                  * 0.14_DP * crater%frad**(0.74_DP) * (lrad / crater%frad)**(-3.0_DP)
+                  ebh =  0.14_DP * crater%frad**(0.74_DP) * (lrad / crater%frad)**(-3.0_DP) * ejdistribution(i,j)
+                  surf(xpi,ypi)%ejcov = surf(xpi,ypi)%ejcov + ebh
+                  kdiff(xpi,ypi) = kdiff(xpi,ypi) +  1.5_DP * ebh**2 * ejdistribution(i,j)
                end do
             end do
             !$OMP END PARALLEL DO
@@ -208,12 +228,12 @@ subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval,kdiff
    kdiff(user%gridsize + 1,:) = kdiff(1,:)
    kdiff(:,user%gridsize + 1) = kdiff(:,1)  
 
-   write(*,*)
-   write(*,*) 'avgkdiff = ',sum(kdiff) / (user%gridsize + 2)**2 / finterval
-   write(*,*)
-   open(unit=55,file='avgkdiff.dat',status='unknown',position='append')
-   write(55,*) sum(kdiff) / (user%gridsize + 2)**2 / finterval
-   close(55)
+   !write(*,*)
+   !write(*,*) 'avgkdiff = ',sum(kdiff) / (user%gridsize + 2)**2 / finterval
+   !write(*,*)
+   !open(unit=55,file='avgkdiff.dat',status='unknown',position='append')
+   !write(55,*) sum(kdiff) / (user%gridsize + 2)**2 / finterval
+   !close(55)
 
 
    call util_diffusion_solver(user,surf,user%gridsize + 2,indarray,kdiff,cumulative_elchange,maxhits)

diff --git a/src/crater/module_crater.f90 b/src/crater/module_crater.f90
@@ -229,12 +229,12 @@ subroutine crater_soften_accumulate(user,surf,crater,domain,kdiff)
       end subroutine crater_soften_accumulate
    end interface
    interface
-      subroutine crater_subpixel_diffusion(user,surf,prod,crtscl,domain,finterval,kdiffin)
+      subroutine crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiffin)
       use module_globals
       implicit none
       type(usertype),intent(in) :: user
       type(surftype),dimension(:,:),intent(inout) :: surf
-      real(DP),dimension(:,:),intent(in) :: prod,crtscl
+      real(DP),dimension(:,:),intent(in) :: nflux
       type(domaintype),intent(in)    :: domain
       real(DP),intent(in) :: finterval
       real(DP),dimension(:,:),intent(inout) :: kdiffin
@@ -262,4 +262,15 @@ function crater_critical_slope(user,crater,lrad) result(critical)
       end function crater_critical_slope
    end interface
 
+   interface
+      function crater_degradation_function(user,crater) result(Kd)
+      use module_globals
+      type(usertype),intent(in) :: user
+      type(cratertype),intent(in) :: crater
+      real(DP) :: Kd
+      end function crater_degradation_function
+   end interface
+
+
+
 end module
diff --git a/src/ejecta/ejecta_emplace.f90 b/src/ejecta/ejecta_emplace.f90
@@ -134,13 +134,14 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot)
 
    ! Increase the box a bit to take into account possible ejecta pattern distortion due to topography
    inc = ceiling(inc * 1.5_DP)
+   krad = user%ejecta_truncation * crater%frad
+   dradsq = int(krad / user%pix) + 3
+   inc = max(inc,dradsq)
+   dradsq = dradsq**2
 
    if (user%dosoftening) then
-      krad = user%ejecta_truncation * crater%frad
-      kdiffmax = user%Kd1 * crater%frad**(user%psi)
-      dradsq = int(krad / user%pix) + 3
-      inc = max(inc,dradsq)
-      dradsq = dradsq**2
+      !kdiffmax = user%Kd1 * crater%frad**(user%psi)
+      kdiffmax = crater_degradation_function(user,crater)
    end if
 
    crater%maxinc = max(crater%maxinc,inc)
@@ -259,6 +260,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot)
          if (user%dosoftening) then
             ! Do extra diffusive degradation over ejecta region
             areafrac =  (1.0_DP - util_area_intersection(crater%frad,xbar,ybar,user%pix)) 
+            areafrac = areafrac * util_area_intersection(crater%fe * crater%frad,xbar,ybar,user%pix)
             kdiff(i,j) = areafrac * diffdistribution(idistorted,jdistorted) * kdiffmax
          end if