attempt to move function pattern() to its own f90 file called ejecta_…

…ray_pattern_function

src/Makefile.am

@@ -65,6 +65,7 @@ io/io_write_pindex_map.f90\
 io/io_write_age_depth.f90\
 ejecta/ejecta_emplace.f90\
 ejecta/ejecta_ray_pattern.f90\
+ejecta/ejecta_ray_pattern_function.f90\
 ejecta/ejecta_blanket.f90\
 ejecta/ejecta_blanket_func.f90\
 ejecta/ejecta_thickness.f90\

src/ejecta/ejecta_emplace.f90

@@ -189,6 +189,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_r
    !!$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
@@ -276,11 +277,13 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_r
          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,age,age_resolution,volm)
             vmelt = vmelt + volm
+            !write(74,*) erad, surf(xpi,ypi)%regolayer%regodata%meltfrac
          end if
 
 
       end do
    end do
+   !close(74)
    !!$OMP END PARALLEL DO
    if(user%doregotrack .and. user%testflag) then
       write(*,*) 'Ejected Melt: ', vmelt

src/ejecta/ejecta_ray_pattern.f90

@@ -58,6 +58,7 @@ subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,diffdistribution,
    real(DP),dimension(xi:xf,yi:yf) :: isray
    real(DP),dimension(:),allocatable :: numinray,totnum
    real(DP),dimension(:),allocatable :: mefarray
+   real(DP) :: ans
 
 
    real(DP) :: C1,C2,p ! Fits to fe vs fd equation for the new ray pattern
@@ -105,8 +106,8 @@ subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,diffdistribution,
             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 * pattern(theta,r,rmin,rmax,thetari,.false.) 
-            ejdistribution(i,j) = areafrac * pattern(theta,r,rmin,rmax,thetari,.true.) 
+            diffdistribution(i,j) = areafrac * ejecta_ray_pattern_function(theta,r,rmin,rmax,thetari,.false.) 
+            ejdistribution(i,j) = areafrac * ejecta_ray_pattern_function(theta,r,rmin,rmax,thetari,.true.) 
          end do
       end do
       !!$OMP END PARALLEL DO
@@ -147,73 +148,76 @@ subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,diffdistribution,
    return  
    contains
 
