diff --git a/src/USER-NETCDF/dump_netcdf.cpp b/src/USER-NETCDF/dump_netcdf.cpp
index c09a131e0a..4dc242b39c 100644
--- a/src/USER-NETCDF/dump_netcdf.cpp
+++ b/src/USER-NETCDF/dump_netcdf.cpp
@@ -55,11 +55,13 @@
 #include "universe.h"
 #include "variable.h"
 #include "force.h"
+#include "output.h"
+#include "thermo.h"
 
 using namespace LAMMPS_NS;
 using namespace MathConst;
 
-enum{INT,DOUBLE};  // same as in dump_custom.cpp
+enum{INT,FLOAT,BIGINT}; // same as in thermo.cpp
 
 const char NC_FRAME_STR[]         = "frame";
 const char NC_SPATIAL_STR[]       = "spatial";
@@ -208,15 +210,15 @@ DumpNetCDF::DumpNetCDF(LAMMPS *lmp, int narg, char **arg) :
     perat[inc].field[idim] = i;
   }
 
-  n_perframe = 0;
-  perframe = NULL;
-
   n_buffer = 0;
   int_buffer = NULL;
   double_buffer = NULL;
 
   double_precision = false;
 
+  thermo = false;
+  thermovar = NULL;
+
   framei = 0;
 }
 
@@ -227,8 +229,7 @@ DumpNetCDF::~DumpNetCDF()
   closefile();
 
   delete [] perat;
-  if (n_perframe > 0)
-    delete [] perframe;
+  if (thermovar)  delete [] thermovar;
 
   if (int_buffer) memory->sfree(int_buffer);
   if (double_buffer) memory->sfree(double_buffer);
@@ -238,6 +239,11 @@ DumpNetCDF::~DumpNetCDF()
 
 void DumpNetCDF::openfile()
 {
+  if (thermo && !singlefile_opened) {
+    if (thermovar)  delete [] thermovar;
+    thermovar = new int[output->thermo->nfield];
+  }
+
   // now the computes and fixes have been initialized, so we can query
   // for the size of vector quantities
   for (int i = 0; i < n_perat; i++) {
@@ -289,30 +295,30 @@ void DumpNetCDF::openfile()
       // dimensions
       NCERRX( nc_inq_dimid(ncid, NC_FRAME_STR, &frame_dim), NC_FRAME_STR );
       NCERRX( nc_inq_dimid(ncid, NC_SPATIAL_STR, &spatial_dim),
-	      NC_SPATIAL_STR );
+          NC_SPATIAL_STR );
       NCERRX( nc_inq_dimid(ncid, NC_VOIGT_STR, &Voigt_dim), NC_VOIGT_STR );
       NCERRX( nc_inq_dimid(ncid, NC_ATOM_STR, &atom_dim), NC_ATOM_STR );
       NCERRX( nc_inq_dimid(ncid, NC_CELL_SPATIAL_STR, &cell_spatial_dim),
-	      NC_CELL_SPATIAL_STR );
+          NC_CELL_SPATIAL_STR );
       NCERRX( nc_inq_dimid(ncid, NC_CELL_ANGULAR_STR, &cell_angular_dim),
-	      NC_CELL_ANGULAR_STR );
+          NC_CELL_ANGULAR_STR );
       NCERRX( nc_inq_dimid(ncid, NC_LABEL_STR, &label_dim), NC_LABEL_STR );
 
       // default variables
       NCERRX( nc_inq_varid(ncid, NC_SPATIAL_STR, &spatial_var),
-	      NC_SPATIAL_STR );
+          NC_SPATIAL_STR );
       NCERRX( nc_inq_varid(ncid, NC_CELL_SPATIAL_STR, &cell_spatial_var),
-	      NC_CELL_SPATIAL_STR);
+          NC_CELL_SPATIAL_STR);
       NCERRX( nc_inq_varid(ncid, NC_CELL_ANGULAR_STR, &cell_angular_var),
-	      NC_CELL_ANGULAR_STR);
+          NC_CELL_ANGULAR_STR);
 
       NCERRX( nc_inq_varid(ncid, NC_TIME_STR, &time_var), NC_TIME_STR );
       NCERRX( nc_inq_varid(ncid, NC_CELL_ORIGIN_STR, &cell_origin_var),
-	      NC_CELL_ORIGIN_STR );
+          NC_CELL_ORIGIN_STR );
       NCERRX( nc_inq_varid(ncid, NC_CELL_LENGTHS_STR, &cell_lengths_var),
-	      NC_CELL_LENGTHS_STR);
+          NC_CELL_LENGTHS_STR);
       NCERRX( nc_inq_varid(ncid, NC_CELL_ANGLES_STR, &cell_angles_var),
-	      NC_CELL_ANGLES_STR);
+          NC_CELL_ANGLES_STR);
 
       // variables specified in the input file
       for (int i = 0; i < n_perat; i++) {
@@ -334,9 +340,12 @@ void DumpNetCDF::openfile()
       }
 
       // perframe variables
