fixed a bug with the melt index arrays

MintonGroup · Oct 26, 2023 · 1a1f03c · 1a1f03c
1 parent 924c89c
commit 1a1f03c
Show file tree

Hide file tree

Showing 8 changed files with 244 additions and 241 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) :: 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
@@ -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
@@ -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