diff --git a/src/linad99/fvar.hpp b/src/linad99/fvar.hpp
index c15c5c94f..4fa3b5476 100755
--- a/src/linad99/fvar.hpp
+++ b/src/linad99/fvar.hpp
@@ -3485,6 +3485,7 @@ int count_factor(const ivector & v);
 
  //void gradcalc( int , double *);
 void gradcalc(int nvar, const dvector & g);
+double gradcalc(int nvar, const dvector& g, dvariable& f);
 void slave_gradcalc(void);
 
 /**
diff --git a/src/linad99/sgradclc.cpp b/src/linad99/sgradclc.cpp
index be4e8fb15..ef38661ac 100755
--- a/src/linad99/sgradclc.cpp
+++ b/src/linad99/sgradclc.cpp
@@ -75,10 +75,11 @@
 #endif
 
 /**
-  Compute the gradient from the data stored in the global \ref gradient_structure.
-  \param nvar Number of variables in the gradient.
-  \param _g Vector from 1 to nvar. On return contains the gradient.
- */
+Compute the gradient from the data stored in the global \ref gradient_structure.
+
+\param nvar Number of variables in the gradient.
+\param _g Vector from 1 to nvar. On return contains the gradient.
+*/
 void gradcalc(int nvar, const dvector& _g)
 {
   if (nvar!=0)
@@ -227,6 +228,20 @@ void gradcalc(int nvar, const dvector& _g)
     gradient_structure::restore_variables();
   }
 }
+/**
+Compute the gradient from the data stored in the global \ref gradient_structure.
+
+\param nvar Number of variables in the gradient.
+\param _g Vector from 1 to nvar. On return contains the gradient.
+\param f objective function
+\returns likelihood value
+*/
+double gradcalc(int nvar, const dvector& _g, dvariable& f)
+{
+  double v = value(f);
+  gradcalc(nvar, _g);
+  return v;
+}
 /**
  */
 void gradient_structure::save_arrays()
diff --git a/src/nh99/mod_hess.cpp b/src/nh99/mod_hess.cpp
index 6069b4930..08aab2520 100755
--- a/src/nh99/mod_hess.cpp
+++ b/src/nh99/mod_hess.cpp
@@ -70,7 +70,6 @@ void function_minimizer::hess_routine_noparallel(void)
   //if (adjm_ptr) set_labels_for_hess(nvar);
   independent_variables x(1,nvar);
   initial_params::xinit(x);        // get the initial values into the x vector
-  double f;
   double delta=1.e-5;
   dvector g1(1,nvar);
   dvector g2(1,nvar);
@@ -95,8 +94,7 @@ void function_minimizer::hess_routine_noparallel(void)
       *objective_function_value::pobjfun=0.0;
       pre_userfunction();
       vf+=*objective_function_value::pobjfun;
-      f=value(vf);
-      gradcalc(nvar,g1);
+      gradcalc(nvar, g1, vf);
     }
     double sdelta1;
     double sdelta2;
@@ -114,8 +112,7 @@ void function_minimizer::hess_routine_noparallel(void)
       *objective_function_value::pobjfun=0.0;
       pre_userfunction();
       vf+=*objective_function_value::pobjfun;
-      double f = value(vf);
-      gradcalc(nvar,g1);
+      gradcalc(nvar, g1, vf);
 
       sdelta2=x(i)-delta;
       useless(sdelta2);
@@ -126,8 +123,7 @@ void function_minimizer::hess_routine_noparallel(void)
       *objective_function_value::pobjfun=0.0;
       pre_userfunction();
       vf+=*objective_function_value::pobjfun;
-      f=value(vf);
-      gradcalc(nvar,g2);
+      gradcalc(nvar, g2, vf);
       x(i)=xsave;
       hess1=(g1-g2)/(sdelta1-sdelta2);
 
@@ -140,8 +136,7 @@ void function_minimizer::hess_routine_noparallel(void)
       *objective_function_value::pobjfun=0.0;
       pre_userfunction();
       vf+=*objective_function_value::pobjfun;
-      f=value(vf);
-      gradcalc(nvar,g1);
+      gradcalc(nvar, g1, vf);
 
       x(i)=xsave-eps*delta;
       sdelta2=x(i)-eps*delta;
@@ -153,8 +148,7 @@ void function_minimizer::hess_routine_noparallel(void)
       *objective_function_value::pobjfun=0.0;
       pre_userfunction();
       vf+=*objective_function_value::pobjfun;
-      f=value(vf);
-      gradcalc(nvar,g2);
+      gradcalc(nvar, g2, vf);
       x(i)=xsave;
 
       vf=initial_params::reset(dvar_vector(x));
@@ -178,7 +172,6 @@ void function_minimizer::hess_routine_and_constraint(int iprof,
   int nvar=initial_params::nvarcalc(); // get the number of active parameters
   independent_variables x(1,nvar);
   initial_params::xinit(x);        // get the initial values into the x vector
-  double f;
   double delta=1.e-6;
   dvector g1(1,nvar);
   dvector g2(1,nvar);
@@ -203,8 +196,7 @@ void function_minimizer::hess_routine_and_constraint(int iprof,
       pre_userfunction();
       vf+=*objective_function_value::pobjfun;
       vf-=lambda*likeprof_params::likeprofptr[iprof]->variable();
-      f=value(vf);
-      gradcalc(nvar,g1);
+      gradcalc(nvar, g1, vf);
     }
     double sdelta1;
     double sdelta2;
@@ -213,7 +205,6 @@ void function_minimizer::hess_routine_and_constraint(int iprof,
     {
       hess_calcreport(i,nvar);
 
-      double f=0.0;
       double xsave=x(i);
       sdelta1=x(i)+delta;
       useless(sdelta1);
@@ -225,8 +216,7 @@ void function_minimizer::hess_routine_and_constraint(int iprof,
       pre_userfunction();
       vf+=*objective_function_value::pobjfun;
       vf-=lambda*likeprof_params::likeprofptr[iprof]->variable();
-      f=value(vf);
-      gradcalc(nvar,g1);
+      gradcalc(nvar, g1, vf);
 
       sdelta2=x(i)-delta;
       useless(sdelta2);
@@ -238,8 +228,7 @@ void function_minimizer::hess_routine_and_constraint(int iprof,
       pre_userfunction();
       vf+=*objective_function_value::pobjfun;
       vf-=lambda*likeprof_params::likeprofptr[iprof]->variable();
-      f=value(vf);
-      gradcalc(nvar,g2);
+      gradcalc(nvar, g2, vf);
       x(i)=xsave;
       hess1=(g1-g2)/(sdelta1-sdelta2);
   /*