diff --git a/src/Makefile.am b/src/Makefile.am index 2e8c1902..efe3c311 100755 --- a/src/Makefile.am +++ b/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 diff --git a/src/crater/crater_subpixel_diffusion.f90 b/src/crater/crater_subpixel_diffusion.f90 index f374f251..3da7b696 100644 --- a/src/crater/crater_subpixel_diffusion.f90 +++ b/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 diff --git a/src/crater/module_crater.f90 b/src/crater/module_crater.f90 index c0dd5dfe..d5108399 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(in) :: user + type(usertype),intent(inout) :: user type(surftype),dimension(:,:),intent(inout) :: surf type(cratertype),intent(inout) :: crater type(domaintype),intent(inout) :: domain diff --git a/src/ejecta/ejecta_emplace.f90 b/src/ejecta/ejecta_emplace.f90 index d2497bd2..fdb05563 100644 --- a/src/ejecta/ejecta_emplace.f90 +++ b/src/ejecta/ejecta_emplace.f90 @@ -93,27 +93,25 @@ 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 @@ -121,10 +119,6 @@ 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) @@ -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 diff --git a/src/ejecta/ejecta_ray_pattern.f90 b/src/ejecta/ejecta_ray_pattern.f90 index 99d9a515..663a216d 100644 --- a/src/ejecta/ejecta_ray_pattern.f90 +++ b/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 diff --git a/src/regolith/module_regolith.f90 b/src/regolith/module_regolith.f90 index 2ab24cd9..0c77a42f 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(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 diff --git a/src/regolith/regolith_subpixel_streamtube.f90 b/src/regolith/regolith_subpixel_streamtube.f90 index acf2b08d..694b23f5 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(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