anchor.diff:              node_xyz(1,np1))**2 + (node_xyz(2,np2) - &
               node_xyz(2,np1))**2 + (node_xyz(3,np2) - &
               node_xyz(3,np1))**2)
  -       !write(*,*),"avrad",ne,elem_field(ne_radius,ne)
  +       !print *,"avrad",ne,elem_field(ne_radius,ne)
          elem_field(ne_radius,ne)=(elem_field(ne_radius_in,ne)+elem_field(ne_radius_out,ne))/2
          ! element Poiseuille (laminar) resistance in units of Pa.s.mm-3        
          resistance = 8.d0*viscosity*elem_field(ne_length,ne)/ &
  
positive.code:        ! element turbulent resistance (flow in bifurcating tubes)
         !reynolds=DABS(elem_field(ne_flow,ne)*2.d0*density/ &
           !   (PI*elem_field(ne_radius,ne)*viscosity))
         zeta = MAX(1.d0,dsqrt(2.d0*elem_field(ne_radius,ne)* &
              reynolds/elem_field(ne_length,ne))*gamma)
         elem_field(ne_resist,ne) = resistance * zeta
         !write(*,*),"TESTING RESISTANCE",ne,elem_field(ne_resist,ne),elem_field(ne_radius,ne),elem_ordrs(2,ne)
      enddo 
  
      call enter_exit(sub_name,2)
  
negative.code: beta=1.0_dp
  Go_artery=6.67_dp
  Ptm_max=32.0_dp
  pleural_density=0.25_dp*0.1e-5_dp !kg/mm**3
  
  !!---------DESCRIPTION OF IMPORTANT PARAMETERS-----------
  !viscosity: fluid viscosity
  