added melt sheet emplacement. Compiles but not tested (want to make mixing optional for test craters first)

src/crater/crater_emplace.f90

@@ -51,7 +51,7 @@
 !  Notes       :  
 !
 !**********************************************************************************************************************************
-subroutine crater_emplace(user,surf,crater,domain,deltaMtot,incval)
+subroutine crater_emplace(user,surf,crater,domain,deltaMtot,incval,nmeltsheet)
    use module_globals
    use module_util
    use module_porosity   
@@ -64,7 +64,7 @@ subroutine crater_emplace(user,surf,crater,domain,deltaMtot,incval)
    type(cratertype),intent(inout) :: crater
    type(domaintype),intent(inout) :: domain
    real(DP),intent(out) :: deltaMtot
-   integer(I4B),intent(out) :: incval
+   integer(I4B),intent(out) :: incval,nmeltsheet
 
    ! Internal variables
    real(DP) :: lradsq,newelev, x_relative, y_relative 
@@ -83,6 +83,7 @@ subroutine crater_emplace(user,surf,crater,domain,deltaMtot,incval)
    deltaMtot = 0.0_DP !ejbmass
    incsq = inc**2
    incval = inc
+   nmeltsheet = 0
 
    ! This loop may not be parallelizable because of the linked list operation inside crater_form_interior
    do j=-inc,inc  ! Do the loop in pixel space
@@ -108,6 +109,7 @@ subroutine crater_emplace(user,surf,crater,domain,deltaMtot,incval)
          if (lradsq > crater%frad**2) cycle
          call crater_form_interior(user,surf(xpi,ypi),crater,x_relative,y_relative,newelev,deltaMi)
          deltaMtot = deltaMtot + deltaMi
+         nmeltsheet = nmeltsheet + 1
 
          ! do porosity computation if (user%doporosity)
          ! It is still important to consider the physical meaning of frad and rad. 

src/crater/crater_populate.f90

@@ -78,7 +78,8 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
    integer(I4B)            :: oldpbarpos
    real(DP),dimension(:,:),allocatable   :: ejecta_dem
    real(DP)                :: hmax, hmin
-   integer(I4B)            :: nmixingtimes, incval
+   integer(I4B)            :: nmixingtimes, incval, nmeltsheet
+   real(DP)                :: vmeltsheet
 
    ! ejecta blanket array
    type(ejbtype),dimension(EJBTABSIZE) :: ejb       ! Ejecta blanket lookup table
@@ -249,7 +250,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
          call crater_averages(user,surf,crater)
 
          ! Place crater onto the surface
-         call crater_emplace(user,surf,crater,domain,ejbmass,incval)
+         call crater_emplace(user,surf,crater,domain,ejbmass,incval,nmeltsheet)
          if (abs(ejbmass) < 2*tiny(1.0_DP)) cycle
 
          call ejecta_distance_estimate(user,crater,domain,crater%ejdis) ! Fast but imprecise estimate of the total ejecta distance
@@ -267,12 +268,12 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
                call ejecta_table_define(user,crater,domain,ejb,ejtble)
                !call ejecta_interpolate(crater,domain,crater%frad,ejb(1:ejtble),ejtble,crater%ejrim)
             end if
-            call ejecta_emplace(user,surf,crater,domain,ejb(1:ejtble),ejtble,ejbmass,age,age_resolution,ejecta_dem)
+            call ejecta_emplace(user,surf,crater,domain,ejb(1:ejtble),ejtble,ejbmass,age,age_resolution,ejecta_dem,vmeltsheet)
          else
             ejtble = 0
          end if
 
-         if (user%doregotrack) call regolith_interior(user,surf,crater,incval)
+         if (user%doregotrack) call regolith_interior(user,surf,crater,domain,incval,nmeltsheet,vmeltsheet)
 
          if (user%dorealistic) call crater_realistic_topography(user,surf,crater,domain,ejecta_dem) 
          deallocate(ejecta_dem)

src/crater/module_crater.f90

@@ -97,15 +97,15 @@ end subroutine crater_averages
    end interface
 
    interface
-      subroutine crater_emplace(user,surf,crater,domain,deltaMtot,incval)
+      subroutine crater_emplace(user,surf,crater,domain,deltaMtot,incval,nmeltsheet)
       use module_globals
       implicit none
       type(usertype),intent(in) :: user
       type(surftype),dimension(:,:),intent(inout) :: surf
       type(cratertype),intent(inout) :: crater
       type(domaintype),intent(inout) :: domain
       real(DP),intent(out) :: deltaMtot
-      integer(I4B),intent(out) :: incval
+      integer(I4B),intent(out) :: incval,nmeltsheet
       end subroutine crater_emplace
    end interface
 

src/ejecta/ejecta_emplace.f90

@@ -75,7 +75,7 @@
 !                The cutoff of ejecta thickness is still buggy.  
 !
 !**********************************************************************************************************************************
-subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_resolution,cumulative_elchange)
+subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_resolution,cumulative_elchange,vmeltsheet)
    use module_globals
    use module_util
    use module_io
@@ -95,6 +95,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_r
    real(DP),intent(in)  :: age
    real(DP),intent(in)  :: age_resolution
    real(DP),dimension(:,:),allocatable,intent(out) :: cumulative_elchange
