diff --git a/src/bindings/c/repro_exports.c b/src/bindings/c/repro_exports.c
index a4fa26f..c94fae2 100644
--- a/src/bindings/c/repro_exports.c
+++ b/src/bindings/c/repro_exports.c
@@ -11,6 +11,9 @@ 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_terminal_perfusion_c(const char *EXNODEFILE, int *EXNODEFILE_LEN, const char *name, int *name_len);
 void export_node_geometry_c(const char *EXNODEFILE, int *EXNODEFILE_LEN, const char *name, int *name_len);
+void export_1d_elem_geometry_grpd_c(const char *EXELEMFILE, int *EXELEMFILE_LEN, const char *name, int *name_len, const char *mesh_code, int *mesh_code_len);
+void export_1d_elem_field_grouped_c(int *ne_field, const char *EXELEMFILE, int *EXELEMFILE_LEN,
+                            const char *group_name, int *group_name_len, const char *field_name, int *field_name_len , const char *mesh_code, int *mesh_code_len);
 
 void export_1d_elem_field(int ne_field, const char *EXELEMFILE, const char *group_name, const char *field_name )
 {
@@ -20,7 +23,14 @@ void export_1d_elem_field(int ne_field, const char *EXELEMFILE, const char *grou
 
   export_1d_elem_field_c(&ne_field, EXELEMFILE, &filename_len, group_name, &group_name_len, field_name, &field_name_len);
 }
-
+void export_1d_elem_field_grouped(int ne_field, const char *EXELEMFILE, const char *group_name, const char *field_name , const char *mesh_code)
+{
+  int filename_len = strlen(EXELEMFILE);
+  int group_name_len = strlen(group_name);
+  int field_name_len = strlen(field_name);
+  int mesh_code_len = strlen(mesh_code);
+  export_1d_elem_field_grouped_c(&ne_field, EXELEMFILE, &filename_len, group_name, &group_name_len, field_name, &field_name_len, mesh_code, &mesh_code_len);
+}
 void export_1d_elem_geometry(const char *EXELEMFILE, const char *name)
 {
   int filename_len = strlen(EXELEMFILE);
@@ -28,7 +38,14 @@ void export_1d_elem_geometry(const char *EXELEMFILE, const char *name)
 
   export_1d_elem_geometry_c(EXELEMFILE, &filename_len, name, &name_len);
 }
+void export_1d_elem_geometry_grpd(const char *EXELEMFILE, const char *name, const char *mesh_code)
+{
+  int filename_len = strlen(EXELEMFILE);
+  int name_len = strlen(name);
+  int mesh_code_len = strlen(mesh_code);
 
+  export_1d_elem_geometry_grpd_c(EXELEMFILE, &filename_len, name, &name_len, mesh_code, &mesh_code_len);
+}
 void export_node_field(int nj_field, const char *EXNODEFIELD, const char *name, const char *field_name)
 {
   int filename_len = strlen(EXNODEFIELD);
diff --git a/src/bindings/c/repro_exports.f90 b/src/bindings/c/repro_exports.f90
index bd3817e..34b9246 100644
--- a/src/bindings/c/repro_exports.f90
+++ b/src/bindings/c/repro_exports.f90
@@ -135,5 +135,57 @@ subroutine export_node_field_c(nj_field, EXNODEFIELD, filename_len, name, name_l
 #endif
 
   end subroutine export_node_field_c
+!!!################################################################
+
+  subroutine export_1d_elem_geometry_grpd_c(EXELEMFILE, filename_len, name, name_len, mesh_code, mesh_code_len) &
+    bind(C, name="export_1d_elem_geometry_grpd_c")
+    use iso_c_binding, only: c_ptr
+    use utils_c, only: strncpy
+    use repro_exports, only: export_1d_elem_geometry_grpd
+    use other_consts, only: MAX_STRING_LEN, MAX_FILENAME_LEN
+    implicit none
+    integer,intent(in) :: filename_len, name_len, mesh_code_len
+    type(c_ptr), value, intent(in) :: EXELEMFILE, name, mesh_code
+    character(len=MAX_FILENAME_LEN) :: filename_f
+    character(len=MAX_STRING_LEN) :: name_f, mesh_code_f
+
+    call strncpy(filename_f, EXELEMFILE, filename_len)
+    call strncpy(name_f, name, name_len)
+    call strncpy(mesh_code_f, mesh_code, mesh_code_len)
+
+#if defined _WIN32 && defined __INTEL_COMPILER
+    call so_export_1d_elem_geometry_grpd(filename_f, name_f, mesh_code_f)
+#else
+    call export_1d_elem_geometry_grpd(filename_f, name_f, mesh_code_f)
+#endif
+
+  end subroutine export_1d_elem_geometry_grpd_c
+
+  subroutine export_1d_elem_field_grouped_c(ne_field, EXELEMFILE, filename_len, group_name, group_name_len, &
+    field_name, field_name_len, mesh_code, mesh_code_len) bind(C, name="export_1d_elem_field_grouped_c")
+    use iso_c_binding, only: c_ptr
+    use utils_c, only: strncpy
+    use repro_exports, only: export_1d_elem_field_grouped
+    use other_consts, only: MAX_STRING_LEN, MAX_FILENAME_LEN
+    implicit none
+    integer,intent(in) :: ne_field, filename_len, group_name_len, field_name_len, mesh_code_len
+    type(c_ptr), value, intent(in) :: EXELEMFILE, group_name, field_name, mesh_code
+    character(len=MAX_FILENAME_LEN) :: filename_f
+    character(len=MAX_STRING_LEN) :: group_name_f, field_name_f, mesh_code_f
+
+    call strncpy(filename_f, EXELEMFILE, filename_len)
+    call strncpy(group_name_f, group_name, group_name_len)
+    call strncpy(field_name_f, field_name, field_name_len)
+    call strncpy(mesh_code_f, mesh_code, mesh_code_len)
+
+#if defined _WIN32 && defined __INTEL_COMPILER
+    call so_export_1d_elem_field_grouped(ne_field, filename_f, group_name_f, field_name_f, mesh_code_f)
+#else
+    call export_1d_elem_field_grouped(ne_field, filename_f, group_name_f, field_name_f, mesh_code_f)
+#endif
+
+  end subroutine export_1d_elem_field_grouped_c
+!!!############################################################################
+
 
 end module repro_exports_c
diff --git a/src/bindings/c/repro_exports.h b/src/bindings/c/repro_exports.h
index 2c2d2bc..bde4b92 100644
--- a/src/bindings/c/repro_exports.h
+++ b/src/bindings/c/repro_exports.h
@@ -9,5 +9,6 @@ SHO_PUBLIC void export_1d_elem_geometry(const char *EXELEMFILE, const char *name
 SHO_PUBLIC void export_node_field(int nj_field, const char *EXNODEFIELD, const char *name, const char *field_name);
 SHO_PUBLIC void export_terminal_perfusion(const char *EXNODEFILE, const char *name);
 SHO_PUBLIC void export_node_geometry(const char *EXNODEFILE, const char *name);
-
+SHO_PUBLIC void export_1d_elem_geometry_grpd(const char *EXELEMFILE, const char *name, const char *mesh_code);
+SHO_PUBLIC void export_1d_elem_field_grouped(int ne_field, const char *EXELEMFILE, const char *group_name, const char *field_name, const char *mesh_code );
 #endif /* REPROSIM_REPRO_EXPORTS_H */
diff --git a/src/lib/geometry.f90 b/src/lib/geometry.f90
index f7728fb..608c889 100644
--- a/src/lib/geometry.f90
+++ b/src/lib/geometry.f90
@@ -831,7 +831,7 @@ subroutine calc_capillary_unit_length
     seg_length=int_length/dble(num_convolutes) !lengh of each intermediate villous segment
     viscosity=0.33600e-02_dp !Pa.s !viscosity: fluid viscosity
     cap_unit_radius = 0.03_dp
-    cap_resistance=(8.0_dp*viscosity*cap_length)/(PI*cap_radius**4.0_dp)/dble(num_parallel) !resistance of each capillary convolute segment (6 capillaries in parallel)
+    cap_resistance=(8.0_dp*"viscosity"*cap_length)/(PI*cap_radius**4.0_dp)/dble(num_parallel) !resistance of each capillary convolute segment (6 capillaries in parallel)
 
     terminal_resistance = 0
 
@@ -1504,7 +1504,7 @@ subroutine define_rad_from_file(FIELDFILE,order_system, s_ratio)
           read(unit=10, fmt="(a)", iostat=ierror) ctemp1
     	  if(index(ctemp1, "value")> 0) then
                 call get_final_real(ctemp1,radius)
-                node_radius(np)=radius   
+                node_radius(np)=radius
           endif
        endif !index
        if(np.ge.radii_num_nodes) exit read_a_node
@@ -1801,7 +1801,7 @@ subroutine element_connectivity_1d()
     character(len=60) :: sub_name
     integer :: orphan_nodes(num_nodes)
     integer :: diagnostics_level
-
+    logical :: TTTS_anast
     sub_name = 'element_connectivity_1d'
     call enter_exit(sub_name,1)
     call get_diagnostics_level(diagnostics_level)
@@ -1810,7 +1810,7 @@ subroutine element_connectivity_1d()
     ! elems_at_node(node np,0)= total number of elements connected to this node
     ! elems_at_node(node np, index of each connected element starting at 1) = connected element
     elems_at_node = 0 !initialise
-
+    TTTS_anast = .TRUE.
     DO ne=1,num_elems
        DO nn=1,2
           np=elem_nodes(nn,ne)
@@ -1853,7 +1853,10 @@ subroutine element_connectivity_1d()
        IF(NNT == 2) THEN !1d
           np1=elem_nodes(1,ne) !first local node
           np2=elem_nodes(2,ne) !second local node
+          IF ((elems_at_node(np2,0)).EQ.3)THEN
 
+
+          END IF
           DO noelem=1,elems_at_node(np2,0) !for each element connected to node np2
              ne2=elems_at_node(np2,noelem) !get the element number connected to node np2
              IF(ne2 /= ne)THEN !if element connected to node np2 is not the current element ne
@@ -2494,7 +2497,7 @@ subroutine get_final_integer(string,num)
     integer :: ibeg,iend,nsign,ntemp
     character :: sub_string*(40)
 
-    iend=len(string) !get the length of the string
+    iend=len(string) !get the length of the stringelem_ordrs
     ibeg=index(string,":")+1 !get location of integer in string, follows ":"
     sub_string = adjustl(string(ibeg:iend)) ! get the characters beyond ":"
     iend=len(sub_string) !length of the sub-string
diff --git a/src/lib/repro_exports.f90 b/src/lib/repro_exports.f90
index 49763b4..4855f08 100644
--- a/src/lib/repro_exports.f90
+++ b/src/lib/repro_exports.f90
@@ -5,9 +5,8 @@ module repro_exports
   implicit none
 
   private
-  public export_1d_elem_geometry,export_node_geometry,export_node_field,&
-       export_terminal_perfusion,&
-       export_1d_elem_field
+  public export_1d_elem_geometry,export_1d_elem_geometry_grpd,export_1d_elem_field_grouped,export_node_geometry,&
+    export_node_field,export_terminal_perfusion,export_1d_elem_field
 
 contains
 !!!################################################################
@@ -168,6 +167,7 @@ subroutine export_node_geometry(EXNODEFILE, name)
              enddo
           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))
@@ -295,6 +295,158 @@ subroutine export_node_field(nj_field, EXNODEFIELD, name, field_name)
 
   end subroutine export_node_field
 
-!!! ###########################################################
+!!!################################################################
+
+  subroutine export_1d_elem_field_grouped(ne_field, EXELEMFILE, group_name, field_name , mesh_type)
+
+    use other_consts, only: MAX_FILENAME_LEN, MAX_STRING_LEN
+    use arrays,only: dp,elem_field,num_elems
+    use indices
+    implicit none
+  !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_EXPORT_1D_ELEM_FIELD_GROUPED" :: EXPORT_1D_ELEM_FIELD_GROUPED
+
+!!! Parameters
+    integer, intent(in) :: ne_field
+    real(dp) :: mesh_code
+    character(len=MAX_FILENAME_LEN), intent(in) :: EXELEMFILE
+    character(len=MAX_STRING_LEN), intent(in) :: field_name
+    character(len=MAX_STRING_LEN), intent(in) :: group_name
+    character(len=MAX_STRING_LEN), intent(in) :: mesh_type
+
+
+!!! Local Variables
+    integer :: len_end,ne
+    logical :: CHANGED
+
+    if(mesh_type.eq.'art') then
+       mesh_code = 0.0_dp
+    elseif(mesh_type.eq.'vein')then
+       mesh_code = 1.0_dp
+    elseif(mesh_type.eq.'cap')then
+       mesh_code = 2.0_dp
+    elseif(mesh_type.eq.'anast')then
+       mesh_code = 3.0_dp
+
+    endif
+    open(10, file=EXELEMFILE, status='replace')
+
+    len_end=len_trim(group_name)
+    !**     write the group name
+    write(10,'( '' Group name: '',A)') group_name(:len_end)
+    !**         write the elements
+    write(10,'( '' Shape.  Dimension=1'' )')
+    CHANGED=.TRUE. !initialise to force output of element information
+    len_end=len_trim(field_name)
+    do ne=1,num_elems
+       if(ne>1) THEN
+          CHANGED=.FALSE.
+       endif
+       if(CHANGED)THEN
+          write(10,'( '' #Scale factor sets=0'' )')
+          write(10,'( '' #Nodes= 0'' )')
+          write(10,'( '' #Fields= 1'' )')
+          write(10,'( '' 1)'',A,'', field, rectangular cartesian, #Components=1'')')&
+               field_name(:len_end)
+          write(10,'( ''  '',A,''.  l.Lagrange, no modify, grid based.'')') &
+               field_name(:len_end)
+          write(10,'( ''  #xi1=1'')')
+       endif
+       if(elem_field(ne_group,ne).eq.mesh_code)then
+         write(10,'(1X,''Element: '',I12,'' 0 0'' )') ne
+         write(10,'(3X,''Values:'' )')
+         write(10,'(4X,2(1X,E12.5))') elem_field(ne_field,ne),elem_field(ne_field,ne)
+       endif
+    enddo !no_nelist (ne)
+    close(10)
+
+  end subroutine export_1d_elem_field_grouped
+
+
+!!!################################################################
+
+
+  subroutine export_1d_elem_geometry_grpd(EXELEMFILE, name, mesh_type)
+
+    use arrays,only: dp,elem_field,elem_nodes,num_elems
+    use other_consts, only: MAX_FILENAME_LEN, MAX_STRING_LEN
+    use indices
+    use diagnostics, only: enter_exit
+
+    implicit none
+  !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_EXPORT_1D_ELEM_GEOMETRY_GRPD" :: EXPORT_1D_ELEM_GEOMETRY_GRPD
+
+!!! Parameters
+    character(len=MAX_FILENAME_LEN), intent(in) :: EXELEMFILE
+    character(len=MAX_STRING_LEN), intent(in) :: name
+    character(len=MAX_STRING_LEN), intent(in) :: mesh_type
+
+!!! Local Variables
+    integer :: len_end,ne,nj,nn
+    real(dp) :: mesh_code
+    character(len=1) :: char1
+    logical :: CHANGED
+    character(len=60) :: sub_name
+
+    sub_name = 'export_1d_elem_geometry'
+    call enter_exit(sub_name,1)
+
+
+
+    if(mesh_type.eq.'art') then
+       mesh_code = 0.0_dp
+    elseif(mesh_type.eq.'vein')then
+       mesh_code = 1.0_dp
+    elseif(mesh_type.eq.'cap')then
+       mesh_code = 2.0_dp
+    elseif(mesh_type.eq.'anast')then
+       mesh_code = 3.0_dp
+
+    endif
+    open(10, file=EXELEMFILE, status='replace')
+    len_end=len_trim(name)
+    !**     write the group name
+    write(10,'( '' Group name: '',A)') name(:len_end)
+    !**         write the elements
+    write(10,'( '' Shape.  Dimension=1'' )')
+    CHANGED=.TRUE. !initialise to force output of element information
+    do ne=1,num_elems
+       if(ne>1) THEN
+          CHANGED=.FALSE.
+       endif
+       if(CHANGED)THEN
+          write(10,'( '' #Scale factor sets=1'' )')
+          write(10,'( ''   l.Lagrange, #Scale factors= 2'' )')
+          write(10,'( '' #Nodes= 2'' )')
+          write(10,'( '' #Fields= 1'' )')
+          write(10,'( '' 1) coordinates, coordinate, rectangular cartesian, #Components=3'')')
+          do nj=1,3
+             if(nj==1) char1='x'; if(nj==2) char1='y'; if(nj==3) char1='z';
+             write(10,'(''  '',A2,''.  l.Lagrange, no modify, standard node based.'')') char1
+             write(10,'( ''     #Nodes= 2'')')
+             do nn=1,2
+                write(10,'(''      '',I1,''.  #Values=1'')') nn
+                write(10,'(''       Value indices:      1 '')')
+                write(10,'(''       Scale factor indices:'',I4)') nn
+             enddo !nn
+          enddo !nj
+       endif
+       if(elem_field(ne_group,ne).eq.mesh_code)then
+
+           write(10,'(1X,''Element: '',I12,'' 0 0'' )') ne
+           !**               write the nodes
+           write(10,'(3X,''Nodes:'' )')
+           write(10,'(4X,2(1X,I12))') elem_nodes(1,ne),elem_nodes(2,ne)
+           !**                 write the scale factors
+           write(10,'(3X,''Scale factors:'' )')
+           write(10,'(4X,2(1X,E12.5))') 1.d0,1.d0
+
+       endif
+    enddo !no_nelist (ne)
+    close(10)
 
+    call enter_exit(sub_name,2)
+  end subroutine export_1d_elem_geometry_grpd
+
+
+!!!##########################################################################
 end module repro_exports
diff --git a/src/lib/repro_exportsSJ.f90 b/src/lib/repro_exportsSJ.f90
new file mode 100644
index 0000000..49763b4
--- /dev/null
+++ b/src/lib/repro_exportsSJ.f90
@@ -0,0 +1,300 @@
+!
+!*Description:* This module handles all export functions
+!
+module repro_exports
+  implicit none
+
+  private
+  public export_1d_elem_geometry,export_node_geometry,export_node_field,&
+       export_terminal_perfusion,&
+       export_1d_elem_field
+
+contains
+!!!################################################################
+
+  subroutine export_1d_elem_field(ne_field, EXELEMFILE, group_name, field_name )
+ 
+    use other_consts, only: MAX_FILENAME_LEN, MAX_STRING_LEN
+    use arrays,only: elem_field,num_elems
+    implicit none
+  !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_EXPORT_1D_ELEM_FIELD" :: EXPORT_1D_ELEM_FIELD
+  
+!!! Parameters
+    integer, intent(in) :: ne_field
+    character(len=MAX_FILENAME_LEN), intent(in) :: EXELEMFILE
+    character(len=MAX_STRING_LEN), intent(in) :: field_name
+    character(len=MAX_STRING_LEN), intent(in) :: group_name
+
+!!! Local Variables
+    integer :: len_end,ne
+    logical :: CHANGED
+
+    open(10, file=EXELEMFILE, status='replace')
+
+    len_end=len_trim(group_name)
+    !**     write the group name
+    write(10,'( '' Group name: '',A)') group_name(:len_end)
+    !**         write the elements
+    write(10,'( '' Shape.  Dimension=1'' )')
+    CHANGED=.TRUE. !initialise to force output of element information
+    len_end=len_trim(field_name)
+    do ne=1,num_elems
+       if(ne>1) THEN
+          CHANGED=.FALSE.
+       endif
+       if(CHANGED)THEN
+          write(10,'( '' #Scale factor sets=0'' )')
+          write(10,'( '' #Nodes= 0'' )')
+          write(10,'( '' #Fields= 1'' )')
+          write(10,'( '' 1)'',A,'', field, rectangular cartesian, #Components=1'')')&
+               field_name(:len_end)
+          write(10,'( ''  '',A,''.  l.Lagrange, no modify, grid based.'')') &
+               field_name(:len_end)
+          write(10,'( ''  #xi1=1'')')
+       endif
+
+       write(10,'(1X,''Element: '',I12,'' 0 0'' )') ne
+       write(10,'(3X,''Values:'' )')
+       write(10,'(4X,2(1X,E12.5))') elem_field(ne_field,ne),elem_field(ne_field,ne)
+    enddo !no_nelist (ne)
+    close(10)
+
+  end subroutine export_1d_elem_field
+
+!!!############################################################################
+
+  subroutine export_1d_elem_geometry(EXELEMFILE, name)
+
+    use arrays,only: elem_nodes,num_elems
+    use other_consts, only: MAX_FILENAME_LEN, MAX_STRING_LEN
+    use diagnostics, only: enter_exit
+
+    implicit none
+  !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_EXPORT_1D_ELEM_GEOMETRY" :: EXPORT_1D_ELEM_GEOMETRY
+
+!!! Parameters
+    character(len=MAX_FILENAME_LEN), intent(in) :: EXELEMFILE
+    character(len=MAX_STRING_LEN), intent(in) :: name
+
+!!! Local Variables
+    integer :: len_end,ne,nj,nn
+    character(len=1) :: char1
+    logical :: CHANGED
+    character(len=60) :: sub_name
+
+    sub_name = 'export_1d_elem_geometry'
+    call enter_exit(sub_name,1)
+
+    open(10, file=EXELEMFILE, status='replace')
+    len_end=len_trim(name)
+    !**     write the group name
+    write(10,'( '' Group name: '',A)') name(:len_end)
+    !**         write the elements
+    write(10,'( '' Shape.  Dimension=1'' )')
+    CHANGED=.TRUE. !initialise to force output of element information
+    do ne=1,num_elems
+       if(ne>1) THEN
+          CHANGED=.FALSE.
+       endif
+       if(CHANGED)THEN
+          write(10,'( '' #Scale factor sets=1'' )')
+          write(10,'( ''   l.Lagrange, #Scale factors= 2'' )')
+          write(10,'( '' #Nodes= 2'' )')
+          write(10,'( '' #Fields= 1'' )')
+          write(10,'( '' 1) coordinates, coordinate, rectangular cartesian, #Components=3'')')
+          do nj=1,3
+             if(nj==1) char1='x'; if(nj==2) char1='y'; if(nj==3) char1='z';
+             write(10,'(''  '',A2,''.  l.Lagrange, no modify, standard node based.'')') char1
+             write(10,'( ''     #Nodes= 2'')')
+             do nn=1,2
+                write(10,'(''      '',I1,''.  #Values=1'')') nn
+                write(10,'(''       Value indices:      1 '')')
+                write(10,'(''       Scale factor indices:'',I4)') nn
+             enddo !nn
+          enddo !nj
+       endif
+       write(10,'(1X,''Element: '',I12,'' 0 0'' )') ne
+       !**               write the nodes
+       write(10,'(3X,''Nodes:'' )')
+       write(10,'(4X,2(1X,I12))') elem_nodes(1,ne),elem_nodes(2,ne)
+       !**                 write the scale factors
+       write(10,'(3X,''Scale factors:'' )')
+       write(10,'(4X,2(1X,E12.5))') 1.d0,1.d0
+    enddo !no_nelist (ne)
+    close(10)
+
+    call enter_exit(sub_name,2)
+  end subroutine export_1d_elem_geometry
+
+
+!!!##########################################################################
+
+  subroutine export_node_geometry(EXNODEFILE, name)
+  
+    use arrays,only: node_xyz,num_nodes
+    use other_consts, only: MAX_FILENAME_LEN, MAX_STRING_LEN
+    implicit none
+  !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_EXPORT_NODE_GEOMETRY" :: EXPORT_NODE_GEOMETRY
+
+!!! Parameters
+    character(len=MAX_FILENAME_LEN),intent(in) :: EXNODEFILE
+    character(len=MAX_STRING_LEN),intent(in) :: name
+
+!!! Local Variables
+    integer :: len_end,nj,np,np_last,VALUE_INDEX
+    logical :: FIRST_NODE
+
+    len_end=len_trim(name)
+    if(num_nodes.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 Geometry
+       do np=1,num_nodes
+          if(np.gt.1) np_last = np
+          !*** Write the field information
+          VALUE_INDEX=1
+          if(FIRST_NODE)THEN
+             write(10,'( '' #Fields=1'' )')
+             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="no") 1,0
+                write(10,'()')
+             enddo
+          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))
+          enddo !njj2
+          FIRST_NODE=.FALSE.
+          np_last=np
+       enddo !nolist (np)
+    endif !num_nodes
+    close(10)
+
+  end subroutine export_node_geometry
+
+!!!########################################################################
+
+  subroutine export_terminal_perfusion(EXNODEFILE, name)
+  
+    use arrays,only: elem_nodes,&
+         node_xyz,num_units,units,unit_field
+    use indices
+    use other_consts, only: MAX_FILENAME_LEN, MAX_STRING_LEN
+    implicit none
+  !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_EXPORT_TERMINAL_PERFUSION" :: EXPORT_TERMINAL_PERFUSION
+
+!!! 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=3'' )')
+             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
+             !perfusion
+             write(10,'('' 2) flow, field, rectangular cartesian, #Components=1'')')
+             write(10,'(2X,''1.  '')',advance="no")
+             write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") VALUE_INDEX,0
+             !Pressure
+             VALUE_INDEX=VALUE_INDEX+1
+             write(10,'('' 3) pressure, 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,F20.6))') (unit_field(nu_perf,NOLIST)) !flow
+           write(10,'(2X,4(1X,F12.6))') (unit_field(nu_blood_press,NOLIST)) !pressure
+          FIRST_NODE=.FALSE.
+          np_last=np
+       enddo !nolist (np)
+    endif !num_nodes
+    close(10)
+
+  end subroutine export_terminal_perfusion
+
+
+
+
+!!! #################################################################
+
+  subroutine export_node_field(nj_field, EXNODEFIELD, name, field_name)
+  
+    use arrays,only: node_field,num_nodes
+    use other_consts, only: MAX_FILENAME_LEN, MAX_STRING_LEN
+    implicit none
+  !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_EXPORT_NODE_FIELD" :: EXPORT_NODE_FIELD
+
+!!! Parameters
+    integer,intent(in) :: nj_field
+    character(len=MAX_FILENAME_LEN),intent(in) :: EXNODEFIELD
+    character(len=MAX_STRING_LEN),intent(in) :: field_name
+    character(len=MAX_STRING_LEN),intent(in) :: name
+
+!!! Local Variables
+    integer :: len_end,np
+    logical :: FIRST_NODE
+
+    open(10, file=EXNODEFIELD, status='replace')
+    !**     write the group name
+    len_end=len_trim(name)
+    write(10,'( '' Group name: '',A)') name(:len_end)
+    len_end=len_trim(field_name)
+    FIRST_NODE=.TRUE.
+    !*** the field as specified by user
+    do np=1,num_nodes
+       !*** Write the field information
+       if(FIRST_NODE)THEN
+          write(10,'( '' #Fields=1'' )')
+          write(10,'('' 1) '',A,'', field, rectangular cartesian, #Components=1'')') &
+               field_name(:len_end)
+          write(10,'(2X,''1.  '')',advance="no")
+          write(10,'(''Value index='',I1,'', #Derivatives='',I1)',advance="yes") 1,0
+       endif !FIRST_NODE
+       !***      write the node
+       write(10,'(1X,''Node: '',I12)') np
+       write(10,'(2X,2(1X,F12.6))') (node_field(nj_field,np))
+       FIRST_NODE=.FALSE.
+    enddo !num_nodes
+    close(10)
+
+  end subroutine export_node_field
+
+!!! ###########################################################
+
+end module repro_exports
diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt
new file mode 100644
index 0000000..1522fb0
--- /dev/null
+++ b/src/tests/CMakeLists.txt
@@ -0,0 +1,2 @@
+add_subdirectory(lib)
+
diff --git a/src/tests/lib/CMakeLists.txt b/src/tests/lib/CMakeLists.txt
new file mode 100644
index 0000000..e13934f
--- /dev/null
+++ b/src/tests/lib/CMakeLists.txt
@@ -0,0 +1,14 @@
+# Add all the files that make a single test, we could have multiple files testing
+# the same module.  Don't add test files into the same test that test different modules.
+# These are all .pf files.
+set(DIAGNOSTICS_TEST_SRCS
+    test_diagnostics.pf)
+
+# Make use of the pFUnit helper function to create a test.
+# Arguments    : - test_package_name: Name of the test package
+#                - test_sources     : List of pf-files to be compiled
+#                - extra_sources    : List of extra Fortran source code used for testing (if none, input empty string "")
+#                - extra_sources_c  : List of extra C/C++ source code used for testing (if none, input empty string "")
+add_pfunit_test(diagnostics_test ${DIAGNOSTICS_TEST_SRCS} "" "")
+# Link the test to the reprosim library target.
+target_link_libraries (diagnostics_test reprosim)
diff --git a/src/tests/lib/test_diagnostics.pf b/src/tests/lib/test_diagnostics.pf
new file mode 100644
index 0000000..55dbf56
--- /dev/null
+++ b/src/tests/lib/test_diagnostics.pf
@@ -0,0 +1,16 @@
+@test
+subroutine testSetDiagnostics()
+   use pfunit_mod
+   use diagnostics, only: get_diagnostics_on, set_diagnostics_on
+   implicit none
+
+   logical :: state
+
+   call get_diagnostics_on(state)
+   @assertFalse(state)
+   call set_diagnostics_on(.true.)
+   call get_diagnostics_on(state)
+   @assertTrue(state)
+
+end subroutine testSetDiagnostics
+