Restructured the ejecta blanket functions to prevent temporary arrays…

… from being generated that are larger than the domain.
MintonGroup · Sep 1, 2023 · 26d1b24 · 26d1b24
1 parent b6c7c5f
commit 26d1b24
Show file tree

Hide file tree

Showing 5 changed files with 162 additions and 212 deletions.
diff --git 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),dimension(:,:),allocatable :: diffdistribution,ejdistribution
+   real(DP) :: eji, diffi
    integer(I4B),dimension(:,:),allocatable :: ejisray
    real(DP) :: xbar,ybar,dD,xp,yp,areafrac,krad,supersize,lrad
    integer(I8B),dimension(user%gridsize,user%gridsize) :: Ngrid
@@ -176,7 +176,6 @@ subroutine crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiffin)
 
             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) &
@@ -185,12 +184,13 @@ subroutine crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiffin)
                do i = imin,imax
                   xpi = crater%xlpx + i
                   ypi = crater%ylpx + j
-                  kdiff(xpi,ypi) = kdiff(xpi,ypi) + dKdN * diffdistribution(i,j)
+                  call ejecta_ray_pattern(user,crater,i,j,diffi,eji)
+                  kdiff(xpi,ypi) = kdiff(xpi,ypi) + dKdN * diffi
                   !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) &
+                  surf(xpi,ypi)%ejcov = surf(xpi,ypi)%ejcov + eji &
                   * 0.14_DP * crater%frad**(0.74_DP) * (lrad / crater%frad)**(-3)
                end do
             end do

diff --git 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),dimension(:,:),allocatable :: ejdistribution, diffdistribution
-   integer(I4B),dimension(:,:),allocatable :: ejisray
+   real(DP) :: diffi, eji
+   logical :: ejisray
 
    ! Melt or glassy ray test
    real(DP) :: xp, yp, lradsq, lrad, erad
@@ -144,40 +144,26 @@ 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
-                  if (ejisray(xpi,ypi) == 0) cycle  
-                  call regolith_superdomain(user,crater,domain,surf(i,j)%regolayer,ejdistribution(xpi,ypi),&
+                  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,&
                        i,j,rm,depthb)
                end do
             end do
 
          end if
 
-         deallocate(ejdistribution,diffdistribution,ejisray)
 
        end do