diff --git a/src/bindings/c/CMakeLists.txt b/src/bindings/c/CMakeLists.txt
index 193548f5..e336d49d 100644
--- a/src/bindings/c/CMakeLists.txt
+++ b/src/bindings/c/CMakeLists.txt
@@ -13,6 +13,7 @@ set(C_FORTRAN_LIB_SRCS
     growtree.f90
     indices.f90
     imports.f90
+    lymphatics.f90
     pressure_resistance_flow.f90
     species_transport.f90
     surface_fitting.f90
@@ -33,6 +34,7 @@ set(C_C_LIB_SRCS
     growtree.c
     indices.c
     imports.c
+    lymphatics.c
     pressure_resistance_flow.c
     species_transport.c
     surface_fitting.c
@@ -54,6 +56,7 @@ set(C_LIB_HDRS
   growtree.h
   indices.h
   imports.h
+  lymphatics.h
   pressure_resistance_flow.h
   species_transport.h
   surface_fitting.h
diff --git a/src/bindings/c/arrays.c b/src/bindings/c/arrays.c
index 7e7b8a14..9e325a99 100644
--- a/src/bindings/c/arrays.c
+++ b/src/bindings/c/arrays.c
@@ -3,10 +3,15 @@
 #include "string.h"
 
 void set_node_field_value_c(int *row, int *col, double *value);
+void update_parameter_c(const char *parameter_name, int *parameter_name_len, double *parameter_value);
 
 void set_node_field_value(int row, int col, double value)
 {
   set_node_field_value_c(&row, &col, &value);
 }
 
-
+void update_parameter(const char *parameter_name, double parameter_value)
+{
+  int parameter_name_len = (int)strlen(parameter_name);
+  update_parameter_c(parameter_name, &parameter_name_len, &parameter_value);
+}
diff --git a/src/bindings/c/arrays.f90 b/src/bindings/c/arrays.f90
index 8efb0e71..ed73ec35 100644
--- a/src/bindings/c/arrays.f90
+++ b/src/bindings/c/arrays.f90
@@ -26,6 +26,24 @@ subroutine set_node_field_value_c(row, col, value) bind(C, name="set_node_field_
 
   end subroutine set_node_field_value_c
 
+  subroutine update_parameter_c(parameter_name, parameter_name_len, parameter_value) bind(C, name="update_parameter_c")
+    use iso_c_binding, only: c_ptr
+    use utils_c, only: strncpy
+    use other_consts, only: max_filename_len
+    use arrays, only: update_parameter
+    implicit none
     
-
+    integer, intent(in) :: parameter_name_len
+    real(dp), intent(in) :: parameter_value
+    type(c_ptr),value, intent(in) :: parameter_name
+    character(len=max_filename_len) :: parameter_name_f
+    
+    call strncpy(parameter_name_f, parameter_name, parameter_name_len)
+#if defined _WIN32 && defined __INTEL_COMPILER
+    call so_update_parameter(parameter_name_f, parameter_value)
+#else
+    call update_parameter(parameter_name_f, parameter_value)
+#endif
+  end subroutine update_parameter_c
+  
 end module arrays_c
diff --git a/src/bindings/c/arrays.h b/src/bindings/c/arrays.h
index 1468b841..15495724 100644
--- a/src/bindings/c/arrays.h
+++ b/src/bindings/c/arrays.h
@@ -5,5 +5,6 @@
 #include "symbol_export.h"
 
 SHO_PUBLIC void set_node_field_value(int row, int col, double value);
+SHO_PUBLIC void update_parameter(const char *parameter_name, double parameter_value);
 
 #endif /* AETHER_ARRAYS_H */
diff --git a/src/bindings/c/exports.c b/src/bindings/c/exports.c
index e3c506ab..bf530142 100644
--- a/src/bindings/c/exports.c
+++ b/src/bindings/c/exports.c
@@ -17,6 +17,7 @@ void export_node_field_c(int *nj_field, const char *EXNODEFIELD, int *EXNODEFIEL
                          const char *name, int *name_len, const char *field_name, int *field_name_len);
 void export_elem_geometry_2d_c(const char *EXELEMFILE, int *EXELEMFILE_LEN,
 			       const char *name, int *name_len, int *offset_elem, int *offset_node);
+void export_terminal_lymphatic_c(const char *EXNODEFILE, int *EXNODEFILE_LEN, const char *name, int *name_len);
 void export_terminal_solution_c(const char *EXNODEFILE, int *EXNODEFILE_LEN,
 				const char *name, int *name_len);
 void export_terminal_perfusion_c(const char *EXNODEFILE, int *EXNODEFILE_LEN,
@@ -87,6 +88,14 @@ void export_node_field(int nj_field, const char *EXNODEFIELD,
 		      &name_len, field_name, &field_name_len);
 }
 
+void export_terminal_lymphatic(const char *EXNODEFILE, const char *name)
+{
+  int filename_len = strlen(EXNODEFILE);
+  int name_len = strlen(name);
+
+  export_terminal_lymphatic_c(EXNODEFILE, &filename_len, name, &name_len);
+}
+
 void export_terminal_solution(const char *EXNODEFILE, const char *name)
 {
   int filename_len = strlen(EXNODEFILE);
diff --git a/src/bindings/c/exports.f90 b/src/bindings/c/exports.f90
index 5d7416d3..6a9c535e 100644
--- a/src/bindings/c/exports.f90
+++ b/src/bindings/c/exports.f90
@@ -164,6 +164,31 @@ subroutine export_data_geometry_c(EXDATAFILE, filename_len, name, name_len, offs
 
   end subroutine export_data_geometry_c
 
+!!!########################################################################
+
+  subroutine export_terminal_lymphatic_c(EXNODEFILE, filename_len, name, name_len) bind(C, name="export_terminal_lymphatic_c")
+
+    use iso_c_binding, only: c_ptr
+    use utils_c, only: strncpy
+    use exports, only: export_terminal_lymphatic
+    use other_consts, only: MAX_STRING_LEN, MAX_FILENAME_LEN
+    implicit none
+    integer,intent(in) :: filename_len, name_len
+    type(c_ptr), value, intent(in) :: EXNODEFILE, name
+    character(len=MAX_FILENAME_LEN) :: filename_f
+    character(len=MAX_STRING_LEN) :: name_f
+
+    call strncpy(filename_f, EXNODEFILE, filename_len)
+    call strncpy(name_f, name, name_len)
+
+#if defined _WIN32 && defined __INTEL_COMPILER
+    call so_export_terminal_lymphatic(filename_f, name_f)
+#else
+    call export_terminal_lymphatic(filename_f, name_f)
+#endif
+
+  end subroutine export_terminal_lymphatic_c
+
 !!!########################################################################
 
   subroutine export_terminal_solution_c(EXNODEFILE, filename_len, name, name_len) bind(C, name="export_terminal_solution_c")
diff --git a/src/bindings/c/exports.h b/src/bindings/c/exports.h
index 17d7ce89..4aaaa0e3 100644
--- a/src/bindings/c/exports.h
+++ b/src/bindings/c/exports.h
@@ -9,6 +9,7 @@ SHO_PUBLIC void export_1d_elem_field(int ne_field, const char *EXELEMFILE, const
 SHO_PUBLIC void export_1d_elem_geometry(const char *EXELEMFILE, const char *name);
 SHO_PUBLIC void export_elem_geometry_2d(const char *EXELEMFILE, const char *name, int offset_elem, int offset_node);
 SHO_PUBLIC void export_node_field(int nj_field, const char *EXNODEFIELD, const char *name, const char *field_name);
+SHO_PUBLIC void export_terminal_lymphatic(const char *EXNODEFILE, const char *name);
 SHO_PUBLIC void export_terminal_solution(const char *EXNODEFILE, const char *name);
 SHO_PUBLIC void export_terminal_perfusion(const char *EXNODEFILE, const char *name);
 SHO_PUBLIC void export_terminal_ssgexch(const char *EXNODEFILE, const char *name);
diff --git a/src/bindings/c/imports.c b/src/bindings/c/imports.c
index 1c7dc2b2..873f51fb 100644
--- a/src/bindings/c/imports.c
+++ b/src/bindings/c/imports.c
@@ -4,9 +4,21 @@
 
 #include <string.h>
 
+void import_capillary_c(const char *FLOWFILE,int *FLOWFILE_LEN);
+void import_terminal_c(const char *FLOWFILE,int *FLOWFILE_LEN);
 void import_ventilation_c(const char *FLOWFILE,int *FLOWFILE_LEN);
 void import_perfusion_c(const char *FLOWFILE,int *FLOWFILE_LEN);
 
+void import_capillary(const char *FLOWFILE)
+{
+	int filename_len = strlen(FLOWFILE);
+	import_capillary_c(FLOWFILE, &filename_len);
+}
+void import_terminal(const char *FLOWFILE)
+{
+	int filename_len = strlen(FLOWFILE);
+	import_terminal_c(FLOWFILE, &filename_len);
+}
 void import_ventilation(const char *FLOWFILE)
 {
 	int filename_len = strlen(FLOWFILE);
diff --git a/src/bindings/c/imports.f90 b/src/bindings/c/imports.f90
index fd63276c..b9e0aec1 100644
--- a/src/bindings/c/imports.f90
+++ b/src/bindings/c/imports.f90
@@ -3,6 +3,54 @@ module imports_c
   private
 
 contains
+!
+!###################################################################################
+!
+  subroutine import_capillary_c(FLOWFILE, filename_len) bind(C, name="import_capillary_c")
+  
+    use iso_c_binding, only: c_ptr
+    use utils_c, only: strncpy
+    use imports, only: import_capillary
+    use other_consts, only: MAX_STRING_LEN, MAX_FILENAME_LEN
+    implicit none
+    integer,intent(in) :: filename_len
+    type(c_ptr), value, intent(in) :: FLOWFILE
+    character(len=MAX_FILENAME_LEN) :: filename_f
+
+    call strncpy(filename_f, FLOWFILE, filename_len)
+
+#if defined _WIN32 && defined __INTEL_COMPILER
+    call so_import_capillary(filename_f)
+#else
+    call import_capillary(filename_f)
+#endif
+
+  end subroutine import_capillary_c
+  
+!
+!###################################################################################
+!
+  subroutine import_terminal_c(FLOWFILE, filename_len) bind(C, name="import_terminal_c")
+  
+    use iso_c_binding, only: c_ptr
+    use utils_c, only: strncpy
+    use imports, only: import_terminal
+    use other_consts, only: MAX_STRING_LEN, MAX_FILENAME_LEN
+    implicit none
+    integer,intent(in) :: filename_len
+    type(c_ptr), value, intent(in) :: FLOWFILE
+    character(len=MAX_FILENAME_LEN) :: filename_f
+
+    call strncpy(filename_f, FLOWFILE, filename_len)
+
+#if defined _WIN32 && defined __INTEL_COMPILER
+    call so_import_terminal(filename_f)
+#else
+    call import_terminal(filename_f)
+#endif
+
+  end subroutine import_terminal_c
+  
 !
 !###################################################################################
 !
diff --git a/src/bindings/c/imports.h b/src/bindings/c/imports.h
index 9a5c6124..3232df27 100644
--- a/src/bindings/c/imports.h
+++ b/src/bindings/c/imports.h
@@ -3,6 +3,8 @@
 
 #include "symbol_export.h"
 
+SHO_PUBLIC void import_capillary(const char *FLOWFILE);
+SHO_PUBLIC void import_terminal(const char *FLOWFILE);
 SHO_PUBLIC void import_ventilation(const char *FLOWFILE);
 SHO_PUBLIC void import_perfusion(const char *FLOWFILE);
 
diff --git a/src/bindings/c/indices.c b/src/bindings/c/indices.c
index 3ab18822..20174baa 100644
--- a/src/bindings/c/indices.c
+++ b/src/bindings/c/indices.c
@@ -4,6 +4,7 @@
 #include <string.h>
 
 void define_problem_type_c(const char *PROBLEMTYPE,int *PROBLEMTYPE_LEN);
+void lymphatic_indices_c();
 void ventilation_indices_c();
 void perfusion_indices_c();
 int get_ne_radius_c();
@@ -18,6 +19,11 @@ void define_problem_type(const char *PROBLEMTYPE)
 }
 
 
+void lymphatic_indices()
+{
+  lymphatic_indices_c();
+}
+
 void ventilation_indices()
 {
   ventilation_indices_c();
diff --git a/src/bindings/c/indices.f90 b/src/bindings/c/indices.f90
index df308ba6..6aa51b2d 100644
--- a/src/bindings/c/indices.f90
+++ b/src/bindings/c/indices.f90
@@ -35,6 +35,23 @@ end subroutine ventilation_indices_c
 !
 !######################################################################
 !
+!> Lymphatic indices
+  subroutine lymphatic_indices_c() bind(C, name="lymphatic_indices_c")
+
+    use indices, only: lymphatic_indices
+    implicit none
+
+#if defined _WIN32 && defined __INTEL_COMPILER
+    call so_lymphatic_indices()
+#else
+    call lymphatic_indices()
+#endif
+
+  end subroutine lymphatic_indices_c
+
+!
+!######################################################################
+!
 !> Perfusion indices
   subroutine perfusion_indices_c() bind(C, name="perfusion_indices_c")
 
diff --git a/src/bindings/c/indices.h b/src/bindings/c/indices.h
index f1ec8f87..2c9c6907 100644
--- a/src/bindings/c/indices.h
+++ b/src/bindings/c/indices.h
@@ -3,6 +3,7 @@
 
 #include "symbol_export.h"
 SHO_PUBLIC void define_problem_type(const char *PROBLEMTYPE);
+SHO_PUBLIC void lymphatic_indices();
 SHO_PUBLIC void ventilation_indices();
 SHO_PUBLIC void perfusion_indices();
 SHO_PUBLIC int get_ne_radius();
diff --git a/src/bindings/c/lymphatics.c b/src/bindings/c/lymphatics.c
new file mode 100644
index 00000000..e0d64424
--- /dev/null
+++ b/src/bindings/c/lymphatics.c
@@ -0,0 +1,17 @@
+#include "lymphatics.h"
+#include "string.h"
+
+void alveolar_capillary_flux_c(int *num_nodes, int *write_out);
+void lymphatic_transport_c(const char *filename, int *filename_len);
+
+
+void alveolar_capillary_flux(int num_nodes, int write_out)
+{
+  alveolar_capillary_flux_c(&num_nodes, &write_out);
+}
+
+void lymphatic_transport(const char *filename)
+{
+  int filename_len = (int)strlen(filename);
+  lymphatic_transport_c(filename, &filename_len);
+}
diff --git a/src/bindings/c/lymphatics.f90 b/src/bindings/c/lymphatics.f90
new file mode 100644
index 00000000..3bcf0054
--- /dev/null
+++ b/src/bindings/c/lymphatics.f90
@@ -0,0 +1,53 @@
+module lymphatics_c
+  implicit none
+
+  private
+
+contains
+!
+!###################################################################################
+!
+!*alveolar_capillary_flux:* 
+  subroutine alveolar_capillary_flux_c(num_nodes,write_out) bind(C, name="alveolar_capillary_flux_c")
+    use lymphatics,only: alveolar_capillary_flux
+    implicit none
+
+    integer,intent(in) :: num_nodes
+    logical,intent(in) :: write_out
+
+#if defined _WIN32 && defined __INTEL_COMPILER
+    call so_alveolar_capillary_flux(num_nodes,write_out)
+#else
+    call alveolar_capillary_flux(num_nodes,write_out)
+#endif
+    
+  end subroutine alveolar_capillary_flux_c
+!
+!###################################################################################
+!
+!*lymphatic_transport* 
+  subroutine lymphatic_transport_c(filename, filename_len) bind(C, name="lymphatic_transport_c")
+    use iso_c_binding, only: c_ptr
+    use utils_c, only: strncpy
+    use other_consts, only: MAX_FILENAME_LEN
+    use lymphatics,only: lymphatic_transport
+    implicit none
+    
+    integer,intent(in) :: filename_len
+    type(c_ptr), value, intent(in) :: filename
+    character(len=MAX_FILENAME_LEN) :: filename_f
+
+    call strncpy(filename_f, filename, filename_len)
+#if defined _WIN32 && defined __INTEL_COMPILER
+    call so_lymphatic_transport(filename_f)
+#else
+    call lymphatic_transport(filename_f)
+#endif
+    
+  end subroutine lymphatic_transport_c
+
+!
+!###################################################################################
+!
+
+end module lymphatics_c
diff --git a/src/bindings/c/lymphatics.h b/src/bindings/c/lymphatics.h
new file mode 100644
index 00000000..83d0e77b
--- /dev/null
+++ b/src/bindings/c/lymphatics.h
@@ -0,0 +1,9 @@
+#ifndef AETHER_LYMPHATICS_H
+#define AETHER_LYMPHATICS_H
+
+#include "symbol_export.h"
+
+SHO_PUBLIC void alveolar_capillary_flux(int num_nodes, int write_out);
+SHO_PUBLIC void lymphatic_transport(const char *filename);
+
+#endif /* AETHER_LYMPHATICS_H */
diff --git a/src/bindings/interface/lymphatics.i b/src/bindings/interface/lymphatics.i
new file mode 100644
index 00000000..99363875
--- /dev/null
+++ b/src/bindings/interface/lymphatics.i
@@ -0,0 +1,10 @@
+
+%module(package="aether") lymphatics
+%include symbol_export.h
+%include lymphatics.h
+
+%{
+#include "lymphatics.h"
+%}
+
+
diff --git a/src/bindings/python/CMakeLists.txt b/src/bindings/python/CMakeLists.txt
index ce3301ba..2cebb1fd 100644
--- a/src/bindings/python/CMakeLists.txt
+++ b/src/bindings/python/CMakeLists.txt
@@ -29,6 +29,7 @@ set(INTERFACE_SRCS
   ../interface/growtree.i
   ../interface/indices.i
   ../interface/imports.i
+  ../interface/lymphatics.i
   ../interface/pressure_resistance_flow.i
   ../interface/species_transport.i
   ../interface/surface_fitting.i
diff --git a/src/lib/.lymphatics.f90.swp b/src/lib/.lymphatics.f90.swp
new file mode 100644
index 00000000..cde3c770
Binary files /dev/null and b/src/lib/.lymphatics.f90.swp differ
diff --git a/src/lib/CMakeLists.txt b/src/lib/CMakeLists.txt
index 0fed5c73..6bd78576 100644
--- a/src/lib/CMakeLists.txt
+++ b/src/lib/CMakeLists.txt
@@ -13,6 +13,7 @@ set(LIB_SRCS
   growtree.f90
   indices.f90
   imports.f90
+  lymphatics.f90
   math_utilities.f90
   mesh_utilities.f90
   other_consts.f90
diff --git a/src/lib/arrays.f90 b/src/lib/arrays.f90
index e642461d..fcf5e533 100644
--- a/src/lib/arrays.f90
+++ b/src/lib/arrays.f90
@@ -126,6 +126,17 @@ module arrays
      real(dp) :: air_viscosity = 1.8e-5_dp        ! Pa.s
      real(dp) :: air_density = 1.146e-6_dp        ! g.mm^-3
   end type fluid_properties
+
+  type default_lymphatic_properties
+     ! default values for lymphatic model parameters
+     real(dp) :: lymphatic_density = 1.0_dp !to be calculated from CT??
+     real(dp) :: lymphatic_integrity = 1.0_dp ! a measure of how 'leaky' the lymphatic vessels are and prone to backflow
+     real(dp) :: reflection_coefficient = 0.0_dp
+     real(dp) :: test_time = 86400.0_dp !number of seconds for test to be run
+  end type default_lymphatic_properties
+
+!!! arrays that start with default values, updated during simulations/by user
+  type(default_lymphatic_properties) :: lymphatic_properties
   
 ! temporary, for debugging:
   real(dp) :: unit_before
@@ -141,8 +152,8 @@ module arrays
        elem_units_below, maxgen,capillary_bf_parameters, zero_tol,loose_tol,gasex_field, &
        num_lines_2d, lines_2d, line_versn_2d, lines_in_elem, nodes_in_line, elems_2d, &
        elem_cnct_2d, elem_nodes_2d, elem_versn_2d, elem_lines_2d, elems_at_node_2d, arclength, &
-       scale_factors_2d, fluid_properties, elasticity_vessels, admittance_param, &
-       elasticity_param, two_parameter, three_parameter, four_parameter, all_admit_param, &
+       scale_factors_2d, fluid_properties, lymphatic_properties, elasticity_vessels, admittance_param, &
+       elasticity_param, two_parameter, three_parameter, four_parameter, all_admit_param, update_parameter, &
        mesh_from_depvar, depvar_at_node, depvar_at_elem, SparseCol, SparseRow, triangle, &
        update_resistance_entries, vertex_xyz, &
        SparseVal, RHS, prq_solution, solver_solution, FIX
@@ -158,4 +169,26 @@ subroutine set_node_field_value(row, col, value)
 
   end subroutine set_node_field_value
 
+  subroutine update_parameter(parameter_name, parameter_value)
+    !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_UPDATE_PARAMETER" :: UPDATE_PARAMETER
+    implicit none
+    real(dp), intent(in) :: parameter_value
+    character(len=*), intent(in) :: parameter_name
+    
+    select case(parameter_name)
+       
+!!! lymphatic_properties
+    case('lymphatic_density')
+       lymphatic_properties%lymphatic_density = parameter_value
+    case('lymphatic_integrity')
+       lymphatic_properties%lymphatic_integrity = parameter_value
+    case('reflection_coefficient')
+       lymphatic_properties%reflection_coefficient = parameter_value
+    case('test_time')
+       lymphatic_properties%test_time = parameter_value
+
+    end select
+    
+  end subroutine update_parameter
+
 end module arrays
diff --git a/src/lib/capillaryflow.f90 b/src/lib/capillaryflow.f90
index 09f68fbe..ffdbe6ab 100644
--- a/src/lib/capillaryflow.f90
+++ b/src/lib/capillaryflow.f90
@@ -136,7 +136,7 @@ subroutine evaluate_ladder(ne,NonZeros,MatrixSize,submatrixsize,ngen,&
     real(dp) :: area_scale,length_scale,alpha_c
 
     ! Local variables
-    integer :: i,iter,j,gen,zone,num_sheet
+    integer :: i,iter,j,gen,zone,num_sheet,nunit
     integer, allocatable :: SparseCol(:)
     integer, allocatable ::SparseRow(:)
     real(dp) :: area_new,ErrorEstimate,Hart,Hven,Pin_SHEET,Pout_SHEET
@@ -332,9 +332,17 @@ subroutine evaluate_ladder(ne,NonZeros,MatrixSize,submatrixsize,ngen,&
 ! Rtot=9 Pa/mm^3 | Blood_vol=10 mm^3| sheet_area= 11 mm^2 | ave_TT=12 s |ave_H=13 um |Ppl=14 Pa
           WRITE(20,&
         '(I6,X,5(F9.2,X),2(F8.5,X),F10.2,X,F8.4,X,F10.4,X,F10.3,X,F8.4,X,F9.4,X)') &
-         ne,x,y,z,Pin,Pout,Q01_mthrees*1.d9,Qtot*1.d9,Rtot/1000.d0**3.d0,&
+         ne,x,y,z,Pin,Pout,Q01_mthrees*1.0e+9_dp,Qtot*1.0e+9_dp,Rtot/1000.0_dp**3.0_dp, &
          TOTAL_CAP_VOL,TOTAL_SHEET_SA,TT_TOTAL,TOTAL_SHEET_H,Ppl
-        ENDIF
+
+          ! record the unit values for mean pressure, transit time, surface area.
+          ! ne is the 'linker' element, so nunit is for the parent element
+          nunit = int(elem_field(ne_unit,elem_cnct(-1,1,ne))) 
+          unit_field(nu_blood_press,nunit) = (Pin+Pout)/2.0_dp
+          unit_field(nu_tt,nunit) = TT_TOTAL
+          unit_field(nu_sa,nunit) = TOTAL_SHEET_SA
+
+       ENDIF
 !
     deallocate (SparseCol, STAT = AllocateStatus)
     deallocate (SparseRow, STAT = AllocateStatus)
diff --git a/src/lib/exports.f90 b/src/lib/exports.f90
index 3606b677..d3c587a1 100644
--- a/src/lib/exports.f90
+++ b/src/lib/exports.f90
@@ -27,6 +27,7 @@ module exports
        export_node_geometry_2d,&
        export_node_field, &
        export_elem_field, &
+       export_terminal_lymphatic, &
        export_terminal_solution, &
        export_terminal_perfusion,&
        export_terminal_ssgexch, &
@@ -766,6 +767,95 @@ subroutine export_data_geometry(EXDATAFILE, name, offset)
 
   end subroutine export_data_geometry
 
+!
+!##############################################################################
+!
+
+  subroutine export_terminal_lymphatic(EXNODEFILE, name)
+  !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_EXPORT_TERMINAL_LYMPHATIC" :: EXPORT_TERMINAL_LYMPHATIC
+
+!!! Parameters
+    character(len=MAX_FILENAME_LEN),intent(in) :: EXNODEFILE
+    character(len=MAX_STRING_LEN),intent(in) :: name
+
+!!! Local Variables
+    integer :: len_end,ne,nj,NOLIST,np,np_last,VALUE_INDEX
+    logical :: FIRST_NODE
+
+    len_end=len_trim(name)
+    if(num_units.GT.0) THEN
+       open(10, file=EXNODEFILE, status='replace')
+       !**     write the group name
+       write(10,'( '' Group name: '',A)') name(:len_end)
+       FIRST_NODE=.TRUE.
+       np_last=1
+       !*** Exporting Terminal Solution
+       do nolist=1,num_units
+          if(nolist.GT.1) np_last = np
+          ne=units(nolist)
+          np=elem_nodes(2,ne)
+          !*** Write the field information
+          VALUE_INDEX=1
+          if(FIRST_NODE)THEN
+             write(10,'( '' #Fields=7'' )')
+             write(10,'('' 1) coordinates, coordinate, rectangular cartesian, #Components=3'')')
+             do nj=1,3
+                if(nj.eq.1) write(10,'(2X,''x.  '')',advance="no")
+                if(nj.eq.2) write(10,'(2X,''y.  '')',advance="no")
+                if(nj.eq.3) write(10,'(2X,''z.  '')',advance="no")
+                write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0
+                VALUE_INDEX=VALUE_INDEX+1
+             enddo
+
+             write(10,'('' 2) flux, field, rectangular cartesian, #Components=1'')')
+             write(10,'(2X,''1.  '')',advance="no")
+             write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0
+             VALUE_INDEX=VALUE_INDEX+1
+
+             write(10,'('' 3) int_sat, field, rectangular cartesian, #Components=1'')')
+             write(10,'(2X,''1.  '')',advance="no")
+             write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0
+             VALUE_INDEX=VALUE_INDEX+1
+
+             write(10,'('' 4) transit_time, field, rectangular cartesian, #Components=1'')')
+             write(10,'(2X,''1.  '')',advance="no")
+             write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0
+             VALUE_INDEX=VALUE_INDEX+1
+
+             write(10,'('' 5) cap_pressure, field, rectangular cartesian, #Components=1'')')
+             write(10,'(2X,''1.  '')',advance="no")
+             write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0
+             VALUE_INDEX=VALUE_INDEX+1
+
+             write(10,'('' 6) cap_SA, field, rectangular cartesian, #Components=1'')')
+             write(10,'(2X,''1.  '')',advance="no")
+             write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0
+             VALUE_INDEX=VALUE_INDEX+1
+
+             write(10,'('' 7) lymf_flow, field, rectangular cartesian, #Components=1'')')
+             write(10,'(2X,''1.  '')',advance="no")
+             write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0
+          endif !FIRST_NODE
+          !***      write the node
+          write(10,'(1X,''Node: '',I12)') np
+          do nj=1,3
+             write(10,'(2X,4(1X,F12.6))') (node_xyz(nj,np))      !Coordinates
+          enddo !njj2
+          write(10,'(2X,4(1X,e12.6))') (unit_field(nu_av_flux,NOLIST))
+          write(10,'(2X,4(1X,e12.6))') (unit_field(nu_intsat,nolist))
+          write(10,'(2X,4(1X,e12.6))') (unit_field(nu_tt,nolist))
+          write(10,'(2X,4(1X,e12.6))') (unit_field(nu_blood_press,nolist))
+          write(10,'(2X,4(1X,e12.6))') (unit_field(nu_sa,nolist))
+          write(10,'(2X,4(1X,e12.6))') (unit_field(nu_lymphflow,nolist))
+          
+          FIRST_NODE=.FALSE.
+          np_last=np
+       enddo !nolist (np)
+    endif !num_nodes
+    close(10)
+
+  end subroutine export_terminal_lymphatic
+  
 !
 !##############################################################################
 !
@@ -803,7 +893,7 @@ subroutine export_terminal_solution(EXNODEFILE, name)
           !*** Write the field information
           VALUE_INDEX=1
           if(FIRST_NODE)THEN
-             write(10,'( '' #Fields=5'' )')
+             write(10,'( '' #Fields=8'' )')
              write(10,'('' 1) coordinates, coordinate, rectangular cartesian, #Components=3'')')
              do nj=1,3
                 if(nj.eq.1) write(10,'(2X,''x.  '')',advance="no")
@@ -812,32 +902,38 @@ subroutine export_terminal_solution(EXNODEFILE, name)
                 write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0
                 VALUE_INDEX=VALUE_INDEX+1
              enddo
+             !element number
+             write(10,'('' 2) terminal_element, field, rectangular cartesian, #Components=1'')')
+             write(10,'(2X,''1.  '')',advance="no")
+             write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0
+             VALUE_INDEX=VALUE_INDEX+1
              !Ventilation (tidal volume/insp time)
-             write(10,'('' 2) flow, field, rectangular cartesian, #Components=1'')')
+             write(10,'('' 3) flow, field, rectangular cartesian, #Components=1'')')
              write(10,'(2X,''1.  '')',advance="no")
              write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0
-             !VALUE_INDEX=VALUE_INDEX+1
-             !Volume
-             !write(10,'('' 3) volume, field, rectangular cartesian, #Components=1'')')
-             !write(10,'(2X,''1.  '')',advance="no")
-             !write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0
-             !VALUE_INDEX=VALUE_INDEX+1
-             !!Pressure
-             !write(10,'('' 4) pressure, field, rectangular cartesian, #Components=1'')')
-             !write(10,'(2X,''1.  '')',advance="no")
-             !write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0
+             VALUE_INDEX=VALUE_INDEX+1
              !Compliance
-             write(10,'('' 5) compliance, field, rectangular cartesian, #Components=1'')')
+             write(10,'('' 4) compliance, field, rectangular cartesian, #Components=1'')')
              write(10,'(2X,''1.  '')',advance="no")
              write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0
              VALUE_INDEX=VALUE_INDEX+1
              !Pleural pressure
-             write(10,'('' 6) pleural pressure, field, rectangular cartesian, #Components=1'')')
+             write(10,'('' 5) pleural_pressure, field, rectangular cartesian, #Components=1'')')
              write(10,'(2X,''1.  '')',advance="no")
              write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0
              VALUE_INDEX=VALUE_INDEX+1
              !Tidal volume
-             write(10,'('' 7) tidal volume, field, rectangular cartesian, #Components=1'')')
+             write(10,'('' 6) tidal_volume, field, rectangular cartesian, #Components=1'')')
+             write(10,'(2X,''1.  '')',advance="no")
+             write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0
+             VALUE_INDEX=VALUE_INDEX+1
+             !maximum Pe at a unit
+             write(10,'('' 7) max_Pe, field, rectangular cartesian, #Components=1'')')
+             write(10,'(2X,''1.  '')',advance="no")
+             write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0
+             VALUE_INDEX=VALUE_INDEX+1
+             !minimum Pe at a unit
+             write(10,'('' 8) min_Pe, field, rectangular cartesian, #Components=1'')')
              write(10,'(2X,''1.  '')',advance="no")
              write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0
           endif !FIRST_NODE
@@ -846,12 +942,14 @@ subroutine export_terminal_solution(EXNODEFILE, name)
           do nj=1,3
              write(10,'(2X,4(1X,F12.6))') (node_xyz(nj,np))      !Coordinates
           enddo !njj2
+          write(10,'(2X,4(1X,F12.6))')  real(ne) !element number
           write(10,'(2X,4(1X,F12.6))') (unit_field(nu_vent,NOLIST)) !Ventilation
-          !write(10,'(2X,4(1X,F12.6))') (unit_field(nu_vol,nolist))   !Volume (end expiration)
-          !write(10,'(2X,4(1X,F12.6))') (unit_field(nu_press,nolist)) !Pressure
           write(10,'(2X,4(1X,F12.6))') (unit_field(nu_comp,nolist))  !Compliance (end exp)
           write(10,'(2X,4(1X,F12.6))') (unit_field(nu_pe,nolist))    !Recoil pressure
           write(10,'(2X,4(1X,F12.6))') (unit_field(nu_vt,nolist))    !Tidal volume
+          write(10,'(2X,4(1X,F12.6))') (unit_field(nu_Pe_max,nolist))    !maximum elastic recoil
+          write(10,'(2X,4(1X,F12.6))') (unit_field(nu_Pe_min,nolist))    !minimum elastic recoil
+          
           FIRST_NODE=.FALSE.
           np_last=np
        enddo !nolist (np)
diff --git a/src/lib/geometry.f90 b/src/lib/geometry.f90
index 558b121d..0f4d021d 100644
--- a/src/lib/geometry.f90
+++ b/src/lib/geometry.f90
@@ -416,13 +416,15 @@ subroutine append_units()
     unit_field=0.0_dp
     units=0
     elem_units_below(1:num_elems) = 0 !initialise the number of terminal units below a branch
-
+    elem_field(ne_unit,:) = 0.0_dp
+    
     nu=0
     do ne=1,num_elems
        if(elem_cnct(1,0,ne).eq.0)THEN
           nu=nu+1
           units(nu)=ne     !Set up units array containing terminals
           elem_units_below(ne)=1
+          elem_field(ne_unit,ne) = real(nu)
        endif
     enddo
 
diff --git a/src/lib/imports.f90 b/src/lib/imports.f90
index 4acf6d16..040b8e58 100644
--- a/src/lib/imports.f90
+++ b/src/lib/imports.f90
@@ -24,11 +24,62 @@ module imports
   !Module variables
 
   !Interfaces
-  private 
+  private
+  public import_capillary
   public import_ventilation
   public import_perfusion
-
+  public import_terminal
+  
 contains
+  !
+  !###########################################################################################
+  !
+  !>*import_capillary:* This subroutine reads in the results for the micro-circulatory
+  ! components (capillary bed within 'units') of a perfusion model .
+  subroutine import_capillary(micro_unit_file)
+    !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_IMPORT_CAPILLARY" :: IMPORT_CAPILLARY
+    
+    character(len=MAX_FILENAME_LEN),intent(in) :: micro_unit_file
+    !local variables
+    integer :: count_units,ios,ne,nunit
+    real(dp) :: Pin,Pout,Qtot,temp(7),TOTAL_CAP_VOL,TOTAL_SHEET_SA,TT_TOTAL
+    
+    character(len=60) :: sub_name
+    
+    sub_name = 'import_capillary'
+    call enter_exit(sub_name,1)
+    
+    open(10, file = micro_unit_file, status='old')
+    
+    ! ios is negative if an end of record condition is encountered or if
+    ! an endfile condition was detected.  It is positive if an error was
+    ! detected.  ios is zero otherwise.
+    ios = 0
+    count_units = 0
+    do while (ios == 0)
+       read(10, '(I6,X,5(F9.2,X),2(F8.5,X),F10.2,X,F8.4,X,F10.4,X,F10.3,X,F8.4,X,F9.4,X)', iostat=ios)  &
+            ne,temp(1:3),Pin,Pout,temp(4),Qtot,temp(5),TOTAL_CAP_VOL,TOTAL_SHEET_SA,TT_TOTAL,temp(6:7)
+       if(ios == 0)then
+          ! record the unit values for mean pressure, transit time, surface area.
+          ! ne is the 'linker' element, so nunit is for its parent element
+          nunit = int(elem_field(ne_unit,elem_cnct(-1,1,ne))) 
+          unit_field(nu_blood_press,nunit) = (Pin+Pout)/2.0_dp
+          unit_field(nu_tt,nunit) = TT_TOTAL
+          unit_field(nu_sa,nunit) = TOTAL_SHEET_SA
+          count_units = count_units + 1
+       endif
+    enddo
+
+    if(count_units.ne.num_units)then
+       write(*,'(''WARNING: the number of capillary units ('',i6,'') does not match the geometric model ('',i6,'')'')') &
+            count_units,num_units
+    endif
+
+    close(10)
+    
+   call enter_exit(sub_name,2)
+   
+ end subroutine import_capillary
 !
 !##############################################################################
 !
@@ -136,4 +187,89 @@ subroutine import_exelemfield(FLOWFILE,field_no)
     call enter_exit(sub_name,2)
  end subroutine import_exelemfield
 
+!
+!##############################################################################
+!
+!>*import_terminal:* This subroutine reads in the content of an exnode field file
+ subroutine import_terminal(EXFILE)
+ !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_IMPORT_TERMINAL" :: IMPORT_TERMINAL
+
+   character(len=MAX_FILENAME_LEN),intent(in) :: EXFILE
+   !local variables
+   integer :: count_units,field_label(10),i,ibeg,iend,ierror,ne,nunit,n_fields
+   real(dp) :: rtemp
+   character(LEN=132) :: ctemp,label
+   character(len=300) :: readfile
+
+   character(len=60) :: sub_name
+
+   sub_name = 'import_terminal'
+   call enter_exit(sub_name,1)
+
+   if(index(EXFILE, ".exnode")> 0) then !full filename is given
+      readfile = EXFILE
+   else ! need to append the correct filename extension
+      readfile = trim(EXFILE)//'.exnode'
+   endif
+   
+   open(10, file=readfile, status='old')
+
+   n_fields = 0
+   ierror = 0
+
+   read_field_labels : do
+      read(10, fmt="(a)", iostat=ierror) ctemp
+      if(index(ctemp, "Node:")> 0) exit read_field_labels
+      if(index(ctemp, ") ")> 0) then
+         ibeg = index(ctemp, ") ")+1 ! beginning of label
+         iend = index(ctemp, ",")-1  ! end of label
+         label = adjustl(ctemp(ibeg:iend))
+         n_fields = n_fields + 1
+         field_label(n_fields) = 0
+         
+         select case(label)
+         case('terminal_element')
+         case('pleural_pressure')
+            field_label(n_fields) = nu_pe
+         case('tidal_volume')
+            field_label(n_fields) = nu_vt
+         case('max_Pe')
+            field_label(n_fields) = nu_Pe_max
+         case('min_Pe')
+            field_label(n_fields) = nu_Pe_min
+         end select
+      endif
+         
+   end do read_field_labels
+
+   nunit = 0
+   ierror = 0
+   
+   do while (ierror == 0)
+      if(index(ctemp, "Node:")> 0) then
+         
+         do i = 1,3  ! read coordinates; not used 
+            read(10, *, iostat=ierror) ctemp
+         enddo
+         read(10, *, iostat=ierror) ctemp ! read element: not used but could be
+
+         nunit = nunit + 1
+         
+         do i = 3,n_fields
+            read(10, *, iostat=ierror) rtemp
+            if(field_label(i).gt.0)then
+               unit_field(field_label(i),nunit) = rtemp
+            endif
+         enddo
+         
+         read(10, *, iostat=ierror) ctemp
+      endif
+   enddo
+      
+   close(10)
+   
+   call enter_exit(sub_name,2)
+   
+ end subroutine import_terminal
+ 
 end module imports
diff --git a/src/lib/indices.f90 b/src/lib/indices.f90
index 1cf75ba8..ed8a31fd 100644
--- a/src/lib/indices.f90
+++ b/src/lib/indices.f90
@@ -12,8 +12,10 @@
 
 module indices
 
+  use arrays
   use diagnostics
   use other_consts
+  use precision
   
   implicit none
   
@@ -28,11 +30,12 @@ module indices
        ne_resist=0,ne_t_resist=0,ne_Vdot=0,ne_Vdot0=0,ne_a_A=0,&
        ne_dvdt=0,ne_radius_in=0,ne_radius_in0=0,&
        ne_radius_out=0,ne_radius_out0=0,ne_group=0,ne_Qdot=0, &
-       ne_vd_bel=0, ne_vol_bel=0
+       ne_vd_bel=0, ne_vol_bel=0, ne_unit=0
   ! indices for unit_field
   integer :: num_nu,nu_vol=0,nu_comp=0,nu_conc2=0,nu_Vdot0=0,nu_Vdot1=0, &
        nu_Vdot2=0,nu_dpdt=0,nu_pe=0,nu_vt=0,nu_air_press=0,nu_conc1=0,nu_vent=0,&
-       nu_vd=0,nu_perf=0,nu_blood_press=0
+       nu_vd=0,nu_perf=0,nu_blood_press=0,nu_sa = 0,nu_tt = 0,nu_Pe_max=0, &
+       nu_Pe_min=0,nu_flux=0,nu_intsat=0,nu_av_flux=0,nu_lymphflow=0,nu_time=0
   !indices for gas exchange field
   ! indices for gasex_field
   integer,parameter :: num_gx = 12
@@ -60,12 +63,13 @@ module indices
        ne_resist,ne_t_resist,ne_Vdot,ne_Vdot0,ne_a_A,&
        ne_dvdt,ne_radius_in,ne_radius_in0,ne_radius_out,&
        ne_radius_out0,ne_group,ne_Qdot, &
-       ne_vd_bel, ne_vol_bel
+       ne_vd_bel, ne_vol_bel,ne_unit
   
   public num_nu,nu_vol,nu_comp, nu_conc2,nu_Vdot0,nu_Vdot1, &
        nu_Vdot2,nu_dpdt,nu_pe,nu_vt,nu_air_press,&
-       nu_conc1,nu_vent,nu_vd,&
-       nu_perf,nu_blood_press
+       nu_conc1,nu_vent,nu_vd,nu_flux,nu_intsat, &
+       nu_perf,nu_blood_press,nu_sa,nu_tt,nu_Pe_max,nu_Pe_min, &
+       nu_av_flux,nu_lymphflow,nu_time
   
   public num_gx, ng_p_alv_o2,ng_p_alv_co2,ng_p_ven_o2,ng_p_ven_co2, &
        ng_p_cap_o2, ng_p_cap_co2,ng_source_o2,ng_source_co2, &
@@ -76,8 +80,8 @@ module indices
   
   !Interfaces
   private
-  public define_problem_type,ventilation_indices, perfusion_indices, get_ne_radius, get_nj_conc1, &
-       growing_indices
+  public define_problem_type, lymphatic_indices,ventilation_indices, perfusion_indices, &
+       get_ne_radius, get_nj_conc1, growing_indices
   
 contains
   
@@ -100,6 +104,10 @@ subroutine define_problem_type(PROBLEM_TYPE)
     case ('gas_transfer')
        print *, 'You are solving a gas transfer model, setting up indices'
        call exchange_indices
+    case ('lymphatic_transport')
+       print *, 'You are solving a lymphatic transport model, setting up indices'
+       call lymphatic_indices
+       call update_units
     case ('perfusion')
        print *, 'You are solving a static perfusion model, setting up indices'
        call perfusion_indices
@@ -228,7 +236,7 @@ subroutine ventilation_indices
     ne_vd_bel = 9
     ne_vol_bel = 10
     ! indices for unit_field
-    num_nu=10
+    num_nu=12
     nu_vol=1
     nu_comp=2
     nu_Vdot0=3
@@ -239,9 +247,39 @@ subroutine ventilation_indices
     nu_vt=8
     nu_air_press=9
     nu_vent=10
+    nu_Pe_max = 11
+    nu_Pe_min = 12
     call enter_exit(sub_name,2)
   end subroutine ventilation_indices
 
+!!!#############################################################################
+  !> Lymphatic indices
+  subroutine lymphatic_indices
+    !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_LYMPHATIC_INDICES" :: LYMPHATIC_INDICES
+    
+    character(len=60) :: sub_name
+    
+    sub_name = 'lymphatic_indices'
+    call enter_exit(sub_name,1)
+    
+    ! indices for unit_field
+    num_nu = 11
+    nu_perf = 1
+    nu_blood_press = 2
+    nu_sa = 3
+    nu_tt = 4
+    nu_Pe_max = 5
+    nu_Pe_min = 6
+    nu_flux = 7
+    nu_intsat = 8
+    nu_av_flux = 9
+    nu_lymphflow = 10
+    nu_time = 11
+    
+    call enter_exit(sub_name,2)
+    
+  end subroutine lymphatic_indices
+
 !!!#############################################################################
 
   subroutine growing_indices
@@ -288,7 +326,7 @@ subroutine perfusion_indices
     num_nj=1
     nj_bv_press=1 !pressure in blood vessel
     ! indices for elem_field
-    num_ne=9
+    num_ne=10
     ne_radius=1 !strained average radius over whole element
     ne_radius_in=2 !strained radius into an element
     ne_radius_out=3 !strained radius out of an element
@@ -298,14 +336,40 @@ subroutine perfusion_indices
     ne_Qdot=7 !flow in an element
     ne_resist=8 !resistance of a blood vessel
     ne_group=9!Groups vessels into arteries (field=0), capillaries (field=1) and veins(field=2)
+    ne_unit=10 !store the unit number for terminal element
     !indices for units
-    num_nu=2
+    num_nu=4
     nu_perf=1
     nu_blood_press=2
+    nu_sa = 3
+    nu_tt = 4
     
     call enter_exit(sub_name,2)
   end subroutine perfusion_indices
+
   
+  subroutine update_units
+    !*update_units:* reallocates unit_field following a change in problem type
+
+    integer :: nu,num_nu_old,nunits
+    real(dp),allocatable :: unit_temp(:,:)
+
+    if(allocated(unit_field))then
+       num_nu_old = size(unit_field,1)
+       allocate(unit_temp(num_nu_old,nunits))
+       unit_temp = unit_field
+       deallocate(unit_field)
+       allocate(unit_field(num_nu,num_units))
+       unit_field = 0.0_dp
+       do nu = 1,num_nu_old
+          unit_field(nu,:) = unit_temp(nu,:)
+       enddo
+       deallocate(unit_temp)
+    endif
+
+  end subroutine update_units
+    
+    
   function get_ne_radius() result(res)
     
     implicit none
diff --git a/src/lib/lymphatics.f90 b/src/lib/lymphatics.f90
new file mode 100644
index 00000000..960733fb
--- /dev/null
+++ b/src/lib/lymphatics.f90
@@ -0,0 +1,419 @@
+module lymphatics
+  !*Brief Description:* This module contains all lymphatic-specific subroutines
+  !
+  !*LICENSE:*
+  !
+  !Test Test Test
+  !
+  !*Full Description:*
+  !
+  !This module contains code for pulmonary fluid flux within the alveolo-capillary network,
+  !and lymph transport through lymphatic collecting vessels.
+
+  use arrays
+  use diagnostics
+  use indices
+  use other_consts
+  use precision ! sets dp for precision
+
+  implicit none
+
+  !Module parameters
+
+  !Module types
+
+  !Module variables
+
+  !Interfaces
+!  private
+  public alveolar_capillary_flux
+  public lymphatic_transport
+contains
+
+!!!#############################################################################
+
+  subroutine alveolar_capillary_flux(ne,write_out)
+    !*alveolar_capillary_flux:* calculate fluid flux from blood to interstitium
+    !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_ALVEOLAR_CAPILLARY_FLUX" :: ALVEOLAR_CAPILLARY_FLUX
+
+    use other_consts,only: pi
+
+    integer,intent(in) :: ne
+    logical,intent(in) :: write_out
+
+    ! Baseline value parameters (eventually will be user-defined?)
+    integer,parameter :: sex = 1 ! 0 = male, 1 = female
+    !sex only determines the weight and therefore size of the lung. Should be updated based on CT
+
+    ! Capillary parameters
+    real(dp),parameter :: capillary_conductivity = 4.41335e-8 !ml/s/Pa  obtained from Parker (6e-8 cm H2O)
+    !real(dp),parameter :: capillary_conductivity = 9.26e-8 !ml/s/Pa  obtained from Parker (6e-8 cm H2O)
+
+    ! Osmotic pressure parameters
+    real(dp),parameter :: capillary_molar_conc = 0.0010250_dp ! this number is currenlty g/L, likely needs to change for osmolar to work...
+    real(dp),parameter :: IGC = 62.3630_dp ! ideal gas constant in mmHg.L.mol-1.K-1
+    real(dp),parameter :: T = 310.0_dp ! temperature in Kelvin - based on 37C blood temp
+
+    ! Simulation parameters
+    real(dp),parameter :: breathing_rate = 15.0_dp !constant but should be imported directly from ventilation model
+
+
+    ! Local variables
+    integer :: i,liflowcount,nunit,n_timesteps,printcount
+    real(dp) :: alveolar_volume,breathing_function,capillary_flow,capillary_vps,capillary_osmotic,capillary_osm_n, &
+         capillary_volume,capillary_volume_raw,cap_osm_conc,diffusion,dt,excess,flux_a,flux_b,flux_c,gas_diffusion_restriction, &
+         initial_lymphatic_flow,initial_lymphatic_pressure,initial_lymphatic_surface_area,initial_lymphatic_volume, &
+         initial_lymph_conc,initial_osm_n,interstitial_capacity,interstitial_capacity_a,interstitial_capacity_b, &
+         interstitial_osmotic,interstitial_pressure_a,interstitial_pressure_b,interstitial_saturation,interstitial_volume, &
+         int_osm_conc,int_osm_n,interstitial_volume_a,interstitial_volume_b,lung_mass,lymphatic_conductivity,max_Pe, &
+         min_Pe,net_flux,fluctuation,mx_pe,mn_pe,intPmax,intPmin,lymphPmax,lymphPmin, &
+         open_capillaries,osm_flux,osm_n_flux,overflow,sumflux,sumuptake,test_time,time,time_period,time_sum, &
+         time_variable,total_flux,total_hydro_flux,total_osm_flux,capillary_pressure,transit_time,capillary_SA
+    real(dp) :: sat1,sat2,sat3,sat4,sat5
+
+    logical :: continue
+
+
+    character(len=60) :: sub_name
+
+    ! --------------------------------------------------------------------------
+
+    sub_name = 'alveolar_capillary_flux'
+    call enter_exit(sub_name,1)
+
+    ! get information for the unit fron unit_field
+    ! ne is the 'linker' element in the artery-capillary-vein model, so nunit is for the parent element
+    nunit = int(elem_field(ne_unit,elem_cnct(-1,1,ne)))
+    capillary_pressure = unit_field(nu_blood_press,nunit)/133.32239_dp !converted to mmHg
+    transit_time = unit_field(nu_tt,nunit)
+    capillary_SA = unit_field(nu_sa,nunit)
+    max_Pe = unit_field(nu_Pe_max,nunit)
+    min_Pe = unit_field(nu_Pe_min,nunit)
+
+
+    if(write_out)then
+       write(*,'('' Unit'',i8,'': Pblood='',f7.2,'' mmHg; TT='',f7.2,'' s; SA='',f8.2,'' mm^2; Pe range='',f6.2,'' mmHg'')') &
+            nunit,capillary_pressure,transit_time,capillary_SA,(max_Pe-min_Pe)/133.32239_dp
+    endif
+
+    ! Calculated values
+    lung_mass = abs(real((1-sex)*840.0_dp))+real(sex)*639.0_dp  ! g;female lung weight of 639g and male of 840g - should be updated from CT
+    open_capillaries = 1.0_dp/6.0_dp !based on open capillaries at rest. Should be solved for by perfusion model?
+    capillary_volume_raw = 213.00_dp  !Gehr 1978 based on having a body mass of 74 kg - should be unique to each person
+    capillary_volume = (capillary_volume_raw*open_capillaries)/real(num_units)
+    !Only used for the osmotic model at the moment (which isn't operational) can volume be obtained elsewhere?
+
+    ! interstitial values
+    ! interstitial_capacity == maximal volume before spillover into alveolar in mm^3 - based on 30ml.100g of fluid (Drake 2002)
+    interstitial_capacity = ((30.0_dp*(lung_mass/100.0_dp))/real(num_units))*1000.0_dp !based on lung mass which should be obtained from CT
+    interstitial_capacity_a = 0.005_dp*interstitial_capacity !arbitrarily sized - needs further studies on the capillary-lymph interface
+    interstitial_capacity_b = 0.995_dp*interstitial_capacity !arbitrarily sized - needs further studies on the capillary-lymph interface
+    interstitial_volume_a = 0.0_dp
+    interstitial_volume_b = 0.48_dp*interstitial_capacity !assumed to be around 48% saturated at rest
+    alveolar_volume = 0.0_dp !alveolar volume likely greater at rest, but is lost to respiration - further information needed to put in model
+    liflowcount = 0
+    initial_lymphatic_surface_area = capillary_SA !initial lymphatic SA assumed to be the same as capillary SA as no other indication.
+    !both initial lymphatic SA and initial lymphatic hydraulic conductivity make up the filtration coefficient and currently it is the
+    !hydraulic conductivity that is adjusted to compensate
+
+    ! initial lymphatic values
+    initial_lymphatic_volume = 0.0_dp
+
+    ! Osmotic pressures
+    i = 1
+    capillary_osmotic = real(i)*capillary_molar_conc*IGC*T
+    capillary_osm_n = 0.0_dp
+    interstitial_osmotic = 0.0_dp
+    int_osm_n = 0.0_dp
+    initial_osm_n = 0.0_dp
+    total_osm_flux = 0.0_dp
+
+    time = 0.0_dp
+    printcount = 0
+    total_hydro_flux = 0.0_dp
+
+    intPmax = -1.00_dp
+    intPmin = -8.00_dp
+    lymphPmax = 1.00_dp
+    lymphPmin = -8.00_dp
+    mx_pe = max_Pe/133.32239_dp
+    mn_pe = min_Pe/133.32239_dp
+    fluctuation = ((mx_pe-mn_pe)/2.0_dp)
+
+    if(write_out) write(*,'(''fluctuation '',f8.4)')fluctuation
+    ! dt or n_timesteps should be controlled by the user
+    n_timesteps = 96
+    dt = transit_time/real(n_timesteps)
+
+    breathing_function = (2.0_dp*pi)/(60.0_dp/breathing_rate)
+
+    continue = .true.
+
+    sat1 = 1.0_dp
+    sat2 = 2.0_dp
+    sat3 = 3.0_dp
+    sat4 = 4.0_dp
+    sat5 = 5.0_dp
+
+    !do while(time < lymphatic_properties%test_time)
+    do while (continue)
+       time_sum = dt
+       if(transit_time.le.0.001)then
+          continue =.false.
+       endif
+       do while(time_sum < transit_time)
+          interstitial_volume = interstitial_volume_a + interstitial_volume_b
+          interstitial_saturation = interstitial_volume / interstitial_capacity  ! saturation as a proportion of 0-100%
+          time_variable = time + time_sum
+
+          ! calculating flux from capillary into interstitium
+          interstitial_pressure_a = fluctuation * sin(time_variable * breathing_function) + &
+               (((intPmin-intPmax+(fluctuation*2.0_dp)) * (interstitial_volume_a / interstitial_capacity_a)**2.0_dp) + &
+               ((intPmin-intPmax+(fluctuation*2.0_dp))*(-2.0_dp)) * (interstitial_volume_a / interstitial_capacity_a) + &
+               (intPmin + fluctuation))
+          !arbitrarily defined mathematical relationship between interstitial volume and pressure: (same for a and b)
+          !interstitial pressure changes a lot at low volumes with a small volume change, but at high volumes
+          !a large volume change is needed to cause a small change in pressure
+          interstitial_pressure_b = fluctuation * sin(time_variable * breathing_function) + &
+               (((intPmin-intPmax+(fluctuation*2.0_dp)) * (interstitial_volume_b / interstitial_capacity_b)**2.0_dp) + &
+               ((intPmin-intPmax+(fluctuation*2.0_dp))*(-2.0_dp)) * (interstitial_volume_b / interstitial_capacity_b) + &
+               (intPmin + fluctuation))
+          !write(*,'(''Pint: '',f8.4)')interstitial_pressure_b
+          ! pressure determined from saturation equation based of literature (currently linear, but likely not)
+          if(capillary_pressure > interstitial_pressure_a)then
+             flux_a = 0.5_dp * (capillary_conductivity * capillary_SA * (capillary_pressure - &
+                  interstitial_pressure_a)) * dt
+          else
+             flux_a = 0.0_dp
+          endif
+          if(capillary_pressure > interstitial_pressure_b)then
+             flux_b = 0.5_dp * (capillary_conductivity * capillary_SA * (capillary_pressure - &
+                  interstitial_pressure_b)) * dt
+          else
+             flux_b = 0.0_dp
+          endif
+          flux_c = flux_a + flux_b
+          total_hydro_flux = total_hydro_flux + flux_c
+
+          if(interstitial_volume_a + flux_a > interstitial_capacity_a)then
+             excess = flux_a - (interstitial_capacity_a - interstitial_volume_a)
+             interstitial_volume_a = interstitial_capacity_a
+             alveolar_volume = alveolar_volume + 0.5_dp*excess
+
+             if((interstitial_volume_b + 0.5_dp*excess) > interstitial_capacity_b)then
+                overflow = 0.5_dp * excess - (interstitial_capacity_b - interstitial_volume_b)
+                alveolar_volume = alveolar_volume + overflow
+                interstitial_volume_b = interstitial_capacity_b
+             else
+                interstitial_volume_b = interstitial_volume_b + 0.5_dp*excess
+             endif
+          else
+             interstitial_volume_a = interstitial_volume_a + flux_a
+          endif
+
+          interstitial_volume_b = interstitial_volume_b + flux_b
+          diffusion = (((interstitial_volume_a/interstitial_capacity_a)-(interstitial_volume_b/interstitial_capacity_b))/ &
+               (200_dp)) * dt
+          interstitial_volume_b = interstitial_volume_b + diffusion
+          interstitial_volume_a = interstitial_volume_a - diffusion
+
+          ! Osmotic
+          int_osm_conc = int_osm_n/interstitial_volume
+          interstitial_osmotic = real(i)*int_osm_conc*IGC*T
+          osm_flux = (lymphatic_properties%reflection_coefficient * capillary_SA * &
+               (capillary_osmotic - interstitial_osmotic))* dt
+          cap_osm_conc = capillary_osm_n / capillary_volume
+          int_osm_conc = int_osm_n / interstitial_volume
+
+          if(capillary_osmotic > interstitial_osmotic)then
+             capillary_volume = capillary_volume - osm_flux
+             interstitial_volume_b = interstitial_volume_b + osm_flux
+          else
+             capillary_volume = capillary_volume + osm_flux
+             interstitial_volume_b = interstitial_volume_b - osm_flux
+          endif
+
+          net_flux = osm_flux + flux_a + flux_b
+          if(net_flux > 0.0_dp)then
+             osm_n_flux = net_flux * cap_osm_conc
+             int_osm_n = int_osm_n + osm_n_flux
+             !assuming capillary is constant and does not need updating
+          else
+             osm_n_flux = net_flux * int_osm_conc
+             int_osm_n = int_osm_n - osm_n_flux
+          endif
+
+          total_osm_flux = total_osm_flux + osm_flux
+
+          !calculating flux from interstitium to initial lymphatics
+          if(interstitial_volume_b/interstitial_capacity_b < 0.3_dp)then
+             lymphatic_conductivity = 1.48_dp * 4.41335e-8 !all calculated does as a function of capillary_conductivity
+          !no information on the size of pores or similar for lympatic conductivity so assumed to be similar to capillary.
+          else
+
+             lymphatic_conductivity = ((845.87_dp * (interstitial_volume_b / interstitial_capacity_b)**5.0_dp) + &
+                  (-2416.7_dp * (interstitial_volume_b / interstitial_capacity_b)**4.0_dp) + (2388.5_dp * &
+                  (interstitial_volume_b / interstitial_capacity_b)**3.0_dp) + (-922.24_dp * (interstitial_volume_b / &
+                  interstitial_capacity_b)**2.0_dp) + (125.85_dp * (interstitial_volume_b / interstitial_capacity_b)) &
+                  - 0.0067_dp)* 4.41335e-8 !(capillary_conductivity)
+          endif
+
+          initial_lymphatic_pressure = fluctuation * sin((time_variable * breathing_function) + pi/2.0_dp) + &
+               ((((lymphPmax-lymphPmin-(fluctuation*2.0_dp))* ((interstitial_volume_b / interstitial_capacity_b)**2.0_dp)) + &
+               (lymphPmin + fluctuation)))
+          !arbitrarily defined mathematical relationship to show that lymphatic pressure does not change much at low volumes with a
+          !large volume change, but at high volumes only a small volume change is needed to cause a large change in pressure
+          !write(*,'(''Plym: '',f8.4)')initial_lymphatic_pressure
+          if(interstitial_volume.le.0.0_dp)then
+             initial_lymphatic_flow = 0.0_dp
+             interstitial_volume = 0.0_dp
+          elseif (interstitial_pressure_b > initial_lymphatic_pressure)then
+             initial_lymphatic_flow = (lymphatic_conductivity * initial_lymphatic_surface_area * (interstitial_pressure_b - &
+                     initial_lymphatic_pressure)) * dt
+          else
+             initial_lymphatic_flow = 0.0_dp
+          endif
+
+          interstitial_volume_b = interstitial_volume_b - initial_lymphatic_flow
+          initial_lymphatic_volume = initial_lymphatic_volume + initial_lymphatic_flow
+
+          int_osm_conc = int_osm_n/interstitial_volume_b
+          liflowcount = liflowcount + initial_lymphatic_flow
+          initial_osm_n = initial_osm_n + (initial_lymphatic_flow*int_osm_conc)
+          int_osm_n = int_osm_n - (initial_lymphatic_flow*int_osm_conc)
+          if (initial_lymphatic_volume > 0.0_dp)then
+             initial_lymph_conc = initial_osm_n/initial_lymphatic_volume
+          else
+             initial_lymph_conc = 0.0_dp
+          endif
+
+          time_sum = time_sum + dt
+
+          total_flux = total_hydro_flux ! +total_osm_flux
+          sumuptake = sumuptake + initial_lymphatic_flow
+
+       enddo
+
+       time = time + transit_time
+
+       printcount = printcount + 1
+       if (printcount.eq.100)then
+          sat5 = sat4
+          sat4 = sat3
+          sat3 = sat2
+          sat2 = sat1
+          sat1 = interstitial_saturation
+          if(write_out)then
+             write(*,'(f12.2, e12.4, f9.4, e12.4, f11.4, f9.4, e12.4, 2(f11.4), f12.8,f12.4)') time_variable, &
+                  flux_c/transit_time*1000.0_dp,interstitial_volume,interstitial_volume_a,interstitial_volume_b, &
+                  100.0_dp*interstitial_saturation,initial_lymphatic_flow*1000.0_dp, &
+                  initial_lymphatic_volume,total_flux,alveolar_volume,interstitial_pressure_b
+          endif
+          printcount = 0
+          if(time.gt.200.0_dp*transit_time)then
+             if((abs(((sat1 + sat2 + sat3 + sat4 +sat5)/5.0_dp)-sat1)).le.0.000005_dp)then
+                continue = .false.
+             endif
+          endif
+       endif
+
+       if(time.gt.10000.0_dp*transit_time) continue = .false.
+
+    enddo !while
+
+    unit_field(nu_intsat,nunit) = interstitial_saturation
+    unit_field(nu_time,nunit) = time
+    unit_field(nu_av_flux,nunit) = total_flux/time
+    unit_field(nu_lymphflow,nunit) = initial_lymphatic_volume/time
+    !write(*,'('' T='',e12.3,'': intsat='',e12.6,'' %; flux='',e12.3,'' ul/s; avFlux='',e12.3,'' ul/s; lyFlo='',e12.6,'' ul/s'')') &
+    !     unit_field(nu_time,nunit),unit_field(nu_intsat,nunit),&
+    !     unit_field(nu_flux,nunit),unit_field(nu_av_flux,nunit),unit_field(nu_lymphflow,nunit)
+
+    call enter_exit(sub_name,2)
+
+  end subroutine alveolar_capillary_flux
+
+!!!#############################################################################
+
+  subroutine lymphatic_transport(filename)
+    !*lymphatic_transport:* whole system transport
+    !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_LYMPHATIC_TRANSPORT" :: LYMPHATIC_TRANSPORT
+
+    character(len=MAX_FILENAME_LEN), intent(in) :: filename
+    ! Local parameters
+    integer :: i,ne,ne_child,np,nunit
+    real(dp) :: time_to_run,time_0
+    character(len=300) :: writefile
+    character(len=60) :: sub_name
+    !real(dp) :: interstitial_saturation,interstitial_pressure_b,nu_av_flux,nu_lymphflow,nu_time
+    ! --------------------------------------------------------------------------
+
+    sub_name = 'lymphatic_transport'
+    call enter_exit(sub_name,1)
+
+    if(index(filename, ".oplymph")> 0) then !full filename is given
+       writefile = filename
+    else ! need to append the correct filename extension
+       writefile = trim(filename)//'.oplymph'
+    endif
+
+    open(10, file=writefile, status='replace')
+
+    call cpu_time(time_0)
+    do ne = 1,num_elems
+       if(elem_field(ne_group,ne).eq.1.0_dp)then!(elem_field(ne_group,ne)-1.0_dp).lt.TOLERANCE)then
+          nunit = int(elem_field(ne_unit,elem_cnct(-1,1,ne)))
+          call alveolar_capillary_flux(ne,.false.)
+          np = elem_nodes(2,ne)
+          write(10,'(i8,11(e14.5))') ne,node_xyz(1:3,np),unit_field(nu_av_flux,nunit),unit_field(nu_intsat,nunit), &
+               unit_field(nu_lymphflow,nunit),unit_field(nu_blood_press,nunit),unit_field(nu_tt,nunit), &
+               unit_field(nu_sa,nunit),unit_field(nu_Pe_max,nunit),unit_field(nu_Pe_min,nunit)
+       endif
+    enddo
+    time_to_run = time_0
+    call cpu_time(time_to_run)
+    time_to_run = time_to_run - time_0
+    write(*,*) 'run time=',time_to_run
+    close(10)
+
+    do ne = num_elems,1,-1
+       if(elem_field(ne_group,ne).eq.1.0_dp)then
+          nunit = int(elem_field(ne_unit,elem_cnct(-1,1,ne)))
+          elem_field(ne_radius_in0,ne) = unit_field(nu_flux,nunit)
+          elem_field(ne_radius_out0,ne) = unit_field(nu_intsat,nunit)
+       else if(elem_field(ne_group,ne).eq.0.0_dp)then ! artery
+          elem_field(ne_radius_in0,ne) = 0.0_dp
+          elem_field(ne_radius_out0,ne) = 0.0_dp
+          do i = 1,elem_cnct(1,0,ne) ! each child branch
+             ne_child = elem_cnct(1,i,ne)
+             elem_field(ne_radius_in0,ne) = elem_field(ne_radius_in0,ne) + elem_field(ne_radius_in0,ne_child)
+             elem_field(ne_radius_out0,ne) = elem_field(ne_radius_out0,ne) + elem_field(ne_radius_out0,ne_child)
+          enddo
+          elem_field(ne_radius_in0,ne) = elem_field(ne_radius_in0,ne)/real(elem_cnct(1,0,ne))
+          elem_field(ne_radius_out0,ne) = elem_field(ne_radius_out0,ne)/real(elem_cnct(1,0,ne))
+       endif
+    enddo
+
+    call enter_exit(sub_name,2)
+
+  end subroutine lymphatic_transport
+
+!!!#############################################################################
+
+end module lymphatics
+
+
+!FUTURE DIRECTIONS
+!input a constant to account for difference between current values and expected values
+     !Model appeared to be working within the range of the literature but is likely off by a factor of 1000 due to nl to ul conversion error.
+     !Need to check that outputted units are correct - most things are in ml and mmHg
+!Lymphatic network tree
+     !currently all lymph is returned to the circulation immediately, in reality it moves up a tree of lymphatics against a pressure gradient
+     !would require excessive modelling perhaps
+!impairment of gas diffusion caused by high interstitial saturation
+     !unclear at what level this would occur
+!alveolar flooding changes
+     !alveolar flooding should be able to move between adjacent compartments
+     !alveolar fluid should be removed via respiration naturally and therefore should always occur naturally at some low level
+!individuality needs to be added in line with the other modules
+     !currently operates only on preset male/female values
diff --git a/src/lib/ventilation.f90 b/src/lib/ventilation.f90
index 97e09793..ae040ace 100644
--- a/src/lib/ventilation.f90
+++ b/src/lib/ventilation.f90
@@ -94,6 +94,8 @@ subroutine evaluate_vent
     sum_tidal = 0.0_dp ! initialise the inspired and expired volumes
     sum_expid = 0.0_dp
     last_vol = 0.0_dp
+    unit_field(nu_Pe_max,:) = -1.0e6_dp
+    unit_field(nu_Pe_min,:) = 1.0e6_dp
 
 !!! set default values for the parameters that control the breathing simulation
 !!! these should be controlled by user input (showing hard-coded for now)
@@ -554,11 +556,17 @@ subroutine tissue_compliance(chest_wall_compliance,undef)
             *(lambda**2+1.0_dp)/lambda**4)
        unit_field(nu_comp,nunit) = undef/unit_field(nu_comp,nunit) ! V/P
        ! add the chest wall (proportionately) in parallel
-       unit_field(nu_comp,nunit) = 1.0_dp/(1.0_dp/unit_field(nu_comp,nunit)&
-            +1.0_dp/(chest_wall_compliance/dble(num_units)))
+!       unit_field(nu_comp,nunit) = 1.0_dp/(1.0_dp/unit_field(nu_comp,nunit)&
+!            +1.0_dp/(chest_wall_compliance/dble(num_units)))
        !estimate an elastic recoil pressure for the unit
        unit_field(nu_pe,nunit) = cc/2.0_dp*(3.0_dp*a+b)*(lambda**2.0_dp &
             -1.0_dp)*exp_term/lambda
+       ! set minimum and maximum Pe values
+       unit_field(nu_Pe_max,nunit) = max(unit_field(nu_Pe_max,nunit), &
+            unit_field(nu_pe,nunit))
+       unit_field(nu_Pe_min,nunit) = min(unit_field(nu_Pe_min,nunit), &
+            unit_field(nu_pe,nunit))
+       
     enddo !nunit
 
     call enter_exit(sub_name,2)