void 
wr_elem_result_exo(Exo_DB *exo, const char *filename, double ***vector, 
		   const int variable_index, const int time_step,
		   const double time_value, 
		   struct Results_Description *rd)
{
  int error, i;
  double local_time_value=time_value;
  /* static char *yo = "wr_elem_result_exo"; */

  /*
   * This file must already exist.
   */

  exo->cmode = EX_WRITE;

  exo->io_wordsize = 0;		/* query */
  exo->exoid = ex_open(filename, exo->cmode, &exo->comp_wordsize, 
		       &exo->io_wordsize, &exo->version);
  EH(exo->exoid, "ex_open");

#ifdef DEBUG
  fprintf(stderr, "\t\tfilename    = \"%s\"\n", filename);
  fprintf(stderr, "\t\tcomp_ws     = %d\n", exo->comp_wordsize);
  fprintf(stderr, "\t\tio_wordsize = %d\n", exo->io_wordsize);
#endif

  error = ex_put_time (exo->exoid, time_step, &local_time_value );
  EH(error, "ex_put_time");

  /* If the truth table has NOT been set up, this will be really slow... */

  for (i = 0; i < exo->num_elem_blocks; i++) {
    if (exo->elem_var_tab_exists == TRUE) {
      /* Only write out vals if this variable exists for the block */
      if (exo->elem_var_tab[i*rd->nev + variable_index] == 1) {
	error = ex_put_var(exo->exoid, time_step, EX_ELEM_BLOCK, variable_index+1,
				exo->eb_id[i], exo->eb_num_elems[i],
				vector[i][variable_index]);
	EH(error, "ex_put_var elem");
      }
    }
    else {
      /* write it anyway (not really recommended from a performance viewpoint) */
      error      = ex_put_var ( exo->exoid,
				time_step, EX_ELEM_BLOCK,
				variable_index+1, /* Convert to 1 based for
						     exodus */
				exo->eb_id[i],
				exo->eb_num_elems[i],
				vector[i][variable_index] );
      EH(error, "ex_put_var elem");
    }
  }

  error      = ex_close ( exo->exoid );
  EH(error, "ex_close");

  return;
}

static void
check_discontinuous_interp_type(PROBLEM_DESCRIPTION_STRUCT *curr_pd,
			        int var_type)
     
    /********************************************************************
     *
     * check_discontinuous_interp_type
     *
     * -> Test whether an interpolation specified in the input deck
     *    is consistent with a discontinuous interpolation at
     *    the interface
     ********************************************************************/
{
  int *v_ptr = curr_pd->v;
  int interp_type = curr_pd->i[var_type];
  if (v_ptr[var_type] & V_MATSPECIFIC) {
    switch (interp_type) {
    case I_NOTHING:
    case I_Q1:
    case I_Q2:
    case I_Q2_LSA:
    case I_H3:
    case I_S2:
    case I_B3:
    case I_Q3:
    case I_Q4:
    case I_SP:
	fprintf(stderr,"check_interpolation_discontinuous ERROR");
	fprintf(stderr," var type to be discontinuous in mat %s\n,",
		curr_pd->MaterialName);
	fprintf(stderr,"\tbut incompatible interp type specified: %d\n",
		interp_type);
	EH(-1,"incompatible interp type");
	break;
	    case I_G0:
    case I_G1:
    case I_P0:
    case I_P1:
    case I_Q1_D:
    case I_Q2_D:
    case I_Q2_D_LSA:
    case I_PQ1:
    case I_PQ2:
         break;
    default:
     	fprintf(stderr,"check_interpolation_discontinuous ERROR");
	fprintf(stderr,"unknown interpolation type\n");
	EH(-1,"unknown interpolation type");
    }
  }
}

