anchor.diff:              endif
             else!notfix depvar2
               grav=0.d0
  -            if(grav_type.eq.'on') then
  -              do nj=1,3                                  
  -                 grav=grav+density*grav_vect(nj)*9810.0_dp*(node_xyz(nj,elem_nodes(1,ne))-node_xyz(nj,elem_nodes(2,ne)))
  -              enddo
  -            endif
  +            do nj=1,3
  +              grav=grav+density*grav_vect(nj)*9810.0_dp*(node_xyz(nj,elem_nodes(1,ne))-node_xyz(nj,elem_nodes(2,ne)))
  +            enddo
                    if(nhs.EQ.1) THEN
                       RHS(nzz_row)=-prq_solution(depvar2,1)+grav
                    elseif(nhs.EQ.2) THEN
  
positive.code: !need to write solution to element/nodal fields for export
      call map_solution_to_mesh(prq_solution,depvar_at_elem,depvar_at_node,mesh_dof)
  
      deallocate (mesh_from_depvar, STAT = AllocateStatus)
      deallocate (depvar_at_elem, STAT = AllocateStatus)
      deallocate (depvar_at_node, STAT = AllocateStatus)
  
negative.code:         CONVERGED=.TRUE.
        else
  !Update vessel radii based on predicted pressures and then update resistance through tree
          call calc_press_area(grav_vect,KOUNT,depvar_at_node,prq_solution,&
             mesh_dof,vessel_type,elasticity_parameters,mechanics_type,mechanics_parameters)
          call calculate_resistance(viscosity,density,gamma)
  
  