anchor.diff:  C ======================================================================
   
         Integer Function Dbc(x, y, z, label, DATA, iSYS, Coef)
  -      implicit none
  -      Include 'assemble.fd'
  -
  -      Real*8  x, y, z, DATA(*), Coef(*)
  -      Integer label, iSYS(*)
  -
  -      iSYS(1) = 1
  -      iSYS(2) = 1
  -
  -      Dbc = BC_NEUMANN
  +          implicit none
  +          Include 'assemble.fd'
  +
  +          Real*8  x, y, z, DATA(*), Coef(*)
  +          Integer label, iSYS(*)
  +
  +          iSYS(1) = 1
  +          iSYS(2) = 1
  +
  +          Dbc = BC_NEUMANN
   
   c === \int\limits_{\partial\Omega} g_N d S = 1.
   c === where \Omega = [0, 1]^3 
   c === and g_N = const
   
  -      Coef(1) = 1d0 / 6d0
  -
  -      Return
  +          Coef(1) = 1d0 / 6d0
  +
  +          Return
         End
   
   C ======================================================================
  
positive.code:      &                 sle_right_hand_side,                              ! b
       &                 sle_solution)                                     ! x
  
        write (*,*) "LSolver has finished!"
  
  c ===   write results to VTK file
          call GMVscalarTet(sle_solution, 
  
negative.code:      &                    face_nodes,
       &                    face_material_labels) 
  
          write (*,*) "GMVscalarTet has finished!"
  
  c ===   I may be wrong about that      
        end program thermoconductivity
  