-      for (int i = 0; i < n_perframe; i++) {
-        NCERRX( nc_inq_varid(ncid, perframe[i].name, &perframe[i].var),
-                perframe[i].name );
+      if (thermo) {
+        Thermo *th = output->thermo;
+        for (int i = 0; i < th->nfield; i++) {
+          NCERRX( nc_inq_varid(ncid, th->keyword[i], &thermovar[i]),
+                  th->keyword[i] );
+        }
       }
 
       size_t nframes;
@@ -355,48 +364,48 @@ void DumpNetCDF::openfile()
       singlefile_opened = 1;
 
       NCERRX( nc_create(filename, NC_64BIT_DATA, &ncid),
-	      filename );
+          filename );
 
       // dimensions
       NCERRX( nc_def_dim(ncid, NC_FRAME_STR, NC_UNLIMITED, &frame_dim),
-	      NC_FRAME_STR );
+          NC_FRAME_STR );
       NCERRX( nc_def_dim(ncid, NC_SPATIAL_STR, 3, &spatial_dim),
-	      NC_SPATIAL_STR );
+          NC_SPATIAL_STR );
       NCERRX( nc_def_dim(ncid, NC_VOIGT_STR, 6, &Voigt_dim),
-	      NC_VOIGT_STR );
+          NC_VOIGT_STR );
       NCERRX( nc_def_dim(ncid, NC_ATOM_STR, ntotalgr, &atom_dim),
-	      NC_ATOM_STR );
+          NC_ATOM_STR );
       NCERRX( nc_def_dim(ncid, NC_CELL_SPATIAL_STR, 3, &cell_spatial_dim),
-	      NC_CELL_SPATIAL_STR );
+          NC_CELL_SPATIAL_STR );
       NCERRX( nc_def_dim(ncid, NC_CELL_ANGULAR_STR, 3, &cell_angular_dim),
-	      NC_CELL_ANGULAR_STR );
+          NC_CELL_ANGULAR_STR );
       NCERRX( nc_def_dim(ncid, NC_LABEL_STR, 10, &label_dim),
-	      NC_LABEL_STR );
+          NC_LABEL_STR );
 
       // default variables
       dims[0] = spatial_dim;
       NCERRX( nc_def_var(ncid, NC_SPATIAL_STR, NC_CHAR, 1, dims, &spatial_var),
-	      NC_SPATIAL_STR );
+          NC_SPATIAL_STR );
       NCERRX( nc_def_var(ncid, NC_CELL_SPATIAL_STR, NC_CHAR, 1, dims,
-			 &cell_spatial_var), NC_CELL_SPATIAL_STR );
+             &cell_spatial_var), NC_CELL_SPATIAL_STR );
       dims[0] = spatial_dim;
       dims[1] = label_dim;
       NCERRX( nc_def_var(ncid, NC_CELL_ANGULAR_STR, NC_CHAR, 2, dims,
-			 &cell_angular_var), NC_CELL_ANGULAR_STR );
+             &cell_angular_var), NC_CELL_ANGULAR_STR );
 
       dims[0] = frame_dim;
       NCERRX( nc_def_var(ncid, NC_TIME_STR, NC_DOUBLE, 1, dims, &time_var),
-	      NC_TIME_STR);
+          NC_TIME_STR);
       dims[0] = frame_dim;
       dims[1] = cell_spatial_dim;
       NCERRX( nc_def_var(ncid, NC_CELL_ORIGIN_STR, NC_DOUBLE, 2, dims,
-			 &cell_origin_var), NC_CELL_ORIGIN_STR );
+             &cell_origin_var), NC_CELL_ORIGIN_STR );
       NCERRX( nc_def_var(ncid, NC_CELL_LENGTHS_STR, NC_DOUBLE, 2, dims,
-			 &cell_lengths_var), NC_CELL_LENGTHS_STR );
+             &cell_lengths_var), NC_CELL_LENGTHS_STR );
       dims[0] = frame_dim;
       dims[1] = cell_angular_dim;
       NCERRX( nc_def_var(ncid, NC_CELL_ANGLES_STR, NC_DOUBLE, 2, dims,
-			 &cell_angles_var), NC_CELL_ANGLES_STR );
+             &cell_angles_var), NC_CELL_ANGLES_STR );
 
       // variables specified in the input file
       dims[0] = frame_dim;
