Updated the porosity function.

src/porosity/porosity_form_interior.f90

@@ -5,34 +5,35 @@
 !  Project     : CTEM
 !  Language    : Fortran 2003
 !
-!  Description : Compute the thickness of ejecta coverage due to a transient crater at each pixel. 
+!  Description : Compute the thickness of porosity development inside the crater. 
 !       
-!
 !  Input
 !    Arguments : user, surfi, crater, elchange, lradsq, lradsq 
-!              : elchange = the elevation change from the previous run to the latest one (defined in crater_form_interior.f90). 
 !              : lradsq   = the distance squared from the crater center (defined in crater_emplace.f90).
-!              : newelev  = the new elevation in the previous run at a given pixel (defined in crater_emplace.f90). 
 !              
 !  Output
 !    Arguments : surfi
 !           
-! 
 !  Notes       : The thickness of ejecta coverage is added to surfi%ejcov 
 !              : The current version only considers a binary mode: one being porosity, the other being non-porosity
 !
-!  Developer   : Toshi Hirabayashi (01/16/16)
+!  Developer   : Toshi Hirabayashi (12/08/16)
 !**********************************************************************************************************************************
-subroutine porosity_form_interior(user,surfi,crater,elchange,lradsq,newelev)
+subroutine porosity_form_interior(user, surfi, crater, lradsq)
    use module_globals
+   use module_util
    use module_porosity, EXCEPT_THIS_ONE => porosity_form_interior
    implicit none
 
    ! Arguments
-   type(usertype),intent(in) :: user
+   type(usertype),intent(in)    :: user
    type(surftype),intent(inout) :: surfi
-   type(cratertype),intent(in) :: crater
-   real(DP),intent(in) :: lradsq,elchange,newelev
+   type(cratertype),intent(in)  :: crater
+   real(DP),intent(in)          :: lradsq
+
+   ! Internal variables
+   type(regodatatype)          :: nlayer  ! Dummy parameters to search for the layers.  
+   type(regolisttype), pointer :: current => null() ! Dummy parameters to search for the layers. 
 
    ! Internal variables
    real(DP) :: newdem ! Same as surfi%dem
@@ -43,23 +44,29 @@ subroutine porosity_form_interior(user,surfi,crater,elchange,lradsq,newelev)
 
    ! Executable code
    ! Interface
-   newdem=surfi%dem
+   newdem = surfi%dem
 
-   !compute Transition crater depth
-   ht=0.3*crater%fcrat               !The depth of a transient crater   
+   ! compute Transient crater depth
+   ht = TRNRATIO * crater%fcrat      !The depth of a transient crater   
    parabht = ht / ((crater%frad)**2) !The factor of a parabolic profile  
    depthht = ht - (parabht * lradsq) !The depth at a given pixel.  
-   newtrn  = newelev - depthht       !The elevation of the bottom of the transient crater at a given pixel.
+   newtrn  = crater%melev - depthht  !The elevation of the bottom of the transient crater at a given pixel.
 
-   ! Impact makes regolith
-   if (surfi%ejcov < max(newdem - newtrn,0.0_DP)) then !This condition was made because depthht is the magnitude
-                                            !while newdem is the elevation. Thus, to have the difference between them, 
-                                            !depthht + newdem. 
-   		surfi%ejcov = max(surfi%ejcov + newdem - newtrn,0.0_DP)
-   else 
-   		surfi%ejcov = max(surfi%ejcov + min(elchange,0.0_DP),0.0_DP) 
-   end if 	
+	nlayer%porosity = 0.2;            !The porosity of the first value. This value is an assumed value. 
+	nlayer%depth    = newtrn; 
 
+	! Currently consider only the first layer. 
+	if (surfi%porolayer%regodata%porosity == nlayer%porosity) then
+		! just update if the newly calculated layer is deeper than the stored one. 
+		if (surfi%porolayer%regodata%depth > nlayer%depth) then
+			surfi%porolayer%regodata%depth    = nlayer%depth
+			surfi%porolayer%regodata%porosity = nlayer%porosity
+		end if
+	else
+		! Linked list if there is only one porosity layer
+		call util_push(surfi%porolayer, nlayer)
+	end if 
+
    return
 end subroutine porosity_form_interior