anchor.diff:  
     subroutine tree_resistance(resistance)
       use indices
  -    use arrays,only: dp,num_elems,elem_cnct,elem_field,elems
  +    use arrays,only: dp,num_elems,elem_cnct,elem_field
       use diagnostics, only: enter_exit
       character(len=60) :: sub_name
   !local variables
       real(dp), intent(out) :: resistance
       real(dp) :: invres,elem_res(num_elems)
  -    integer :: num,num2,ne,ne2
  +    integer :: num2,ne,ne2
   
       sub_name = 'tree_resistance'
       call enter_exit(sub_name,1)
  
positive.code: !
  !*calc_press_area:* Calculates new radii based on pressure area relnships
  
  subroutine calc_press_area(KOUNT,variables_at_node,prq_solution,mesh_dof,ptrans,beta,Go_artery,Ptm_max)
  use indices
  use arrays,only: dp,num_nodes,num_elems,elem_field,elem_nodes,elem_cnct!,num_elems,elem_ordrs,elems,elem_nodes,num_nodes,elems_at_node,elem_field
  use diagnostics, only: enter_exit
  integer,intent(in) :: KOUNT,mesh_dof
  integer,intent(in) :: variables_at_node(num_nodes,0:2,2)
  real(dp) ::  prq_solution(mesh_dof,2)
  real(dp) :: ptrans,beta,Go_artery,Ptm_max
  
  !local variables
  integer :: np,de,ne,ny,numelem,nn
  double precision :: HEIGHT(3),G_PLEURAL,Ptm,R0,Pblood,Ppl
  
  character(len=60) :: sub_name
  sub_name = 'calc_press_area'
  call enter_exit(sub_name,1)
  
  IF(KOUNT.EQ.1)THEN !store initial, unstressed radius values
    DO ne=1,num_elems
      elem_field(ne_radius_in0,ne)=elem_field(ne_radius_in,ne)
      elem_field(ne_radius_out0,ne)=elem_field(ne_radius_in,ne)
  ENDDO !elems
  ENDIF
  
  DO ne=1,num_elems
    do nn=1,2
    np=elem_nodes(1,ne)
    ny=variables_at_node(np,0,1)
    Ppl=(ptrans*98.07d0/1000.0d0)!-G_PLEURAL !cmH2O->kPa
    Pblood=prq_solution(ny,1)/1000.0d0 ! Pa->kPa
    Ptm=Pblood+Ppl     ! kPa
    R0=elem_field(ne_radius_in0,ne)
  !...ARC: giving a maximum distension
    IF(Ptm.LT.Ptm_max.and.Go_artery.gt.0.d0)THEN
       if(nn.eq.1) elem_field(ne_radius_in,ne)=R0*((Ptm/Go_artery)+1.d0)**(1.d0/beta)
       if(nn.eq.2) elem_field(ne_radius_out,ne)=R0*((Ptm/Go_artery)+1.d0)**(1.d0/beta)
    ELSEIF(Ptm.lt.0.or.Go_artery.eq.0.d0)THEN
      IF(Ptm.lt.0)write(*,*) 'Transmural pressure < zero',ne,Ptm,Pblood,Ppl
       if(nn.eq.1) elem_field(ne_radius_in,ne)=R0
       if(nn.eq.2) elem_field(ne_radius_out,ne)=R0
    ELSE!ptm>ptmmax
      if(nn.eq.1) elem_field(ne_radius_in,ne)=R0*((Ptm_max/Go_artery)+1.d0)**(1.d0/beta)
      if(nn.eq.2) elem_field(ne_radius_in,ne)=R0*((Ptm_max/Go_artery)+1.d0)**(1.d0/beta)
    ENDIF
    enddo!nn
  enddo!ne
  !ny=variables_at_node(np,0,1)
  
  !...  KSB: Adding gravity-dependent pleural pressure
  !G_PLEURAL=0.d0    !gravitational force
  !DO nj=1,3
  !HEIGHT(nj)=node_xyz(nj,np)-node_xyz(nj,np_in)
  !G_PLEURAL=G_PLEURAL+PLEURAL_DENSITY*grav_vect(nj)*9810.d0*HEIGHT(nj) !kg
  !ENDDO
  !Ppl=(ptrans*98.07d0/1000.0d0)!-G_PLEURAL !cmH2O->kPa
  !Pblood=solution_prq(ny,1)/1000.0d0 ! Pa->kPa
  !Ptm=Pblood+Ppl     ! kPa
  
  
  !R0=elem_field(nj_radius0,np)
  !...ARC: giving a maximum distension
  !IF(Ptm.LT.Ptm_max.and.Go_vein.gt.0.d0)THEN
  !node_field(nj_radius,np)=R0*((Ptm/Go)+1.d0)**(1.d0/beta)
  !ELSEIF(Ptm.lt.0.or.Go.eq.0.d0)THEN
  !IF(Ptm.lt.0)write(*,*) 'Transmural pressure < zero',ne,Ptm,Pblood,Ppl
  !node_field(nj_radius,np)=R0
  !ELSE
  !node_field(nj_radius,np)=R0*((Ptm_max/Go)+1.d0)**(1.d0/beta)
  !write(*,*) 'Transmural pressure < zero',ne,Ptm,Pblood,Ppl,Go,Go_vein
  !ENDIF
  !Hypoxia removed
  !ENDDO     !nonode
  
  
  
  
negative.code: !!--TRANSFER SOLVER SOLUTIONS TO FULL SOLUTIONS
        ERR=0.0_dp
        no=0
        do depvar=1,mesh_dof
           if(.NOT.FIX(depvar)) THEN
              no=no+1
              prq_solution(depvar,2)=prq_solution(depvar,1) !temp storage of previous solution
              prq_solution(depvar,1)=solver_solution(no) !new pressure & flow solutions
              if(DABS(prq_solution(depvar,1)).GT.0.d-6)THEN
                 ERR=ERR+(prq_solution(depvar,2)-prq_solution(depvar,1))**2.d0/prq_solution(depvar,1)**2
              endif
           endif
        enddo !no2
  