-   pure function ray(theta,thetar,r,n,w) result(ans)
-   implicit none
-   real(DP) :: ans
-   real(DP),intent(in) :: theta,thetar,r,w
-   integer(I4B),intent(in) :: n
-   real(DP) :: thetap,thetapp,a,b,c,dtheta
-
-   c = w / r
-   b = thetar 
-   dtheta = min(2*pi - abs(theta - b),abs(theta - b))
-   a = sqrt(2 * pi) / (n * c * erf(pi / (2 *sqrt(2._DP) * c)))
-   ans = a * exp(-dtheta**2 / (2 * c**2))
-
-   return
-   end function ray
-
-   pure function pattern(theta,r,rmin,rmax,thetari,ej) result(ans)
-   implicit none
-   real(DP) :: ans
-   real(DP),intent(in) :: r,rmin,rmax,theta
-   real(DP),dimension(:),intent(in) :: thetari
-   logical,intent(in) :: ej
-   real(DP) :: a,c
-   real(DP) :: thetar,rw,rw0,rw1
-   real(DP) :: f,rtrans,length,rpeak,minray,FF
-   integer(I4B) :: n,i
-
-
-   minray = rmin * 3
-
-   if (r > rmax) then
-      ans = 0._DP
-   else if (r < 1.0_DP) then
-      if (ej) then
-         ans = 1.0_DP
-      else
-         ans = 0.0_DP
-      end if
-   else
-      rw0 = rmin * pi / Nraymax / 2
-      rw1 = 2 * pi / Nraymax
-      rw = rw0 * (1._DP - (1.0_DP - rw1 / rw0) * exp(1._DP - (r / rmin)**2))
-      n = max(min(floor((Nraymax**rayp - (Nraymax**rayp - 1) * log(r/minray) / log(rray/minray))**(1._DP/rayp)),Nraymax),1) ! Exponential decay of ray number with distance
-      ans = 0._DP
-      rtrans = r - 1.0_DP
-      c = rw / r
-      a = sqrt(2 * pi) / (n * c * erf(pi / (2 *sqrt(2._DP) * c)))
-      do i = 1,Nraymax
-         length = minray * exp(log(rray/minray) * ((Nraymax - i + 1)**rayp - 1_DP) / ((Nraymax**rayp - 1)))
-         rpeak = (length - 1_DP) * 0.5_DP
-         if (ej) then
-            FF = 1.0_DP
-            if (r > length) then
-               f = 0.0_DP
-            else
-               f = a 
-            end if
-         else
-            FF = rayfmult * (20 / rmax)**(0.5_DP) * 0.25_DP 
-            f = FF * fpeak * (rtrans / rpeak)**rayq * exp(1._DP / rayq * (1.0_DP - (rtrans / rpeak)**rayq)) 
-         end if
-         ans = ans + ray(theta,thetari(i),r,n,rw) * f / a 
-      end do
-   end if 
+   ! pure function ray(theta,thetar,r,n,w) result(ans)
+   ! implicit none
+   ! real(DP) :: ans
+   ! real(DP),intent(in) :: theta,thetar,r,w
+   ! integer(I4B),intent(in) :: n
+   ! real(DP) :: thetap,thetapp,a,b,c,dtheta
+
+   ! c = w / r
+   ! b = thetar 
+   ! dtheta = min(2*pi - abs(theta - b),abs(theta - b))
+   ! a = sqrt(2 * pi) / (n * c * erf(pi / (2 *sqrt(2._DP) * c)))
+   ! ans = a * exp(-dtheta**2 / (2 * c**2))
+
+   ! return
+   ! end function ray
+
+   ! !pure function pattern(theta,r,rmin,rmax,thetari,ej) result(ans)
+   ! !function pattern(theta,r,rmin,rmax,thetari,ej) result(ans)
+   ! subroutine pattern(theta,r,rmin,rmax,thetari,ej,ans)
+   ! implicit none
+   ! real(DP) :: ans
+   ! real(DP),intent(in) :: r,rmin,rmax,theta
+   ! real(DP),dimension(:),intent(in) :: thetari
+   ! logical,intent(in) :: ej
+   ! real(DP) :: a,c
+   ! real(DP) :: thetar,rw,rw0,rw1
+   ! real(DP) :: f,rtrans,length,rpeak,minray,FF
+   ! integer(I4B) :: n,i
+
+
+   ! minray = rmin * 3
+
+   ! if (r > rmax) then
+   !    ans = 0._DP
+   ! else if (r < 1.0_DP) then
+   !    if (ej) then
+   !       ans = 1.0_DP
+   !    else
+   !       ans = 0.0_DP
+   !    end if
+   ! else
+   !    rw0 = rmin * pi / Nraymax / 2
+   !    rw1 = 2 * pi / Nraymax
+   !    rw = rw0 * (1._DP - (1.0_DP - rw1 / rw0) * exp(1._DP - (r / rmin)**2))
+   !    n = max(min(floor((Nraymax**rayp - (Nraymax**rayp - 1) * log(r/minray) / log(rray/minray))**(1._DP/rayp)),Nraymax),1) ! Exponential decay of ray number with distance
+   !    ans = 0._DP
+   !    rtrans = r - 1.0_DP
+   !    c = rw / r
+   !    a = sqrt(2 * pi) / (n * c * erf(pi / (2 *sqrt(2._DP) * c)))
+   !    do i = 1,Nraymax
+   !       length = minray * exp(log(rray/minray) * ((Nraymax - i + 1)**rayp - 1_DP) / ((Nraymax**rayp - 1)))
+   !       rpeak = (length - 1_DP) * 0.5_DP
+   !       if (ej) then
+   !          FF = 1.0_DP
+   !          if (r > length) then
+   !             f = 0.0_DP
+   !          else
+   !             f = a 
+   !          end if
+   !       else
+   !          FF = rayfmult * (20 / rmax)**(0.5_DP) * 0.25_DP 
+   !          f = FF * fpeak * (rtrans / rpeak)**rayq * exp(1._DP / rayq * (1.0_DP - (rtrans / rpeak)**rayq)) 
+   !       end if
+   !       ans = ans + ray(theta,thetari(i),r,n,rw) * f / a 
+   !    end do
+   ! end if 
 
 
-   end function pattern
+   ! !end function pattern
+   ! end subroutine pattern
 
    subroutine shuffle(a)
    real(DP), intent(inout) :: a(:)

src/ejecta/ejecta_ray_pattern_function.f90