@@ -423,17 +432,17 @@ void DumpNetCDF::openfile()
             // this is a tensor in Voigt notation
             dims[2] = Voigt_dim;
             NCERRX( nc_def_var(ncid, perat[i].name, xtype, 2, dims+1,
-			       &perat[i].var), perat[i].name );
+                   &perat[i].var), perat[i].name );
           }
           else if (perat[i].dims == 3) {
             // this is a vector, we need to store x-, y- and z-coordinates
             dims[2] = spatial_dim;
             NCERRX( nc_def_var(ncid, perat[i].name, xtype, 2, dims+1,
-			       &perat[i].var), perat[i].name );
+                   &perat[i].var), perat[i].name );
           }
           else if (perat[i].dims == 1) {
             NCERRX( nc_def_var(ncid, perat[i].name, xtype, 1, dims+1,
-			       &perat[i].var), perat[i].name );
+                   &perat[i].var), perat[i].name );
           }
           else {
             char errstr[1024];
@@ -448,17 +457,17 @@ void DumpNetCDF::openfile()
             // this is a tensor in Voigt notation
             dims[2] = Voigt_dim;
             NCERRX( nc_def_var(ncid, perat[i].name, xtype, 3, dims,
-			       &perat[i].var), perat[i].name );
+                   &perat[i].var), perat[i].name );
           }
           else if (perat[i].dims == 3) {
             // this is a vector, we need to store x-, y- and z-coordinates
             dims[2] = spatial_dim;
             NCERRX( nc_def_var(ncid, perat[i].name, xtype, 3, dims,
-			       &perat[i].var), perat[i].name );
+                   &perat[i].var), perat[i].name );
           }
           else if (perat[i].dims == 1) {
             NCERRX( nc_def_var(ncid, perat[i].name, xtype, 2, dims,
-			       &perat[i].var), perat[i].name );
+                   &perat[i].var), perat[i].name );
           }
           else {
             char errstr[1024];
@@ -471,14 +480,21 @@ void DumpNetCDF::openfile()
       }
 
       // perframe variables