static void
rearrange_mat_increasing(int var_type, int subvarIndex,
			 NODAL_VARS_STRUCT *nv_old,
			 NODAL_VARS_STRUCT *nv_new, int *rec_used)
    
      /*****************************************************************
       *
       * rearrange_mat_increasing():
       *
       *   Search for the record with the matching variable type and
       *   subvar index that  has the lowest MatID field, that is
       *   currently unused. It then copies that variable type from
       *   the old nodal variable structure to the new.
       *****************************************************************/
{
  int j, mn_min, var_rec, i, num;
  VARIABLE_DESCRIPTION_STRUCT *vd;

  /*
   *  Loop over the number of records in the old structure that
   *  have the given variable type.
   */
  num =  nv_old->Num_Var_Desc_Per_Type[var_type];
  if (var_type == MASS_FRACTION) {
    num =  nv_old->Num_Var_Desc_Per_Type[var_type] /
	   upd->Max_Num_Species_Eqn;
  }
  
  for (j = 0; j < num; j++) {
 
    /*
     * Search for the record with the matching variable type and
     * subvar index that
     * has the lowest MatID field, that is currently unused.
     */
    mn_min = INT_MAX;
    var_rec = -1;
    for (i = 0; i < nv_old->Num_Var_Desc; i++) {
      if (! rec_used[i]) {
	vd = nv_old->Var_Desc_List[i];
	if ((var_type == vd->Variable_Type) &&
	    (subvarIndex == (int) vd->Subvar_Index) &&
	    (mn_min > vd->MatID)) {
	  var_rec = i;
	  mn_min = vd->MatID;
	}
      } 
    }
    if (var_rec == -1) {
      EH(-1,"AUGH! we shouln't be here");
    }
    rec_used[var_rec] = 1;
      
    /*
     * Now fill in fields in the nodal vars struct
     */
    add_var_to_nv_struct(nv_old->Var_Desc_List[var_rec],
			 nv_new);
  }
}

int
gsum_Int(const int value)
    
    /********************************************************************
     *
     * gsum_Int
     *
     *   This routine will return the sum of a single integer
     *   distributed across the processors.
     *
     *
     *  Return
     *  -------
     *  Routine returns the sum.
     ********************************************************************/
{
#ifdef PARALLEL
  int out_buf, err;
  err = MPI_Allreduce((void *) &value, (void *) &out_buf, 1, MPI_INT,
		      MPI_SUM, MPI_COMM_WORLD);
  if (err != MPI_SUCCESS) {
    EH(-1, "gsum_Int: MPI_Allreduce returned an error");
  }
  return out_buf;
#else
  return value;
#endif
}

double
gavg_double(const double value)

    /********************************************************************
     *
     * gavg_double
     *
     *   This routine will find the average value of an input
     *   across all of the processors
     *
     *
     *  Return
     *  -------
     *  Routine returns the average value
     ********************************************************************/
{
#ifdef PARALLEL
  int err;
  double out_buf;
  err = MPI_Allreduce((void *) &value, (void *) &out_buf, 1, MPI_DOUBLE,
	     	      MPI_SUM, MPI_COMM_WORLD);
  if (err != MPI_SUCCESS) {
    EH(-1, "gavg_double: MPI_Allreduce returned an error");
  }
  return out_buf / Num_Proc;
#else
  return value;
#endif
}

void
retrieve_parameterS(double *lambda, /* PARAMETER VALUE */
		  double *x, 	 /* UNKNOWN VECTOR */
		  double *xdot,  /* UNKNOWN_DOT VECTOR */
		  int sens_type,/*  type of sensitivity variable(BC or MT) */
		  int sens_id,		/*  variable id BCID or MT# */
		  int sens_flt,		/* data float id or matl prop id */
		  int sens_flt2,	/* data float id for UM */
		  Comm_Ex *cx,	 /* array of communications structures */
		  Exo_DB *exo,	 /* ptr to the finite element mesh database */
		  Dpi *dpi)	 /* distributed processing information */
{

  int ic;
  int mn,mpr;
  int ibc, idf;
  
#ifdef DEBUG
  static const char yo[]="retrieve_parameterS";

  fprintf(stderr, "%s() begins...\n", yo);
#endif

  if (sens_type == 1) {

    ibc = sens_id;
    idf = sens_flt;

    retrieve_BC_parameter(lambda, ibc, idf, cx, exo, dpi);

  	}
  else
  if (sens_type == 2) {

    mn = sens_id;
    mpr = sens_flt;

      retrieve_MT_parameter(lambda, mn, mpr, cx, exo, dpi);
	}
  else
  if (sens_type == 3) {
 
     ibc = sens_id;
     idf = sens_flt;
 
      retrieve_AC_parameter(lambda, ibc, idf, cx, exo, dpi);
 	}
  else 
  if (sens_type == 4) {

    mn = sens_id;
    mpr = sens_flt;
    ic = sens_flt2;

      retrieve_UM_parameter(lambda, mn, mpr, ic, cx, exo, dpi);
  }

  else EH(-1, "Bad sens. parameter type!");
 
}/* END of routine retrieve_parameterS  */