+   real(DP),intent(out) :: vmeltsheet
 
    ! Internal variables
    real(DP) :: lrad,lradsq
@@ -129,6 +130,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_r
 
    if (user%doregotrack) call regolith_melt_zone(user,crater,crater%imp,crater%impvel,rm,dm,totmelt)
    vmelt = 0.0_DP
+   vmeltsheet = 0.0_DP
 
    crater%vdepth = crater%rimheight + crater%floordepth
    crater%vrim   = crater%rimheight
@@ -290,6 +292,8 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_r
       write(*,*) 'Total Melt: ', totmelt
       write(*,*) 'ejected / total melt:', vmelt/totmelt
    end if
+   vmeltsheet = totmelt - vmelt   
+
    ejbmass = sum(cumulative_elchange)
 
    ! Create buffer to prevent infinite hole bug

src/ejecta/module_ejecta.f90

@@ -25,7 +25,7 @@ module module_ejecta
 save
 
    interface
-      subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_resolution,cumulative_elchange)
+      subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_resolution,cumulative_elchange,vmeltsheet)
       use module_globals
       implicit none
       type(usertype),intent(in) :: user
@@ -38,6 +38,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_r
       real(DP),intent(in)  :: age
       real(DP),intent(in)  :: age_resolution
       real(DP),dimension(:,:),allocatable,intent(inout) :: cumulative_elchange
+      real(DP),intent(out) :: vmeltsheet
       end subroutine ejecta_emplace
    end interface
 

src/regolith/module_regolith.f90

@@ -326,12 +326,14 @@ end function regolith_shock_damage
    end interface
 
    interface
-      subroutine regolith_interior(user,surf,crater,incval)
+      subroutine regolith_interior(user,surf,crater,domain,incval,nmeltsheet,vmeltsheet)
       use module_globals
       type(usertype),intent(in) :: user
       type(surftype),dimension(:,:),intent(inout) :: surf
+      type(domaintype),intent(in) :: domain
       type(cratertype),intent(in) :: crater
-      integer(I4B),intent(in)     :: incval
+      integer(I4B),intent(in)     :: incval, nmeltsheet
+      real(DP),intent(in)         :: vmeltsheet
       end subroutine regolith_interior
    end interface
 

src/regolith/regolith_interior.f90

@@ -9,16 +9,16 @@
 !  
 !
 !  Input
-!    Arguments : user,surf,crater,inc
+!    Arguments : user,surf,domain,crater,inc,nmeltsheet,vmeltsheet
 !
 !  Output
 !    Arguments : surf
 !           
 ! 
-!  Notes       :
+!  Notes       : nmeltsheet is the number of pixels that will contain the melt sheet; vmeltsheet is the melt sheet volume
 !
 !**********************************************************************************************************************************
-subroutine regolith_interior(user,surf,crater,incval)
+subroutine regolith_interior(user,surf,crater,domain,incval,nmeltsheet,vmeltsheet)
     use module_globals 
     use module_util
     use module_regolith, EXCEPT_THIS_ONE => regolith_interior
@@ -27,16 +27,22 @@ subroutine regolith_interior(user,surf,crater,incval)
     !Arguments
     type(usertype),intent(in) :: user
     type(surftype),dimension(:,:),intent(inout) :: surf
+    type(domaintype),intent(in) :: domain
     type(cratertype),intent(in) :: crater
-    integer(I4B),intent(in)     :: incval
+    integer(I4B),intent(in)     :: incval, nmeltsheet
+    real(DP),intent(in)         :: vmeltsheet
 
     !internal variables
     integer(I4B) xpi,ypi,i,j,inc,incsq,iradsq
-    real(DP) :: lradsq, x_relative, y_relative, xp, yp 
+    real(DP) :: lradsq, x_relative, y_relative, xp, yp, hmeltsheet 
     type(regodatatype),dimension(:),allocatable :: poppedarray
+    type(regodatatype) :: newlayer
 
     !Executable code
 
+    hmeltsheet = vmeltsheet / (nmeltsheet*user%gridsize*user%gridsize)
+    allocate(newlayer%meltdist(domain%rcnum))
+
     inc = incval
 
     do j=-inc,inc
@@ -58,10 +64,23 @@ subroutine regolith_interior(user,surf,crater,incval)
             call util_traverse_pop_array(surf(xpi,ypi)%regolayer,surf(xpi,ypi)%abselc,poppedarray)
             deallocate(poppedarray)
 
-            !Adding the melt sheet goes here! :)
+            !fill top layer with melt sheet of given thickness hmeltsheet
+            newlayer%meltfrac = 1.0_DP
+            newlayer%ejm = 0.0_DP
+            newlayer%ejmf = 0.0_DP
+            newlayer%thickness = hmeltsheet
+            newlayer%meltvolume = vmeltsheet / nmeltsheet
+            newlayer%totvolume = newlayer%meltvolume
+            newlayer%meltdist(:) = 0.0_SP
+            if(domain%currentqmc) then
+                newlayer%meltdist(domain%nqmc) = newlayer%meltfrac
+            end if
+            call util_push_array(surf(xpi,ypi)%regolayer,newlayer)
         end do
     end do
 
+    deallocate(newlayer%meltdist)
+
     return
 end subroutine regolith_interior