Fixed bugs in the ejecta restructuring and also fixed some inconsiste…

…ncy in intent statements between interface modules and implementations.

src/Makefile.am

@@ -11,7 +11,7 @@ AM_FCFLAGS =  $(IPRODUCTION)
 #ifort debug flags
 
 #gfortran optimized flags
-#AM_FCFLAGS =  -O3 -fopenmp -ffree-form  -g -fbounds-check -fbacktrace
+AM_FCFLAGS =  -O3 -fopenmp -ffree-form  -g -fbounds-check -fbacktrace
 #gfortran debug flags
 #AM_FCFLAGS =  -O0 -g -fopenmp -fbounds-check -Wall -Warray-bounds -Warray-temporaries -Wimplicit-interface -ffree-form -fsanitize-address-use-after-scope -fstack-check  -fsanitize=bounds-strict -fsanitize=undefined -fsanitize=signed-integer-overflow -fsanitize=object-size -fstack-protector-all
 

src/crater/crater_subpixel_diffusion.f90

@@ -174,8 +174,6 @@ 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))
             ! Loop over affected matrix area
             !!$OMP PARALLEL DO DEFAULT(SHARED) IF(inc > INCPAR) &
             !!$OMP FIRSTPRIVATE(jmin,jmax,imin,imax) &
@@ -195,7 +193,6 @@ subroutine crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiffin)
                end do
             end do
             !!$OMP END PARALLEL DO
-            deallocate(diffdistribution,ejdistribution)
          end do
       end if
 

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(in) :: user
+      type(usertype),intent(inout) :: user
       type(surftype),dimension(:,:),intent(inout)  :: surf
       type(cratertype),intent(inout)               :: crater
       type(domaintype),intent(inout)               :: domain

src/ejecta/ejecta_emplace.f90

@@ -93,38 +93,32 @@ 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(out) :: cumulative_elchange
+   real(DP),dimension(:,:),allocatable,intent(inout) :: 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,k,n,inc,incsq,iradsq,idistorted,jdistorted
+   integer(I4B) :: xpi,ypi,i,j,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 :: big_cumulative_elchange,kdiff,big_kdiff,cel,big_cel
-   integer(I4B),dimension(:,:,:),allocatable :: indarray,big_indarray
-   integer(I4B) :: bigi,bigj,maxhits,nin,nnot,dradsq
+   real(DP),dimension(:,:),allocatable :: kdiff,cel
+   integer(I4B),dimension(:,:,:),allocatable :: indarray
+   integer(I4B) :: 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, eradc
+   real(DP) :: rm, dm, melt
 
    ! Ejecta pattern distortion parameters
    real(DP) :: distance,erad,craterslope,landslope,baseline,lrange,frac,ejheight,ebh0,maxdistance
    real(DP) :: maxslp
    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)
@@ -368,7 +362,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ
       cel = 0.0_DP
 
       if (bigej) then
-         call util_diffusion_solver(user,user%gridsize + 2,indarray,kdiff,cel,maxhits)
+         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)
@@ -382,7 +376,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ
          indarray(1:2,:,0) = user%gridsize
          indarray(1:2,:,user%gridsize+1) = 1
 
-         call ejecta_soften(user,user%gridsize + 2,indarray,cumulative_elchange)
+         call ejecta_soften(user,surf,user%gridsize + 2,indarray,cumulative_elchange)
 
          ! Add the ejecta back to the DEM
          do ypi = 1,user%gridsize

src/ejecta/ejecta_ray_pattern.f90

@@ -53,7 +53,6 @@ subroutine ejecta_ray_pattern(user,crater,i,j,diffi,eji)
    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
@@ -76,8 +75,8 @@ subroutine ejecta_ray_pattern(user,crater,i,j,diffi,eji)
    diffi = 0.0_DP
 
    if (user%dorays) then
-      do i = 1,Nraymax
-         thetari(i) = 2 * pi * i / Nraymax
+      do k = 1,Nraymax
+         thetari(k) = 2 * pi * k / Nraymax
       end do
       call shuffle(thetari) ! randomize the ray pattern
 
@@ -102,26 +101,19 @@ subroutine ejecta_ray_pattern(user,crater,i,j,diffi,eji)
       eji = areafrac * ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,.true.) 
 
    else
-      !Do simple circular region
-      incsq = inc**2
-      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
+      xpi = crater%xlpx + i
+      ypi = crater%ylpx + j
 
-         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
+      ! 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  

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(out)            :: meltinejecta,totvol
+      real(DP),intent(inout)            :: meltinejecta,totvol
       real(DP),intent(out) :: vmare,totseb
       real(SP),dimension(:),intent(inout) :: age_collector
       real(DP),intent(in)             :: xmints

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(out)              :: meltinejecta, totvol
+   real(DP),intent(inout)              :: meltinejecta, totvol
    real(DP),intent(out)                :: vmare,totseb
    real(SP),dimension(:),intent(inout) :: age_collector
    real(DP),intent(in)                 :: xmints