int
variable_description_init(VARIABLE_DESCRIPTION_STRUCT **vd_hdl)

    /*******************************************************************
     *
     * variable_description_init:
     *
     *  Mallocs and initializes a variable description structure
     *  to its default state
     *
     *  Input
     *   *el_hdl => If NULL will allocate a structure. If nonnull it will
     *              assume that the structure has already been malloced.
     *  Return:
     *    Success:           CPC_SUCCESS
     *    Interface Failure: CPC_PUB_BAD
     *    *el_hdl   => address of the new structure that has just been
     *                 malloced.
     ******************************************************************/
{
  if (*vd_hdl == NULL) {
    *vd_hdl = smalloc(sizeof(VARIABLE_DESCRIPTION_STRUCT));
    if (*vd_hdl == NULL) EH(-1, "variable_description_init");
  }
  return variable_description_default(*vd_hdl);
}

/*
 * Checksum for SRECORD.  Add all the bytes from the record length field through the data, and take
 * the ones complement.  This includes the address field, which for SRECORDs can be 2 bytes, 3 bytes or
 * 4 bytes.  In this case, since the upper bytes of the address will be 0 for records of the smaller sizes,
 * just add all 4 bytes of the address in all cases.
 * 
 * The arguments are:
 * 
 *  buf		a pointer to the beginning of the data for the record
 *  length	the length of the data portion of the record
 *  chunk_len	the length of the record starting at the address and including the checksum
 *  chunk_addr	starting address for the data in the record
 *
 * Returns an unsigned char with the checksum
 */
static unsigned char srec_csum(unsigned char *buf, unsigned int length, int chunk_len, int chunk_addr) {
  int sum = chunk_len + (LO(chunk_addr)) + (HI(chunk_addr)) + (EX(chunk_addr)) + (EH(chunk_addr));
  unsigned int i;
  for(i = 0; i < length; i++) {
    sum += buf[i];
  }
  return ~sum & 0xff;
}

int 
get_nv_offset_idof(NODAL_VARS_STRUCT *nv, const int varType, 
		   const int idof, int subVarType,
		   VARIABLE_DESCRIPTION_STRUCT **vd_ptr)

    /**************************************************************
     *
     * get_nv_offset_idof():
     *
     *      This function returns the offset into the local solution
     * vector of the idof'th occurrence of variable type, varType.
     * MASS_FRACTION variables types also distinguish between 
     * subvartypes. It also returns the address of the variable
     * description structure pertaining to this unknown.
     **************************************************************/
{
  int i, index, offset = -1, nfound;
  int num = nv->Num_Var_Desc_Per_Type[varType];
  VARIABLE_DESCRIPTION_STRUCT *vd;
#ifdef DEBUG_HKM
  if (idof >= num) {
    EH(-1, "ERROR");
  }
#endif
  if (varType == MASS_FRACTION) {
    for (i = 0, nfound = 0; i < num; i++) {
      index = nv->Var_Type_Index[varType][i];
      vd = nv->Var_Desc_List[index];
      if (vd->Subvar_Index == subVarType) {
	if (nfound == idof) {
	  if (vd_ptr) *vd_ptr = vd;
	  offset = nv->Nodal_Offset[index];
	  return offset;
	} else {
	  nfound++;
	}
      }
    }
    EH(-1,"NOT FOUND");
  } else {
    index = nv->Var_Type_Index[varType][idof];
    if (vd_ptr) *vd_ptr = nv->Var_Desc_List[index];    
    offset = nv->Nodal_Offset[index] + subVarType;
  }
  return offset;
}

