From 101937203060be326cb9c97e87cdfc72082302d6 Mon Sep 17 00:00:00 2001 From: David Minton Date: Sat, 14 Sep 2019 23:14:08 +0200 Subject: [PATCH] removed variable fe code --- src/Makefile.am | 1 + src/Makefile.in | 4 ++ src/crater/crater_degradation_function.f90 | 46 ++++++++++++++ src/crater/crater_generate.f90 | 18 ++++++ src/crater/crater_populate.f90 | 2 +- src/crater/crater_soften.f90 | 13 ++-- src/crater/crater_subpixel_diffusion.f90 | 74 ++++++++++++++-------- src/crater/module_crater.f90 | 15 ++++- src/ejecta/ejecta_emplace.f90 | 12 ++-- 9 files changed, 141 insertions(+), 44 deletions(-) create mode 100644 src/crater/crater_degradation_function.f90 diff --git a/src/Makefile.am b/src/Makefile.am index ca1acb36..8b28057f 100755 --- 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 index 572d9828..e8c77274 100644 --- 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 new file mode 100644 index 00000000..da5ed994 --- /dev/null +++ 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 index 9b5f6704..37b2c968 100644 --- 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 index 3ec949bd..e0f60bb9 100644 --- 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 index f9d81235..d9d7e9ec 100644 --- 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 index 7f8c27a3..271a777f 100644 --- 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,7 +29,7 @@ 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 @@ -37,18 +37,19 @@ subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval,kdiff ! 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 index d6cddf87..07be6ffa 100644 --- 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 index 16c4a9c1..431612eb 100644 --- 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