-      for (int i = 0; i < n_perframe; i++) {
-        if (perframe[i].type == THIS_IS_A_BIGINT) {
-          NCERRX( nc_def_var(ncid, perframe[i].name, NC_LONG, 1, dims,
-			     &perframe[i].var), perframe[i].name );
-        }
-        else {
-          NCERRX( nc_def_var(ncid, perframe[i].name, NC_DOUBLE, 1, dims,
-			     &perframe[i].var), perframe[i].name );
+      if (thermo) {
+        Thermo *th = output->thermo;
+        for (int i = 0; i < th->nfield; i++) {
+          if (th->vtype[i] == FLOAT) {
+            NCERRX( nc_def_var(ncid, th->keyword[i], NC_DOUBLE, 1, dims,
+                               &thermovar[i]), th->keyword[i] );
+          }
+          else if (th->vtype[i] == INT) {
+            NCERRX( nc_def_var(ncid, th->keyword[i], NC_INT, 1, dims,
+                               &thermovar[i]), th->keyword[i] );
+          }
+          else if (th->vtype[i] == BIGINT) {
+            NCERRX( nc_def_var(ncid, th->keyword[i], NC_LONG, 1, dims,
+                               &thermovar[i]), th->keyword[i] );
+          }
         }
       }
 
@@ -622,46 +638,30 @@ void DumpNetCDF::write()
   start[0] = framei-1;
   start[1] = 0;
 
-  for (int i = 0; i < n_perframe; i++) {
-
-    if (perframe[i].type == THIS_IS_A_BIGINT) {
-      bigint data;
-      (this->*perframe[i].compute)((void*) &data);
-
-      if (filewriter)
+  if (thermo) {
+    Thermo *th = output->thermo;
+    for (int i = 0; i < th->nfield; i++) {
+      th->call_vfunc(i);
+      if (filewriter) {
+        if (th->vtype[i] == FLOAT) {
+          NCERRX( nc_put_var1_double(ncid, thermovar[i], start,
+                                     &th->dvalue),
+                  th->keyword[i] );
+        }
+        else if (th->vtype[i] == INT) {
+          NCERRX( nc_put_var1_int(ncid, thermovar[i], start, &th->ivalue),
+                  th->keyword[i] );
+        }
+        else if (th->vtype[i] == BIGINT) {
 #if defined(LAMMPS_SMALLBIG) || defined(LAMMPS_BIGBIG)
-        NCERR( nc_put_var1_long(ncid, perframe[i].var, start, &data) );
+          NCERRX( nc_put_var1_long(ncid, thermovar[i], start, &th->bivalue),
+                  th->keyword[i] );
 #else
-        NCERR( nc_put_var1_int(ncid, perframe[i].var, start, &data) );
+          NCERRX( nc_put_var1_int(ncid, thermovar[i], start, &th->bivalue),
+                  th->keyword[i] );
 #endif
-    }
-    else {
-      double data;
-      int j = perframe[i].index;
-      int idim = perframe[i].dim;
-
-      if (perframe[i].type == THIS_IS_A_COMPUTE) {
-        if (idim >= 0) {
-          modify->compute[j]->compute_vector();
-          data = modify->compute[j]->vector[idim];
         }
-        else
-          data = modify->compute[j]->compute_scalar();
       }
-      else if (perframe[i].type == THIS_IS_A_FIX) {
-        if (idim >= 0) {
-          data = modify->fix[j]->compute_vector(idim);
-        }
-        else
-          data = modify->fix[j]->compute_scalar();
-      }
-      else if (perframe[i].type == THIS_IS_A_VARIABLE) {
-        j = input->variable->find(perframe[i].id);
-        data = input->variable->compute_equal(j);
-      }
-
-      if (filewriter)
-        NCERR( nc_put_var1_double(ncid, perframe[i].var, start, &data) );
     }
   }
 
@@ -908,126 +908,19 @@ int DumpNetCDF::modify_param(int narg, char **arg)
     iarg++;
     return 2;
   }
-  else if (strcmp(arg[iarg],"global") == 0) {
-    // "perframe" quantities, i.e. not per-atom stuff
-
+  else if (strcmp(arg[iarg],"thermo") == 0) {
     iarg++;
-
-    n_perframe = narg-iarg;
-    perframe = new nc_perframe_t[n_perframe];
-
-    for (int i = 0; iarg < narg; iarg++, i++) {
-      int n;
-      char *suffix=NULL;
-
-      if (!strcmp(arg[iarg],"step")) {
-        perframe[i].type = THIS_IS_A_BIGINT;
-        perframe[i].compute = &DumpNetCDF::compute_step;
-        strcpy(perframe[i].name, arg[iarg]);
-      }
-      else if (!strcmp(arg[iarg],"elapsed")) {
-        perframe[i].type = THIS_IS_A_BIGINT;
-        perframe[i].compute = &DumpNetCDF::compute_elapsed;
-        strcpy(perframe[i].name, arg[iarg]);
-      }
-      else if (!strcmp(arg[iarg],"elaplong")) {
-        perframe[i].type = THIS_IS_A_BIGINT;
-        perframe[i].compute = &DumpNetCDF::compute_elapsed_long;
-        strcpy(perframe[i].name, arg[iarg]);
-      }
-      else {
-
-        n = strlen(arg[iarg]);
-
-        if (n > 2) {
-          suffix = new char[n-1];
-          strcpy(suffix, arg[iarg]+2);
-        }
-        else {
-          char errstr[1024];
-          sprintf(errstr, "perframe quantity '%s' must thermo quantity or "
-                  "compute, fix or variable", arg[iarg]);
-          error->all(FLERR,errstr);
-        }
-
-        if (!strncmp(arg[iarg], "c_", 2)) {
-          int idim = -1;
-          char *ptr = strchr(suffix, '[');
-
-          if (ptr) {
-            if (suffix[strlen(suffix)-1] != ']')
-              error->all(FLERR,"Missing ']' in dump modify command");
-            *ptr = '\0';
-            idim = ptr[1] - '1';
-          }
-
-          n = modify->find_compute(suffix);
-          if (n < 0)
-            error->all(FLERR,"Could not find dump modify compute ID");
-          if (modify->compute[n]->peratom_flag != 0)
-            error->all(FLERR,"Dump modify compute ID computes per-atom info");
-          if (idim >= 0 && modify->compute[n]->vector_flag == 0)
-            error->all(FLERR,"Dump modify compute ID does not compute vector");
-          if (idim < 0 && modify->compute[n]->scalar_flag == 0)
-            error->all(FLERR,"Dump modify compute ID does not compute scalar");
-
-          perframe[i].type = THIS_IS_A_COMPUTE;
-          perframe[i].dim = idim;
-          perframe[i].index = n;
-          strcpy(perframe[i].name, arg[iarg]);
-        }
-        else if (!strncmp(arg[iarg], "f_", 2)) {
-          int idim = -1;
-          char *ptr = strchr(suffix, '[');
-
-          if (ptr) {
-            if (suffix[strlen(suffix)-1] != ']')
-              error->all(FLERR,"Missing ']' in dump modify command");
-            *ptr = '\0';
-            idim = ptr[1] - '1';
-          }
-
-          n = modify->find_fix(suffix);
-          if (n < 0)
-            error->all(FLERR,"Could not find dump modify fix ID");
-          if (modify->fix[n]->peratom_flag != 0)
-            error->all(FLERR,"Dump modify fix ID computes per-atom info");
-          if (idim >= 0 && modify->fix[n]->vector_flag == 0)
-            error->all(FLERR,"Dump modify fix ID does not compute vector");
-          if (idim < 0 && modify->fix[n]->scalar_flag == 0)
-            error->all(FLERR,"Dump modify fix ID does not compute vector");
-
-          perframe[i].type = THIS_IS_A_FIX;
-          perframe[i].dim = idim;
-          perframe[i].index = n;
-          strcpy(perframe[i].name, arg[iarg]);
-        }
-        else if (!strncmp(arg[iarg], "v_", 2)) {
-          n = input->variable->find(suffix);
-          if (n < 0)
-            error->all(FLERR,"Could not find dump modify variable ID");
-          if (!input->variable->equalstyle(n))
-            error->all(FLERR,"Dump modify variable must be of style equal");
-
-          perframe[i].type = THIS_IS_A_VARIABLE;
-          perframe[i].dim = 1;
-          perframe[i].index = n;
-          strcpy(perframe[i].name, arg[iarg]);
-          strcpy(perframe[i].id, suffix);
-        }
-        else {
-          char errstr[1024];
-          sprintf(errstr, "perframe quantity '%s' must be compute, fix or "
-                  "variable", arg[iarg]);
-          error->all(FLERR,errstr);
-        }
-
-        delete [] suffix;
-
-      }
+    if (iarg >= narg)
+      error->all(FLERR,"expected 'yes' or 'no' after 'thermo' keyword.");
+    if (strcmp(arg[iarg],"yes") == 0) {
+      thermo = true;
     }
-
-    return narg;
+    else if (strcmp(arg[iarg],"no") == 0) {
+      thermo = false;
+    }
+    else error->all(FLERR,"expected 'yes' or 'no' after 'thermo' keyword.");
+    iarg++;
+    return 2;
   } else return 0;
 }
 
@@ -1101,41 +994,14 @@ void DumpNetCDF::ncerr(int err, const char *descr, int line)
     char errstr[1024];
     if (descr) {
       sprintf(errstr, "NetCDF failed with error '%s' (while accessing '%s') "
-	      " in line %i of %s.", nc_strerror(err), descr, line, __FILE__);
+          " in line %i of %s.", nc_strerror(err), descr, line, __FILE__);
     }
     else {
       sprintf(errstr, "NetCDF failed with error '%s' in line %i of %s.",
-	      nc_strerror(err), line, __FILE__);
+          nc_strerror(err), line, __FILE__);
     }
     error->one(FLERR,errstr);
   }
 }
 
-/* ----------------------------------------------------------------------
-   one method for every keyword thermo can output
-   called by compute() or evaluate_keyword()
-   compute will have already been called
-   set ivalue/dvalue/bivalue if value is int/double/bigint
-   customize a new keyword by adding a method
-------------------------------------------------------------------------- */
-
-void DumpNetCDF::compute_step(void *r)
-{
-  *((bigint *) r) = update->ntimestep;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void DumpNetCDF::compute_elapsed(void *r)
-{
-  *((bigint *) r) = update->ntimestep - update->firststep;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void DumpNetCDF::compute_elapsed_long(void *r)
-{
-  *((bigint *) r) = update->ntimestep - update->beginstep;
-}
-
 #endif /* defined(LMP_HAS_NETCDF) */
diff --git a/src/USER-NETCDF/dump_netcdf.h b/src/USER-NETCDF/dump_netcdf.h
index daf4e9d0de..049d8a91e0 100644
--- a/src/USER-NETCDF/dump_netcdf.h
+++ b/src/USER-NETCDF/dump_netcdf.h
@@ -69,22 +69,6 @@ class DumpNetCDF : public DumpCustom {
     int ndumped;                  // number of enties written for this prop.
   };
 
-  typedef void (DumpNetCDF::*funcptr_t)(void *);
-
-  // per-frame quantities (variables, fixes or computes)
-  struct nc_perframe_t {
-    char name[NC_FIELD_NAME_MAX]; // field name
-    int var;                      // NetCDF variable
-    int type;                     // variable, fix, compute or callback
-    int index;                    // index in fix/compute list
-    funcptr_t compute;            // compute function
-    int dim;                      // dimension
-    char id[NC_FIELD_NAME_MAX];   // variable id
-
-    bigint bigint_data;           // actual data
-    double double_data;           // actual data
-  };
-
   int framei;                  // current frame index
   int blocki;                  // current block index
   int ndata;                   // number of data blocks to expect
@@ -94,10 +78,10 @@ class DumpNetCDF : public DumpCustom {
   int n_perat;                 // # of netcdf per-atom properties
   nc_perat_t *perat;           // per-atom properties
 
-  int n_perframe;              // # of global netcdf (not per-atom) fix props
-  nc_perframe_t *perframe;     // global properties
+  int *thermovar;              // NetCDF variables for thermo output
 
   bool double_precision;       // write everything as double precision
+  bool thermo;                 // write thermo output to netcdf file
 
   bigint n_buffer;             // size of buffer
   int *int_buffer;             // buffer for passing data to netcdf
@@ -131,10 +115,6 @@ class DumpNetCDF : public DumpCustom {
   virtual int modify_param(int, char **);
 
   void ncerr(int, const char *, int);
-
-  void compute_step(void *);
-  void compute_elapsed(void *);
-  void compute_elapsed_long(void *);
 };
 
 }