void
wr_global_result_exo( Exo_DB *exo, 
		      const char *filename, 
		      const int time_step, 
		      const int ngv, 
		      double u[] )
{
     /*****************************************************************
      * write_global_result_exo() 
      *     -- open/write/close EXODUS II db for all global values
      *
      * The output EXODUS II database contains the original model
      * information with some minor QA and info additions, with new 
      * global data written.
      * 
      ******************************************************************/
  int error;


  /* 
   * This capability is deactivated for parallel processing.
   * brkfix doesn't support global variables, when this
   * changes this restriction should be removed. TAB 3/2002 */

  if( u == NULL ) return ; /* Do nothing if this is NULL */

  exo->cmode = EX_WRITE;
  exo->io_wordsize = 0;		/* query */
  exo->exoid = ex_open(filename, exo->cmode, &exo->comp_wordsize, 
		       &exo->io_wordsize, &exo->version);
  if (exo->exoid < 0) {
    EH(-1,"wr_nodal_result_exo: could not open the output file");
  }

  error = ex_put_var( exo->exoid, time_step, EX_GLOBAL, 1, 0, ngv, u );

  EH(error, "ex_put_var glob_vars");

  error = ex_close( exo->exoid);


  return;
}

void
update_user_TP_parameter(double lambda, double *x, double *xdot,
                         double *x_AC, Comm_Ex *cx, Exo_DB *exo, Dpi *dpi)
/*
 * This function handles updates to the second (TP) parameter
 * when LOCA bifurcation tracking algorithms (turning point,
 * pitchfork, or Hopf are used. It works the same as the above
 * function, but here lambda is the TP parameter.
 *
 * The examples in the above function also apply identically to this
 * one, and for brevity are not repeated. The standard call to
 * do_user_update is also the same as above. The main difference is
 * the variables first_cp and first_tp are set up to indicate which
 * parameter is being updated and whether this is the first update
 * call for that parameter. The user needn't worry about this detail.
 */

{
    /* Actual function begins here (refer to previous example) */
    static int first_tp = 1;
    // int first_cp = -1;
    // int n = 0;
    // int Type, BCID, DFID, MTID, MPID, MDID;
    // double value;
    /* Declare any additional variables here */


    /* If using this function, comment out this line. */
    EH(-1, "No user TP continuation conditions entered!");


    /* Evaluate any intermediate quantities and/or assign constants here */


    /* Enter each continuation condition in sequence in this space */


    /* ID's */


    /* Value */


    /* Update call - copy from example */


    /* Done */
    first_tp = 0;
    return;

} /* END of routine update_user_TP_parameter */

void 
assign_species_desc_suffix(const int species_var_name, char * retn_string)

   /***************************************************************************
    *
    * assign_species_desc_suffix:
    *
    *  This function assigns a string to the return string, associated with the
    *  type of species variable being considered.
    *
    *  Input
    * --------
    *  species_var_name: Valid species type
    *                   (The valid species types are listed in rf_fem_const.h)
    * Output
    * --------
    *  retn_string:  Has to have an input length of at least 24
    ***************************************************************************/
{
  if (retn_string == NULL) {
   EH(-1, "assign_species_desc_suffix: bad interface\n");
  }
  switch (species_var_name) {
  case SPECIES_MASS_FRACTION:
      (void) strcpy(retn_string, "Mass Fraction");
      break;
  case SPECIES_MOLE_FRACTION:
      (void) strcpy(retn_string, "Mole Fraction");
      break;
  case SPECIES_VOL_FRACTION:
      (void) strcpy(retn_string, "Volume Fraction");
      break;
  case SPECIES_DENSITY:
      (void) strcpy(retn_string, "Density");
      break;
  case SPECIES_CONCENTRATION:
      (void) strcpy(retn_string, "Concentration");
      break;
  case SPECIES_CAP_PRESSURE:
      (void) strcpy(retn_string, "Capillary Pressure");
      break;
  case SPECIES_UNDEFINED_FORM:
      (void) strcpy(retn_string, "Species Unknown");
      break;
  default:
      (void) strcpy(retn_string, "Mass Fraction(default)");
      printf("assign_species_desc_suffix: WARNING unknown form of "
	     "species vars: %d\n", species_var_name);
  }
}