@@ -0,0 +1,99 @@
+!****f* regolith/regolith_streamtube_volume_func
+! Name
+!   ejecta_ray_pattern_func -- Calculate ejecta ray pattern
+! SYNOPSIS
+!   This uses
+!   * module_globals
+!   * module_ejecta
+!
+!   ans = regolith_streamtube_volume_func()
+!
+! DESCRIPTION
+!   
+!   This function was separated into an independent function for debugging purposes.
+!  
+! ARGUMENTS
+!   Input
+!   * theta      -- 
+!   * r     -- 
+!   * rmin    -- 
+!   * rmax        -- 
+!   * thetari      -- 
+!   * ej     -- 
+!
+!   Output
+!   * ans  -- THe result of this function, used in ejecta_ray_pattern.f90
+! 
+!***
+
+!**********************************************************************************************************
+function pattern(theta,r,rmin,rmax,thetari,ej) result(ans)
+    use module_globals
+    use module_ejecta, EXCEPT_THIS_ONE => ejecta_ray_pattern_function
+    implicit none
+    real(DP) :: ans
+    real(DP),intent(in) :: r,rmin,rmax,theta
+    real(DP),dimension(:),intent(in) :: thetari
+    logical,intent(in) :: ej
+    real(DP) :: a,c
+    real(DP) :: thetar,rw,rw0,rw1
+    real(DP) :: f,rtrans,length,rpeak,minray,FF
+    integer(I4B) :: n,i
+
+
+    minray = rmin * 3
+
+    if (r > rmax) then
+       ans = 0._DP
+    else if (r < 1.0_DP) then
+       if (ej) then
+          ans = 1.0_DP
+       else
+          ans = 0.0_DP
+       end if
+    else
+       rw0 = rmin * pi / Nraymax / 2
+       rw1 = 2 * pi / Nraymax
+       rw = rw0 * (1._DP - (1.0_DP - rw1 / rw0) * exp(1._DP - (r / rmin)**2))
+       n = max(min(floor((Nraymax**rayp - (Nraymax**rayp - 1) * log(r/minray) / log(rray/minray))**(1._DP/rayp)),Nraymax),1) ! Exponential decay of ray number with distance
+       ans = 0._DP
+       rtrans = r - 1.0_DP
+       c = rw / r
+       a = sqrt(2 * pi) / (n * c * erf(pi / (2 *sqrt(2._DP) * c)))
+       do i = 1,Nraymax
+          length = minray * exp(log(rray/minray) * ((Nraymax - i + 1)**rayp - 1_DP) / ((Nraymax**rayp - 1)))
+          rpeak = (length - 1_DP) * 0.5_DP
+          if (ej) then
+             FF = 1.0_DP
+             if (r > length) then
+                f = 0.0_DP
+             else
+                f = a 
+             end if
+          else
+             FF = rayfmult * (20 / rmax)**(0.5_DP) * 0.25_DP 
+             f = FF * fpeak * (rtrans / rpeak)**rayq * exp(1._DP / rayq * (1.0_DP - (rtrans / rpeak)**rayq)) 
+          end if
+          ans = ans + ray(theta,thetari(i),r,n,rw) * f / a 
+       end do
+    end if 
+    !return
+    contains
+
+    pure function ray(theta,thetar,r,n,w) result(ans)
+    implicit none
+    real(DP) :: ans
+    real(DP),intent(in) :: theta,thetar,r,w
+    integer(I4B),intent(in) :: n
+    real(DP) :: thetap,thetapp,a,b,c,dtheta
+
+    c = w / r
+    b = thetar 
+    dtheta = min(2*pi - abs(theta - b),abs(theta - b))
+    a = sqrt(2 * pi) / (n * c * erf(pi / (2 *sqrt(2._DP) * c)))
+    ans = a * exp(-dtheta**2 / (2 * c**2))
+
+    return
+    end function ray
+
+end function pattern

src/ejecta/module_ejecta.f90

@@ -54,6 +54,17 @@ subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,diffdistribution,
       end subroutine ejecta_ray_pattern
    end interface
 
+   interface
+      function ejecta_ray_pattern_function(theta,r,rmin,rmax,thetari,ej) result(ans)
+      use module_globals
+      implicit none
+      real(DP) :: ans
+      real(DP),intent(in) :: r,rmin,rmax,theta
+      real(DP),dimension(:),intent(in) :: thetari
+      logical,intent(in) :: ej
+      end function ejecta_ray_pattern_function
+   end interface
+
    interface
       subroutine ejecta_rootfind(user,crater,domain,erad,lrad,vejsq,ejang,firstrun)
       use module_globals

src/regolith/regolith_melt_zone.f90

@@ -74,6 +74,8 @@ subroutine regolith_melt_zone(user,crater,dimp,vimp,rmelt,depthb,volm)
    depthb = 0.5 * dimp
    vtc = PI/3.0 * (crater%rad)**3
    volm = 2.9 * vtc * Em**(-0.85) * (dimp)**0.66 * (user%gaccel) **0.66 * vimp**(0.37) 
+   !increase or decrease volm by and arbitrary amount to see what happens-- ***NOT PHYSICAL!!!***
+   !volm = volm / 3.0
    ! Calculate radius of melt zone (shifted melt zone model Pierazzo et al. 1997)
    ! rmelt = (rvapor**3 + 1.5/PI*volm)**(1.0/3.0) (hemisphere model)
    ! R_m ^3 + 1.5 * R_imp * R_m ^2 - 2.5 * R_imp ^3 - 3/(2 * PI) * volm = 0