diff --git a/src/AnisotropicallyAdapted3DDistributedDiscreteModels.jl b/src/AnisotropicallyAdapted3DDistributedDiscreteModels.jl
index 363a0b5..10177bc 100644
--- a/src/AnisotropicallyAdapted3DDistributedDiscreteModels.jl
+++ b/src/AnisotropicallyAdapted3DDistributedDiscreteModels.jl
@@ -125,6 +125,81 @@ function vertically_adapt(model::OctreeDistributedDiscreteModel{Dc,Dp},
   return ref_model, adaptivity_glue
 end
 
+function _vertically_uniformly_refine!(model::OctreeDistributedDiscreteModel{Dc,Dp}) where {Dc,Dp}
+  pXest_type = model.pXest_type
+  function refine_layer_callback(p6est::Ptr{p6est_t},
+    which_tree::p4est_topidx_t,
+    column::Ptr{p4est_quadrant_t},
+    layer::Ptr{p2est_quadrant_t})
+    Cint(1)
+  end
+  
+  refine_layer_fn_callback_c = @cfunction($refine_layer_callback, Cint, 
+                                    (Ptr{p6est_t}, p4est_topidx_t, Ptr{p4est_quadrant_t}, Ptr{p2est_quadrant_t}))
+  
+  # # Copy input p4est, refine and balance
+  ptr_new_pXest = pXest_copy(pXest_type, model.ptr_pXest)
+  p6est_vertically_refine!(ptr_new_pXest,
+                           refine_layer_fn_callback_c,
+                           C_NULL)
+  ptr_new_pXest
+end 
+
+function vertically_uniformly_refine(model::OctreeDistributedDiscreteModel)  
+  ptr_new_pXest = _vertically_uniformly_refine!(model)
+
+  # Extract ghost and lnodes
+  ptr_pXest_ghost  = setup_pXest_ghost(model.pXest_type, ptr_new_pXest)
+  ptr_pXest_lnodes = setup_pXest_lnodes_nonconforming(model.pXest_type, ptr_new_pXest, ptr_pXest_ghost)
+
+  # Build fine-grid mesh
+  fmodel,non_conforming_glue = setup_non_conforming_distributed_discrete_model(model.pXest_type,
+                                                                               model.parts,
+                                                                               model.coarse_model,
+                                                                               model.ptr_pXest_connectivity,
+                                                                               ptr_new_pXest,
+                                                                               ptr_pXest_ghost,
+                                                                               ptr_pXest_lnodes)
+    
+  pXest_ghost_destroy(model.pXest_type,ptr_pXest_ghost)
+  pXest_lnodes_destroy(model.pXest_type,ptr_pXest_lnodes)
+
+  pXest_refinement_rule_type = PXestVerticalRefinementRuleType()
+  _refinement_and_coarsening_flags = map(partition(get_cell_gids(model))) do indices
+    flags  = Vector{Cint}(undef,length(local_to_global(indices)))
+    flags .= refine_flag
+  end
+
+  stride = pXest_stride_among_children(model.pXest_type,pXest_refinement_rule_type,model.ptr_pXest)
+  adaptivity_glue = _compute_fine_to_coarse_model_glue(model.pXest_type,
+                                                       pXest_refinement_rule_type,
+                                                       model.parts,
+                                                       model.dmodel,
+                                                       fmodel,
+                                                       _refinement_and_coarsening_flags,
+                                                       stride)
+  adaptive_models = map(local_views(model),
+                        local_views(fmodel),
+                        adaptivity_glue) do model, fmodel, glue 
+      Gridap.Adaptivity.AdaptedDiscreteModel(fmodel,model,glue)
+  end
+  fmodel = GridapDistributed.GenericDistributedDiscreteModel(adaptive_models,get_cell_gids(fmodel))
+  ref_model = OctreeDistributedDiscreteModel(3,3,
+                                             model.parts,
+                                             fmodel,
+                                             non_conforming_glue,
+                                             model.coarse_model,
+                                             model.ptr_pXest_connectivity,
+                                             ptr_new_pXest,
+                                             model.pXest_type,
+                                             model.pXest_refinement_rule_type,
+                                             false,
+                                             model)
+  return ref_model, adaptivity_glue
+end
+
+
+
 function setup_non_conforming_distributed_discrete_model(pXest_type::P6estType,
                                                          parts,
                                                          coarse_discrete_model,
diff --git a/src/GridapP4est.jl b/src/GridapP4est.jl
index a56bcca..bcb6c15 100644
--- a/src/GridapP4est.jl
+++ b/src/GridapP4est.jl
@@ -26,6 +26,7 @@ module GridapP4est
   export AnisotropicallyAdapted3DDistributedDiscreteModel
   export vertically_adapt 
   export horizontally_adapt
+  export vertically_uniformly_refine
   export nothing_flag, refine_flag, coarsen_flag
   export FixedFractionAdaptiveFlagsMarkingStrategy
   export update_adaptivity_flags!
diff --git a/src/OctreeDistributedDiscreteModels.jl b/src/OctreeDistributedDiscreteModels.jl
index 3f950bf..c839f10 100644
--- a/src/OctreeDistributedDiscreteModels.jl
+++ b/src/OctreeDistributedDiscreteModels.jl
@@ -241,7 +241,7 @@ function Init(a::OctreeDistributedDiscreteModel)
 end
 
 function Finalize(a::OctreeDistributedDiscreteModel)
-  # octree_distributed_discrete_model_free!(a)
+  octree_distributed_discrete_model_free!(a)
   return nothing
 end
 
diff --git a/test/PoissonAnisotropicOctreeModelsTests.jl b/test/PoissonAnisotropicOctreeModelsTests.jl
index 16d1c79..df0c5c0 100644
--- a/test/PoissonAnisotropicOctreeModelsTests.jl
+++ b/test/PoissonAnisotropicOctreeModelsTests.jl
@@ -245,6 +245,9 @@ module PoissonAnisotropicOctreeModelsTests
     rdmodel=dmodel
     for i=1:num_amr_steps
      rdmodel=test_transfer_ops_and_redistribute(ranks,rdmodel,order,T)
+     if i==1
+      rdmodel,_=vertically_uniformly_refine(dmodel)
+     end
     end
   end 
 
@@ -271,11 +274,8 @@ module PoissonAnisotropicOctreeModelsTests
     end
     for perm in (1,2), order in (1,4), scalar_or_vector in (:vector,)
       test(ranks,perm,order,_field_type(Val{3}(),scalar_or_vector))
-    end 
-    for order=1:1, scalar_or_vector in (:scalar,:vector)
-      test_3d(ranks,order,_field_type(Val{3}(),scalar_or_vector), num_amr_steps=4)
     end
-    for order=1:1, scalar_or_vector in (:scalar,)
+    for order=1:1, scalar_or_vector in (:scalar,:vector)
       test_3d(ranks,order,_field_type(Val{3}(),scalar_or_vector), num_amr_steps=4)
     end
   end