static void
define_dimension(const int unit,
		 const char *string,
		 const int value,
		 int *identifier)
{
  int err;
  char err_msg[MAX_CHAR_ERR_MSG];

  /*
   * Dimensions with value zero are problematic, so filter them out.
   */

  if ( value <= 0 )
    {
      *identifier = -1;
      return;
    }
#ifdef NETCDF_3  
  err  = nc_def_dim(unit, string, value, identifier);
  if ( err != NC_NOERR )
    {
      sprintf(err_msg, "nc_def_dim() on %s [<%d] id=%d", string, 
		   value, *identifier);
      EH(-1, err_msg);
    }
#endif
#ifdef NETCDF_2
  err  = ncdimdef(unit, string, value);
  sprintf(err_msg, "ncdimdef() on %s [<%d] rtn %d", string, 
		 value, err);
  EH(err, err_msg);
  *identifier = err;
#endif  

  return;
}

int
dof_lnode_var_type(const int n, const int Element_Type,
		   const int proc_node_num, const int var_type,
		   PROBLEM_DESCRIPTION_STRUCT *pd_ptr)

    /**********************************************************************
     *
     * dof_lnode_var_type:
     *
     *  This is a wrapper around dof_lnode_interp_type(). It calculates
     *  the number of degrees of freedom located at this local node number,
     *  that is actually interpolated on this element (i.e., active on this
     *  element)
     *  given the interpolation type, given the variable type and the current
     *  problem description structure. Within that structure it defines
     *  the interpolation to be used for that variable type on that
     *  element type.
     *  See the dof_lnode_interp_type() description for more information.
     *
     *  proc_node_num is needed to look up whether this node is on the
     *  edge of a domain.
     *
     *  So, if a variable is located at that node, but is not active for
     *  that element, the answer is zero.
     *
     *  The return value is equal to the number of degrees of freedom at
     *  a local node for a variable type corresponding to the current
     *  element, given the problem description structure.
     *********************************************************************/
{
  int retn, interp_type, edge;
  /*
   * Store the interpolation type for the input variable
   */
  interp_type =  pd_ptr->i[var_type];

  /*
   * Store whether this node is on an edge of the grid or not
   */
  edge = (int) Nodes[proc_node_num]->EDGE;

  /*
   * Now, given the local node number, n, the element type, and the
   * interpolation type, return the number of degrees of freedom
   */
  retn = dof_lnode_interp_type(n, Element_Type, interp_type, edge);
  EH(retn, "dof_lnode_var_type: ERROR");  
  return retn;
}

void
ReduceBcast_BOR(int *ivec, int length)

   /********************************************************************
    *
    * ReduceBcast_BOR()
    *
    *  Does a logical bitwize OR operation on a vector of ints
    *  distibuted across different processors.
    ********************************************************************/
{
#ifdef PARALLEL
  int k, err, *ivec_recv;
  if (length <= 0) {
    printf(" ReduceBcast_BOR Warning, length = %d\n", length);
    return;
  }
  if (ivec == NULL) {
    EH(-1, " ReduceBcast_BOR fatal Interface error");
  }
  ivec_recv = alloc_int_1(length, INT_NOINIT);
  err = MPI_Reduce(ivec, ivec_recv, length, MPI_INT, MPI_BOR, 0,
	           MPI_COMM_WORLD);
  if (err != MPI_SUCCESS) {
    EH(-1, " ReduceBcast_BOR fatal MPI Error");
  }
  err = MPI_Bcast(ivec_recv, V_LAST, MPI_INT, 0, MPI_COMM_WORLD);
  if (err != MPI_SUCCESS) {
    EH(-1, " ReduceBcast_BOR fatal MPI Error");
  }  
  for (k = 0; k < V_LAST; k++) {
    ivec[k] = ivec_recv[k];
  }
  safer_free((void **) &ivec_recv);
#endif
}

