diff --git a/src/Makefile.am b/src/Makefile.am index d6766570..8d53b108 100755 --- a/src/Makefile.am +++ b/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\ diff --git a/src/ejecta/ejecta_emplace.f90 b/src/ejecta/ejecta_emplace.f90 index a49ddacd..1539bba7 100644 --- a/src/ejecta/ejecta_emplace.f90 +++ b/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 diff --git a/src/ejecta/ejecta_ray_pattern.f90 b/src/ejecta/ejecta_ray_pattern.f90 index 32982add..bc280982 100644 --- a/src/ejecta/ejecta_ray_pattern.f90 +++ b/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(:) diff --git a/src/ejecta/ejecta_ray_pattern_function.f90 b/src/ejecta/ejecta_ray_pattern_function.f90 new file mode 100644 index 00000000..2f0e9d4b --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/ejecta/module_ejecta.f90 b/src/ejecta/module_ejecta.f90 index de3e5c3a..578af4fb 100644 --- a/src/ejecta/module_ejecta.f90 +++ b/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 diff --git a/src/regolith/regolith_melt_zone.f90 b/src/regolith/regolith_melt_zone.f90 index 58ea3566..a62901a5 100644 --- a/src/regolith/regolith_melt_zone.f90 +++ b/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