From 1a1f03c73df4faf759227e3f2a771b396eb1b245 Mon Sep 17 00:00:00 2001 From: Austin Blevins Date: Thu, 26 Oct 2023 15:11:44 -0400 Subject: [PATCH] fixed a bug with the melt index arrays --- src/crater/crater_subpixel_diffusion.f90 | 11 +- src/crater/crater_superdomain.f90 | 26 +- src/crater/module_crater.f90 | 60 +--- src/ejecta/ejecta_emplace.f90 | 270 +++++++++--------- src/ejecta/ejecta_ray_pattern.f90 | 105 ++++--- src/ejecta/module_ejecta.f90 | 9 +- src/regolith/module_regolith.f90 | 2 +- src/regolith/regolith_subpixel_streamtube.f90 | 2 +- 8 files changed, 244 insertions(+), 241 deletions(-) diff --git a/src/crater/crater_subpixel_diffusion.f90 b/src/crater/crater_subpixel_diffusion.f90 index 3da7b696..630d4401 100644 --- a/src/crater/crater_subpixel_diffusion.f90 +++ b/src/crater/crater_subpixel_diffusion.f90 @@ -45,7 +45,7 @@ subroutine crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiffin) real(DP),dimension(3) :: rn real(DP) :: superlen,rayfrac,cutout,fe,fd type(cratertype) :: crater - real(DP) :: eji, diffi + real(DP),dimension(:,:),allocatable :: diffdistribution,ejdistribution integer(I4B),dimension(:,:),allocatable :: ejisray real(DP) :: xbar,ybar,dD,xp,yp,areafrac,krad,supersize,lrad integer(I8B),dimension(user%gridsize,user%gridsize) :: Ngrid @@ -174,6 +174,9 @@ subroutine crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiffin) jmax = ypi - crater%ylpx + allocate(diffdistribution(imin:imax,jmin:jmax)) + allocate(ejdistribution(imin:imax,jmin:jmax)) + call ejecta_ray_pattern(user,surf,crater,inc,imin,imax,jmin,jmax,diffdistribution,ejdistribution) ! Loop over affected matrix area !!$OMP PARALLEL DO DEFAULT(SHARED) IF(inc > INCPAR) & !!$OMP FIRSTPRIVATE(jmin,jmax,imin,imax) & @@ -182,17 +185,17 @@ subroutine crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiffin) do i = imin,imax xpi = crater%xlpx + i ypi = crater%ylpx + j - call ejecta_ray_pattern(user,crater,i,j,diffi,eji) - kdiff(xpi,ypi) = kdiff(xpi,ypi) + dKdN * diffi + kdiff(xpi,ypi) = kdiff(xpi,ypi) + dKdN * diffdistribution(i,j) !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 + eji & + surf(xpi,ypi)%ejcov = surf(xpi,ypi)%ejcov + ejdistribution(i,j) & * 0.14_DP * crater%frad**(0.74_DP) * (lrad / crater%frad)**(-3) end do end do !!$OMP END PARALLEL DO + deallocate(diffdistribution,ejdistribution) end do end if diff --git a/src/crater/crater_superdomain.f90 b/src/crater/crater_superdomain.f90 index 997bd473..1d1f6092 100644 --- a/src/crater/crater_superdomain.f90 +++ b/src/crater/crater_superdomain.f90 @@ -42,8 +42,8 @@ subroutine crater_superdomain(user,surf,prod,nflux,domain,finterval) real(DP),dimension(2) :: rn real(DP) :: superlen,rayfrac type(cratertype) :: crater - real(DP) :: diffi, eji - logical :: ejisray + real(DP),dimension(:,:),allocatable :: ejdistribution, diffdistribution + integer(I4B),dimension(:,:),allocatable :: ejisray ! Melt or glassy ray test real(DP) :: xp, yp, lradsq, lrad, erad @@ -144,26 +144,40 @@ subroutine crater_superdomain(user,surf,prod,nflux,domain,finterval) lrad = minval(lradif) if (lrad > crater%ejdis) cycle + allocate(ejdistribution(xi:xf,yi:yf)) + allocate(diffdistribution(xi:xf,yi:yf)) + allocate(ejisray(xi:xf,yi:yf)) ! Now generate ray pattern + call ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,diffdistribution,ejdistribution) ! Now if doregotrack is on, do melt zone calculation if (user%doregotrack) call regolith_melt_zone_superdomain(user,crater,domain,rm,depthb) + do j = yi,yf + do i = xi,xf + ! Ejecta ray distribution + if (ejdistribution(i,j) > 0.0001_DP) then + ejisray(i,j) = 1 + else + ejisray(i,j) = 0 + end if + end do + end do + if (user%doregotrack) then do j = 1, user%gridsize do i = 1, user%gridsize xpi = i - crater%xlpx ypi = j - crater%ylpx if ((abs(xpi) > inc) .or. (abs(ypi) > inc)) cycle - call ejecta_ray_pattern(user,crater,i,j,diffi,eji) - ejisray = (eji > 0.0001_DP) - if (.not.ejisray) cycle - call regolith_superdomain(user,crater,domain,surf(i,j)%regolayer,eji,& + if (ejisray(xpi,ypi) == 0) cycle + call regolith_superdomain(user,crater,domain,surf(i,j)%regolayer,ejdistribution(xpi,ypi),& i,j,rm,depthb) end do end do end if + deallocate(ejdistribution,diffdistribution,ejisray) end do diff --git a/src/crater/module_crater.f90 b/src/crater/module_crater.f90 index 1c8e9253..c0dd5dfe 100644 --- a/src/crater/module_crater.f90 +++ b/src/crater/module_crater.f90 @@ -33,7 +33,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt mass,fracdone,nflux,ntotcrat,curyear,rclist) use module_globals implicit none - type(usertype),intent(inout) :: user + type(usertype),intent(in) :: user type(surftype),dimension(:,:),intent(inout) :: surf type(cratertype),intent(inout) :: crater type(domaintype),intent(inout) :: domain @@ -338,62 +338,4 @@ subroutine crater_superdomain(user,surf,prod,nflux,domain,finterval) end subroutine crater_superdomain end interface - - - -! new subroutines added by jundu on 10/15/2022 - - -! Rim_crest: - - ! interface - - ! subroutine Calculate_am_wl_phase_from_diameter(psd_1D,amplitude,wavelength,phase) - ! use module_globals - ! implicit none - ! ! in and out - ! type(psdtype),intent(inout) :: psd_1D - ! real(DP),dimension(:), allocatable,intent(out) :: amplitude,wavelength,phase - ! end subroutine Calculate_am_wl_phase_from_diameter - - ! subroutine Calculate_breakpoint_slope_from_diameter(psd_1D) - ! use module_globals - ! implicit none - ! ! in and out - ! type(psdtype),intent(inout) :: psd_1D - ! end subroutine Calculate_breakpoint_slope_from_diameter - - ! subroutine Calculate_targetPSD_from_breakpoint_slope(psd_1D,wavelength,psd) - ! use module_globals - ! implicit none - ! !in and out - ! type(psdtype), intent(in) :: psd_1D - ! real(DP) ,dimension(:),allocatable,intent(out) :: wavelength,psd - ! end subroutine Calculate_targetPSD_from_breakpoint_slope - - ! subroutine Calculate_am_wl_phase_from_targetPSD(psd_1D,wavelength,psd,amplitude,phase) - ! use module_globals - ! implicit none - ! !in and out - ! type(psdtype), intent(in) :: psd_1D - ! real(DP) ,dimension(:),intent(in) :: wavelength,psd - ! real(DP) ,dimension(:),allocatable,intent(out) :: amplitude,phase - ! end subroutine Calculate_am_wl_phase_from_targetPSD - - - ! subroutine Create_rim(arc_length,psd_1D,amplitude,wavelength,phase,rim_parameter) - ! use module_globals - ! implicit none - ! ! in and out - ! real(DP),intent(in) :: arc_length - ! type(psdtype), intent(in) :: psd_1D - ! real(DP),dimension(:),intent(in) :: amplitude - ! real(DP),dimension(:),intent(in) :: wavelength - ! real(DP),dimension(:),intent(in) :: phase - ! real(DP),intent(out) :: rim_parameter - ! end subroutine Create_rim - - ! end interface - - end module diff --git a/src/ejecta/ejecta_emplace.f90 b/src/ejecta/ejecta_emplace.f90 index 8f97fd41..20706f1c 100644 --- a/src/ejecta/ejecta_emplace.f90 +++ b/src/ejecta/ejecta_emplace.f90 @@ -1,4 +1,4 @@ -!***** ejecta/ejecta_emplace +!****f* ejecta/ejecta_emplace ! Name ! ejecta_emplace -- Calculate ejecta mass during excavation stage. ! SYNOPSIS @@ -93,25 +93,28 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ integer(I4B),intent(in) :: ejtble type(ejbtype),dimension(:),intent(inout) :: ejb real(DP),intent(in) :: deltaMtot - real(DP),dimension(:,:),allocatable,intent(inout) :: cumulative_elchange + real(DP),dimension(:,:),allocatable,intent(out) :: cumulative_elchange integer(I4B),intent(in) :: nmeltsheet real(DP),intent(out) :: vmeltsheet ! Internal variables real(DP) :: lrad,lradsq integer(I4B),parameter :: MAXLOOP = 100 ! Maximum number of times to loop the ejecta angle correction calculation - integer(I4B) :: xpi,ypi,i,j,n,inc,incsq,iradsq,idistorted,jdistorted + integer(I4B) :: xpi,ypi,i,j,k,n,inc,incsq,iradsq,idistorted,jdistorted real(DP) :: xp,yp,fradsq,fradpxsq,radsq,ebh,ejdissq,ejbmass,fmasscons,areafrac,xbar,ybar,krad,kdiffmax - real(DP),dimension(:,:),allocatable :: kdiff,cel - integer(I4B),dimension(:,:,:),allocatable :: indarray - integer(I4B) :: maxhits,nin,nnot,dradsq + real(DP),dimension(:,:),allocatable :: big_cumulative_elchange,kdiff,big_kdiff,cel,big_cel + integer(I4B),dimension(:,:,:),allocatable :: indarray,big_indarray + real(DP),dimension(:,:),allocatable :: ejdistribution,diffdistribution + integer(I4B) :: bigi,bigj,maxhits,nin,nnot,dradsq character(len=MESSAGESIZE) :: message ! message for the progress bar real(DP) :: vmelt, totmelt, volm - real(DP) :: diffi,eji - logical :: bigej + + + ! Ray mixing model variables + real(DP) :: dsc ! Melt zone's radius - real(DP) :: rm, dm, melt + real(DP) :: rm, dm, melt, eradc ! Ejecta pattern distortion parameters real(DP) :: distance,erad,craterslope,landslope,baseline,lrange,frac,ejheight,ebh0,maxdistance @@ -119,6 +122,10 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ real(DP) :: vsq, ejtheta integer(I4B) :: ind,klo + ! Age + real(SP) :: age_mean + + ! Executable code if (user%doregotrack) call regolith_melt_zone(user,crater,crater%imp,crater%impvel,rm,dm,totmelt) @@ -152,39 +159,39 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ incsq = inc**2 - bigej = (inc >= user%gridsize / 2) - if (bigej) then - allocate(cumulative_elchange(0:user%gridsize+1,0:user%gridsize+1)) - allocate(cel(0:user%gridsize+1,0:user%gridsize+1)) - allocate(kdiff(0:user%gridsize+1,0:user%gridsize+1)) - allocate(indarray(2,0:user%gridsize+1,0:user%gridsize+1)) - cumulative_elchange(:,:) = 0.0_DP - cel(:,:) = 0.0_DP - kdiff(:,:) = 0.0_DP - + if (inc >= user%gridsize / 2) then if (user%testflag) then - write(*,*) 'Big ejecta: fcrat =',crater%fcrat, ' Ej/S =',(crater%ejdispx*user%pix)/domain%side, ' Ejrim =', crater%ejrim - else + write(*,*) 'Big ejecta: fcrat =',crater%fcrat, ' Ej/S =',(crater%ejdispx*user%pix)/domain%side, ' Ejrim =', crater%ejrim + else write(message,'("Ejb: Dc=",ES9.2," Ej/S=",F0.3)') crater%fcrat,(crater%ejdispx*user%pix)/domain%side call io_updatePbar(message) - end if - else - allocate(cumulative_elchange(-inc:inc,-inc:inc)) - allocate(cel(-inc:inc,-inc:inc)) - allocate(kdiff(-inc:inc,-inc:inc)) - allocate(indarray(2,-inc:inc,-inc:inc)) - end if + end if + endif + + allocate(cumulative_elchange(-inc:inc,-inc:inc)) + allocate(cel(-inc:inc,-inc:inc)) + allocate(kdiff(-inc:inc,-inc:inc)) + allocate(indarray(2,-inc:inc,-inc:inc)) + allocate(ejdistribution(-inc:inc,-inc:inc)) + allocate(diffdistribution(-inc:inc,-inc:inc)) cumulative_elchange = 0.0_DP kdiff = 0.0_DP indarray = inc - 1 ! initialize this array to point to a corner (this should have 0 elevation change since we're only doing work ! within a circle of radius irad + ! *************************** Continuous Ejecta Formula *****************************! + call ejecta_ray_pattern(user,surf,crater,inc,-inc,inc,-inc,inc,diffdistribution,ejdistribution) ejbmass = 0.0_DP nin = 0 nnot = 0 + !!$OMP PARALLEL DO DEFAULT(PRIVATE) IF(inc > INCPAR) & + !!$OMP SHARED(user,domain,crater,surf,ejb,ejtble) & + !!$OMP SHARED(inc,incsq) & + !!$OMP SHARED(cumulative_elchange,kdiff,kdiffmax,indarray,ejdistribution,diffdistribution) + !open(74, file='meltvserad.csv', status='replace') do j = -inc,inc do i = -inc,inc ! find distance from crater center @@ -199,10 +206,8 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ ! periodic boundary conditions call util_periodic(xpi,ypi,user%gridsize) - if (.not.bigej) then - indarray(1,i,j) = xpi - indarray(2,i,j) = ypi - end if + indarray(1,i,j) = xpi + indarray(2,i,j) = ypi lradsq = (crater%xl - xp)**2 + (crater%yl - yp)**2 lrad = sqrt(lradsq) @@ -246,68 +251,52 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ if (vsq < 0.0_DP) cycle if (distance /= distance) cycle - idistorted = int(i * distance / lrad) if (abs(idistorted) > inc) cycle jdistorted = int(j * distance / lrad) if (abs(jdistorted) > inc) cycle - + iradsq = idistorted**2 + jdistorted**2 if ((iradsq > incsq).or.(distance <= crater%ejrad)) cycle - call ejecta_ray_pattern(user,crater,idistorted,jdistorted,diffi,eji) ! we need to cut a hole out from the inside of the crater xbar = xpi * user%pix - crater%xl ybar = ypi * user%pix - crater%yl areafrac = (1.0_DP - util_area_intersection(crater%ejrad,xbar,ybar,user%pix)) - ebh = areafrac * eji * ebh - if (bigej) then - cumulative_elchange(xpi,ypi) = cumulative_elchange(xpi,ypi) + ebh + crater_profile(user, crater, lrad) - else - cumulative_elchange(i,j) = ebh + crater_profile(user, crater, lrad) - end if + ebh = areafrac * ejdistribution(idistorted,jdistorted) * ebh + cumulative_elchange(i,j) = ebh + crater_profile(user, crater, lrad) 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) - if (bigej) then - kdiff(xpi,ypi) = kdiff(xpi,ypi) + areafrac * diffi * kdiffmax - else - kdiff(i,j) = areafrac * diffi * kdiffmax - end if - + kdiff(i,j) = areafrac * diffdistribution(idistorted,jdistorted) * kdiffmax end if end do end do + !close(74) + !!$OMP END PARALLEL DO + ! if(user%doregotrack .and. user%testflag) then + ! write(*,*) 'Ejected Melt: ', vmelt + ! write(*,*) 'Total Melt: ', totmelt + ! write(*,*) 'ejected / total melt:', vmelt/totmelt + ! end if ejbmass = sum(cumulative_elchange) ! Create buffer to prevent infinite hole bug - if (bigej) then - kdiff(0,:) = 0.0_DP - kdiff(user%gridsize+1,:) = 0.0_DP - kdiff(:,0) = 0.0_DP - kdiff(:,user%gridsize+1) = 0.0_DP - - cumulative_elchange(0,:) = 0.0_DP - cumulative_elchange(user%gridsize+1,:) = 0.0_DP - cumulative_elchange(:,0) = 0.0_DP - cumulative_elchange(:,user%gridsize+1) = 0.0_DP - else - kdiff(-inc,:) = 0.0_DP - kdiff(inc,:) = 0.0_DP - kdiff(:,-inc) = 0.0_DP - kdiff(:,inc) = 0.0_DP - - cumulative_elchange(-inc,:) = 0.0_DP - cumulative_elchange(inc,:) = 0.0_DP - cumulative_elchange(:,-inc) = 0.0_DP - cumulative_elchange(:,inc) = 0.0_DP - end if + kdiff(-inc,:) = 0.0_DP + kdiff(inc,:) = 0.0_DP + kdiff(:,-inc) = 0.0_DP + kdiff(:,inc) = 0.0_DP + + cumulative_elchange(-inc,:) = 0.0_DP + cumulative_elchange(inc,:) = 0.0_DP + cumulative_elchange(:,-inc) = 0.0_DP + cumulative_elchange(:,inc) = 0.0_DP ! Do mass conservation by adjusting ejecta thickness fmasscons = (-deltaMtot)/ ejbmass @@ -331,68 +320,40 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ ! periodic boundary conditions call util_periodic(xpi,ypi,user%gridsize) - - if (.not.bigej) then - indarray(1,i,j) = xpi - indarray(2,i,j) = ypi - end if + + indarray(1,i,j) = xpi + indarray(2,i,j) = ypi lradsq = (crater%xl - xp)**2 + (crater%yl - yp)**2 lrad = sqrt(lradsq) if (lrad < crater%ejrad) cycle - - if (bigej) then - ebh = cumulative_elchange(xpi,ypi) - crater_profile(user, crater, lrad) - else - ebh = cumulative_elchange(i,j) - crater_profile(user, crater, lrad) - end if - - if (user%doregotrack .and. ebh>1.0e-8_DP) then - call regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,lrad,ebh,rm,vsq,volm) - vmelt = vmelt + volm - end if + + + + ebh = cumulative_elchange(i,j) - crater_profile(user, crater, lrad) + + + if (user%doregotrack .and. ebh>1.0e-8_DP) then + call regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,lrad,ebh,rm,vsq,volm) + vmelt = vmelt + volm + !write(74,*) erad, surf(xpi,ypi)%regolayer%regodata%meltfrac + end if end do end do end if - 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 + 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 - ! extra soften calculation - if (user%dosoftening) then - cel = 0.0_DP - - if (bigej) then - call util_diffusion_solver(user,surf,user%gridsize + 2,indarray,kdiff,cel,maxhits) - do ypi = 1,user%gridsize - do xpi = 1,user%gridsize - surf(xpi,ypi)%dem = surf(xpi,ypi)%dem + cel(xpi,ypi) - surf(xpi,ypi)%ejcov = max(surf(xpi,ypi)%ejcov + cel(xpi,ypi),0.0_DP) - indarray(1,xpi,ypi) = xpi - indarray(2,xpi,ypi) = ypi - end do - end do - indarray(1:2,0,:) = user%gridsize - indarray(1:2,user%gridsize+1,:) = 1 - indarray(1:2,:,0) = user%gridsize - indarray(1:2,:,user%gridsize+1) = 1 - - call ejecta_soften(user,surf,user%gridsize + 2,indarray,cumulative_elchange) - - ! Add the ejecta back to the DEM - do ypi = 1,user%gridsize - do xpi = 1,user%gridsize - surf(xpi,ypi)%dem = surf(xpi,ypi)%dem + cumulative_elchange(xpi,ypi) - surf(xpi,ypi)%ejcov = max(surf(xpi,ypi)%ejcov + cumulative_elchange(xpi,ypi), 0.0_DP) - end do - end do + ! Create box for soften calculation (will be no bigger than the grid itself) + if (2 * inc + 1 < user%gridsize) then - else + ! extra soften calculation + if (user%dosoftening) then + cel = 0.0_DP call util_diffusion_solver(user,surf,2 * inc + 1,indarray,kdiff,cel,maxhits) do j = -inc,inc do i = -inc,inc @@ -402,22 +363,75 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ surf(xpi,ypi)%ejcov = max(surf(xpi,ypi)%ejcov + cel(i,j),0.0_DP) end do end do + end if - call ejecta_soften(user,surf,2 * inc + 1,indarray,cumulative_elchange) + call ejecta_soften(user,surf,2 * inc + 1,indarray,cumulative_elchange) - ! Add the ejecta back to the DEM - do j = -inc,inc - do i = -inc,inc - xpi = indarray(1,i,j) - ypi = indarray(2,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) + ! Add the ejecta back to the DEM + do j = -inc,inc + do i = -inc,inc + xpi = indarray(1,i,j) + ypi = indarray(2,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 + end do + + + else ! Ejecta wraps around the grid. + ! We will therefore send in the whole grid with the total ejecta thickness added to each pixel + allocate(big_cumulative_elchange(0:user%gridsize+1,0:user%gridsize+1)) + allocate(big_cel(0:user%gridsize+1,0:user%gridsize+1)) + allocate(big_indarray(2,0:user%gridsize+1,0:user%gridsize+1)) + allocate(big_kdiff(0:user%gridsize+1,0:user%gridsize+1)) + + do bigj = 0,user%gridsize + 1 + do bigi = 0,user%gridsize + 1 + xpi = bigi + ypi = bigj + call util_periodic(xpi,ypi,user%gridsize) + big_indarray(1,bigi,bigj) = xpi + big_indarray(2,bigi,bigj) = ypi + big_cumulative_elchange(bigi,bigj) = 0.0_DP + end do + end do + + big_kdiff = 0.0_DP + + do j = -inc,inc + do i = -inc,inc + xpi = indarray(1,i,j) + ypi = indarray(2,i,j) + big_cumulative_elchange(xpi,ypi) = big_cumulative_elchange(xpi,ypi) + cumulative_elchange(i,j) + big_kdiff(xpi,ypi) = big_kdiff(xpi,ypi) + kdiff(i,j) + end do + end do + + if (user%dosoftening) then + big_cel = 0.0_DP + call util_diffusion_solver(user,surf,user%gridsize + 2,big_indarray,big_kdiff,big_cel,maxhits) + + do j = 1,user%gridsize + do i = 1,user%gridsize + surf(i,j)%dem = surf(i,j)%dem + big_cel(i,j) + surf(i,j)%ejcov = max(surf(i,j)%ejcov + big_cel(i,j),0.0_DP) end do end do end if + + call ejecta_soften(user,surf,user%gridsize + 2,big_indarray,big_cumulative_elchange) + + do i = 1,user%gridsize + do j = 1,user%gridsize + surf(i,j)%dem = surf(i,j)%dem + big_cumulative_elchange(i,j) + surf(i,j)%ejcov = max(surf(i,j)%ejcov + big_cumulative_elchange(i,j),0.0_DP) + end do + end do + + deallocate(big_cumulative_elchange,big_indarray,big_kdiff,big_cel) end if - deallocate(indarray,kdiff,cel) + deallocate(indarray,diffdistribution,ejdistribution,kdiff,cel) return end subroutine ejecta_emplace diff --git a/src/ejecta/ejecta_ray_pattern.f90 b/src/ejecta/ejecta_ray_pattern.f90 index 663a216d..352cea4a 100644 --- a/src/ejecta/ejecta_ray_pattern.f90 +++ b/src/ejecta/ejecta_ray_pattern.f90 @@ -26,11 +26,12 @@ ! * domain -- Simulation domain variable container ! ! Output +! * ejdist -- logical array containing the pixels that contain ejecta ! !*** !********************************************************************************************************************************** -subroutine ejecta_ray_pattern(user,crater,i,j,diffi,eji) +subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,diffdistribution,ejdistribution) use module_globals use module_util use module_io @@ -41,22 +42,24 @@ subroutine ejecta_ray_pattern(user,crater,i,j,diffi,eji) ! Arguments type(usertype),intent(in) :: user + type(surftype),dimension(:,:),intent(in) :: surf type(cratertype),intent(inout) :: crater - integer(I4B),intent(in) :: i,j - real(DP),intent(out) :: diffi - real(DP),intent(out) :: eji + integer(I4B),intent(in) :: inc,xi,xf,yi,yf + real(DP),dimension(xi:xf,yi:yf),intent(out) :: diffdistribution + real(DP),dimension(xi:xf,yi:yf),intent(out) :: ejdistribution ! Internal variables - integer(I4B) :: nrays,k,n,nef,incsq,iradsq,xpi,ypi,ejpxsq + integer(I4B) :: nrays,i,j,k,n,nef,incsq,iradsq,xpi,ypi,ejpxsq real(DP) :: frac,mef,lrad,lradsq,xp,yp,binres,areafrac,xbar,ybar real(DP) :: rn real(DP) :: theta, lradp, maxdistance real(DP), parameter :: n1 = 4.0_DP real(DP) :: n2, mag + real(DP),dimension(xi:xf,yi:yf) :: isray real(DP),dimension(:),allocatable :: numinray,totnum real(DP),dimension(:),allocatable :: mefarray real(DP) :: ans - logical :: bigej + real(DP) :: C1,C2,p ! Fits to fe vs fd equation for the new ray pattern real(DP) :: rmin,rmax,r @@ -71,12 +74,11 @@ subroutine ejecta_ray_pattern(user,crater,i,j,diffi,eji) logical :: ej real(DP),dimension(Nraymax) :: thetari - eji = 0.0_DP - diffi = 0.0_DP + !TEMPORARY if (user%dorays) then - do k = 1,Nraymax - thetari(k) = 2 * pi * k / Nraymax + do i = 1,Nraymax + thetari(i) = 2 * pi * i / Nraymax end do call shuffle(thetari) ! randomize the ray pattern @@ -84,38 +86,65 @@ subroutine ejecta_ray_pattern(user,crater,i,j,diffi,eji) rmax = user%ejecta_truncation rmin = crater%continuous / crater%frad crater%fe = 10.0_DP ! Estimate the equivalent degradation radius - xpi = crater%xlpx + i - ypi = crater%ylpx + j - - ! Find distance from crater center to current pixel center in real space - xp = xpi * user%pix - yp = ypi * user%pix - - xbar = xp - crater%xl - ybar = yp - crater%yl - - areafrac = util_area_intersection(user%ejecta_truncation * crater%frad,xbar,ybar,user%pix) - r = sqrt(xbar**2 + ybar**2) / crater%frad - theta = mod(atan2(ybar,xbar) + pi + rn * 2 * pi,2 * pi) - diffi = areafrac * ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,.false.) - eji = areafrac * ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,.true.) + !ejdistribution = 0.0_DP + !diffdistribution = 0.0_DP + !!$OMP PARALLEL DO DEFAULT(PRIVATE) & + !!$OMP SHARED(user,crater) & + !!$OMP SHARED(xi,xf,yi,yf,rn,diffdistribution,ejdistribution,thetari,rmin) + do j = yi,yf + do i = xi,xf + xpi = crater%xlpx + i + ypi = crater%ylpx + j + + ! Find distance from crater center to current pixel center in real space + xp = xpi * user%pix + yp = ypi * user%pix + + xbar = xp - crater%xl + ybar = yp - crater%yl + + areafrac = util_area_intersection(user%ejecta_truncation * crater%frad,xbar,ybar,user%pix) + r = sqrt(xbar**2 + ybar**2) / crater%frad + theta = mod(atan2(ybar,xbar) + pi + rn * 2 * pi,2 * pi) + diffdistribution(i,j) = areafrac * ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,.false.) + ejdistribution(i,j) = areafrac * ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,.true.) + end do + end do + !!$OMP END PARALLEL DO + else + !Do simple circular region + incsq = inc**2 + ejdistribution = 0.0_DP + diffdistribution = 0.0_DP + !$OMP PARALLEL DO DEFAULT(PRIVATE) & + !$OMP SHARED(user,crater) & + !$OMP SHARED(inc,incsq,xi,xf,yi,yf,diffdistribution,ejdistribution) + do j = yi,yf + do i = xi,xf + iradsq = i*i + j*j + + if (iradsq < incsq) then + + xpi = crater%xlpx + i + ypi = crater%ylpx + j + + ! Find distance from crater center to current pixel center in real space + xp = xpi * user%pix + yp = ypi * user%pix + + xbar = xp - crater%xl + ybar = yp - crater%yl + areafrac = util_area_intersection(user%ejecta_truncation * crater%frad,xbar,ybar,user%pix) ! uniform circular + diffdistribution(i,j) = areafrac + ejdistribution(i,j) = areafrac + end if + end do + end do + !$OMP END PARALLEL DO - xpi = crater%xlpx + i - ypi = crater%ylpx + j - - ! Find distance from crater center to current pixel center in real space - xp = xpi * user%pix - yp = ypi * user%pix - - xbar = xp - crater%xl - ybar = yp - crater%yl - areafrac = util_area_intersection(user%ejecta_truncation * crater%frad,xbar,ybar,user%pix) ! uniform circular - diffi = areafrac - eji = areafrac end if - return contains diff --git a/src/ejecta/module_ejecta.f90 b/src/ejecta/module_ejecta.f90 index 9e75eb4d..06f8d9e9 100644 --- a/src/ejecta/module_ejecta.f90 +++ b/src/ejecta/module_ejecta.f90 @@ -43,14 +43,15 @@ end subroutine ejecta_emplace end interface interface - subroutine ejecta_ray_pattern(user,crater,i,j,diffi,eji) + subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,diffdistribution,ejdistribution) use module_globals implicit none type(usertype),intent(in) :: user + type(surftype),dimension(:,:),intent(in) :: surf type(cratertype),intent(inout) :: crater - integer(I4B),intent(in) :: i,j - real(DP),intent(out) :: diffi - real(DP),intent(out) :: eji + integer(I4B),intent(in) :: inc,xi,xf,yi,yf + real(DP),dimension(xi:xf,yi:yf),intent(out) :: diffdistribution + real(DP),dimension(xi:xf,yi:yf),intent(out) :: ejdistribution end subroutine ejecta_ray_pattern end interface diff --git a/src/regolith/module_regolith.f90 b/src/regolith/module_regolith.f90 index 0c77a42f..2ab24cd9 100644 --- a/src/regolith/module_regolith.f90 +++ b/src/regolith/module_regolith.f90 @@ -102,7 +102,7 @@ subroutine regolith_traverse_streamtube(user,surfi,deltar,ri,rip1,eradi,erado,ne type(surftype),intent(inout) :: surfi real(DP),intent(in) :: deltar,ri,rip1,eradi,erado type(regodatatype),intent(inout) :: newlayer - real(DP),intent(inout) :: meltinejecta,totvol + real(DP),intent(out) :: meltinejecta,totvol real(DP),intent(out) :: vmare,totseb real(SP),dimension(:),intent(inout) :: age_collector real(DP),intent(in) :: xmints diff --git a/src/regolith/regolith_subpixel_streamtube.f90 b/src/regolith/regolith_subpixel_streamtube.f90 index 694b23f5..acf2b08d 100644 --- a/src/regolith/regolith_subpixel_streamtube.f90 +++ b/src/regolith/regolith_subpixel_streamtube.f90 @@ -61,7 +61,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer type(surftype),intent(inout) :: surfi real(DP),intent(in) :: deltar,ri,rip1,eradi type(regodatatype),intent(inout) :: newlayer - real(DP),intent(inout) :: meltinejecta, totvol + real(DP),intent(out) :: meltinejecta, totvol real(DP),intent(out) :: vmare,totseb real(SP),dimension(:),intent(inout) :: age_collector real(DP),intent(in) :: xmints