void 
wr_nodal_result_exo(Exo_DB *exo, char *filename, double vector[],
		    int variable_index, int time_step, double time_value)

     /*****************************************************************
      * write_nodal_result_exo() 
      *     -- open/write/close EXODUS II db for 1 nodal var at one
      *        time step.
      *
      * The output EXODUS II database contains the original model
      * information with some minor QA and info additions, with new 
      * nodal value solution data written.
      * 
      ******************************************************************/
{
  char err_msg[MAX_CHAR_IN_INPUT];
  int error;
  exo->cmode = EX_WRITE;
  exo->io_wordsize = 0;		/* query */
  exo->exoid = ex_open(filename, exo->cmode, &exo->comp_wordsize, 
		       &exo->io_wordsize, &exo->version);
  if (exo->exoid < 0)
    {
      sr = sprintf(err_msg, 
		   "ex_open() = %d on \"%s\" failure @ step %d, time = %g",
		   exo->exoid, filename, time_step, time_value);
      EH(-1, err_msg);
    }
  error      = ex_put_time(exo->exoid, time_step, &time_value);
  EH(error, "ex_put_time");
  error      = ex_put_var(exo->exoid, time_step, EX_NODAL, variable_index, 1,
				exo->num_nodes, vector);
  EH(error, "ex_put_var nodal");
  error      = ex_close(exo->exoid);
  return;
}

void
AF_restore_Jacobian_Flag(void)

    /*************************************************************************
     *
     *  AF_restore_Jacobian_Flag
     *
     *
     *************************************************************************/
{
  if (originalChoice_Jacobian == -23) {
    EH(-1,"restore_Jacobian_Action_Flag ERROR: no original choice");
  }
  af->Assemble_Jacobian = originalChoice_Jacobian;
  originalChoice_Jacobian = -23;
}

void
get_nv_vd_from_offset(NODAL_VARS_STRUCT *nv, int offset,
		      VARIABLE_DESCRIPTION_STRUCT **vd_ptr,  int *idof)

    /*
     *
     */
{
  int i, sum = 0;
  VARIABLE_DESCRIPTION_STRUCT *vd;
  for (i = 0; i < nv->Num_Var_Desc; i++) {
    vd = nv->Var_Desc_List[i];
    if (vd->Ndof + sum > offset) {
      *vd_ptr = vd;
      *idof = offset - sum;
      return;
    }
    sum += vd->Ndof; 
  }
  EH(-1,"RAN OUT OF vds");
}

int
nsid2nn ( int psid )

{
  int nsp;			/* node set pointer for this pointset */

  nsp     = match_nsid(psid);	/* node set pointer this node set */

  if( nsp == -1 )
    {
      if( Num_Proc == 1 )
       {
        EH(nsp, "Failed to match nodeset.");
       }
      else
        return(-1);
    }

  return( Proc_NS_List[ Proc_NS_Pointers[nsp] ] ); /* global node num */

}/* END of routine nsid2nn */

int
gminloc_int(const int value, const int local_loc, int *global_minloc)

    /********************************************************************
     *
     * gminloc_int():
     *
     *   This routine will find the minimum integer value input from
     *   all processors. It will also return the location from which
     *   the minimum occurred. In case of ties, it will return the
     *   location from the lowest numbered processor ID with the
     *   min value.
     *
     *  Output
     *  -------
     *  *global_minloc = Returns the value of local_loc from the
     *                   processor containing the min value.
     *
     *  Return
     *  -------
     *  Routine returns the minimum value of "value" input from
     *  any of the processors.
     ********************************************************************/
{
#ifdef PARALLEL
  int  in_buf[2], out_buf[2], err;
  in_buf[0] = value;
  in_buf[1] = local_loc;
  err = MPI_Allreduce((void *)in_buf, (void *) out_buf, 1, MPI_2INT,
		      MPI_MINLOC, MPI_COMM_WORLD);
  if (err != MPI_SUCCESS) {
    EH(-1, "gminloc_int: MPI_Allreduce returned an error");
  }
  *global_minloc = out_buf[1];
  return out_buf[0];
#else
  *global_minloc = 0;
  return value;
#endif
}