diff --git a/ccx/src/CalculiXstep.c b/ccx/src/CalculiXstep.c
index 8e83d3e..cd77276 100644
--- a/ccx/src/CalculiXstep.c
+++ b/ccx/src/CalculiXstep.c
@@ -272,11 +272,11 @@ void CalculiXstep(int argc,char argv[][133],ITG **nelemloadp,double **xloadp,
       if(jin==0){strcpy(jobnamec,argv[1]);strcpy1(jobnamef,argv[1],132);}
 
       /* next lines deactivated on March 18, 2020 */
-    
-      /*    for(i=1;i<argc;i++){
-	    if(strcmp1(argv[i],"-o")==0) {
-	    strcpy(output,argv[i+1]);break;}
-	    }*/
+      /* Reactived by Pete to allow exodus output */
+      for(i=1;i<argc;i++){
+	if(strcmp1(argv[i],"-o")==0) {
+	  strcpy(output,argv[i+1]);break;}
+      }
     }
 
     putenv("CCX_JOBNAME_GETJOBNAME=jobnamec");
diff --git a/ccx/src/Makefile b/ccx/src/Makefile
index 055fe6e..461838a 100755
--- a/ccx/src/Makefile
+++ b/ccx/src/Makefile
@@ -150,8 +150,8 @@ LDFLAGS += `pkg-config --libs pardiso`
 ## EXODUSII
 ## CFLAGS  += -DEXODUSII
 ## LDFLAGS += -lexoIIv2c -lnetcdf
-## CFLAGS  += `pkg-config --cflags exodusii` -DEXODUSII
-## LDFLAGS += `pkg-config --libs exodusii`
+CFLAGS  += `pkg-config --cflags exodusii` -DEXODUSII
+LDFLAGS += `pkg-config --libs exodusii`
 
 
 
diff --git a/ccx/src/Makefile.inc b/ccx/src/Makefile.inc
index cfb5727..00fad78 100644
--- a/ccx/src/Makefile.inc
+++ b/ccx/src/Makefile.inc
@@ -1029,7 +1029,15 @@ dynacont.c \
 dynboun.c \
 electromagnetics.c \
 elemChecker.c \
+exo.c \
+exoselect.c \
+exosetfind.c \
+exovector.c \
 elementcpuload.c \
+exo.c \
+exoselect.c \
+exosetfind.c \
+exovector.c \
 expand.c \
 external.c \
 extrapol_dpelmain.c \
diff --git a/ccx/src/calinput.f b/ccx/src/calinput.f
index fd54645..9997c30 100644
--- a/ccx/src/calinput.f
+++ b/ccx/src/calinput.f
@@ -425,10 +425,12 @@ c     write(*,*) textpart(1)
      &       ipoinp,inp,ipoinpc,imat,ier)
       elseif((textpart(1)(1:12).eq.'*CONTACTFILE').or.
      &       (textpart(1)(1:14).eq.'*CONTACTOUTPUT')) then
-        if(textpart(1)(1:12).eq.'*CONTACTFILE') then
-          output='asc '
-        else
-          output='bin '
+         if(output(1:3).ne.'exo') then
+            if(textpart(1)(1:12).eq.'*CONTACTFILE') then
+               output='asc '
+            else
+               output='bin '
+            endif
         endif
         ifile_output=3
         call noelfiles(inpc,textpart,jout,filab,nmethod,
@@ -624,10 +626,12 @@ c
 !     
       elseif((textpart(1)(1:7).eq.'*ELFILE').or.
      &       (textpart(1)(1:14).eq.'*ELEMENTOUTPUT')) then
-        if(textpart(1)(1:7).eq.'*ELFILE') then
-          output='asc '
-        else
-          output='bin '
+         if(output(1:3).ne.'exo') then
+            if(textpart(1)(1:7).eq.'*ELFILE') then
+               output='asc '
+            else
+               output='bin '
+            endif
         endif
         ifile_output=2
         call noelfiles(inpc,textpart,jout,filab,nmethod,
@@ -858,10 +862,12 @@ c
 !     
       elseif((textpart(1)(1:9).eq.'*NODEFILE').or.
      &       (textpart(1)(1:11).eq.'*NODEOUTPUT')) then
-        if(textpart(1)(1:9).eq.'*NODEFILE') then
-          output='asc '
-        else
-          output='bin '
+         if(output(1:3).ne.'exo') then
+            if(textpart(1)(1:9).eq.'*NODEFILE') then
+               output='asc '
+            else
+               output='bin '
+            endif
         endif
         ifile_output=1
         call noelfiles(inpc,textpart,jout,filab,nmethod,
diff --git a/ccx/src/ccx_2.17.c b/ccx/src/ccx_2.17.c
index 8f5ce61..3c1d3ed 100644
--- a/ccx/src/ccx_2.17.c
+++ b/ccx/src/ccx_2.17.c
@@ -126,11 +126,11 @@ int main(int argc,char *argv[])
     if(jin==0){strcpy(jobnamec,argv[1]);strcpy1(jobnamef,argv[1],132);}
 
     /* next lines deactivated on March 18, 2020 */
-    
-    /*    for(i=1;i<argc;i++){
+    /* Reactived by Pete to allow exodus output */
+    for(i=1;i<argc;i++){
       if(strcmp1(argv[i],"-o")==0) {
 	strcpy(output,argv[i+1]);break;}
-	}*/
+    }
   }
 
   putenv("CCX_JOBNAME_GETJOBNAME=jobnamec");
@@ -145,6 +145,7 @@ int main(int argc,char *argv[])
 
   printf("\n************************************************************\n\n");
   printf("CalculiX Version 2.17, Copyright(C) 1998-2020 Guido Dhondt\n");
+  printf("CalculiX Extras version 2.17, Copyright(C) 2013-2021 Peter A. Gustafson\n");
   printf("CalculiX comes with ABSOLUTELY NO WARRANTY. This is free\n");
   printf("software, and you are welcome to redistribute it under\n");
   printf("certain conditions, see gpl.htm\n\n");
@@ -1742,8 +1743,14 @@ int main(int argc,char *argv[])
     printf("*ERROR in frd: cannot open frd file for writing...");
     exit(0);
   }
+#ifdef EXODUSII
+ if(strcmp1(output,"exo")==1){
+#endif
   fprintf(f1," 9999\n");
   fclose(f1);
+#ifdef EXODUSII
+ }
+#endif
 
   /* deallocating the fields
      this section is addressed immediately after leaving calinput */
diff --git a/ccx/src/exo.c b/ccx/src/exo.c
new file mode 100755
index 0000000..9b62690
--- /dev/null
+++ b/ccx/src/exo.c
@@ -0,0 +1,2009 @@
+/*     Calculix - A 3-dimensional finite element program                 */
+/*              Copyright (C) 1998-2018 Guido Dhondt                     */
+/*     This subroutine                                                   */
+/*              Copyright (C) 2013-2021 Peter A. Gustafson               */
+/*                                                                       */
+/*     This program is free software; you can redistribute it and/or     */
+/*     modify it under the terms of the GNU General Public License as    */
+/*     published by the Free Software Foundation(version 2);    */
+/*                                                                       */
+
+/*     This program is distributed in the hope that it will be useful,   */
+/*     but WITHOUT ANY WARRANTY; without even the implied warranty of    */
+/*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the      */
+/*     GNU General Public License for more details.                      */
+
+/*     You should have received a copy of the GNU General Public License */
+/*     along with this program; if not, write to the Free Software       */
+/*     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.         */
+
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+#include "CalculiX.h"
+#ifdef EXODUSII
+#include "exodusII.h"
+#include "exo.h"
+#endif
+
+#define min(a,b) ((a) <= (b) ? (a) : (b))
+#define max(a,b) ((a) >= (b) ? (a) : (b))
+
+void exo(double *co,ITG *nk,ITG *kon,ITG *ipkon,char *lakon,ITG *ne0,
+	 double *v,double *stn,ITG *inum,ITG *nmethod,ITG *kode,
+	 char *filab,double *een,double *t1,double *fn,double *time,
+	 double *epn,ITG *ielmat,char *matname,double *enern,
+	 double *xstaten,ITG *nstate_,ITG *istep,ITG *iinc,
+	 ITG *ithermal,double *qfn,ITG *mode,ITG *noddiam,
+	 double *trab,ITG *inotr,ITG *ntrans,double *orab,
+	 ITG *ielorien,ITG *norien,char *description,ITG *ipneigh,
+	 ITG *neigh,ITG *mi,double *stx,double *vr,double *vi,
+	 double *stnr,double *stni,double *vmax,double *stnmax,
+	 ITG *ngraph,double *veold,double *ener,ITG *ne,double *cs,
+	 char *set,ITG *nset,ITG *istartset,ITG *iendset,ITG *ialset,
+	 double *eenmax,double *fnr,double *fni,double *emn,
+	 double *thicke,char *jobnamec,char *output,double *qfx,
+         double *cdn,ITG *mortar,double *cdnr,double *cdni,ITG *nmat,
+	 ITG *ielprop,double *prop){
+
+#ifdef EXODUSII
+  /* stores the results in exo format
+
+     iselect selects which nodes are to be stored:
+     iselect=-1 means only those nodes for which inum negative
+     ist, i.e. network nodes
+     iselect=+1 means only those nodes for which inum positive
+     ist, i.e. structural nodes
+     iselect=0  means both of the above */
+
+  /* Note for frd strcmp1(output,"asc")==0 defines ascii file, otherwise binary */
+
+  char fneig[132]="",
+    material[59]="                                                          ",
+    text[2]=" ";
+
+  /* For call compatibility with frd */
+  char m1[4]="";
+  char m3[4]="";
+  FILE *f1=NULL;
+
+  static ITG nkcoords,iaxial;
+  ITG nout,noutmin,noutplus;
+  ITG nterms;
+
+  ITG i,j,k,l,m,n,o,indexe,nemax,nlayer,noutloc,iset,iselect,ncomp,nope,
+    nodes,ifield[7],nfield[2],icomp[7],ifieldstate[*nstate_],icompstate[*nstate_],nelout;
+
+  ITG mt=mi[1]+1;
+
+  ITG ncompscalar=1,ifieldscalar[1]={1},icompscalar[1]={0},nfieldscalar[2]={1,0};
+  ITG ncompvector=3,ifieldvector[3]={1,1,1},icompvector[3]={0,1,2},nfieldvector1[2]={3,0},nfieldvector0[2]={mi[1]+1,0},icompvectorlast[3]={3,4,5};
+  ITG ncomptensor=6,ifieldtensor[6]={1,1,1,1,1,1},icomptensor[6]={0,1,2,3,5,4},nfieldtensor[2]={6,0};
+  ITG ncompscalph=2,ifieldscalph[2]={1,2},icompscalph[2]={0,0},nfieldscalph[2]={0,0};
+  ITG ncompvectph=6,ifieldvectph[6]={1,1,1,2,2,2},icompvectph[6]={1,2,3,1,2,3},nfieldvectph[2]={mi[1]+1,mi[1]+1};
+  ITG ncomptensph=12,ifieldtensph[12]={1,1,1,1,1,1,2,2,2,2,2,2},icomptensph[12]={0,1,2,3,5,4,0,1,2,3,5,4},nfieldtensph[2]={6,6};
+
+  double *vold=NULL;
+ 
+  int errr=0, exoid=0;
+  ITG num_dim, num_elem;
+  ITG num_elem_blk; /* Node element blocks.  One eltype per block*/
+  ITG num_ns, num_ss, num_es, num_fs; /* Node sets, side sets, element sets, face sets */
+  int CPU_word_size = sizeof(float);
+  int IO_word_size = sizeof(float);
+
+  /* Filename */
+  strcpy (fneig, jobnamec);
+  strcat (fneig, ".exo");
+
+  /* nkcoords is the number of nodes at the time when
+     the nodal coordinates are stored in the exo file */
+  nkcoords = *nk;
+  ITG num_nodes = nkcoords;
+
+  /* determining nout, noutplus and noutmin
+     nout: number of structural and network nodes
+     noutplus: number of structural nodes
+     noutmin: number of network nodes */
+  if(*nmethod!=0){
+    nout=0;
+    noutplus=0;
+    noutmin=0;
+    for(i=0;i<*nk;i++){
+      if(inum[i]==0) continue;
+      if(inum[i]>0) noutplus++;
+      if(inum[i]<0) noutmin++;
+      nout++;
+    }
+  }else{
+    nout=*nk;
+  }
+
+  // Allocate memory for node positions
+  float *x, *y, *z;
+  x = (float *) calloc(nout, sizeof(float));
+  y = (float *) calloc(nout, sizeof(float));
+  z = (float *) calloc(nout, sizeof(float));
+
+  // Write optional node map
+  j = 0; // Counter for the exo order of the nodes
+  ITG *node_map,*node_map_inv;
+  // Note that the node_map_inv changes in subsequent passes if
+  // element deletion or incomplete results storage is used.  However,
+  // the map is recreated in exovector.c and exoselect.c so it only
+  // needs to be correct here once.
+  node_map = (ITG *)     calloc(nout, sizeof(ITG));
+  node_map_inv = (ITG *) calloc(nkcoords, sizeof(ITG));
+
+  /* storing the coordinates of the nodes */
+  if(*nmethod!=0){
+    // for(i=0;i<*nk;i++){
+    for(i=0;i<nkcoords;i++){
+      if(inum[i]==0){continue;}
+      // The difference between i and j is that not all values of i
+      // increment j.
+      node_map[j] = i+1;
+      node_map_inv[i] = j+1;
+      x[j]   = co[3*i];
+      y[j]   = co[3*i+1];
+      z[j++] = co[3*i+2];
+    }
+  }else{
+    // for(i=0;i<*nk;i++){
+    for(i=0;i<nkcoords;i++){
+      node_map[j] = i+1;
+      node_map_inv[i] = j+1;
+      x[j]   = co[3*i];
+      y[j]   = co[3*i+1];
+      z[j++] = co[3*i+2];
+    }
+  }
+
+  /* first time something is written in the exo-file: store
+     computational metadata, the nodal coordinates and the
+     topology */
+  if((*kode==1)&&((*nmethod!=5)||(*mode!=0))){
+    iaxial=0.;
+    num_dim = 3;
+    num_elem_blk = 17;
+
+    // Find the number of sets
+    exosetfind(set, nset, ialset, istartset, iendset,
+	       &num_ns, &num_ss, &num_es, &num_fs, NULL, exoid, (int) 0, nk);
+    printf ("Side sets to exo file not implemented.\n");
+    num_ss=0;
+
+#ifdef LONGLONG
+    // This handles LONGLONG transparently
+    remove(fneig);
+    exoid = ex_create (fneig, /*Filename*/
+    		       EX_ALL_INT64_API,	/* create mode */
+    		       &CPU_word_size,  /* CPU float word size in bytes */
+    		       &IO_word_size);  /* I/O float word size in bytes */
+#else
+    exoid = ex_create (fneig, /*Filename*/
+		       EX_CLOBBER,	/* create mode */
+		       &CPU_word_size,  /* CPU float word size in bytes */
+		       &IO_word_size);  /* I/O float word size in bytes */
+#endif
+
+    /* determining the number of elements */
+    if(*nmethod!=0){
+      nelout=0;
+      for(i=0;i<*ne0;i++){
+	if(ipkon[i]<=-1){
+	  continue;
+	}else if(strcmp1(&lakon[8*i],"ESPRNGC")==0){
+	  continue;
+	}else if(strcmp1(&lakon[8*i],"ESPRNGF")==0){
+	  continue;
+	}else if(strcmp1(&lakon[8*i],"DCOUP3D")==0){
+	  continue;
+	}else if(strcmp2(&lakon[8*i+6],"LC",2)==0){
+	  // Count the number of layers
+	  nlayer=0;
+	  for(k=0;k<mi[2];k++){
+	    if(ielmat[i*mi[2]+k]==0) break;
+	    nlayer++;
+	  }
+	  // Allow an element for each layer
+	  for(k=0;k<nlayer;k++){
+	    nelout++;
+	  }
+	}else{
+	  nelout++;
+	}
+      }
+    }else{
+      nelout=*ne;
+    }
+    num_elem = nelout;
+
+    /* initialize file with parameters */
+    printf("\nData writen to the .exo file\n");
+    num_nodes=nout;
+    printf("Number of nodes: %" ITGFORMAT "\n", num_nodes);
+    printf("Number of elements %" ITGFORMAT "\n", num_elem);
+    printf("Number of element blocks %" ITGFORMAT "\n", num_elem_blk);
+    printf("Number of node sets %" ITGFORMAT "\n", num_ns);
+    printf("Number of side sets %" ITGFORMAT "\n", num_ss);
+
+    errr = ex_put_init (exoid, "CalculiX EXO File",
+			num_dim, num_nodes,
+			num_elem, num_elem_blk,
+			num_ns, num_ss);
+
+    // Write values to database
+    errr = ex_put_coord (exoid, x, y, z);
+    if(errr)printf("*ERROR in exo: failed node positions");
+    errr = ex_put_id_map (exoid, EX_NODE_MAP, node_map);
+    if(errr)printf("*ERROR in exo: failed node map");
+
+    // Deallocate
+    free (node_map);
+    free (x);
+    free (y);
+    free (z);
+
+    // Write coordinate names
+    char *coord_names[3];
+    coord_names[0] = "x";
+    coord_names[1] = "y";
+    coord_names[2] = "z";
+    errr = ex_put_coord_names (exoid, coord_names);
+    if(errr){printf("*ERROR in exo: failed coordinate names");}
+
+    // Initialize enough memory to store the element numbers
+    ITG *elem_map;
+    elem_map = (ITG *) calloc(num_elem, sizeof(ITG));
+    char curblk[6]="     ";
+
+    ITG *blkassign;
+    blkassign = (ITG *) calloc(num_elem, sizeof(ITG));
+
+    l=0;
+    for(i=0;i<*ne0;i++){ // For each element.  Composite elements are
+			 // one increment in this loop and all layers
+			 // have the same element number.
+      if(ipkon[i]<=-1){
+	continue;
+      }else if(strcmp1(&lakon[8*i],"F")==0){
+	continue;
+      }else if(strcmp1(&lakon[8*i],"ESPRNGC")==0){
+	continue;
+      }else if(strcmp1(&lakon[8*i],"ESPRNGF")==0){
+	continue;
+      }else if(strcmp1(&lakon[8*i],"DCOUP3D")==0){
+	continue;
+      }else{
+	indexe=ipkon[i];
+      }
+
+      elem_map[l] = i+1;
+
+      strcpy1(curblk,&lakon[8*i],5);
+      // strcpy1(curblk,&lakon[8*i],10);
+      // printf ("%s\n", curblk);
+      strcpy1(material,&matname[80*(ielmat[i*mi[2]]-1)],5);
+      // printf ("TODO store material identifier and name.\n");
+
+      // Identify element type
+      if(strcmp1(&lakon[8*i+3],"2")==0){
+	/* 20-node brick element */
+	if(((strcmp1(&lakon[8*i+6]," ")==0)||
+	    (strcmp1(&filab[4],"E")==0)||
+	    (strcmp1(&lakon[8*i+6],"I")==0))&&
+	   (strcmp2(&lakon[8*i+6],"LC",2)!=0)){
+	  blkassign[l++]=1;
+	}else if(strcmp2(&lakon[8*i+6],"LC",2)==0){
+	  /* composite material */
+	  nlayer=0;
+	  for(k=0;k<mi[2];k++){
+	    if(ielmat[i*mi[2]+k]==0) break;
+	    nlayer++;
+	  }
+	  for(k=0;k<nlayer;k++){
+	    nemax++;
+	    elem_map[l] = i+1;
+	    blkassign[l++]=2;
+	  }
+	}else if(strcmp1(&lakon[8*i+6],"B")==0){
+	  /* 3-node beam element */
+	  blkassign[l++]=3;
+	}else{
+	  /* 8-node 2d element */
+	  blkassign[l++]=4;
+	}
+      }else if(strcmp1(&lakon[8*i+3],"8")==0){
+	if((strcmp1(&lakon[8*i+6]," ")==0)||
+	   (strcmp1(&filab[4],"E")==0)){
+	  /* 8-node brick element */
+	  blkassign[l++]=5;
+	}else if(strcmp1(&lakon[8*i+6],"B")==0){
+	  /* 2-node 1d element */
+	  blkassign[l++]=6;
+	  // if(strcmp1(&lakon[8*i+4],"R")==0){
+	  //   blkassign[l++]=6;
+	  // }else if(strcmp1(&lakon[8*i+4],"I")==0){
+	  //   blkassign[l++]=6;
+	  // }
+	}else{
+	  /* 4-node 2d element */
+	  blkassign[l++]=7;
+	  /* not sure exactly what this does, probably axisymmetric? */
+	  // if(strcmp1(&lakon[8*i+6],"A")==0) iaxial=1;
+	  // 
+	  // if(strcmp1(&lakon[8*i+4],"R")==0){
+	  //   blkassign[l++]=7;
+	  // }else if(strcmp1(&lakon[8*i+4],"I")==0){
+	  //   blkassign[l++]=7;
+	  // }	    
+	}
+      }else if((strcmp1(&lakon[8*i+3],"10")==0)||
+	       (strcmp1(&lakon[8*i+3],"14")==0)){
+	/* 10-node tetrahedral element */
+	blkassign[l++]=8;
+      }else if(strcmp1(&lakon[8*i+3],"4")==0){
+	/* 4-node tetrahedral element */
+	blkassign[l++]=9;
+      }else if(strcmp1(&lakon[8*i+3],"15")==0){
+	if(((strcmp1(&lakon[8*i+6]," ")==0)||
+	    (strcmp1(&filab[4],"E")==0))&&
+           (strcmp2(&lakon[8*i+6],"LC",2)!=0)){
+	  /* 15-node wedge element */
+	  blkassign[l++]=10;
+	}else{
+	  /* 6-node 2d element */
+	  blkassign[l++]=11;
+	  /* not sure exactly what this does */
+	  if(strcmp1(&lakon[8*i+6],"A")==0) iaxial=1;
+	}
+      }else if(strcmp1(&lakon[8*i+3],"6")==0){
+	if((strcmp1(&lakon[8*i+6]," ")==0)||
+	   (strcmp1(&filab[4],"E")==0)){
+	  /* 6-node wedge element */
+	  blkassign[l++]=12;
+	}else{
+	  /* 3-node 2d element */ /* Shells and triangles */
+	  /* not sure exactly what this does */
+	  if(strcmp1(&lakon[8*i+6],"A")==0) iaxial=1;
+	  blkassign[l++]=13;
+	}
+	//      }else if((strcmp1(&lakon[8*i],"D")==0)&&
+	//	       (strcmp1(&lakon[8*i],"DCOUP3D")!=0)){
+      }else if(strcmp1(&lakon[8*i],"D")==0){
+	if(kon[indexe]==0){
+	  /* 2-node 1d element (network entry element) */
+	  blkassign[l++]=14;
+	}else if(kon[indexe+2]==0){
+	  /* 2-node 1d element (network exit element) */
+	  blkassign[l++]=15;
+	}else{
+	  /* 3-node 1d element (genuine network element) */
+	  blkassign[l++]=16;
+	}
+      }else if((strcmp1(&lakon[8*i],"E")==0)&&
+	       (strcmp1(&lakon[8*i+6],"A")==0)){
+	/* Not sure exactly what iaxial does yet */
+	if(strcmp1(&lakon[8*i+6],"A")==0) iaxial=1;
+	/* 2-node 1d element (spring element) */
+	blkassign[l++]=17;
+      }
+    }
+
+
+    int num_nodes_per_elem[num_elem_blk], num_edges_per_elem[num_elem_blk], num_faces_per_elem[num_elem_blk];
+    char *blknames[num_elem_blk];
+    j=0;
+    num_nodes_per_elem[j]=1;  num_edges_per_elem[j]=0;  num_faces_per_elem[j]=0;  blknames[j++]="PNT";
+    num_nodes_per_elem[j]=20; num_edges_per_elem[j]=12; num_faces_per_elem[j]=6;  blknames[j++]="C3D20[R]";
+    num_nodes_per_elem[j]=20; num_edges_per_elem[j]=12; num_faces_per_elem[j]=6;  blknames[j++]="COMPOSITE LAYER C3D20";
+    num_nodes_per_elem[j]=3;  num_edges_per_elem[j]=4;  num_faces_per_elem[j]=4;  blknames[j++]="Beam B32[R]";
+    num_nodes_per_elem[j]=8;  num_edges_per_elem[j]=4;  num_faces_per_elem[j]=4;  blknames[j++]="CPS8 CPE8 CAX8 S8[R]";
+    num_nodes_per_elem[j]=8;  num_edges_per_elem[j]=12; num_faces_per_elem[j]=6;  blknames[j++]="C3D8[R]";
+    num_nodes_per_elem[j]=2;  num_edges_per_elem[j]=0;  num_faces_per_elem[j]=0;  blknames[j++]="TRUSS2";
+    num_nodes_per_elem[j]=4;  num_edges_per_elem[j]=4;  num_faces_per_elem[j]=4;  blknames[j++]="CPS4[RI] CPE4[RI] CAX4 S4[R]";
+    num_nodes_per_elem[j]=10; num_edges_per_elem[j]=6;  num_faces_per_elem[j]=4;  blknames[j++]="C3D10";
+    num_nodes_per_elem[j]=4;  num_edges_per_elem[j]=6;  num_faces_per_elem[j]=4;  blknames[j++]="C3D4";
+    num_nodes_per_elem[j]=15; num_edges_per_elem[j]=9;  num_faces_per_elem[j]=5;  blknames[j++]="C3D15";
+    num_nodes_per_elem[j]=6;  num_edges_per_elem[j]=3;  num_faces_per_elem[j]=3;  blknames[j++]="CPS6 CPE6 S6";
+    num_nodes_per_elem[j]=6;  num_edges_per_elem[j]=3;  num_faces_per_elem[j]=3;  blknames[j++]="C3D6";
+    num_nodes_per_elem[j]=3;  num_edges_per_elem[j]=3;  num_faces_per_elem[j]=3;  blknames[j++]="CPS3 CPE3 S3";
+    num_nodes_per_elem[j]=2;  num_edges_per_elem[j]=0;  num_faces_per_elem[j]=0;  blknames[j++]="2-node 1d network entry elem";
+    num_nodes_per_elem[j]=2;  num_edges_per_elem[j]=0;  num_faces_per_elem[j]=0;  blknames[j++]="2-node 1d network exit elem";
+    num_nodes_per_elem[j]=3;  num_edges_per_elem[j]=0;  num_faces_per_elem[j]=0;  blknames[j++]="2-node 1d genuine network elem";
+    num_nodes_per_elem[j]=2;  num_edges_per_elem[j]=0;  num_faces_per_elem[j]=0;  blknames[j++]="2-node 1d spring elem";
+
+    errr = ex_put_names (exoid, EX_ELEM_BLOCK, blknames);
+    if(errr){printf("*ERROR in exo: cannot write block names");}
+
+
+    /* write element connectivity */
+    ITG *connect;
+    ITG num_elem_in_blk;
+    ITG blksize[num_elem_blk];
+
+    for(l=0;l<num_elem_blk;l++){ 
+      // Go over each element type and figure out how many elements
+      // exist.  Called block size.
+      j=0;
+      m=0;
+      for(i=0;i<*ne0;i++){
+	if(ipkon[i]<0) continue;
+	if(blkassign[m]==l){
+	  if(blkassign[m]==2){//composite
+	    for(k=0;k<mi[2];k++){
+	      if(ielmat[i*mi[2]+k]==0) break;
+	      j++;
+	    }
+	  }else{
+	    j++;
+	  }
+	}
+	m++;
+      }
+
+      blksize[l]=j;
+      num_elem_in_blk=blksize[l];
+
+      connect = (ITG *) calloc (num_elem_in_blk*num_nodes_per_elem[l], sizeof(ITG));
+      // printf ("Size of connect %" ITGFORMAT "\n", num_elem_in_blk*num_nodes_per_elem[l]*sizeof(ITG));
+      k=0; o=0;
+
+      // Now connectivity
+      for(i=0;i<*ne0;i++){
+	if(ipkon[i]<0) continue;
+	indexe=ipkon[i];
+	if (blkassign[o]==l){
+	  // printf ("block assignment %" ITGFORMAT "\n", blkassign[o]);
+	  if(blkassign[o]==1){ // C3D20
+	    for(m=0;m<12;m++){connect[k++] = node_map_inv[kon[indexe+m]-1];}
+	    for(m=16;m<20;m++){connect[k++] = node_map_inv[kon[indexe+m]-1];}
+	    for(m=12;m<16;m++){connect[k++] = node_map_inv[kon[indexe+m]-1];}
+	  }else if(blkassign[o]==10){ // C3D15
+	    for(m=0;m<9;m++){connect[k++] = node_map_inv[kon[indexe+m]-1];}
+	    for(m=12;m<15;m++){connect[k++] = node_map_inv[kon[indexe+m]-1];}
+	    for(m=9;m<12;m++){connect[k++] = node_map_inv[kon[indexe+m]-1];}
+	  }else if (blkassign[o]==2){ // Composite
+	    nlayer=0;
+	    for(m=0;m<mi[2];m++){
+	      if(ielmat[i*mi[2]+m]==0) break;
+	      nlayer++;
+	    }
+	    for(n=0;n<nlayer;n++){
+	      for(m=0;m<12;m++){connect[k++] = node_map_inv[kon[indexe+28+20*n+m]-1];}
+	      for(m=16;m<20;m++){connect[k++] = node_map_inv[kon[indexe+28+20*n+m]-1];}
+	      for(m=12;m<16;m++){connect[k++] = node_map_inv[kon[indexe+28+20*n+m]-1];}
+	    }
+	  }else if(blkassign[o]==4){ // 8 Node 2D elements CAX8 S8 S8R etc
+	    for (j = 0; j <num_nodes_per_elem[l]; j++){
+	      connect[k++] = node_map_inv[kon[indexe+20+j]-1];
+	    }
+	  }else if(blkassign[o]== 5 || // C3D8 or C3D8R
+		   blkassign[o]== 8 || // C3D10
+		   blkassign[o]== 9 || // C3D4
+		   blkassign[o]==12){  // C3D6
+	    for (j = 0; j <num_nodes_per_elem[l]; j++){
+	      connect[k++] = node_map_inv[kon[indexe+j]-1];
+	    }
+	  }else if(blkassign[o]==11){ // 6-node 2D element (S6)
+	    for (j = 0; j <num_nodes_per_elem[l]; j++){
+	      connect[k++] = node_map_inv[kon[indexe+15+j]-1];
+	    }
+	  }else if(blkassign[o]==13){ // 2D triangle or 3-node shell
+	    for (j = 0; j <num_nodes_per_elem[l]; j++){
+	      connect[k++] = node_map_inv[kon[indexe+6+j]-1];
+	    }
+	  }else { // 4 node shell element?
+	    for (j = 0; j <num_nodes_per_elem[l]; j++){
+	      // RETAIN FOR A TIME TO SEE IF ANYTHING BREAKS.
+	      // Introduced new element classifications above in 2.8
+	      connect[k++] = node_map_inv[kon[indexe+8+j]-1];
+	    }
+	  }
+	}
+	o++;
+      }
+
+      int num_attr=0;
+      switch (l)
+	{
+	case 0:
+	  errr = ex_put_block (exoid, EX_ELEM_BLOCK, l, "SPHERE", num_elem_in_blk, num_nodes_per_elem[l], num_edges_per_elem[l], num_faces_per_elem[l], num_attr);
+	  break;
+	case 1:
+	  errr = ex_put_block (exoid, EX_ELEM_BLOCK, l, "HEX", num_elem_in_blk, num_nodes_per_elem[l], num_edges_per_elem[l], num_faces_per_elem[l], num_attr);
+	  break;
+	case 2:
+	  errr = ex_put_block (exoid, EX_ELEM_BLOCK, l, "HEX", num_elem_in_blk, num_nodes_per_elem[l], num_edges_per_elem[l], num_faces_per_elem[l], num_attr);
+	  break;
+	case 3:
+	  errr = ex_put_block (exoid, EX_ELEM_BLOCK, l, "TRUSS", num_elem_in_blk, num_nodes_per_elem[l], num_edges_per_elem[l], num_faces_per_elem[l], num_attr);
+	  break;
+	case 4:
+	  errr = ex_put_block (exoid, EX_ELEM_BLOCK, l, "QUAD", num_elem_in_blk, num_nodes_per_elem[l], num_edges_per_elem[l], num_faces_per_elem[l], num_attr);
+	  break;
+	case 5:
+	  errr = ex_put_block (exoid, EX_ELEM_BLOCK, l, "HEX", num_elem_in_blk, num_nodes_per_elem[l], num_edges_per_elem[l], num_faces_per_elem[l], num_attr);
+	  break;
+	case 6:
+	  errr = ex_put_block (exoid, EX_ELEM_BLOCK, l, "TRUSS", num_elem_in_blk, num_nodes_per_elem[l], num_edges_per_elem[l], num_faces_per_elem[l], num_attr);
+	  break;
+	case 7:
+	  errr = ex_put_block (exoid, EX_ELEM_BLOCK, l, "SHELL", num_elem_in_blk, num_nodes_per_elem[l], num_edges_per_elem[l], num_faces_per_elem[l], num_attr);
+	  break;
+	case 8:
+	  errr = ex_put_block (exoid, EX_ELEM_BLOCK, l, "TETRA", num_elem_in_blk, num_nodes_per_elem[l], num_edges_per_elem[l], num_faces_per_elem[l], num_attr);
+	  break;
+	case 9:
+	  errr = ex_put_block (exoid, EX_ELEM_BLOCK, l, "TETRA", num_elem_in_blk, num_nodes_per_elem[l], num_edges_per_elem[l], num_faces_per_elem[l], num_attr);
+	  break;
+	case 10:
+	  errr = ex_put_block (exoid, EX_ELEM_BLOCK, l, "WEDGE", num_elem_in_blk, num_nodes_per_elem[l], num_edges_per_elem[l], num_faces_per_elem[l], num_attr);
+	  break;
+	case 11:
+	  errr = ex_put_block (exoid, EX_ELEM_BLOCK, l, "HEX", num_elem_in_blk, num_nodes_per_elem[l], num_edges_per_elem[l], num_faces_per_elem[l], num_attr);
+	  break;
+	case 12:
+	  errr = ex_put_block (exoid, EX_ELEM_BLOCK, l, "WEDGE", num_elem_in_blk, num_nodes_per_elem[l], num_edges_per_elem[l], num_faces_per_elem[l], num_attr);
+	  break;
+	case 13:
+	  // KEEP FOR AWHILE CHECKING BREAKAGE
+	  // errr = ex_put_block (exoid, EX_ELEM_BLOCK, l, "WEDGE", num_elem_in_blk, num_nodes_per_elem[l], num_edges_per_elem[l], num_faces_per_elem[l], num_attr);
+	  errr = ex_put_block (exoid, EX_ELEM_BLOCK, l, "TRIANGLE", num_elem_in_blk, num_nodes_per_elem[l], num_edges_per_elem[l], num_faces_per_elem[l], num_attr);
+	  break;
+	default:
+	  // case 14:
+	  // case 15:
+	  // case 16:
+	  // case 17:
+	  errr = ex_put_block (exoid, EX_ELEM_BLOCK, l, "TRUSS", num_elem_in_blk, num_nodes_per_elem[l], num_edges_per_elem[l], num_faces_per_elem[l], num_attr);
+	  break;
+	};
+
+
+
+      if (num_elem_in_blk>0){
+	errr = ex_put_conn (exoid, EX_ELEM_BLOCK, l, connect, NULL, NULL);
+	if (errr)
+	  printf ("ERROR in ex_put_conn %i\n", errr);
+      }
+      free (connect);
+    }
+
+    // Write the element map into the file
+    errr = ex_put_id_map (exoid, EX_ELEM_MAP, elem_map);
+    if (errr)
+      printf ("ERROR in ex_put_id_map %i\n", errr);
+
+    // Write the node sets into the file
+    exosetfind(set, nset, ialset, istartset, iendset,
+	       &num_ns, &num_ss, &num_es, &num_fs, node_map_inv, exoid, (int) 1, nk);
+
+    // Free up memory which is gathering dust
+    free (elem_map);
+    free (blkassign);
+
+    // Close files
+    ex_update (exoid);
+    ex_close (exoid);
+
+    if(*nmethod==0){return;}
+
+    /* End of if(*kode==1) */
+  }
+
+
+
+
+  /* Start the data storage section */
+  float version;
+  exoid = ex_open (fneig, /*Filename*/
+		   EX_WRITE,	/* create mode */
+		   &CPU_word_size,  /* CPU float word size in bytes */
+		   &IO_word_size, /* I/O float word size in bytes */
+		   &version);
+
+  // Define an output variable for general use
+  float *nodal_var_vals;
+
+  int num_time_steps;
+  float fdum;
+  char cdum = 0;
+
+  float timet;
+  errr = ex_inquire (exoid, EX_INQ_TIME, &num_time_steps, &fdum, &cdum);
+  errr = ex_get_time (exoid, num_time_steps, &timet);
+  if (num_time_steps>0){
+    printf ("\t%i Time periods in exo file, most recent at time=%f\n", num_time_steps, timet);
+  } else {
+    printf ("\t0 Time periods in exo file\n");
+  }
+  timet = (float) *time;
+  ++num_time_steps;
+  printf ("\tWriting new time period %" ITGFORMAT " at time=%f\n", num_time_steps, *time);
+  errr = ex_put_time (exoid, num_time_steps, &timet);
+  if (errr) printf ("Error storing time into exo file.\n");
+
+  // Statically allocate an array of pointers.  Can't figure out how to do this dynamically.
+  // 100 should be enough to store the variable names.
+  char *var_names[100];
+
+  int countvars=0;
+  int countbool=3;
+
+  while(countbool>0){
+    // First time count
+    // Second time assign names
+    // Last time save data
+
+    /*  for cyclic symmetry frequency calculations only results for
+	even numbers (= odd modes, numbering starts at 0) are stored */
+    if((*nmethod==2)&&(((*mode/2)*2!=*mode)&&(*noddiam>=0))){ex_close(exoid);return;}
+
+    /* storing the displacements in the nodes */
+    if((*nmethod!=5)||(*mode==-1)){
+      if((strcmp1(filab,"U ")==0)&&(*ithermal!=2)){
+	if (countbool==3){
+	  countvars+=3;
+	}else if(countbool==2){
+	  var_names[countvars++]="Ux";
+	  var_names[countvars++]="Uy";
+	  var_names[countvars++]="Uz";
+	}else{
+	  iselect=1;
+
+	  frdset(filab,set,&iset,istartset,iendset,ialset,
+		 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+		 ngraph);
+
+	  exovector(v,&iset,ntrans,filab,&nkcoords,inum,m1,inotr,
+		    trab,co,istartset,iendset,ialset,mi,ngraph,f1,output,m3,
+		    exoid, num_time_steps,countvars);
+	  countvars+=3;
+	}
+      }
+    }
+    
+    /*     storing the imaginary part of displacements in the nodes
+	   for the odd modes of cyclic symmetry calculations */
+    if(*noddiam>=0){
+      if(strcmp1(filab,"U ")==0){
+	if (countbool==3){
+	  countvars+=3;
+	}else if(countbool==2){
+	  var_names[countvars++]="U-imag-x";
+	  var_names[countvars++]="U-imag-y";
+	  var_names[countvars++]="U-imag-z";
+	}else{
+	  exovector(&v[*nk*mt],&iset,ntrans,filab,&nkcoords,inum,m1,inotr,
+		    trab,co,istartset,iendset,ialset,mi,ngraph,f1,output,m3,
+		    exoid,num_time_steps,countvars);
+
+	  countvars+=3;
+	}
+      }
+    }
+    
+    /*     storing the imaginary part of displacements in the nodes
+	   for steady state calculations */
+    if((*nmethod==5)&&(*mode==0)){
+      if((strcmp1(filab,"U ")==0)&&(*ithermal!=2)){
+	if (countbool==3){
+	  countvars+=3;
+	}else if(countbool==2){
+	  var_names[countvars++]="U-imag-x";
+	  var_names[countvars++]="U-imag-y";
+	  var_names[countvars++]="U-imag-z";
+	}else{
+	  iselect=1;
+	  exovector(v,&iset,ntrans,filab,&nkcoords,inum,m1,inotr,
+		    trab,co,istartset,iendset,ialset,mi,ngraph,f1,output,m3,
+		    exoid,num_time_steps,countvars);
+	  countvars+=3;
+	}
+      }
+    }
+
+    /* storing the velocities in the nodes */
+    if((strcmp1(&filab[1740],"V   ")==0)&&(*ithermal!=2)){
+      if (countbool==3){
+	countvars+=3;
+      }else if(countbool==2){
+      	var_names[countvars++]="Vx";
+	var_names[countvars++]="Vy";
+	var_names[countvars++]="Vz";
+      }else{
+	iselect=1;
+
+	frdset(&filab[1740],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+
+	exovector(veold,&iset,ntrans,&filab[1740],&nkcoords,inum,m1,inotr,
+		  trab,co,istartset,iendset,ialset,mi,ngraph,f1,output,m3,
+		  exoid,num_time_steps,countvars);
+
+	countvars+=3;
+      }
+    }
+
+    /* storing the temperatures in the nodes */
+    if(strcmp1(&filab[87],"NT  ")==0){
+      if (countbool==3){
+	countvars+=1;
+      }else if(countbool==2){
+      	var_names[countvars++]="NT";
+      }else{
+	iselect=0;
+
+	frdset(&filab[87],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+
+	if(*ithermal<=1){
+	  exoselect(t1,t1,&iset,&nkcoords,inum,istartset,iendset,
+		    ialset,ngraph,&ncompscalar,ifieldscalar,icompscalar,
+		    nfieldscalar,&iselect,exoid,num_time_steps,countvars);
+	}else{
+	  exoselect(v,v,&iset,&nkcoords,inum,istartset,iendset,
+		    ialset,ngraph,&ncompscalar,ifieldscalar,icompscalar,
+		    nfieldvector0,&iselect,exoid,num_time_steps,countvars);
+	}
+	countvars+=1;
+      }
+    }
+
+    /* storing the electrical potential in the nodes */
+    if((strcmp1(&filab[3654],"POT ")==0)&&(*ithermal==2)){
+      if (countbool==3){
+	countvars+=1;
+      }else if(countbool==2){
+      	var_names[countvars++]="POT";
+      }else{
+	iselect=0;
+    
+	frdset(&filab[3654],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+    
+	exoselect(v,v,&iset,&nkcoords,inum,istartset,iendset,
+		  ialset,ngraph,&ncompscalar,ifieldscalar,icompscalar,
+		  nfieldvector0,&iselect,exoid,num_time_steps,countvars);
+	printf ("Warning: export POT to exo not tested.\n");
+	countvars+=1;
+      }
+    }
+    
+    /* storing the stresses in the nodes */
+    if((*nmethod!=5)||(*mode==-1)){
+      if((strcmp1(&filab[174],"S   ")==0)&&(*ithermal!=2)){
+	if (countbool==3){
+	  countvars+=6;
+	}else if(countbool==2){
+	  // Note reordered relative to frd file
+	  // Order must be xx  yy  zz  xy  xz  yz
+	  var_names[countvars++]="Sxx";
+	  var_names[countvars++]="Syy";
+	  var_names[countvars++]="Szz";
+	  var_names[countvars++]="Sxy";
+	  var_names[countvars++]="Sxz";
+	  var_names[countvars++]="Syz";
+	}else{
+	  iselect=1;
+
+	  frdset(&filab[174],set,&iset,istartset,iendset,ialset,
+		 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+		 ngraph);
+	  exoselect(stn,stn,&iset,&nkcoords,inum,istartset,iendset,
+		    ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
+		    nfieldtensor,&iselect,exoid,num_time_steps,countvars);
+	  countvars+=6;
+	}
+      }
+    }
+      
+    /* storing the imaginary part of the stresses in the nodes
+       for the odd modes of cyclic symmetry calculations */
+    if(*noddiam>=0){
+      if((strcmp1(&filab[174],"S   ")==0)&&(*ithermal!=2)){
+	if (countbool==3){
+	  countvars+=6;
+	}else if(countbool==2){
+	  // Note reordered relative to frd file
+	  // Order must be xx  yy  zz  xy  xz  yz
+	  var_names[countvars++]="S-imagxx";
+	  var_names[countvars++]="S-imagyy";
+	  var_names[countvars++]="S-imagzz";
+	  var_names[countvars++]="S-imagxy";
+	  var_names[countvars++]="S-imagzx";
+	  var_names[countvars++]="S-imagyz";
+	}else{
+	  exoselect(&stn[6**nk],stn,&iset,&nkcoords,inum,istartset,iendset,
+	  	    ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
+	  	    nfieldtensor,&iselect,exoid,num_time_steps,countvars);
+	  countvars+=6;
+	}
+      }
+    }
+
+    /* storing the imaginary part of the stresses in the nodes
+       for steady state calculations */
+    if((*nmethod==5)&&(*mode==0)){
+      if((strcmp1(&filab[174],"S   ")==0)&&(*ithermal!=2)){
+	if (countbool==3){
+	  countvars+=6;
+	}else if(countbool==2){
+	  // Note reordered relative to frd file
+	  // Order must be xx  yy  zz  xy  xz  yz
+	  var_names[countvars++]="S-imagxx";
+	  var_names[countvars++]="S-imagyy";
+	  var_names[countvars++]="S-imagzz";
+	  var_names[countvars++]="S-imagxy";
+	  var_names[countvars++]="S-imagzx";
+	  var_names[countvars++]="S-imagyz";
+	}else{
+	  iselect=1;
+
+	  frdset(&filab[174],set,&iset,istartset,iendset,ialset,
+		 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+		 ngraph);
+      
+	  exoselect(stn,stn,&iset,&nkcoords,inum,istartset,iendset,
+		    ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
+		    nfieldtensor,&iselect,exoid,num_time_steps,countvars);
+
+	  countvars+=6;
+	}      
+      }
+    }
+
+    /* storing the electromagnetic field E in the nodes */
+    if((strcmp1(&filab[3741],"EMFE")==0)&&(*ithermal!=2)){
+      if (countbool==3){
+	countvars+=3;
+      }else if(countbool==2){
+      	var_names[countvars++]="EMF-Ex";
+	var_names[countvars++]="EMF-Ey";
+	var_names[countvars++]="EMF-Ez";
+      }else{
+	iselect=1;
+
+	frdset(&filab[3741],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+    
+	exoselect(stn,stn,&iset,&nkcoords,inum,istartset,iendset,
+		  ialset,ngraph,&ncompvector,ifieldvector,icompvector,
+		  nfieldtensor,&iselect,exoid,num_time_steps,countvars);
+	printf ("Warning: export EMF-E to exo not tested.\n");
+	countvars+=3;
+      }
+    }
+
+    /* storing the electromagnetic field B in the nodes */
+    if((strcmp1(&filab[3828],"EMFB")==0)&&(*ithermal!=2)){
+      if (countbool==3){
+	countvars+=3;
+      }else if(countbool==2){
+      	var_names[countvars++]="EMF-Bx";
+	var_names[countvars++]="EMF-By";
+	var_names[countvars++]="EMF-Bz";
+      }else{
+	iselect=1;
+
+	frdset(&filab[3828],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+    
+	exoselect(stn,stn,&iset,&nkcoords,inum,istartset,iendset,
+		  ialset,ngraph,&ncompvector,ifieldvector,icompvectorlast,
+		  nfieldtensor,&iselect,exoid,num_time_steps,countvars);
+	printf ("Warning: export EMF-B to exo not tested.\n");
+	countvars+=3;
+      }
+    }
+    
+    /* storing the total strains in the nodes */
+    if((*nmethod!=5)||(*mode==-1)){
+      if((strcmp1(&filab[261],"E   ")==0)&&(*ithermal!=2)){
+	if (countbool==3){
+	  countvars+=6;
+	}else if(countbool==2){
+	  // Note reordered relative to frd file
+	  // Order must be xx  yy  zz  xy  xz  yz
+	  var_names[countvars++]="Exx";
+	  var_names[countvars++]="Eyy";
+	  var_names[countvars++]="Ezz";
+	  var_names[countvars++]="Exy";
+	  var_names[countvars++]="Exz";
+	  var_names[countvars++]="Eyz";
+	}else{
+	  iselect=1;
+
+
+	  frdset(&filab[261],set,&iset,istartset,iendset,ialset,
+		 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+		 ngraph);
+
+	  exoselect(een,een,&iset,&nkcoords,inum,istartset,iendset,
+		    ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
+		    nfieldtensor,&iselect,exoid,num_time_steps,countvars);
+
+	  countvars+=6;
+	}
+      }
+    }
+    
+    /* storing the imaginary part of the total strains in the nodes
+       for the odd modes of cyclic symmetry calculations */
+    if(*noddiam>=0){
+      if(strcmp1(&filab[261],"E   ")==0){
+	if (countbool==3){
+	  countvars+=6;
+	}else if(countbool==2){
+	  // Note reordered relative to frd file
+	  // Order must be xx  yy  zz  xy  xz  yz
+	  var_names[countvars++]="E-imagxx";
+	  var_names[countvars++]="E-imagyy";
+	  var_names[countvars++]="E-imagzz";
+	  var_names[countvars++]="E-imagxy";
+	  var_names[countvars++]="E-imagxz";
+	  var_names[countvars++]="E-imagyz";
+	}else{
+	  exoselect(&een[6**nk],een,&iset,&nkcoords,inum,istartset,iendset,
+		    ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
+		    nfieldtensor,&iselect,exoid,num_time_steps,countvars);
+
+	  countvars+=6;
+	}
+      }
+    }
+    
+    /* storing the imaginary part of the total strains in the nodes
+       for steady state calculations */
+    if((*nmethod==5)&&(*mode==0)){
+      if((strcmp1(&filab[261],"E   ")==0)&&(*ithermal!=2)){
+	if (countbool==3){
+	  countvars+=6;
+	}else if(countbool==2){
+	  // Note reordered relative to frd file
+	  // Order must be xx  yy  zz  xy  xz  yz
+	  var_names[countvars++]="E-imagxx";
+	  var_names[countvars++]="E-imagyy";
+	  var_names[countvars++]="E-imagzz";
+	  var_names[countvars++]="E-imagxy";
+	  var_names[countvars++]="E-imagxz";
+	  var_names[countvars++]="E-imagyz";
+	}else{
+	  iselect=1;
+	  
+	  frdset(&filab[261],set,&iset,istartset,iendset,ialset,
+		 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+		 ngraph);
+	  
+	  exoselect(een,een,&iset,&nkcoords,inum,istartset,iendset,
+		    ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
+		    nfieldtensor,&iselect,exoid,num_time_steps,countvars);
+	  printf ("Warning: export E-imag to exo not tested.\n");
+	  countvars+=6;
+	}
+      }
+    }
+    
+    /* storing the mechanical strains in the nodes */
+    if((*nmethod!=5)||(*mode==-1)){
+      if((strcmp1(&filab[2697],"ME  ")==0)&&(*ithermal!=2)){
+	if (countbool==3){
+	  countvars+=6;
+	}else if(countbool==2){
+	  // Note reordered relative to frd file
+	  // Order must be xx  yy  zz  xy  xz  yz
+	  var_names[countvars++]="MExx";
+	  var_names[countvars++]="MEyy";
+	  var_names[countvars++]="MEzz";
+	  var_names[countvars++]="MExy";
+	  var_names[countvars++]="MExz";
+	  var_names[countvars++]="MEyz";
+	}else{
+	  iselect=1;
+	  frdset(&filab[2697],set,&iset,istartset,iendset,ialset,
+		 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+		 ngraph);
+	  exoselect(emn,emn,&iset,&nkcoords,inum,istartset,iendset,
+		    ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
+		    nfieldtensor,&iselect,exoid,num_time_steps,countvars);
+	  countvars+=6;
+	}
+      }
+    }
+
+    /* storing the imaginary part of the mechanical strains in the nodes
+       for the odd modes of cyclic symmetry calculations */
+    if((*noddiam>=0)||((*nmethod==5)&&(*mode==0))){
+      if((strcmp1(&filab[2697],"ME  ")==0)&&(*ithermal!=2)){
+	if (countbool==3){
+	  countvars+=6;
+	}else if(countbool==2){
+	  // Note reordered relative to frd file
+	  // Order must be xx  yy  zz  xy  xz  yz
+	  var_names[countvars++]="ME-imagxx";
+	  var_names[countvars++]="ME-imagyy";
+	  var_names[countvars++]="ME-imagzz";
+	  var_names[countvars++]="ME-imagxy";
+	  var_names[countvars++]="ME-imagxz";
+	  var_names[countvars++]="ME-imagyz";
+	}else{
+	  exoselect(&emn[6**nk],een,&iset,&nkcoords,inum,istartset,iendset,
+		    ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
+		    nfieldtensor,&iselect,exoid,num_time_steps,countvars);
+	  countvars+=6;
+	}
+      }
+    }
+  
+    /* storing the forces in the nodes */
+    if((*nmethod!=5)||(*mode==-1)){    
+      if((strcmp1(&filab[348],"RF  ")==0)&&(*ithermal!=2)){
+	if (countbool==3){
+	  countvars+=3;
+	}else if(countbool==2){
+	  var_names[countvars++]="RFx";
+	  var_names[countvars++]="RFy";
+	  var_names[countvars++]="RFz";
+	}else{
+	  iselect=1;
+	  frdset(&filab[348],set,&iset,istartset,iendset,ialset,
+		 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+		 ngraph);
+	  if((iaxial==1)&&(strcmp1(&filab[352],"I")==0)){for(i=0;i<*nk;i++){fn[1+i*mt]*=180.;fn[2+i*mt]*=180.;fn[3+i*mt]*=180.;}}
+	  exovector(fn,&iset,ntrans,&filab[348],&nkcoords,inum,m1,inotr,
+		    trab,co,istartset,iendset,ialset,mi,ngraph,f1,output,m3,
+		    exoid,num_time_steps,countvars);
+	  if((iaxial==1)&&(strcmp1(&filab[352],"I")==0)){for(i=0;i<*nk;i++){fn[1+i*mt]/=180.;fn[2+i*mt]/=180.;fn[3+i*mt]/=180.;}}
+	  countvars+=3;
+	}
+      }
+    }
+    
+    /*     storing the imaginary part of the forces in the nodes
+	   for the odd modes of cyclic symmetry calculations */
+    if((*noddiam>=0)||((*nmethod==5)&&(*mode==0))){
+      if((strcmp1(&filab[348],"RF  ")==0)&&(*ithermal!=2)){
+	if (countbool==3){
+	  countvars+=3;
+	}else if(countbool==2){
+	  var_names[countvars++]="RF-imagx";
+	  var_names[countvars++]="RF-imagy";
+	  var_names[countvars++]="RF-imagz";
+	}else{
+	  exovector(&fn[*nk*mt],&iset,ntrans,filab,&nkcoords,inum,m1,inotr,
+		    trab,co,istartset,iendset,ialset,mi,ngraph,f1,output,m3,
+		    exoid,num_time_steps,countvars);
+	  countvars+=3;
+	}
+      }
+    }
+
+    /* storing the equivalent plastic strains in the nodes */
+    if((strcmp1(&filab[435],"PEEQ")==0)&&(*ithermal!=2)){
+      if (countbool==3){
+	countvars+=1;
+      }else if(countbool==2){
+	var_names[countvars++]="PEEQ";
+      }else{
+	iselect=1;
+	frdset(&filab[435],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+	exoselect(epn,epn,&iset,&nkcoords,inum,istartset,iendset,
+		  ialset,ngraph,&ncompscalar,ifieldscalar,icompscalar,
+		  nfieldscalar,&iselect,exoid,num_time_steps,countvars);
+	printf ("Warning: export PEEQ to exo not tested.\n");
+	countvars+=1;
+      }
+    }
+
+
+    /* storing the energy in the nodes */
+    if((*nmethod!=5)||(*mode==-1)){
+      if((strcmp1(&filab[522],"ENER")==0)&&(*ithermal!=2)){
+	if (countbool==3){
+	  countvars+=1;
+	}else if(countbool==2){
+	  var_names[countvars++]="ENER";
+	}else{
+	  iselect=1;
+	  frdset(&filab[522],set,&iset,istartset,iendset,ialset,
+		 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+		 ngraph);
+	  exoselect(enern,enern,&iset,&nkcoords,inum,istartset,iendset,
+		    ialset,ngraph,&ncompscalar,ifieldscalar,icompscalar,
+		    nfieldscalar,&iselect,exoid,num_time_steps,countvars);
+	  printf ("Warning: export ENER to exo not tested.\n");
+	  countvars+=1;
+	}
+      }
+    }
+
+    /* storing the contact displacements and stresses at the slave nodes */
+    /* node-to-face penalty */
+    if((strcmp1(&filab[2175],"CONT")==0)&&(*mortar!=1)&&(*ithermal!=2)&&(*nmethod!=2)){
+      if (countbool==3){
+	countvars+=6;
+      }else if(countbool==2){
+	var_names[countvars++]="COPEN";
+	var_names[countvars++]="CSLIP1";
+	var_names[countvars++]="CSLIP2";
+	var_names[countvars++]="CPRESS";
+	var_names[countvars++]="CSHEAR1";
+	var_names[countvars++]="CSHEAR2";
+      }else{
+
+	nodal_var_vals = (float *) calloc (nout, sizeof(float));
+
+	for(i=*ne-1;i>=0;i--){
+	  if((strcmp1(&lakon[8*i+1],"S")!=0)||(strcmp1(&lakon[8*i+6],"C")!=0))
+	    break;
+	}
+	noutloc=*ne-i-1;
+
+	for(j=0;j<6;j++){
+	  for(i=*ne-1;i>=0;i--){
+	    if((strcmp1(&lakon[8*i+1],"S")!=0)||(strcmp1(&lakon[8*i+6],"C")!=0))
+	      break;
+	    strcpy1(text,&lakon[8*i+7],1);
+	    nope=atoi(text)+1;
+	    nodes=node_map_inv[kon[ipkon[i]+nope-1]-1]-1;
+	    nodal_var_vals[nodes]=stx[6*mi[0]*i+j];
+	  }
+
+	  errr = ex_put_var (exoid, num_time_steps, EX_NODAL, 1+countvars++, 1, nout, nodal_var_vals);
+	  if (errr) printf ("ERROR storing contact data into exo file.\n");
+	}
+
+	free(nodal_var_vals);
+	printf ("Warning: export CONT to exo not tested.\n");
+      }
+    }
+
+    /* face-to-face penalty */
+    if((*nmethod!=5)||(*mode==-1)){
+      if((strcmp1(&filab[2175],"CONT")==0)&&(*mortar==1)&&(*ithermal!=2)){
+	if (countbool==3){
+	  countvars+=6;
+	}else if(countbool==2){
+	  var_names[countvars++]="COPEN";
+	  var_names[countvars++]="CSLIP1";
+	  var_names[countvars++]="CSLIP2";
+	  var_names[countvars++]="CPRESS";
+	  var_names[countvars++]="CSHEAR1";
+	  var_names[countvars++]="CSHEAR2";
+	}else{
+	  iselect=1;
+
+	  frdset(&filab[2175],set,&iset,istartset,iendset,ialset,
+		 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+		 ngraph);
+    
+	  exoselect(cdn,cdn,&iset,&nkcoords,inum,istartset,iendset,
+		    ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
+		    nfieldtensor,&iselect,exoid,num_time_steps,countvars);
+	  printf ("Warning: export CONT to exo not tested.\n");
+	  countvars+=6;
+	}
+      }
+    }
+
+    /* storing imaginary part of the differential contact displacements 
+       and the contact stresses for the odd modes of cyclic symmetry
+       calculations */
+    if((*noddiam>=0)||((*nmethod==5)&&(*mode==0))){
+      if((strcmp1(&filab[2175],"CONT")==0)&&(*mortar==1)&&(*ithermal!=2)){
+	if (countbool==3){
+	  countvars+=6;
+	}else if(countbool==2){
+	  var_names[countvars++]="COPEN-Imag";
+	  var_names[countvars++]="CSLIP1-Imag";
+	  var_names[countvars++]="CSLIP2-Imag";
+	  var_names[countvars++]="CPRESS-Imag";
+	  var_names[countvars++]="CSHEAR1-Imag";
+	  var_names[countvars++]="CSHEAR2-Imag";
+	}else{
+	  iselect=1;
+	  exoselect(&cdn[6**nk],cdn,&iset,&nkcoords,inum,istartset,iendset,
+		    ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
+		    nfieldtensor,&iselect,exoid,num_time_steps,countvars);
+	  countvars+=6;
+	}
+      }
+    }
+  
+    /* storing the contact energy at the slave nodes */
+    if((strcmp1(&filab[2262],"CELS")==0)&&(*ithermal!=2)){
+      if (countbool==3){
+	countvars+=1;
+      }else if(countbool==2){
+	var_names[countvars++]="CELS";
+      }else{
+	for(i=*ne-1;i>=0;i--){
+	  if((strcmp1(&lakon[8*i+1],"S")!=0)||(strcmp1(&lakon[8*i+6],"C")!=0))
+	    break;
+	}
+	noutloc=*ne-i-1;
+
+	nodal_var_vals = (float *) calloc (nkcoords, sizeof(float));
+	for(i=*ne-1;i>=0;i--){
+	  if((strcmp1(&lakon[8*i+1],"S")!=0)||(strcmp1(&lakon[8*i+6],"C")!=0))
+	    break;
+	  nope=atoi(&lakon[8*i+7])+1;
+	  nodes=node_map_inv[kon[ipkon[i]+nope-1]-1]-1;
+	  nodal_var_vals[nodes]=ener[i*mi[0]];
+	}
+
+	countvars+=1;
+	int errr = ex_put_var (exoid, num_time_steps, EX_NODAL, countvars, 1, nout, nodal_var_vals);
+	if (errr) printf ("ERROR storing CELS data into exo file.\n");
+	printf ("Warning: export CELS to exo not tested.\n");
+	free(nodal_var_vals);
+      }
+    }
+
+    /* storing the internal state variables in the nodes */
+    if(strcmp1(&filab[609],"SDV ")==0){
+      if (countbool==3){
+	countvars+=*nstate_;
+      }else if(countbool==2){
+	for(j=1;j<=*nstate_;j++){
+	  sprintf(var_names[countvars++],"SDV%" ITGFORMAT,j);
+	}
+      }else{
+	iselect=1;
+
+	frdset(&filab[609],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+
+	for(i=0;i<*nstate_;i++){
+	  ifieldstate[i]=1;icompstate[i]=i;
+	}
+	nfield[0]=*nstate_;
+
+	exoselect(xstaten,xstaten,&iset,&nkcoords,inum,istartset,iendset,
+		  ialset,ngraph,nstate_,ifieldstate,icompstate,
+		  nfield,&iselect,exoid,num_time_steps,countvars);
+	countvars+=*nstate_;
+      }
+    }
+
+    /* storing the heat flux in the nodes
+       the heat flux has been extrapolated from the integration points
+       in subroutine extropolate.f, taking into account whether the
+       results are requested in the global system or in a local system.
+       Therefore, subroutine exovector cannot be used, since it assumes
+       the values are stored in the global system */
+    if((strcmp1(&filab[696],"HFL ")==0)&&(*ithermal>1)){
+      if (countbool==3){
+	countvars+=3;
+      }else if(countbool==2){
+	var_names[countvars++]="HFLx";
+	var_names[countvars++]="HFLy";
+	var_names[countvars++]="HFLz";
+      }else{
+	iselect=1;
+
+	frdset(&filab[696],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+
+	exoselect(qfn,qfn,&iset,&nkcoords,inum,istartset,iendset,
+		  ialset,ngraph,&ncompvector,ifieldvector,icompvector,
+		  nfieldvector1,&iselect,exoid,num_time_steps,countvars);
+	countvars+=3;
+      }
+    }
+
+    /* storing the electrical current in the nodes
+       (cf. heat flux HFL above)  */
+    
+    if((strcmp1(&filab[3567],"ECD ")==0)&&(*ithermal==2)){
+      if (countbool==3){
+	countvars+=3;
+      }else if(countbool==2){
+	var_names[countvars++]="ECDx";
+	var_names[countvars++]="ECDy";
+	var_names[countvars++]="ECDz";
+      }else{
+	iselect=1;
+    
+	frdset(&filab[3567],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+    
+	exoselect(qfn,qfn,&iset,&nkcoords,inum,istartset,iendset,
+		  ialset,ngraph,&ncompvector,ifieldvector,icompvector,
+		  nfieldvector1,&iselect,exoid,num_time_steps,countvars);
+	printf ("Warning: export ECD to exo not tested.\n");
+	countvars+=3;
+      }
+    }
+    
+    /* storing the heat generation in the nodes */
+    if((strcmp1(&filab[783],"RFL ")==0)&&(*ithermal>1)){
+      if (countbool==3){
+	countvars+=1;
+      }else if(countbool==2){
+	var_names[countvars++]="RFL";
+      }else{
+	iselect=1;
+
+	frdset(&filab[783],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+	exoselect(fn,fn,&iset,&nkcoords,inum,istartset,iendset,
+		  ialset,ngraph,&ncompscalar,ifieldscalar,icompscalar,
+		  nfieldvector0,&iselect,exoid,num_time_steps,countvars);
+	countvars+=1;
+      }
+    }
+
+    /* storing the Zienkiewicz-Zhu improved stresses in the nodes */
+    if((*nmethod!=5)||(*mode==-1)){
+      if((strcmp1(&filab[1044],"ZZS")==0)&&(*ithermal!=2)){
+	if (countbool==3){
+	  countvars+=6;
+	}else if(countbool==2){
+	  // Note reordered relative to frd file
+	  // Order must be xx  yy  zz  xy  xz  yz
+	  var_names[countvars++]="ZZSxx";
+	  var_names[countvars++]="ZZSyy";
+	  var_names[countvars++]="ZZSzz";
+	  var_names[countvars++]="ZZSxy";
+	  var_names[countvars++]="ZZSxz";
+	  var_names[countvars++]="ZZSyz";
+	}else{
+	  FORTRAN(zienzhu,(co,nk,kon,ipkon,lakon,ne0,stn,ipneigh,neigh,
+			   stx,&mi[0]));
+
+	  iselect=1;
+
+	  frdset(&filab[1044],set,&iset,istartset,iendset,ialset,
+		 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+		 ngraph);
+
+	  exoselect(stn,stn,&iset,&nkcoords,inum,istartset,iendset,
+		    ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
+		    nfieldtensor,&iselect,exoid,num_time_steps,countvars);
+	  countvars+=6;
+	}
+      }
+    }
+    
+    /* storing the imaginary part of the Zienkiewicz-Zhu
+       improved stresses in the nodes
+       for the odd modes of cyclic symmetry calculations */
+    
+    if((*noddiam>=0)||((*nmethod==5)&&(*mode==0))){
+      if((strcmp1(&filab[1044],"ZZS")==0)&&(*ithermal!=2)){
+	if (countbool==3){
+	  countvars+=6;
+	}else if(countbool==2){
+	  // Note reordered relative to frd file
+	  // Order must be xx  yy  zz  xy  xz  yz
+	  var_names[countvars++]="ZZS-imagxx";
+	  var_names[countvars++]="ZZS-imagyy";
+	  var_names[countvars++]="ZZS-imagzz";
+	  var_names[countvars++]="ZZS-imagxy";
+	  var_names[countvars++]="ZZS-imagxz";
+	  var_names[countvars++]="ZZS-imagyz";
+	}else{
+
+	  FORTRAN(zienzhu,(co,nk,kon,ipkon,lakon,ne0,stn,ipneigh,neigh,
+			   &stx[6*mi[0]**ne],&mi[0]));
+
+	  exoselect(stn,stn,&iset,&nkcoords,inum,istartset,iendset,
+		    ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
+		    nfieldtensor,&iselect,exoid,num_time_steps,countvars);
+	  printf ("Warning: export of ZZSTR imaginary to exo not tested.\n");
+	  countvars+=6;
+	}
+      }
+    }
+
+    /* storing the error estimator in the nodes */
+    if((*nmethod!=5)||(*mode==-1)){
+      if((strcmp1(&filab[1044],"ERR")==0)&&(*ithermal!=2)){
+	if (countbool==3){
+	  countvars+=2;
+	}else if(countbool==2){
+	  var_names[countvars++]="PSTDERROR";
+	  var_names[countvars++]="VMSTDERROR";
+	}else{
+
+	  nterms=6;
+	  FORTRAN(errorestimator,(stx,stn,ipkon,kon,lakon,nk,ne,
+				  mi,ielmat,&nterms,inum,co,vold,&filab[1048],
+				  ielprop,prop));
+
+	  iselect=1;
+
+	  frdset(&filab[1044],set,&iset,istartset,iendset,ialset,
+		 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+		 ngraph);
+
+	  ncomp=2;
+	  ifield[0]=1;ifield[1]=1;
+	  icomp[0]=0;icomp[1]=1;
+
+	  exoselect(stn,stn,&iset,&nkcoords,inum,istartset,iendset,
+		    ialset,ngraph,&ncomp,ifield,icomp,
+		    nfieldtensor,&iselect,exoid,num_time_steps,countvars);
+
+	  countvars+=2;
+	}
+      }
+    }
+
+    /* storing the imaginary part of the error estimator in the nodes
+       for the odd modes of cyclic symmetry calculations */
+    if((*noddiam>=0)||((*nmethod==5)&&(*mode==0))){
+      if((strcmp1(&filab[1044],"ERR")==0)&&(*ithermal!=2)){
+	if (countbool==3){
+	  countvars+=2;
+	}else if(countbool==2){
+	  var_names[countvars++]="PSTD-imag-ERROR";
+	  var_names[countvars++]="VMSTD-imag-ERROR";
+	}else{
+	  nterms=6;
+	  FORTRAN(errorestimator,(&stx[6*mi[0]**ne],stn,ipkon,kon,lakon,nk,ne,
+				  mi,ielmat,&nterms,inum,co,vold,&filab[1048],
+				  ielprop,prop));
+
+	  ncomp=2;
+	  ifield[0]=1;ifield[1]=1;
+	  icomp[0]=0;icomp[1]=1;
+
+	  exoselect(stn,stn,&iset,&nkcoords,inum,istartset,iendset,
+		    ialset,ngraph,&ncomp,ifield,icomp,
+		    nfieldtensor,&iselect,exoid,num_time_steps,countvars);
+
+	  countvars+=2;
+	}
+      }
+    }
+
+    /* storing the thermal error estimator in the nodes */
+    if((*nmethod!=5)||(*mode==-1)){
+      if((strcmp1(&filab[2784],"HER")==0)&&(*ithermal>1)){
+	if (countbool==3){
+	  countvars+=1;
+	}else if(countbool==2){
+	  var_names[countvars++]="HFLSTD";
+	}else{
+	  nterms=3;
+	  FORTRAN(errorestimator,(qfx,qfn,ipkon,kon,lakon,nk,ne,
+				  mi,ielmat,&nterms,inum,co,vold,&filab[2788],
+				  ielprop,prop));
+
+	  iselect=1;
+	  frdset(&filab[2784],set,&iset,istartset,iendset,ialset,
+		 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+		 ngraph);
+
+	  ncomp=1;
+	  ifield[0]=1;
+	  icomp[0]=1;
+
+	  exoselect(qfn,qfn,&iset,&nkcoords,inum,istartset,iendset,
+		    ialset,ngraph,&ncomp,ifield,icomp,
+		    nfieldvector1,&iselect,exoid,num_time_steps,countvars);
+	  countvars+=1;
+	}
+      }
+    }
+
+    /* storing the imaginary part of the thermal error estimator in the nodes
+       for the odd modes of cyclic symmetry calculations */
+  
+    if((*noddiam>=0)||((*nmethod==5)&&(*mode==0))){
+      if(strcmp1(&filab[2784],"HER")==0){
+	if (countbool==3){
+	  countvars+=1;
+	}else if(countbool==2){
+	  var_names[countvars++]="HFLSTD-IMG";
+	}else{
+	  nterms=3;
+	  FORTRAN(errorestimator,(&qfx[3*mi[0]**ne],qfn,ipkon,kon,lakon,nk,ne,
+				  mi,ielmat,&nterms,inum,co,vold,&filab[2788],
+				  ielprop,prop));
+
+	  ncomp=1;
+	  ifield[0]=1;
+	  icomp[0]=1;
+	  
+	  exoselect(qfn,qfn,&iset,&nkcoords,inum,istartset,iendset,
+		    ialset,ngraph,&ncomp,ifield,icomp,
+		    nfieldvector1,&iselect,exoid,num_time_steps,countvars);
+
+	  countvars+=1;
+	}
+      }
+    }
+
+    /* storing the total temperatures in the network nodes */
+    if((strcmp1(&filab[1131],"TT  ")==0)&&(*ithermal>1)){
+      if (countbool==3){
+	countvars+=1;
+      }else if(countbool==2){
+	var_names[countvars++]="TT";
+      }else{
+
+	iselect=-1;
+	frdset(&filab[1131],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+
+	exoselect(v,v,&iset,&nkcoords,inum,istartset,iendset,
+		  ialset,ngraph,&ncompscalar,ifieldscalar,icompscalar,
+		  nfieldvector0,&iselect,exoid,num_time_steps,countvars);
+	printf ("Warning: export of TT to exo not tested.\n");
+	countvars+=1;
+      }
+    }
+
+    /* storing the total temperatures in the network nodes */
+    if((strcmp1(&filab[1218],"MF  ")==0)&&(*ithermal>1)){
+      if (countbool==3){
+	countvars+=1;
+      }else if(countbool==2){
+	var_names[countvars++]="MF";
+      }else{
+
+	iselect=-1;
+	frdset(&filab[1218],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+
+	icomp[0]=1;
+	if((iaxial==1)&&(strcmp1(&filab[1222],"I")==0)){for(i=0;i<*nk;i++)v[1+i*mt]*=180.;}
+	exoselect(v,v,&iset,&nkcoords,inum,istartset,iendset,
+		  ialset,ngraph,&ncompscalar,ifieldscalar,icomp,
+		  nfieldvector0,&iselect,exoid,num_time_steps,countvars);
+	if((iaxial==1)&&(strcmp1(&filab[1222],"I")==0)){for(i=0;i<*nk;i++)v[1+i*mt]/=180.;}
+	countvars+=1;
+      }
+    }
+
+    /* storing the total pressure in the network nodes */
+    if((strcmp1(&filab[1305],"PT  ")==0)&&(*ithermal>1)){
+      if (countbool==3){
+	countvars+=1;
+      }else if(countbool==2){
+	var_names[countvars++]="PT";
+      }else{
+
+	iselect=-1;
+	frdset(&filab[1305],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+
+	icomp[0]=2;
+	exoselect(v,v,&iset,&nkcoords,inum,istartset,iendset,
+		  ialset,ngraph,&ncompscalar,ifieldscalar,icomp,
+		  nfieldvector0,&iselect,exoid,num_time_steps,countvars);
+	countvars+=1;
+      }
+    }
+
+    /* storing the static pressure in the liquid network nodes */
+    if((strcmp1(&filab[1827],"PS  ")==0)&&(*ithermal>1)){
+      if (countbool==3){
+	countvars+=1;
+      }else if(countbool==2){
+	var_names[countvars++]="STPRESS PS";
+      }else{
+
+	iselect=-1;
+	frdset(&filab[1827],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+
+	icomp[0]=2;
+	exoselect(v,v,&iset,&nkcoords,inum,istartset,iendset,
+		  ialset,ngraph,&ncompscalar,ifieldscalar,icomp,
+		  nfieldvector0,&iselect,exoid,num_time_steps,countvars);
+	printf ("Warning: export of PS to exo not tested.\n");
+	countvars+=1;
+      }
+    }
+
+    /* storing the liquid depth in the channel nodes */
+    if(strcmp1(&filab[2349],"PS  ")==0){
+      if (countbool==3){
+	countvars+=1;
+      }else if(countbool==2){
+	var_names[countvars++]="DEPTH";
+      }else{
+
+	iselect=-1;
+	frdset(&filab[2349],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+
+	icomp[0]=2;
+	exoselect(v,v,&iset,&nkcoords,inum,istartset,iendset,
+		  ialset,ngraph,&ncompscalar,ifieldscalar,icomp,
+		  nfieldvector0,&iselect,exoid,num_time_steps,countvars);
+	printf ("Warning: export of DEPTH to exo not tested.\n");
+	countvars+=1;
+      }
+    }
+
+    /* storing the critical depth in the channel nodes */
+    if((strcmp1(&filab[2436],"HCRI")==0)&&(*ithermal>1)){
+      if (countbool==3){
+	countvars+=1;
+      }else if(countbool==2){
+	var_names[countvars++]="HCRIT";
+      }else{
+
+	iselect=-1;
+	frdset(&filab[2436],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+
+	icomp[0]=3;
+	exoselect(v,v,&iset,&nkcoords,inum,istartset,iendset,
+		  ialset,ngraph,&ncompscalar,ifieldscalar,icomp,
+		  nfieldvector0,&iselect,exoid,num_time_steps,countvars);
+	printf ("Warning: export of HCRIT to exo not tested.\n");
+	countvars+=1;
+      }
+    }
+
+    /* storing the static temperature in the network nodes */
+    if((strcmp1(&filab[1392],"TS  ")==0)&&(*ithermal>1)){
+      if (countbool==3){
+	countvars+=1;
+      }else if(countbool==2){
+	var_names[countvars++]="STTEMP";
+      }else{
+
+	iselect=-1;
+	frdset(&filab[1392],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+
+	icomp[0]=3;
+	exoselect(v,v,&iset,&nkcoords,inum,istartset,iendset,
+		  ialset,ngraph,&ncompscalar,ifieldscalar,icomp,
+		  nfieldvector0,&iselect,exoid,num_time_steps,countvars);
+	printf ("Warning: export of STTEMP to exo not tested.\n");
+	countvars+=1;
+      }
+    }
+
+    /*  the remaining lines only apply to frequency calculations
+	with cyclic symmetry, complex frequency and steady state calculations */
+
+    if((*nmethod!=2)&&(*nmethod!=5)&&(*nmethod!=6)&&(*nmethod!=7)){goto WRITENAMES;}
+    if((*nmethod==5)&&(*mode==-1)){goto WRITENAMES;}
+
+    /* storing the displacements in the nodes (magnitude, phase) */
+    if((strcmp1(&filab[870],"PU  ")==0)&&(*ithermal!=2)){
+      if (countbool==3){
+	countvars+=6;
+      }else if(countbool==2){
+	// Note reordered relative to frd file
+	// Order must be xx  yy  zz  xy  xz  yz
+	var_names[countvars++]="PDISP-MAGx";
+	var_names[countvars++]="PDISP-MAGy";
+	var_names[countvars++]="PDISP-MAGz";
+	var_names[countvars++]="PDISP-PHAx";
+	var_names[countvars++]="PDISP-PHAz";
+	var_names[countvars++]="PDISP-PHAy";
+      }else{
+	iselect=1;
+
+	frdset(&filab[870],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+
+	exoselect(vr,vi,&iset,&nkcoords,inum,istartset,iendset,
+		  ialset,ngraph,&ncompvectph,ifieldvectph,icompvectph,
+		  nfieldvectph,&iselect,exoid,num_time_steps,countvars);
+	countvars+=6;
+      }
+    }
+
+    /* storing the temperatures in the nodes (magnitude, phase) */
+    if((strcmp1(&filab[957],"PNT ")==0)&&(*ithermal>1)){
+      if (countbool==3){
+	countvars+=2;
+      }else if(countbool==2){
+	var_names[countvars++]="PNDTEMP MAG1";
+	var_names[countvars++]="PNDTEMP PHA2";
+      }else{
+	iselect=1;
+
+	frdset(&filab[957],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+
+	exoselect(vr,vi,&iset,&nkcoords,inum,istartset,iendset,
+		  ialset,ngraph,&ncompscalph,ifieldscalph,icompscalph,
+		  nfieldscalph,&iselect,exoid,num_time_steps,countvars);
+	printf ("Warning: export of PNDTEMP to exo not tested.\n");
+	countvars+=2;
+      }
+    }
+
+    /* storing the stresses in the nodes (magnitude, phase) */
+	  
+    if((strcmp1(&filab[1479],"PHS ")==0)&&(*ithermal!=2)){
+      if (countbool==3){
+	countvars+=12;
+      }else if(countbool==2){
+	// Note reordered relative to frd file
+	// Order must be xx  yy  zz  xy  xz  yz
+	var_names[countvars++]="PSTRESS MAGXX";
+	var_names[countvars++]="PSTRESS MAGYY";
+	var_names[countvars++]="PSTRESS MAGZZ";
+	var_names[countvars++]="PSTRESS MAGXY";
+	var_names[countvars++]="PSTRESS MAGXZ";
+	var_names[countvars++]="PSTRESS MAGYZ";
+	var_names[countvars++]="PSTRESS PHAXX";
+	var_names[countvars++]="PSTRESS PHAYY";
+	var_names[countvars++]="PSTRESS PHAZZ";
+	var_names[countvars++]="PSTRESS PHAXY";
+	var_names[countvars++]="PSTRESS PHAXZ";
+	var_names[countvars++]="PSTRESS PHAYZ";
+      }else{
+	iselect=1;
+
+	frdset(&filab[1479],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+
+	exoselect(stnr,stni,&iset,&nkcoords,inum,istartset,iendset,
+		  ialset,ngraph,&ncomptensph,ifieldtensph,icomptensph,
+		  nfieldtensph,&iselect,exoid,num_time_steps,countvars);
+	countvars+=12;
+      }
+    }
+
+    /* storing the differential contact displacements and
+       the contact stresses in the nodes (magnitude, phase)
+       only for face-to-face penalty contact */
+    if((strcmp1(&filab[3915],"PCON")==0)&&(*ithermal!=2)&&(*mortar==1)){
+      if (countbool==3){
+	countvars+=12;
+      }else if(countbool==2){
+	// NOT SURE IF THE ORDER IS CORRECT HERE... ORDERING GETS
+	// CHANGED IN EXOSELECT
+	var_names[countvars++]="MAGO  ";
+	var_names[countvars++]="MAGSL1";
+	var_names[countvars++]="MAGSL2";
+	var_names[countvars++]="MAGP  ";
+	var_names[countvars++]="MAGSH1";
+	var_names[countvars++]="MAGSH2";
+	var_names[countvars++]="PHAO  ";
+	var_names[countvars++]="PHASL1";
+	var_names[countvars++]="PHASL2";
+	var_names[countvars++]="PHAP  ";
+	var_names[countvars++]="PHASH1";
+	var_names[countvars++]="PHASH2";
+      }else{
+	iselect=1;
+
+	frdset(&filab[3915],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+    
+	exoselect(cdnr,cdni,&iset,&nkcoords,inum,istartset,iendset,
+		  ialset,ngraph,&ncomptensph,ifieldtensph,icomptensph,
+		  nfieldtensph,&iselect,exoid,num_time_steps,countvars);
+	countvars+=12;
+	printf ("Warning: export PCON to exo not tested and likely has incorrect ordering.\n");
+      }
+    }
+
+    /* storing the displacements in the nodes (magnitude, phase) */
+    if(strcmp1(&filab[2610],"PRF ")==0){
+      if (countbool==3){
+	countvars+=6;
+      }else if(countbool==2){
+	// Note reordered relative to frd file
+	// Order must be xx  yy  zz  xy  xz  yz
+	var_names[countvars++]="PFORC MAG1";
+	var_names[countvars++]="PFORC MAG2";
+	var_names[countvars++]="PFORC MAG3";
+	var_names[countvars++]="PFORC PHA1";
+	var_names[countvars++]="PFORC PHA3";
+	var_names[countvars++]="PFORC PHA2";
+      }else{
+	iselect=1;
+
+	frdset(&filab[2610],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+
+	exoselect(fnr,fni,&iset,&nkcoords,inum,istartset,iendset,
+		  ialset,ngraph,&ncompvectph,ifieldvectph,icompvectph,
+		  nfieldvectph,&iselect,exoid,num_time_steps,countvars);
+	printf ("Warning: export of PFORC to exo not tested.\n");
+	countvars+=6;
+      }
+    }
+
+    /* the remaining parts are for frequency calculations with cyclic symmetry only */
+    if(*nmethod!=2){goto WRITENAMES;}
+
+    /* storing the maximum displacements of the nodes in the base sector
+       (components, magnitude) */
+	  
+    if((strcmp1(&filab[1566],"MAXU")==0)&&(*ithermal!=2)){
+      if (countbool==3){
+	countvars+=4;
+      }else if(countbool==2){
+	var_names[countvars++]="MDISP DX";
+	var_names[countvars++]="MDISP DY";
+	var_names[countvars++]="MDISP DZ";
+	var_names[countvars++]="MDISP ANG";
+      }else{
+	iselect=1;
+
+	frdset(&filab[1566],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+
+	ncomp=4;
+	ifield[0]=1;icomp[0]=1;
+	ifield[1]=1;icomp[1]=2;
+	ifield[2]=1;icomp[2]=3;
+	ifield[3]=1;icomp[3]=0;
+	nfield[0]=4;nfield[1]=4;
+
+	exoselect(vmax,vmax,&iset,&nkcoords,inum,istartset,iendset,
+		  ialset,ngraph,&ncomp,ifield,icomp,
+		  nfield,&iselect,exoid,num_time_steps,countvars);
+	printf ("Warning: export of MDISP to exo not tested.\n");
+	countvars+=4;
+      }
+    }
+
+    /* storing the worst principal stress at the nodes in the basis
+       sector (components, magnitude)
+    
+       the worst principal stress is the maximum of the absolute value
+       of all principal stresses, times its original sign */
+      
+    if((strcmp1(&filab[1653],"MAXS")==0)&&(*ithermal!=2)){
+      if (countbool==3){
+	countvars+=7;
+      }else if(countbool==2){
+	// Note reordered relative to frd file
+	// Order must be xx  yy  zz  xy  xz  yz
+	var_names[countvars++]="MSTRESS SXX";
+	var_names[countvars++]="MSTRESS SYY";
+	var_names[countvars++]="MSTRESS SZZ";
+	var_names[countvars++]="MSTRESS SXY";
+	var_names[countvars++]="MSTRESS SXZ";
+	var_names[countvars++]="MSTRESS SYZ";
+	var_names[countvars++]="MSTRESS MAG";
+      }else{
+	iselect=1;
+
+	frdset(&filab[1653],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+
+	ncomp=7;
+	ifield[0]=1;icomp[0]=1;
+	ifield[1]=1;icomp[1]=2;
+	ifield[2]=1;icomp[2]=3;
+	ifield[3]=1;icomp[3]=4;
+	ifield[4]=1;icomp[4]=6;
+	ifield[5]=1;icomp[5]=5;
+	ifield[6]=1;icomp[6]=0;
+	nfield[0]=7;nfield[1]=7;
+
+	exoselect(stnmax,stnmax,&iset,&nkcoords,inum,istartset,iendset,
+		  ialset,ngraph,&ncomp,ifield,icomp,
+		  nfield,&iselect,exoid,num_time_steps,countvars);
+	printf ("Warning: export of MSTRESS to exo not tested.\n");
+	countvars+=7;
+      }
+    }
+
+    /* storing the worst principal strain at the nodes
+       in the basis sector (components, magnitude)
+
+       the worst principal strain is the maximum of the
+       absolute value of all principal strains, times
+       its original sign */
+    if((strcmp1(&filab[2523],"MAXE")==0)&&(*ithermal!=2)){
+      if (countbool==3){
+	countvars+=7;
+      }else if(countbool==2){
+	// Note reordered relative to frd file
+	// Order must be xx  yy  zz  xy  xz  yz
+	var_names[countvars++]="MAXE EXX";
+	var_names[countvars++]="MAXE EYY";
+	var_names[countvars++]="MAXE EZZ";
+	var_names[countvars++]="MAXE EXY";
+	var_names[countvars++]="MAXE EXZ";
+	var_names[countvars++]="MAXE EYZ";
+	var_names[countvars++]="MAXE MAG";
+      }else{
+	iselect=1;
+
+	frdset(&filab[2523],set,&iset,istartset,iendset,ialset,
+	       inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
+	       ngraph);
+
+	ncomp=7;
+	ifield[0]=1;icomp[0]=1;
+	ifield[1]=1;icomp[1]=2;
+	ifield[2]=1;icomp[2]=3;
+	ifield[3]=1;icomp[3]=4;
+	ifield[4]=1;icomp[4]=6;
+	ifield[5]=1;icomp[5]=5;
+	ifield[6]=1;icomp[6]=0;
+	nfield[0]=7;nfield[1]=7;
+
+	exoselect(eenmax,eenmax,&iset,&nkcoords,inum,istartset,iendset,
+		  ialset,ngraph,&ncomp,ifield,icomp,
+		  nfield,&iselect,exoid,num_time_steps,countvars);
+	printf ("Warning: export of MAXE to exo not tested.\n");
+	countvars+=7;
+      }
+    }
+
+  WRITENAMES:
+    if (countbool==3){
+      errr = ex_put_variable_param (exoid, EX_NODAL, countvars);
+      ex_update (exoid);
+      // var_names = (char *) calloc (countvars, sizeof (char));
+      // var_names = (char *) calloc (countvars, MAX_STR_LENGTH);
+      int ii=0; // 4 bit integer
+      for (ii=0; ii<countvars; ii++)
+	{
+	  var_names[ii] = (char *) calloc ((MAX_STR_LENGTH+1), sizeof(char));
+	}
+    }else if(countbool==2){
+      errr = ex_put_variable_names (exoid, EX_NODAL, countvars, var_names);
+      if (errr) {
+	printf ("Unable to update variable names.  Was output data requested?\n\n");
+	printf ("  NOTE CalculiX Extras doesn't currently support\n");
+	printf ("  changing output variables between steps.\n");
+	printf ("  Requested output variables in this step are:\n");
+	int ii=0; // 4 bit integer
+	for (ii=0; ii<countvars; ii++){
+	  printf("    %i %s\n", ii+1, var_names[ii]);
+	}
+      }
+
+      //  for (i=0; i<countvars; i++)
+      // 	{
+      // 	  free(var_names[i]);
+      // 	}
+      ex_update (exoid);
+    }
+
+    countvars=0;
+    // printf ("Decrement countbool to %i\n",--countbool);
+    --countbool;
+  };
+
+  // free (var_names);
+  free (node_map_inv);
+  ex_update (exoid);
+  ex_close(exoid);
+
+#else
+  printf ("Warning: requested exodus output support not included at compile time.\n");
+#endif
+  return;
+}
diff --git a/ccx/src/exo.h b/ccx/src/exo.h
new file mode 100644
index 0000000..07ea546
--- /dev/null
+++ b/ccx/src/exo.h
@@ -0,0 +1,52 @@
+/*     Calculix - A 3-dimensional finite element program                 */
+/*              Copyright (C) 1998-2018 Guido Dhondt                     */
+/*     This subroutine                                                   */
+/*              Copyright (C) 2013-2021 Peter A. Gustafson               */
+/*                                                                       */
+/*     This program is free software; you can redistribute it and/or     */
+/*     modify it under the terms of the GNU General Public License as    */
+/*     published by the Free Software Foundation(version 2);    */
+/*                                                                       */
+
+/*     This program is distributed in the hope that it will be useful,   */
+/*     but WITHOUT ANY WARRANTY; without even the implied warranty of    */
+/*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the      */
+/*     GNU General Public License for more details.                      */
+
+/*     You should have received a copy of the GNU General Public License */
+/*     along with this program; if not, write to the Free Software       */
+/*     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.         */
+
+void exo(double *co,ITG *nk,ITG *kon,ITG *ipkon,char *lakon,ITG *ne0,
+	 double *v,double *stn,ITG *inum,ITG *nmethod,ITG *kode,
+	 char *filab,double *een,double *t1,double *fn,double *time,
+	 double *epn,ITG *ielmat,char *matname,double *enern,
+	 double *xstaten,ITG *nstate_,ITG *istep,ITG *iinc,
+	 ITG *ithermal,double *qfn,ITG *mode,ITG *noddiam,
+	 double *trab,ITG *inotr,ITG *ntrans,double *orab,
+	 ITG *ielorien,ITG *norien,char *description,ITG *ipneigh,
+	 ITG *neigh,ITG *mi,double *stx,double *vr,double *vi,
+	 double *stnr,double *stni,double *vmax,double *stnmax,
+	 ITG *ngraph,double *veold,double *ener,ITG *ne,double *cs,
+	 char *set,ITG *nset,ITG *istartset,ITG *iendset,ITG *ialset,
+	 double *eenmax,double *fnr,double *fni,double *emn,
+	 double *thicke,char *jobnamec,char *output,double *qfx,
+         double *cdn,ITG *mortar,double *cdnr,double *cdni,ITG *nmat,
+	 ITG *ielprop,double *prop);
+
+void exosetfind(char *set, ITG *nset, ITG *ialset, ITG *istartset, ITG *iendset,
+		ITG *num_ns, ITG *num_ss, ITG *num_es, ITG *num_fs, ITG *node_map_inv,
+		int exoid, int store, ITG *nk);
+
+ITG exoset_check(ITG n, ITG *node_map_inv, ITG *nk, int *dropped, int *unidentified);
+
+void exovector(double *v,ITG *iset,ITG *ntrans,char * filabl,ITG *nkcoords,
+               ITG *inum,char *m1,ITG *inotr,double *trab,double *co,
+               ITG *istartset,ITG *iendset,ITG *ialset,ITG *mi,ITG *ngraph,
+               FILE *f1,char *output,char *m3, int exoid, ITG time_step,
+	       int countvar);
+
+void exoselect(double *field1,double *field2,ITG *iset,ITG *nkcoords,ITG *inum,
+	       ITG *istartset,ITG *iendset,ITG *ialset,ITG *ngraph,ITG *ncomp,
+	       ITG *ifield,ITG *icomp,ITG *nfield,ITG *iselect,ITG exoid,
+	       ITG time_step, int countvar);
diff --git a/ccx/src/exoselect.c b/ccx/src/exoselect.c
new file mode 100644
index 0000000..2da54db
--- /dev/null
+++ b/ccx/src/exoselect.c
@@ -0,0 +1,192 @@
+/*     CalculiX - A 3-dimensional finite element program                 */
+/*              Copyright (C) 1998-2018 Guido Dhondt                     */
+/*     This subroutine                                                   */
+/*              Copyright (C) 2013-2021 Peter A. Gustafson               */
+/*                                                                       */
+/*     This program is free software; you can redistribute it and/or     */
+/*     modify it under the terms of the GNU General Public License as    */
+/*     published by the Free Software Foundation(version 2);             */
+/*                                                                       */
+
+/*     This program is distributed in the hope that it will be useful,   */
+/*     but WITHOUT ANY WARRANTY; without even the implied warranty of    */
+/*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the      */
+/*     GNU General Public License for more details.                      */
+
+/*     You should have received a copy of the GNU General Public License */
+/*     along with this program; if not, write to the Free Software       */
+/*     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.         */
+
+#ifdef EXODUSII
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+#include "CalculiX.h"
+#include "exodusII.h"
+#include "exo.h"
+
+#define min(a,b) ((a) <= (b) ? (a) : (b))
+#define max(a,b) ((a) >= (b) ? (a) : (b))
+
+void exoselect(double *field1,double *field2,ITG *iset,ITG *nkcoords,ITG *inum,
+	       ITG *istartset,ITG *iendset,ITG *ialset,ITG *ngraph,ITG *ncomp,
+	       ITG *ifield,ITG *icomp,ITG *nfield,ITG *iselect,ITG exoid,
+	       ITG time_step, int countvar){
+
+  /* storing scalars, components of vectors and tensors without additional
+     transformations */
+
+  /* number of components in field1: nfield[0]
+     number of components in field2: nfield[1]
+
+     number of entities to store: ncomp
+     for each entity i, 0<=i<ncomp:
+     - ifield[i]: 1=field1,2=field2
+     - icomp[i]: component: 0...,(nfield[0]-1 or nfield[1]-1) */
+
+  ITG i,j,k,l,m,nksegment;
+
+  /* When initializing parameter values:
+     "g" (or "G")
+     For global variables.
+     "n" (or "N")
+     For nodal variables.
+     "e" (or "E")
+     For element variables.
+     "m" (or "M")
+     For nodeset variables.
+     "s" (or "S")
+     For sideset variables.
+  */
+  
+  // Get required info out of the exodus file.  Most important number of nodes stored.
+  int num_dim, num_nodes, num_elem, num_elem_blk, num_node_sets, num_side_sets, errr;
+  char title [ MAX_LINE_LENGTH +1];
+  errr = ex_get_init (exoid, title, &num_dim, &num_nodes, &num_elem, &num_elem_blk, &num_node_sets, &num_side_sets);
+  
+  float *nodal_var_vals = (float *) calloc(num_nodes, sizeof(float));
+  ITG   *node_map       = (ITG *)   calloc(num_nodes, sizeof(ITG));
+  ITG   *inode_map      = (ITG *)   calloc(*nkcoords, sizeof(ITG));
+  
+  errr = ex_get_id_map (exoid, EX_NODE_MAP, node_map);
+  if(errr)printf("*ERROR in exo: failed to get prior node map");
+
+  for (i=0; i<num_nodes; i++){
+    // Seed with NaN for non-stored output
+    nodal_var_vals[i]=nanf("");
+  }
+  
+  {j=0;
+    // Pull node map and create inverse
+    // WHY DO WE NEED A INVERSE NODE MAP?
+    // FRD doesn't require results to be sequential be we do!
+    // FRD can be out of sequence (based on set order)
+    // FRD tags each result with the associated node number.
+    // Hence we need inverse and FRD does not
+    // Note also that changes in output sets, element deletion, etc
+    // means the initial from exo.c must be recreated each time
+    for(i=0;i<*nkcoords;i++){
+      if(inum[i]<=0) continue;
+      while(node_map[j%*nkcoords]!=i+1){j++;}
+      inode_map[i] = (j%*nkcoords);
+    }
+  }
+
+  for(j=0;j<*ncomp;j++){
+    if(*iset==0){
+      for(i=0;i<*nkcoords;i++){
+
+	/* check whether output is requested for solid nodes or
+	   network nodes */
+
+	if(*iselect==1){
+	  if(inum[i]<=0) continue;
+	}else if(*iselect==-1){
+	  if(inum[i]>=0) continue;
+	}else{
+	  if(inum[i]==0) continue;
+	}
+
+	/* storing the entities */
+	if(ifield[j]==1){
+	  // It is getting here and assigning -nan
+	  nodal_var_vals[inode_map[i]]=field1[i*nfield[0]+icomp[j]];
+	}else{
+	  nodal_var_vals[inode_map[i]]=field2[i*nfield[1]+icomp[j]];
+	}
+      }
+    }else{
+      nksegment=(*nkcoords)/(*ngraph);
+      for(k=istartset[*iset-1]-1;k<iendset[*iset-1];k++){
+	if(ialset[k]>0){
+	
+	  for(l=0;l<*ngraph;l++){
+	    i=ialset[k]+l*nksegment-1;
+	
+	    /* check whether output is requested for solid nodes or
+	       network nodes */
+	
+	    if(*iselect==1){
+	      if(inum[i]<=0) continue;
+	    }else if(*iselect==-1){
+	      if(inum[i]>=0) continue;
+	    }else{
+	      if(inum[i]==0) continue;
+	    }
+	
+	    /* storing the entities */
+	    if(ifield[j]==1){
+	      nodal_var_vals[inode_map[i]]=field1[i*nfield[0]+icomp[j]];
+	    }else{
+	      nodal_var_vals[inode_map[i]]=field2[i*nfield[1]+icomp[j]];
+	    }
+	  }
+	}else{
+	  l=ialset[k-2];
+	  do{
+	    l-=ialset[k];
+	    if(l>=ialset[k-1]) break;
+	    for(m=0;m<*ngraph;m++){
+	      i=l+m*nksegment-1;
+	
+	      /* check whether output is requested for solid nodes or
+		 network nodes */
+	
+	      if(*iselect==1){
+		if(inum[i]<=0) continue;
+	      }else if(*iselect==-1){
+		if(inum[i]>=0) continue;
+	      }else{
+		if(inum[i]==0) continue;
+	      }
+	
+	      /* storing the entities */
+	      if(ifield[j]==1){
+		nodal_var_vals[inode_map[i]]=field1[i*nfield[0]+icomp[j]];
+	      }else{
+		nodal_var_vals[inode_map[i]]=field2[i*nfield[1]+icomp[j]];
+	      }
+	    }
+	  }while(1);
+	}
+      }
+    }
+
+    // Note: the results for exo tensor format are strangely ordered,
+    if (*ncomp==12){
+      k=j;
+    }else{k=j;}
+
+    int errr = ex_put_var (exoid, time_step, EX_NODAL, k+1+countvar, 1, num_nodes, nodal_var_vals);
+    if (errr) printf ("ERROR exoselect data for dim %i record %i.\n", j, k+1+countvar);
+  }
+
+  free(node_map);
+  free(inode_map);
+  free(nodal_var_vals);
+  return;
+
+}
+
+#endif
diff --git a/ccx/src/exosetfind.c b/ccx/src/exosetfind.c
new file mode 100644
index 0000000..dc7639b
--- /dev/null
+++ b/ccx/src/exosetfind.c
@@ -0,0 +1,205 @@
+/*     CalculiX - A 3-dimensional finite element program                 */
+/*              Copyright (C) 1998-2018 Guido Dhondt                     */
+/*     This subroutine                                                   */
+/*              Copyright (C) 2013-2021 Peter A. Gustafson               */
+/*                                                                       */
+/*     This program is free software; you can redistribute it and/or     */
+/*     modify it under the terms of the GNU General Public License as    */
+/*     published by the Free Software Foundation(version 2);             */
+/*                                                                       */
+
+/*     This program is distributed in the hope that it will be useful,   */
+/*     but WITHOUT ANY WARRANTY; without even the implied warranty of    */
+/*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the      */
+/*     GNU General Public License for more details.                      */
+
+/*     You should have received a copy of the GNU General Public License */
+/*     along with this program; if not, write to the Free Software       */
+/*     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.         */
+
+#ifdef EXODUSII
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+#include "CalculiX.h"
+#include "exodusII.h"
+#include "exo.h"
+
+void exosetfind(char *set, ITG *nset, ITG *ialset, ITG *istartset, ITG *iendset,
+		ITG *num_ns, ITG *num_ss, ITG *num_es, ITG *num_fs, ITG *node_map_inv,
+		int exoid, int store, ITG *nk){
+
+  ITG i,j,k,l,n,s,e,gen;
+  char tmpstr[81];
+  char *space = " ";
+  char *pos;
+  ITG *set_nums;
+  int errr;
+  int dropped=0, unidentified=0;
+
+  *num_ns = 0;
+  *num_ss = 0;
+  *num_es = 0;
+  *num_fs = 0;
+
+  for (i=0; i<*nset; i++){
+
+    if (store) {
+      // Find and store the set numbers
+      // The pointer integers are 1 based (fortran)
+      s=istartset[i]-1;
+      e=iendset[i]-1;
+      // Determine if generate was used
+      gen=0; l=0;
+      for (j=s; j<=e; j++){
+	if (ialset[j]<0) {
+	  gen+=(ialset[j-1]-ialset[j-2])/(-ialset[j])+1;
+	  l-=3;
+	}
+      }
+
+      // Now set the length of the set allocation
+      l=e-s+1+gen+l;
+      set_nums = (ITG *) calloc(l, sizeof(ITG));
+
+      /* Only add the generate code if there are at least
+	 three points in the vector */
+      n=0; j=s;
+      if (l>2){
+	while (j<=e-2){
+	  // Account for generated ids
+	  if (ialset[j+2]<0) {
+	    for (k=ialset[j]; k<=ialset[j+1]; k-=ialset[j+2]){
+	      set_nums[n++]=exoset_check(k-1, node_map_inv, nk, &dropped, &unidentified);
+	    }
+	    j+=3;
+	  }else{
+	    // Account for directly added id
+	    gen=ialset[j++]-1;
+	    set_nums[n++]=exoset_check(gen, node_map_inv, nk, &dropped, &unidentified);
+	  }
+	}
+	// Must finish the last two of directly added set
+	if (ialset[e]>0){ // only if the last set is not a generated set
+	  // 1+n++ and -1+n++ to preserve order
+	  set_nums[1+n++]=exoset_check(ialset[e]-2, node_map_inv, nk, &dropped, &unidentified);
+	  if (ialset[e-1]>0){
+	    set_nums[-1+n++]=exoset_check(ialset[e]-3, node_map_inv, nk, &dropped, &unidentified);
+	  }
+	}
+      }else if(l>1){
+	set_nums[n++]=exoset_check(ialset[s]-1, node_map_inv, nk, &dropped, &unidentified);
+	set_nums[n++]=exoset_check(ialset[e]-1, node_map_inv, nk, &dropped, &unidentified);
+      }else{
+	set_nums[n++]=exoset_check(ialset[e]-1, node_map_inv, nk, &dropped, &unidentified);
+      }
+    } // if (store)
+
+    strncpy(tmpstr,set+i*81,81);
+    pos = strpbrk(tmpstr, space)-1;
+
+    if(strcmp1(pos,"N")==0){
+      (*num_ns)++; // printf ("Node set identified\n");
+      if (store){
+	errr = ex_put_set_param (exoid, EX_NODE_SET, i, n, 0); // CURRENTLY NO DISTRIBUTIONS ADDED
+	if (errr) printf ("Error writing set parameters\n");
+	errr = ex_put_set       (exoid, EX_NODE_SET, i, set_nums, NULL);
+	if (errr) printf ("Error writing set numbers\n");
+	// ex_put_set_dist_fact (exoid, EX_NODE_SET, i, set_nums);  //
+      }
+    }
+    if(strcmp1(pos,"E")==0) {
+      (*num_es)++; // printf ("Element set identified\n");}
+      /* No element set storage mechanism?
+      if (store){
+      	errr = ex_put_set_param (exoid, EX_ELEM_SET, i, n, 0);
+      	if (errr) printf ("Error writing set parameters\n");
+      	errr = ex_put_set       (exoid, EX_ELEM_SET, i, set_nums, NULL);
+      	if (errr) printf ("Error writing set numbers\n");
+      	// ex_put_set_dist_fact (exoid, EX_ELEM_SET, i, set_nums);  //
+	} */
+    }
+    if(strcmp1(pos,"S")==0) {
+      (*num_ss)++; // printf ("Node side set surface identified\n");}
+      /* Side sets (node surfaces) not yet implemented
+	 if (store){
+      	errr = ex_put_set_param (exoid, EX_SIDE_SET, i, n, 0);
+      	if (errr) printf ("Error writing set parameters\n");
+      	errr = ex_put_set       (exoid, EX_SIDE_SET, i, set_nums, side_set_side_list);
+      	if (errr) printf ("Error writing set numbers\n");
+      	// ex_put_set_dist_fact (exoid, EX_SIDE_SET, i, set_nums);  //
+	} */
+    }
+    if(strcmp1(pos,"T")==0) {
+      (*num_fs)++; // printf ("Face set surface identified\n");}
+      /* Face sets not yet implemented
+	 if (store){
+	 errr = ex_put_set_param (exoid, EX_FACE_SET, i, n, 0);
+	 if (errr) printf ("Error writing set parameters\n");
+	 errr = ex_put_set       (exoid, EX_FACE_SET, i, set_nums);
+	 if (errr) printf ("Error writing set numbers\n");
+	 // ex_put_set_dist_fact (exoid, EX_FACE_SET, i, set_nums);  //
+	 } */
+    }
+    if (store) {free (set_nums);}
+  }
+
+  if (store){
+    // char *namesnset[*num_ns]; j=0;
+    // char *namessset[*num_ss]; k=0;
+    // char *nameseset[*num_es]; l=0;
+    // char *namesfset[*num_fs]; n=0;
+
+    for (i=0; i<*nset; i++){
+      strncpy(tmpstr,set+i*81,81);
+      pos = strpbrk(tmpstr, space)-1;
+      // This crashed in valgrind
+      // if(strcmp1(pos,"N")==0) {strncpy(namesnset[j++], tmpstr, MAX_STR_LENGTH);}
+      // if(strcmp1(pos,"E")==0) {strcpy(nameseset[l++],tmpstr);}
+      // if(strcmp1(pos,"S")==0) {strcpy(namessset[k++],tmpstr);}
+      // if(strcmp1(pos,"T")==0) {strcpy(namesfset[n++],tmpstr);}
+    }
+
+    // if (*num_ns>0){
+    //   errr = ex_put_names (exoid, EX_NODE_SET, namesnset);
+    //   if (errr) printf ("Error writing node set names\n");
+    // }
+    /* side sets not implemented yet
+       if (*num_ss>0){
+       errr = ex_put_names (exoid, EX_SIDE_SET, namessset);
+       if (errr) printf ("Error writing side set names\n");
+       }
+    */
+  }
+
+  if (dropped){
+    printf ("\nExodus Output WARNING: At least one node or element is dropped from a set.\n");
+    printf ("  This may be due rigid bodies or 3D expansion (beams, shells, OUTPUT=3D).\n\n");
+  }
+
+  if (unidentified){
+    printf ("\nExodus Output WARNING: At least one unidentified node or element is dropped from a set.\n\n");
+  }
+
+  return;
+}
+
+
+ITG exoset_check(ITG n, ITG *node_map_inv, ITG *nk, int *dropped, int *unidentified){
+  ITG val=0;
+  // printf ("%" ITGFORMAT ", %" ITGFORMAT "\n", n, *nk);
+  if (n<=*nk){
+    val = node_map_inv[n]-1;
+    // printf ("val node number %" ITGFORMAT "\n", val+1);
+    if (val==-1) {
+      *dropped = 1;
+    }
+  }else{
+    // printf("UNIDENTIFIED %" ITGFORMAT "\n", n);
+    *unidentified = 1;
+  }
+  return val;
+}
+
+#endif
diff --git a/ccx/src/exovector.c b/ccx/src/exovector.c
new file mode 100644
index 0000000..41260a3
--- /dev/null
+++ b/ccx/src/exovector.c
@@ -0,0 +1,161 @@
+/*     CalculiX - A 3-dimensional finite element program                 */
+/*              Copyright (C) 1998-2020 Guido Dhondt                     */
+/*     This subroutine                                                   */
+/*              Copyright (C) 2013-2021 Peter A. Gustafson               */
+/*                                                                       */
+/*     This program is free software; you can redistribute it and/or     */
+/*     modify it under the terms of the GNU General Public License as    */
+/*     published by the Free Software Foundation(version 2);             */
+/*                                                                       */
+
+/*     This program is distributed in the hope that it will be useful,   */
+/*     but WITHOUT ANY WARRANTY; without even the implied warranty of    */
+/*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the      */
+/*     GNU General Public License for more details.                      */
+
+/*     You should have received a copy of the GNU General Public License */
+/*     along with this program; if not, write to the Free Software       */
+/*     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.         */
+
+#ifdef EXODUSII
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+#include "CalculiX.h"
+#include "exodusII.h"
+#include "exo.h"
+
+void exovector(double *v,ITG *iset,ITG *ntrans,char * filabl,ITG *nkcoords,
+               ITG *inum,char *m1,ITG *inotr,double *trab,double *co,
+               ITG *istartset,ITG *iendset,ITG *ialset,ITG *mi,ITG *ngraph,
+               FILE *f1,char *output,char *m3, int exoid, ITG time_step,
+	       int countvar){
+
+  ITG nksegment;
+  ITG i,j,k,l,m,ii,jj,kk;
+
+  double a[9];
+
+  /* When initializing parameter values:
+     "g" (or "G")
+     For global variables.
+     "n" (or "N")
+     For nodal variables.
+     "e" (or "E")
+     For element variables.
+     "m" (or "M")
+     For nodeset variables.
+     "s" (or "S")
+     For sideset variables.
+  */
+
+  // Get required info out of the exodus file.  Most important number of nodes stored.
+  int num_dim, num_nodes, num_elem, num_elem_blk, num_node_sets, num_side_sets, errr;
+  char title [ MAX_LINE_LENGTH +1];
+  errr = ex_get_init (exoid, title, &num_dim, &num_nodes, &num_elem, &num_elem_blk, &num_node_sets, &num_side_sets);
+  
+  float *nodal_var_vals = (float *) calloc(num_nodes, sizeof(float));
+  ITG   *node_map       = (ITG *)   calloc(num_nodes, sizeof(ITG));
+  ITG   *inode_map      = (ITG *)   calloc(*nkcoords, sizeof(ITG));
+  
+  errr = ex_get_id_map (exoid, EX_NODE_MAP, node_map);
+  if(errr)printf("*ERROR in exo: failed to get prior node map");
+
+  for (i=0; i<num_nodes; i++){
+    // Seed with NaN for non-stored output
+    nodal_var_vals[i]=nanf("");
+  }
+  
+  {j=0;
+    // Pull node map and create inverse
+    // WHY DO WE NEED A INVERSE NODE MAP?
+    // FRD doesn't require results to be sequential be we do!
+    // FRD can be out of sequence (based on set order)
+    // FRD tags each result with the associated node number.
+    // Hence we need inverse and FRD does not.
+    // Note also that changes in output sets, element deletion, etc
+    // means the initial from exo.c must be recreated each time
+    for(i=0;i<*nkcoords;i++){
+      if(inum[i]<=0) continue;
+      while(node_map[j%*nkcoords]!=i+1){j++;}
+      inode_map[i] = (j%*nkcoords);
+    }
+  }
+  
+  int num_nod_vars=3;
+  
+  for (j=1; j<=num_nod_vars; j++){ // For each direction
+    if(*iset==0){
+      if((*ntrans==0)||(strcmp1(&filabl[5],"G")==0)){
+	for(i=0;i<*nkcoords;i++){
+	  if(inum[i]<=0) continue;
+	  nodal_var_vals[inode_map[i]]=v[(mi[1]+1)*(i)+j];
+	}
+      }else{
+	for(i=0;i<*nkcoords;i++){
+	  if(inum[i]<=0) continue;
+
+	  if(inotr[2*i]==0){
+	    nodal_var_vals[inode_map[i]]=v[(mi[1]+1)*i+j];
+	  }else{
+	    ii=(mi[1]+1)*i+1;
+	    jj=(mi[1]+1)*i+2;
+	    kk=(mi[1]+1)*i+3;
+	    FORTRAN(transformatrix,(&trab[7*(inotr[2*i]-1)],&co[3*i],a));
+	    nodal_var_vals[inode_map[i]]=v[ii]*a[0+(j-1)*3]+v[jj]*a[1+(j-1)*3]+v[kk]*a[2+(j-1)*3];
+	  }
+	}
+      }
+    }else{
+      nksegment=(*nkcoords)/(*ngraph);
+      for(k=istartset[*iset-1]-1;k<iendset[*iset-1];k++){
+	if(ialset[k]>0){
+	  for(l=0;l<*ngraph;l++){
+	    i=ialset[k]+l*nksegment-1;
+	    if(inum[i]<=0) continue;
+	    if((*ntrans==0)||(strcmp1(&filabl[5],"G")==0)||(inotr[2*i]==0)){
+	      nodal_var_vals[inode_map[i]]=v[(mi[1]+1)*i+j];
+	    }else{
+	      // Transform into global directions
+	      FORTRAN(transformatrix,(&trab[7*(inotr[2*i]-1)],&co[3*i],a));
+	      ii=(mi[1]+1)*i+1;
+	      jj=(mi[1]+1)*i+2;
+	      kk=(mi[1]+1)*i+3;
+	      nodal_var_vals[inode_map[i]]=v[ii]*a[0+(j-1)*3]+v[jj]*a[1+(j-1)*3]+v[kk]*a[2+(j-1)*3];
+	    }
+	  }
+	}else{
+	  l=ialset[k-2];
+	  do{
+	    l-=ialset[k];
+	    if(l>=ialset[k-1]) break;
+	    for(m=0;m<*ngraph;m++){
+	      i=l+m*nksegment-1;
+	      if(inum[i]<=0) continue;
+	      if((*ntrans==0)||(strcmp1(&filabl[5],"G")==0)||(inotr[2*i]==0)){
+		nodal_var_vals[inode_map[i]]=v[(mi[1]+1)*i+j];
+	      }else{
+		FORTRAN(transformatrix,(&trab[7*(inotr[2*i]-1)],&co[3*i],a));
+		ii=(mi[1]+1)*i+1;
+		jj=(mi[1]+1)*i+2;
+		kk=(mi[1]+1)*i+3;
+		nodal_var_vals[inode_map[i]]=v[ii]*a[0+(j-1)*3]+v[jj]*a[1+(j-1)*3]+v[kk]*a[2+(j-1)*3];
+	      }
+	    }
+	  }while(1);
+	}
+      }
+    }
+
+    errr = ex_put_var (exoid, time_step, EX_NODAL, j+countvar, 1, num_nodes, nodal_var_vals);
+    if (errr) printf ("ERROR storing vector data into exo file.\n");
+  }
+    
+  free(node_map);
+  free(inode_map);
+  free(nodal_var_vals);
+  return;
+
+}
+#endif
diff --git a/ccx/src/extrapolate.f b/ccx/src/extrapolate.f
index 46a3aa0..fba1e14 100644
--- a/ccx/src/extrapolate.f
+++ b/ccx/src/extrapolate.f
@@ -102,11 +102,20 @@
 !
       data a2 /1.1455d0,-0.1455d0,1.1455d0,-0.1455d0,1.1455d0,-0.1455d0,
      &        -0.1455d0,1.1455d0,-0.1455d0,1.1455d0,-0.1455d0,1.1455d0/
-      data a4 /  1.92705d0, -0.30902d0, -0.30902d0, -0.30902d0,
-     &          -0.30902d0,  1.92705d0, -0.30902d0, -0.30902d0,
-     &          -0.30902d0, -0.30902d0,  1.92705d0, -0.30902d0,
-     &          -0.30902d0, -0.30902d0, -0.30902d0,  1.92705d0/
-!
+
+!     Precision increased on these 8/2021 Pete Gustafson (after a course
+!     excercise demonstrating extrapolation). Current precision is
+!     ridiculous, but comes for free.
+      data a4 /
+     & +1.9270509831248423537886083d+0, -3.0901699437494745126286944d-1,
+     & -3.0901699437494745126286944d-1, -3.0901699437494745126286944d-1,
+     & -3.0901699437494745126286944d-1, +1.9270509831248423537886083d+0,
+     & -3.0901699437494745126286944d-1, -3.0901699437494745126286944d-1,
+     & -3.0901699437494745126286944d-1, -3.0901699437494745126286944d-1,
+     & +1.9270509831248423537886083d+0, -3.0901699437494745126286944d-1,
+     & -3.0901699437494745126286944d-1, -3.0901699437494745126286944d-1,
+     & -3.0901699437494745126286944d-1, +1.9270509831248423537886083d+0/
+!     
 !     extrapolation from a 6 integration point scheme in a wedge to
 !     the vertex nodes
 !
@@ -143,19 +152,44 @@
      &          1.63138d0/
 !
 !     extrapolation from a 2x2x2=8 integration point scheme in a hex to
-!     the vertex nodes
-!    
-      data a8 /2.549d0,-.683d0,.183d0,-.683d0,-.683d0,.183d0,
-     &        -.04904d0,.183d0,-.683d0,2.549d0,-.683d0,.183d0,
-     &        .183d0,-.683d0,.183d0,-.04904d0,-.683d0,.183d0,
-     &        -.683d0,2.549d0,.183d0,-.04904d0,.183d0,-.683d0,
-     &        .183d0,-.683d0,2.549d0,-.683d0,-.04904d0,.183d0,
-     &        -.683d0,.183d0,-.683d0,.183d0,-.04904d0,.183d0,
-     &        2.549d0,-.683d0,.183d0,-.683d0,.183d0,-.683d0,
-     &        .183d0,-.04904d0,-.683d0,2.549d0,-.683d0,.183d0,
-     &        .183d0,-.04904d0,.183d0,-.683d0,-.683d0,.183d0,
-     &        -.683d0,2.549d0,-.04904d0,.183d0,-.683d0,.183d0,
-     &        .183d0,-.683d0,2.549d0,-.683d0/  
+!     the vertex nodes.  Precision increased on these 8/2021 Pete
+!     Gustafson (after a course excercise demonstrating
+!     extrapolation). Current precision is ridiculous, but comes for
+!     free.
+!     
+      data a8 /
+     & +2.5490381056766566736371260d+0, -6.8301270189221807704882394d-1,
+     & -6.8301270189221940931645349d-1, +1.8301270189221954809433157d-1,
+     & -6.8301270189221896522724364d-1, +1.8301270189221824358227764d-1,
+     & +1.8301270189221996442796581d-1, -4.9038105676658283460511711d-2,
+     & -6.8301270189221874318263872d-1, +2.5490381056766580059047556d+0,
+     & +1.8301270189221963136105842d-1, -6.8301270189221918727184857d-1,
+     & +1.8301270189221860440476064d-1, -6.8301270189221807704882394d-1,
+     & -4.9038105676659414500218048d-2, +1.8301270189221935380530226d-1,
+     & +1.8301270189222007545026827d-1, -6.8301270189221818807112641d-1,
+     & -6.8301270189221907624954611d-1, +2.5490381056766580059047556d+0,
+     & -4.9038105676658678977464234d-2, +1.8301270189221757744846286d-1,
+     & +1.8301270189221979789451211d-1, -6.8301270189221829909342887d-1,
+     & -6.8301270189221874318263872d-1, +1.8301270189221782724864340d-1,
+     & +2.5490381056766575618155457d+0, -6.8301270189222018647257073d-1,
+     & +1.8301270189222024198372196d-1, -4.9038105676655389941753782d-2,
+     & -6.8301270189222040851717566d-1, +1.8301270189221932604972665d-1,
+     & -6.8301270189221874318263872d-1, +1.8301270189221952033875596d-1,
+     & +1.8301270189221940931645349d-1, -4.9038105676658602649631291d-2,
+     & +2.5490381056766566736371260d+0, -6.8301270189221907624954611d-1,
+     & -6.8301270189221918727184857d-1, +1.8301270189221968687220965d-1,
+     & +1.8301270189221968687220965d-1, -6.8301270189221963136105842d-1,
+     & -4.9038105676657818554620150d-2, +1.8301270189221965911663403d-1,
+     & -6.8301270189221818807112641d-1, +2.5490381056766553413694965d+0,
+     & +1.8301270189222118567329289d-1, -6.8301270189221863216033626d-1,
+     & -4.9038105676657922638028708d-2, +1.8301270189221990891681457d-1,
+     & +1.8301270189221979789451211d-1, -6.8301270189222029749487319d-1,
+     & +1.8301270189221993667239019d-1, -6.8301270189221741091500917d-1,
+     & -6.8301270189222107465099043d-1, +2.5490381056766566736371260d+0,
+     & +1.8301270189221985340566334d-1, -4.9038105676658318154981231d-2,
+     & -6.8301270189221918727184857d-1, +1.8301270189222004769469265d-1,
+     & -6.8301270189221852113803379d-1, +1.8301270189221771622634094d-1,
+     & +2.5490381056766584499939654d+0, -6.8301270189221885420494118d-1/
 !
 !     extrapolation from a 14 integration point scheme in a hex to
 !     the vertex nodes
@@ -191,164 +225,282 @@
      &  0.2069d+00,0.2069d+00,0.2069d+00,0.2069d+00/
 !
 !     extrapolation from a 3x3x3=27 integration point scheme in a hex to
-!     the all nodes in a 20-node element
+!     the all nodes in a 20-node element. Precision increased on these
+!     8/2021 Pete Gustafson (after a course excercise demonstrating
+!     extrapolation). Current precision is ridiculous, but comes for
+!     free.
 !    
       data a27 /
-     &  2.37499d0,-0.12559d0,-0.16145d0,-0.12559d0,-0.12559d0,
-     & -0.16145d0, 0.11575d0,
-     & -0.16145d0, 0.32628d0, 0.11111d0, 0.11111d0, 0.32628d0,
-     &  0.11111d0,-0.10405d0,
-     & -0.10405d0, 0.11111d0, 0.32628d0, 0.11111d0,-0.10405d0,
-     &  0.11111d0,-0.31246d0,
-     & -0.31246d0, 0.31481d0, 0.31481d0, 0.31481d0, 0.31481d0,
-     & -0.16902d0,-0.16902d0,
-     &  1.28439d0,-0.27072d0,-0.19444d0,-0.27072d0,-0.19444d0,
-     &  0.15961d0,-0.00661d0,
-     &  0.15961d0,-0.27072d0,-0.27072d0, 0.15961d0, 0.15961d0,
-     & -0.12559d0, 2.37499d0,
-     & -0.12559d0,-0.16145d0,-0.16145d0,-0.12559d0,-0.16145d0,
-     &  0.11575d0, 0.32628d0,
-     &  0.32628d0, 0.11111d0, 0.11111d0, 0.11111d0, 0.11111d0,
-     & -0.10405d0,-0.10405d0,
-     &  0.11111d0, 0.32628d0, 0.11111d0,-0.10405d0,-0.31246d0,
-     &  0.31481d0, 0.31481d0,
-     & -0.31246d0, 0.31481d0,-0.16902d0,-0.16902d0, 0.31481d0,
-     & -0.27072d0,-0.19444d0,
-     & -0.27072d0, 1.28439d0, 0.15961d0,-0.00661d0, 0.15961d0,
-     & -0.19444d0,-0.27072d0,
-     &  0.15961d0, 0.15961d0,-0.27072d0,-0.48824d0,-0.48824d0,
-     & -0.48824d0,-0.48824d0,
-     &  0.22898d0, 0.22898d0, 0.22898d0, 0.22898d0, 0.05556d0,
-     &  0.05556d0, 0.05556d0,
-     &  0.05556d0, 0.05556d0, 0.05556d0, 0.05556d0, 0.05556d0,
-     & -0.22222d0,-0.22222d0,
-     & -0.22222d0,-0.22222d0, 0.31481d0,-0.31246d0,-0.31246d0,
-     &  0.31481d0,-0.16902d0,
-     &  0.31481d0, 0.31481d0,-0.16902d0,-0.27072d0, 1.28439d0,
-     & -0.27072d0,-0.19444d0,
-     &  0.15961d0,-0.19444d0, 0.15961d0,-0.00661d0, 0.15961d0,
-     & -0.27072d0,-0.27072d0,
-     &  0.15961d0,-0.12559d0,-0.16145d0,-0.12559d0, 2.37499d0,
-     & -0.16145d0, 0.11575d0,
-     & -0.16145d0,-0.12559d0, 0.11111d0, 0.11111d0, 0.32628d0,
-     &  0.32628d0,-0.10405d0,
-     & -0.10405d0, 0.11111d0, 0.11111d0, 0.11111d0,-0.10405d0,
-     &  0.11111d0, 0.32628d0,
-     &  0.31481d0, 0.31481d0,-0.31246d0,-0.31246d0,-0.16902d0,
-     & -0.16902d0, 0.31481d0,
-     &  0.31481d0,-0.19444d0,-0.27072d0, 1.28439d0,-0.27072d0,
-     & -0.00661d0, 0.15961d0,
-     & -0.19444d0, 0.15961d0, 0.15961d0, 0.15961d0,-0.27072d0,
-     & -0.27072d0,-0.16145d0,
-     & -0.12559d0, 2.37499d0,-0.12559d0, 0.11575d0,-0.16145d0,
-     & -0.12559d0,-0.16145d0,
-     &  0.11111d0, 0.32628d0, 0.32628d0, 0.11111d0,-0.10405d0,
-     &  0.11111d0, 0.11111d0,
-     & -0.10405d0,-0.10405d0, 0.11111d0, 0.32628d0, 0.11111d0,
-     & -0.31246d0, 0.31481d0,
-     & -0.16902d0, 0.31481d0,-0.31246d0, 0.31481d0,-0.16902d0,
-     &  0.31481d0,-0.27072d0,
-     &  0.15961d0, 0.15961d0,-0.27072d0,-0.27072d0, 0.15961d0,
-     &  0.15961d0,-0.27072d0,
-     &  1.28439d0,-0.19444d0,-0.00661d0,-0.19444d0,-0.48824d0,
-     & -0.48824d0, 0.22898d0,
-     &  0.22898d0,-0.48824d0,-0.48824d0, 0.22898d0, 0.22898d0,
-     &  0.05556d0,-0.22222d0,
-     &  0.05556d0,-0.22222d0, 0.05556d0,-0.22222d0, 0.05556d0,
-     & -0.22222d0, 0.05556d0,
-     &  0.05556d0, 0.05556d0, 0.05556d0, 0.31481d0,-0.31246d0,
-     &  0.31481d0,-0.16902d0,
-     &  0.31481d0,-0.31246d0, 0.31481d0,-0.16902d0,-0.27072d0,
-     & -0.27072d0, 0.15961d0,
-     &  0.15961d0,-0.27072d0,-0.27072d0, 0.15961d0, 0.15961d0,
-     & -0.19444d0, 1.28439d0,
-     & -0.19444d0,-0.00661d0,-0.48824d0, 0.22898d0, 0.22898d0,
-     & -0.48824d0,-0.48824d0,
-     &  0.22898d0, 0.22898d0,-0.48824d0,-0.22222d0, 0.05556d0,
-     & -0.22222d0, 0.05556d0,
-     & -0.22222d0, 0.05556d0,-0.22222d0, 0.05556d0, 0.05556d0,
-     &  0.05556d0, 0.05556d0,
-     &  0.05556d0,-0.29630d0,-0.29630d0,-0.29630d0,-0.29630d0,
-     & -0.29630d0,-0.29630d0,
-     & -0.29630d0,-0.29630d0,-0.11111d0,-0.11111d0,-0.11111d0,
-     & -0.11111d0,-0.11111d0,
-     & -0.11111d0,-0.11111d0,-0.11111d0,-0.11111d0,-0.11111d0,
-     & -0.11111d0,-0.11111d0,
-     &  0.22898d0,-0.48824d0,-0.48824d0, 0.22898d0, 0.22898d0,
-     & -0.48824d0,-0.48824d0,
-     &  0.22898d0,-0.22222d0, 0.05556d0,-0.22222d0, 0.05556d0,
-     & -0.22222d0, 0.05556d0,
-     & -0.22222d0, 0.05556d0, 0.05556d0, 0.05556d0, 0.05556d0,
-     &  0.05556d0, 0.31481d0,
-     & -0.16902d0, 0.31481d0,-0.31246d0, 0.31481d0,-0.16902d0,
-     &  0.31481d0,-0.31246d0,
-     &  0.15961d0, 0.15961d0,-0.27072d0,-0.27072d0, 0.15961d0,
-     &  0.15961d0,-0.27072d0,
-     & -0.27072d0,-0.19444d0,-0.00661d0,-0.19444d0, 1.28439d0,
-     &  0.22898d0, 0.22898d0,
-     & -0.48824d0,-0.48824d0, 0.22898d0, 0.22898d0,-0.48824d0,
-     & -0.48824d0, 0.05556d0,
-     & -0.22222d0, 0.05556d0,-0.22222d0, 0.05556d0,-0.22222d0,
-     &  0.05556d0,-0.22222d0,
-     &  0.05556d0, 0.05556d0, 0.05556d0, 0.05556d0,-0.16902d0,
-     &  0.31481d0,-0.31246d0,
-     &  0.31481d0,-0.16902d0, 0.31481d0,-0.31246d0, 0.31481d0,
-     &  0.15961d0,-0.27072d0,
-     & -0.27072d0, 0.15961d0, 0.15961d0,-0.27072d0,-0.27072d0
-     & , 0.15961d0,-0.00661d0,
-     & -0.19444d0, 1.28439d0,-0.19444d0,-0.12559d0,-0.16145d0,
-     &  0.11575d0,-0.16145d0,
-     &  2.37499d0,-0.12559d0,-0.16145d0,-0.12559d0, 0.11111d0,
-     & -0.10405d0,-0.10405d0,
-     &  0.11111d0, 0.32628d0, 0.11111d0, 0.11111d0, 0.32628d0,
-     &  0.32628d0, 0.11111d0,
-     & -0.10405d0, 0.11111d0, 0.31481d0, 0.31481d0,-0.16902d0,
-     & -0.16902d0,-0.31246d0,
-     & -0.31246d0, 0.31481d0, 0.31481d0,-0.19444d0, 0.15961d0,
-     & -0.00661d0, 0.15961d0,
-     &  1.28439d0,-0.27072d0,-0.19444d0,-0.27072d0,-0.27072d0,
-     & -0.27072d0, 0.15961d0,
-     &  0.15961d0,-0.16145d0,-0.12559d0,-0.16145d0, 0.11575d0,
-     & -0.12559d0, 2.37499d0,
-     & -0.12559d0,-0.16145d0, 0.11111d0, 0.11111d0,-0.10405d0,
-     & -0.10405d0, 0.32628d0,
-     &  0.32628d0, 0.11111d0, 0.11111d0, 0.11111d0, 0.32628d0,
-     &  0.11111d0,-0.10405d0,
-     &  0.31481d0,-0.16902d0,-0.16902d0, 0.31481d0,-0.31246d0,
-     &  0.31481d0, 0.31481d0,
-     & -0.31246d0, 0.15961d0,-0.00661d0, 0.15961d0,-0.19444d0,
-     & -0.27072d0,-0.19444d0,
-     & -0.27072d0, 1.28439d0,-0.27072d0, 0.15961d0, 0.15961d0,
-     & -0.27072d0, 0.22898d0,
-     &  0.22898d0, 0.22898d0, 0.22898d0,-0.48824d0,-0.48824d0,
-     & -0.48824d0,-0.48824d0,
-     &  0.05556d0, 0.05556d0, 0.05556d0, 0.05556d0, 0.05556d0,
-     &  0.05556d0, 0.05556d0,
-     &  0.05556d0,-0.22222d0,-0.22222d0,-0.22222d0,-0.22222d0,
-     & -0.16902d0, 0.31481d0,
-     &  0.31481d0,-0.16902d0, 0.31481d0,-0.31246d0,-0.31246d0,
-     &  0.31481d0, 0.15961d0,
-     & -0.19444d0, 0.15961d0,-0.00661d0,-0.27072d0, 1.28439d0,
-     & -0.27072d0,-0.19444d0,
-     &  0.15961d0,-0.27072d0,-0.27072d0, 0.15961d0,-0.16145d0,
-     &  0.11575d0,-0.16145d0,
-     & -0.12559d0,-0.12559d0,-0.16145d0,-0.12559d0, 2.37499d0,
-     & -0.10405d0,-0.10405d0,
-     &  0.11111d0, 0.11111d0, 0.11111d0, 0.11111d0, 0.32628d0,
-     &  0.32628d0, 0.11111d0,
-     & -0.10405d0, 0.11111d0, 0.32628d0,-0.16902d0,-0.16902d0,
-     &  0.31481d0, 0.31481d0,
-     &  0.31481d0, 0.31481d0,-0.31246d0,-0.31246d0,-0.00661d0,
-     &  0.15961d0,-0.19444d0,
-     &  0.15961d0,-0.19444d0,-0.27072d0, 1.28439d0,-0.27072d0,
-     &  0.15961d0, 0.15961d0,
-     & -0.27072d0,-0.27072d0, 0.11575d0,-0.16145d0,-0.12559d0,
-     & -0.16145d0,-0.16145d0,
-     & -0.12559d0, 2.37499d0,-0.12559d0,-0.10405d0, 0.11111d0,
-     &  0.11111d0,-0.10405d0,
-     &  0.11111d0, 0.32628d0, 0.32628d0, 0.11111d0,-0.10405d0,
-     &  0.11111d0, 0.32628d0,
-     &  0.11111d0/
+     & +2.3749905130142061615572402d+0, -1.2558804006385479445917497d-1,
+     & -1.6144899697318323972261567d-1, -1.2558804006385501650377989d-1,
+     & -1.2558804006385515528165797d-1, -1.6144899697318335074491813d-1,
+     & +1.1575022772654103286882332d-1, -1.6144899697318332298934251d-1,
+     & +3.2627685256707739114290234d-1, +1.1111111111111149352126404d-1,
+     & +1.1111111111111166005471773d-1, +3.2627685256707633643102895d-1,
+     & +1.1111111111111099392090296d-1, -1.0405463034485473716728166d-1,
+     & -1.0405463034485518125649151d-1, +1.1111111111111171556586896d-1,
+     & +3.2627685256707716909829742d-1, +1.1111111111111129923223473d-1,
+     & -1.0405463034485522288985493d-1, +1.1111111111111163229914212d-1,
+     & -3.1246265455939625788062131d-1, -3.1246265455939359334536221d-1,
+     & +3.1481481481481216055939854d-1, +3.1481481481481499162811133d-1,
+     & +3.1481481481481182749249115d-1, +3.1481481481481465856120394d-1,
+     & -1.6901882692208583525150800d-1, -1.6901882692208175518189250d-1,
+     & +1.2843861132567893079681198d+0, -2.7072129701152342073910972d-1,
+     & -1.9444444444444389241688498d-1, -2.7072129701152392033947081d-1,
+     & -1.9444444444444367037228005d-1, +1.5961018590041298192971908d-1,
+     & -6.6083354790133445183508698d-3, +1.5961018590041253784050923d-1,
+     & -2.7072129701152364278371465d-1, -2.7072129701152358727256342d-1,
+     & +1.5961018590041348153008016d-1, +1.5961018590041123332845530d-1,
+     & -1.2558804006385776430576584d-1, +2.3749905130142034970219811d+0,
+     & -1.2558804006385507201493112d-1, -1.6144899697318287890013266d-1,
+     & -1.6144899697318521036848438d-1, -1.2558804006385462792572127d-1,
+     & -1.6144899697318235154419597d-1, +1.1575022772654162961369906d-1,
+     & +3.2627685256707922301089297d-1, +3.2627685256707933403319544d-1,
+     & +1.1111111111111089677638830d-1, +1.1111111111111129923223473d-1,
+     & +1.1111111111111106330984200d-1, +1.1111111111111049432054187d-1,
+     & -1.0405463034485712414678460d-1, -1.0405463034485641637960640d-1,
+     & +1.1111111111111222904401785d-1, +3.2627685256707822381017081d-1,
+     & +1.1111111111111078575408584d-1, -1.0405463034485666617978694d-1,
+     & -3.1246265455939692401443608d-1, +3.1481481481481471407235517d-1,
+     & +3.1481481481481521367271625d-1, -3.1246265455939792321515824d-1,
+     & +3.1481481481481587980653103d-1, -1.6901882692208511360654200d-1,
+     & -1.6901882692208533565114692d-1, +3.1481481481481415896084286d-1,
+     & -2.7072129701152436442868066d-1, -1.9444444444444536346239261d-1,
+     & -2.7072129701152353176141219d-1, +1.2843861132567899741019346d+0,
+     & +1.5961018590041256559608485d-1, -6.6083354790130886466381632d-3,
+     & +1.5961018590041356479680701d-1, -1.9444444444444441977282167d-1,
+     & -2.7072129701152392033947081d-1, +1.5961018590041231579590431d-1,
+     & +1.5961018590041278764068977d-1, -2.7072129701152264358299249d-1,
+     & -4.8823919872290821864169175d-1, -4.8823919872290794108593559d-1,
+     & -4.8823919872291077215464838d-1, -4.8823919872290899579780898d-1,
+     & +2.2897993946364816886251958d-1, +2.2897993946364941786342229d-1,
+     & +2.2897993946364988970820775d-1, +2.2897993946365155504274469d-1,
+     & +5.5555555555554851643318415d-2, +5.5555555555555712066162499d-2,
+     & +5.5555555555555725943950307d-2, +5.5555555555556031255282079d-2,
+     & +5.5555555555555920232979616d-2, +5.5555555555555365121467304d-2,
+     & +5.5555555555554796132167183d-2, +5.5555555555555392877042920d-2,
+     & -2.2222222222222168253047414d-1, -2.2222222222222254295331823d-1,
+     & -2.2222222222222151599702045d-1, -2.2222222222222257070889384d-1,
+     & +3.1481481481481549122847241d-1, -3.1246265455939775668170455d-1,
+     & -3.1246265455939814525976317d-1, +3.1481481481481582429537980d-1,
+     & -1.6901882692208500258423953d-1, +3.1481481481481621287343842d-1,
+     & +3.1481481481481532469501872d-1, -1.6901882692208539116229815d-1,
+     & -2.7072129701152164438227032d-1, +1.2843861132567930827264036d+0,
+     & -2.7072129701152281011644618d-1, -1.9444444444444328179422143d-1,
+     & +1.5961018590041139986190899d-1, -1.9444444444444586306275369d-1,
+     & +1.5961018590041073372809421d-1, -6.6083354790144668844398268d-3,
+     & +1.5961018590041303744087031d-1, -2.7072129701152297664989987d-1,
+     & -2.7072129701152392033947081d-1, +1.5961018590041164966208953d-1,
+     & -1.2558804006385637652698506d-1, -1.6144899697318712550320186d-1,
+     & -1.2558804006384907681059815d-1, +2.3749905130142003883975121d+0,
+     & -1.6144899697317796616324870d-1, +1.1575022772654076919085497d-1,
+     & -1.6144899697318212949959104d-1, -1.2558804006385926310684908d-1,
+     & +1.1111111111111218741065443d-1, +1.1111111111111010574248326d-1,
+     & +3.2627685256707761318750727d-1, +3.2627685256707888994398559d-1,
+     & -1.0405463034485737394696514d-1, -1.0405463034485638862403079d-1,
+     & +1.1111111111111186822153485d-1, +1.1111111111111141025453719d-1,
+     & +1.1111111111111179883259581d-1, -1.0405463034485576412357943d-1,
+     & +1.1111111111110900939724644d-1, +3.2627685256707994465585898d-1,
+     & +3.1481481481481549122847241d-1, +3.1481481481481282669321331d-1,
+     & -3.1246265455939442601263067d-1, -3.1246265455939720157019224d-1,
+     & -1.6901882692208336500527821d-1, -1.6901882692208519687326884d-1,
+     & +3.1481481481481515816156502d-1, +3.1481481481481332629357439d-1,
+     & -1.9444444444444225483792366d-1, -2.7072129701152342073910972d-1,
+     & +1.2843861132567906402357494d+0, -2.7072129701152181091572402d-1,
+     & -6.6083354790147739304950747d-3, +1.5961018590041103903942599d-1,
+     & -1.9444444444444558550699753d-1, +1.5961018590041214926245061d-1,
+     & +1.5961018590041251008493361d-1, +1.5961018590041275988511416d-1,
+     & -2.7072129701152436442868066d-1, -2.7072129701152253256069002d-1,
+     & -1.6144899697318501607945507d-1, -1.2558804006385732021655599d-1,
+     & +2.3749905130142030529327712d+0, -1.2558804006385482221475058d-1,
+     & +1.1575022772654303127026765d-1, -1.6144899697318307318916197d-1,
+     & -1.2558804006385507201493112d-1, -1.6144899697318496056830384d-1,
+     & +1.1111111111111199312162512d-1, +3.2627685256707950056664913d-1,
+     & +3.2627685256707839034362451d-1, +1.1111111111111152127683965d-1,
+     & -1.0405463034485801232520430d-1, +1.1111111111111041105381503d-1,
+     & +1.1111111111111079963187365d-1, -1.0405463034485795681405307d-1,
+     & -1.0405463034485722129129925d-1, +1.1111111111111175719923239d-1,
+     & +3.2627685256707750216520481d-1, +1.1111111111111160454356650d-1,
+     & -3.1246265455939781219285578d-1, +3.1481481481481216055939854d-1,
+     & -1.6901882692208303193837082d-1, +3.1481481481481571327307734d-1,
+     & -3.1246265455939448152378191d-1, +3.1481481481481587980653103d-1,
+     & -1.6901882692208552994017623d-1, +3.1481481481481204953709607d-1,
+     & -2.7072129701152292113874864d-1, +1.5961018590041298192971908d-1,
+     & +1.5961018590041237130705554d-1, -2.7072129701152497505134420d-1,
+     & -2.7072129701152436442868066d-1, +1.5961018590041345377450455d-1,
+     & +1.5961018590041300968529470d-1, -2.7072129701152336522795849d-1,
+     & +1.2843861132567895300127248d+0, -1.9444444444444416997264113d-1,
+     & -6.6083354790131745154502241d-3, -1.9444444444444400343918744d-1,
+     & -4.8823919872290866273090160d-1, -4.8823919872291138277731193d-1,
+     & +2.2897993946365219342098385d-1, +2.2897993946365047257529568d-1,
+     & -4.8823919872290605370679373d-1, -4.8823919872290855170859913d-1,
+     & +2.2897993946364939010784667d-1, +2.2897993946364658679470949d-1,
+     & +5.5555555555556107583115022d-2, -2.2222222222222276499792315d-1,
+     & +5.5555555555554921032257454d-2, -2.2222222222222173804162537d-1,
+     & +5.5555555555554997360090397d-2, -2.2222222222222245968659138d-1,
+     & +5.5555555555556412894446794d-2, -2.2222222222222173804162537d-1,
+     & +5.5555555555554907154469646d-2, +5.5555555555556489222279737d-2,
+     & +5.5555555555554768376591568d-2, +5.5555555555555864721828385d-2,
+     & +3.1481481481481582429537980d-1, -3.1246265455939853383782179d-1,
+     & +3.1481481481481804474142905d-1, -1.6901882692208644587417155d-1,
+     & +3.1481481481481710105185812d-1, -3.1246265455939759014825086d-1,
+     & +3.1481481481481565776192610d-1, -1.6901882692208705649683509d-1,
+     & -2.7072129701152131131536294d-1, -2.7072129701152486402904174d-1,
+     & +1.5961018590041101128385037d-1, +1.5961018590041312070759716d-1,
+     & -2.7072129701152231051608510d-1, -2.7072129701152375380601711d-1,
+     & +1.5961018590041206599572376d-1, +1.5961018590041239906263115d-1,
+     & -1.9444444444444369812785567d-1, +1.2843861132567919725033789d+0,
+     & -1.9444444444444680675232462d-1, -6.6083354790126783845360947d-3,
+     & -4.8823919872290966193162376d-1, +2.2897993946365174933177400d-1,
+     & +2.2897993946364897377421244d-1, -4.8823919872290882926435529d-1,
+     & -4.8823919872291038357658977d-1, +2.2897993946364955664130036d-1,
+     & +2.2897993946364897377421244d-1, -4.8823919872290783006363313d-1,
+     & -2.2222222222222218213083522d-1, +5.5555555555556114522008926d-2,
+     & -2.2222222222222243193101576d-1, +5.5555555555555927171873520d-2,
+     & -2.2222222222222151599702045d-1, +5.5555555555555295732528265d-2,
+     & -2.2222222222222232090871330d-1, +5.5555555555554927971151358d-2,
+     & +5.5555555555555420632618535d-2, +5.5555555555554879398894030d-2,
+     & +5.5555555555556045133069887d-2, +5.5555555555555143076862379d-2,
+     & -2.9629629629629505860322070d-1, -2.9629629629629572473703547d-1,
+     & -2.9629629629629583575933793d-1, -2.9629629629629633535969901d-1,
+     & -2.9629629629629594678164040d-1, -2.9629629629629655740430394d-1,
+     & -2.9629629629629583575933793d-1, -2.9629629629629583575933793d-1,
+     & -1.1111111111111098004311515d-1, -1.1111111111111106330984200d-1,
+     & -1.1111111111111071636514680d-1, -1.1111111111111100779869076d-1,
+     & -1.1111111111111098004311515d-1, -1.1111111111111082738744926d-1,
+     & -1.1111111111111118820993227d-1, -1.1111111111111088289860049d-1,
+     & -1.1111111111111088289860049d-1, -1.1111111111111095228753953d-1,
+     & -1.1111111111111109106541761d-1, -1.1111111111111099392090296d-1,
+     & +2.2897993946365111095353484d-1, -4.8823919872291154931076562d-1,
+     & -4.8823919872290760801902820d-1, +2.2897993946365025053069076d-1,
+     & +2.2897993946365105544238361d-1, -4.8823919872290849619744790d-1,
+     & -4.8823919872290844068629667d-1, +2.2897993946364875172960751d-1,
+     & -2.2222222222222173804162537d-1, +5.5555555555555559410496613d-2,
+     & -2.2222222222222171028604976d-1, +5.5555555555555351243679496d-2,
+     & -2.2222222222222245968659138d-1, +5.5555555555555767577313730d-2,
+     & -2.2222222222222171028604976d-1, +5.5555555555555767577313730d-2,
+     & +5.5555555555555295732528265d-2, +5.5555555555555705127268595d-2,
+     & +5.5555555555555004298984301d-2, +5.5555555555555892477404001d-2,
+     & +3.1481481481481554673962364d-1, -1.6901882692208303193837082d-1,
+     & +3.1481481481481121686982760d-1, -3.1246265455939586930256269d-1,
+     & +3.1481481481481082829176898d-1, -1.6901882692208441971715160d-1,
+     & +3.1481481481481587980653103d-1, -3.1246265455939370436766467d-1,
+     & +1.5961018590041253784050923d-1, +1.5961018590041117781730406d-1,
+     & -2.7072129701152325420565603d-1, -2.7072129701152192193802648d-1,
+     & +1.5961018590041328724105085d-1, +1.5961018590041095577269914d-1,
+     & -2.7072129701152414238407573d-1, -2.7072129701152169989342156d-1,
+     & -1.9444444444444267117155789d-1, -6.6083354790150705682094667d-3,
+     & -1.9444444444444575204045123d-1, +1.2843861132567913063695642d+0,
+     & +2.2897993946365069461990061d-1, +2.2897993946365294282152547d-1,
+     & -4.8823919872291138277731193d-1, -4.8823919872290899579780898d-1,
+     & +2.2897993946364636475010457d-1, +2.2897993946364950113014913d-1,
+     & -4.8823919872290888477550652d-1, -4.8823919872290616472909619d-1,
+     & +5.5555555555555143076862379d-2, -2.2222222222222237641986453d-1,
+     & +5.5555555555556142277584541d-2, -2.2222222222222215437525961d-1,
+     & +5.5555555555556274116568716d-2, -2.2222222222222207110853276d-1,
+     & +5.5555555555554844704424511d-2, -2.2222222222222157150817168d-1,
+     & +5.5555555555556010438600367d-2, +5.5555555555555087565711148d-2,
+     & +5.5555555555556246360993100d-2, +5.5555555555554782254379376d-2,
+     & -1.6901882692208533565114692d-1, +3.1481481481481687900725319d-1,
+     & -3.1246265455939903343818287d-1, +3.1481481481481426998314532d-1,
+     & -1.6901882692208636260744470d-1, +3.1481481481481493611696010d-1,
+     & -3.1246265455939770117055332d-1, +3.1481481481481687900725319d-1,
+     & +1.5961018590041040066118683d-1, -2.7072129701152281011644618d-1,
+     & -2.7072129701152225500493387d-1, +1.5961018590041073372809421d-1,
+     & +1.5961018590041234355147992d-1, -2.7072129701152303216105111d-1,
+     & -2.7072129701152258807184126d-1, +1.5961018590041131659518214d-1,
+     & -6.6083354790149621479922182d-3, -1.9444444444444450303954852d-1,
+     & +1.2843861132567944149940331d+0, -1.9444444444444569652929999d-1,
+     & -1.2558804006385707041637545d-1, -1.6144899697318029763160041d-1,
+     & +1.1575022772653729974390302d-1, -1.6144899697318318421146444d-1,
+     & +2.3749905130141977238622530d+0, -1.2558804006385462792572127d-1,
+     & -1.6144899697318174092153242d-1, -1.2558804006385054785610578d-1,
+     & +1.1111111111111142413232500d-1, -1.0405463034485587514588190d-1,
+     & -1.0405463034485570861242820d-1, +1.1111111111111218741065443d-1,
+     & +3.2627685256707839034362451d-1, +1.1111111111111081350966145d-1,
+     & +1.1111111111110963389769779d-1, +3.2627685256707772420980973d-1,
+     & +3.2627685256707900096628805d-1, +1.1111111111111046656496626d-1,
+     & -1.0405463034485554207897451d-1, +1.1111111111111003635354422d-1,
+     & +3.1481481481481532469501872d-1, +3.1481481481481626838458965d-1,
+     & -1.6901882692208602954053731d-1, -1.6901882692208566871805431d-1,
+     & -3.1246265455939892241588041d-1, -3.1246265455939670196983116d-1,
+     & +3.1481481481481549122847241d-1, +3.1481481481481587980653103d-1,
+     & -1.9444444444444364261670444d-1, +1.5961018590041220477360184d-1,
+     & -6.6083354790139551410144136d-3, +1.5961018590041317621874839d-1,
+     & +1.2843861132567917504587740d+0, -2.7072129701152458647328558d-1,
+     & -1.9444444444444597408505615d-1, -2.7072129701152358727256342d-1,
+     & -2.7072129701152147784881663d-1, -2.7072129701152325420565603d-1,
+     & +1.5961018590041162190651391d-1, +1.5961018590041162190651391d-1,
+     & -1.6144899697318335074491813d-1, -1.2558804006385107521204247d-1,
+     & -1.6144899697318709774762624d-1, +1.1575022772654164349148687d-1,
+     & -1.2558804006385970719605893d-1, +2.3749905130142017206651417d+0,
+     & -1.2558804006385543283741413d-1, -1.6144899697317885434166840d-1,
+     & +1.1111111111111038329823941d-1, +1.1111111111111066085399557d-1,
+     & -1.0405463034485656903527229d-1, -1.0405463034485608331269901d-1,
+     & +3.2627685256708022221161514d-1, +3.2627685256707772420980973d-1,
+     & +1.1111111111110931470857821d-1, +1.1111111111111240945525935d-1,
+     & +1.1111111111111300620013509d-1, +3.2627685256707733563175111d-1,
+     & +1.1111111111111181271038362d-1, -1.0405463034485783191396280d-1,
+     & +3.1481481481481593531768226d-1, -1.6901882692208522462884446d-1,
+     & -1.6901882692208569647362992d-1, +3.1481481481481443651659902d-1,
+     & -3.1246265455939636890292377d-1, +3.1481481481481526918386749d-1,
+     & +3.1481481481481449202775025d-1, -3.1246265455939570276910899d-1,
+     & +1.5961018590041384235256317d-1, -6.6083354790135110518045636d-3,
+     & +1.5961018590041292641856785d-1, -1.9444444444444311526076774d-1,
+     & -2.7072129701152347625026096d-1, -1.9444444444444591857390492d-1,
+     & -2.7072129701152508607364666d-1, +1.2843861132567906402357494d+0,
+     & -2.7072129701152203296032894d-1, +1.5961018590041201048457253d-1,
+     & +1.5961018590041189946227007d-1, -2.7072129701152380931716834d-1,
+     & +2.2897993946365061135317376d-1, +2.2897993946364780804003658d-1,
+     & +2.2897993946365219342098385d-1, +2.2897993946364972317475406d-1,
+     & -4.8823919872290655330715481d-1, -4.8823919872290932886471637d-1,
+     & -4.8823919872290927335356514d-1, -4.8823919872291082766579962d-1,
+     & +5.5555555555555968805236944d-2, +5.5555555555555864721828385d-2,
+     & +5.5555555555555524716027094d-2, +5.5555555555555170832437994d-2,
+     & +5.5555555555554539393092739d-2, +5.5555555555555899416297905d-2,
+     & +5.5555555555555996560812559d-2, +5.5555555555555434510406343d-2,
+     & -2.2222222222222198784180591d-1, -2.2222222222222154375259606d-1,
+     & -2.2222222222222295928695246d-1, -2.2222222222222151599702045d-1,
+     & -1.6901882692208602954053731d-1, +3.1481481481481576878422857d-1,
+     & +3.1481481481481371487163301d-1, -1.6901882692208530789557130d-1,
+     & +3.1481481481481382589393547d-1, -3.1246265455939703503673854d-1,
+     & -3.1246265455939703503673854d-1, +3.1481481481481582429537980d-1,
+     & +1.5961018590041045617233806d-1, -1.9444444444444489161760714d-1,
+     & +1.5961018590041189946227007d-1, -6.6083354790152527141744443d-3,
+     & -2.7072129701152247704953879d-1, +1.2843861132567935268156134d+0,
+     & -2.7072129701152242153838756d-1, -1.9444444444444400343918744d-1,
+     & +1.5961018590041031739445998d-1, -2.7072129701152364278371465d-1,
+     & -2.7072129701152214398263141d-1, +1.5961018590041220477360184d-1,
+     & -1.6144899697318362830067429d-1, +1.1575022772653784097762752d-1,
+     & -1.6144899697318049192062972d-1, -1.2558804006385601570450206d-1,
+     & -1.2558804006385279605773064d-1, -1.6144899697318265685552774d-1,
+     & -1.2558804006385446139226758d-1, +2.3749905130141986120406727d+0,
+     & -1.0405463034485584739030628d-1, -1.0405463034485686046881625d-1,
+     & +1.1111111111111091065417611d-1, +1.1111111111111157678799088d-1,
+     & +1.1111111111111068860957118d-1, +1.1111111111111046656496626d-1,
+     & +3.2627685256707894545513682d-1, +3.2627685256707839034362451d-1,
+     & +1.1111111111111113269878103d-1, -1.0405463034485616657942586d-1,
+     & +1.1111111111110917593070013d-1, +3.2627685256707922301089297d-1,
+     & -1.6901882692208505809539076d-1, -1.6901882692208772263064986d-1,
+     & +3.1481481481481715656300935d-1, +3.1481481481481510265041379d-1,
+     & +3.1481481481481765616337043d-1, +3.1481481481481515816156502d-1,
+     & -3.1246265455939731259249470d-1, -3.1246265455939992161660257d-1,
+     & -6.6083354790145171914206301d-3, +1.5961018590041184395111884d-1,
+     & -1.9444444444444558550699753d-1, +1.5961018590041017861658190d-1,
+     & -1.9444444444444464181742660d-1, -2.7072129701152175540457279d-1,
+     & +1.2843861132567935268156134d+0, -2.7072129701152303216105111d-1,
+     & +1.5961018590041028963888436d-1, +1.5961018590041234355147992d-1,
+     & -2.7072129701152281011644618d-1, -2.7072129701152203296032894d-1,
+     & +1.1575022772654330882602380d-1, -1.6144899697318310094473759d-1,
+     & -1.2558804006385540508183851d-1, -1.6144899697318235154419597d-1,
+     & -1.6144899697318287890013266d-1, -1.2558804006385546059298974d-1,
+     & +2.3749905130141999443083023d+0, -1.2558804006385557161529221d-1,
+     & -1.0405463034485845641441415d-1, +1.1111111111111122984329569d-1,
+     & +1.1111111111111109106541761d-1, -1.0405463034485783191396280d-1,
+     & +1.1111111111111139637674938d-1, +3.2627685256708005567816144d-1,
+     & +3.2627685256708055527852252d-1, +1.1111111111111132698781034d-1,
+     & -1.0405463034485880335910934d-1, +1.1111111111111135474338596d-1,
+     & +3.2627685256708066630082499d-1, +1.1111111111111106330984200d-1/
 !
       data iflag /1/
 !
diff --git a/ccx/src/extrapolate_se.f b/ccx/src/extrapolate_se.f
index 33a841d..74afad0 100644
--- a/ccx/src/extrapolate_se.f
+++ b/ccx/src/extrapolate_se.f
@@ -104,10 +104,19 @@
 !
       data a2 /1.1455d0,-0.1455d0,1.1455d0,-0.1455d0,1.1455d0,-0.1455d0,
      &        -0.1455d0,1.1455d0,-0.1455d0,1.1455d0,-0.1455d0,1.1455d0/
-      data a4 /  1.92705d0, -0.30902d0, -0.30902d0, -0.30902d0,
-     &          -0.30902d0,  1.92705d0, -0.30902d0, -0.30902d0,
-     &          -0.30902d0, -0.30902d0,  1.92705d0, -0.30902d0,
-     &          -0.30902d0, -0.30902d0, -0.30902d0,  1.92705d0/
+
+!     Precision increased on these 8/2021 Pete Gustafson (after a course
+!     excercise demonstrating extrapolation). Current precision is
+!     ridiculous, but comes for free.
+      data a4 /
+     & +1.9270509831248423537886083d+0, -3.0901699437494745126286944d-1,
+     & -3.0901699437494745126286944d-1, -3.0901699437494745126286944d-1,
+     & -3.0901699437494745126286944d-1, +1.9270509831248423537886083d+0,
+     & -3.0901699437494745126286944d-1, -3.0901699437494745126286944d-1,
+     & -3.0901699437494745126286944d-1, -3.0901699437494745126286944d-1,
+     & +1.9270509831248423537886083d+0, -3.0901699437494745126286944d-1,
+     & -3.0901699437494745126286944d-1, -3.0901699437494745126286944d-1,
+     & -3.0901699437494745126286944d-1, +1.9270509831248423537886083d+0/
 !
 !     extrapolation from a 6 integration point scheme in a wedge to
 !     the vertex nodes
@@ -124,7 +133,7 @@
      &         -0.34151d0,
      &          0.09151d0, 0.09151d0,-0.45753d0,-0.34151d0,-0.34151d0,
      &          1.70753d0/
-!
+
       data a9 / 1.63138d0,-0.32628d0,-0.32628d0,-0.52027d0, 0.10405d0,
      &          0.10405d0,
      &         -0.32628d0, 1.63138d0,-0.32628d0, 0.10405d0,-0.52027d0,
@@ -145,19 +154,44 @@
      &          1.63138d0/
 !
 !     extrapolation from a 2x2x2=8 integration point scheme in a hex to
-!     the vertex nodes
-!    
-      data a8 /2.549d0,-.683d0,.183d0,-.683d0,-.683d0,.183d0,
-     &        -.04904d0,.183d0,-.683d0,2.549d0,-.683d0,.183d0,
-     &        .183d0,-.683d0,.183d0,-.04904d0,-.683d0,.183d0,
-     &        -.683d0,2.549d0,.183d0,-.04904d0,.183d0,-.683d0,
-     &        .183d0,-.683d0,2.549d0,-.683d0,-.04904d0,.183d0,
-     &        -.683d0,.183d0,-.683d0,.183d0,-.04904d0,.183d0,
-     &        2.549d0,-.683d0,.183d0,-.683d0,.183d0,-.683d0,
-     &        .183d0,-.04904d0,-.683d0,2.549d0,-.683d0,.183d0,
-     &        .183d0,-.04904d0,.183d0,-.683d0,-.683d0,.183d0,
-     &        -.683d0,2.549d0,-.04904d0,.183d0,-.683d0,.183d0,
-     &        .183d0,-.683d0,2.549d0,-.683d0/  
+!     the vertex nodes.  Precision increased on these 8/2021 Pete
+!     Gustafson (after a course excercise demonstrating
+!     extrapolation). Current precision is ridiculous, but comes for
+!     free.
+!     
+      data a8 /
+     & +2.5490381056766566736371260d+0, -6.8301270189221807704882394d-1,
+     & -6.8301270189221940931645349d-1, +1.8301270189221954809433157d-1,
+     & -6.8301270189221896522724364d-1, +1.8301270189221824358227764d-1,
+     & +1.8301270189221996442796581d-1, -4.9038105676658283460511711d-2,
+     & -6.8301270189221874318263872d-1, +2.5490381056766580059047556d+0,
+     & +1.8301270189221963136105842d-1, -6.8301270189221918727184857d-1,
+     & +1.8301270189221860440476064d-1, -6.8301270189221807704882394d-1,
+     & -4.9038105676659414500218048d-2, +1.8301270189221935380530226d-1,
+     & +1.8301270189222007545026827d-1, -6.8301270189221818807112641d-1,
+     & -6.8301270189221907624954611d-1, +2.5490381056766580059047556d+0,
+     & -4.9038105676658678977464234d-2, +1.8301270189221757744846286d-1,
+     & +1.8301270189221979789451211d-1, -6.8301270189221829909342887d-1,
+     & -6.8301270189221874318263872d-1, +1.8301270189221782724864340d-1,
+     & +2.5490381056766575618155457d+0, -6.8301270189222018647257073d-1,
+     & +1.8301270189222024198372196d-1, -4.9038105676655389941753782d-2,
+     & -6.8301270189222040851717566d-1, +1.8301270189221932604972665d-1,
+     & -6.8301270189221874318263872d-1, +1.8301270189221952033875596d-1,
+     & +1.8301270189221940931645349d-1, -4.9038105676658602649631291d-2,
+     & +2.5490381056766566736371260d+0, -6.8301270189221907624954611d-1,
+     & -6.8301270189221918727184857d-1, +1.8301270189221968687220965d-1,
+     & +1.8301270189221968687220965d-1, -6.8301270189221963136105842d-1,
+     & -4.9038105676657818554620150d-2, +1.8301270189221965911663403d-1,
+     & -6.8301270189221818807112641d-1, +2.5490381056766553413694965d+0,
+     & +1.8301270189222118567329289d-1, -6.8301270189221863216033626d-1,
+     & -4.9038105676657922638028708d-2, +1.8301270189221990891681457d-1,
+     & +1.8301270189221979789451211d-1, -6.8301270189222029749487319d-1,
+     & +1.8301270189221993667239019d-1, -6.8301270189221741091500917d-1,
+     & -6.8301270189222107465099043d-1, +2.5490381056766566736371260d+0,
+     & +1.8301270189221985340566334d-1, -4.9038105676658318154981231d-2,
+     & -6.8301270189221918727184857d-1, +1.8301270189222004769469265d-1,
+     & -6.8301270189221852113803379d-1, +1.8301270189221771622634094d-1,
+     & +2.5490381056766584499939654d+0, -6.8301270189221885420494118d-1/
 !
 !     extrapolation from a 14 integration point scheme in a hex to
 !     the vertex nodes
@@ -193,164 +227,282 @@
      &  0.2069d+00,0.2069d+00,0.2069d+00,0.2069d+00/
 !
 !     extrapolation from a 3x3x3=27 integration point scheme in a hex to
-!     the all nodes in a 20-node element
+!     the all nodes in a 20-node element. Precision increased on these
+!     8/2021 Pete Gustafson (after a course excercise demonstrating
+!     extrapolation). Current precision is ridiculous, but comes for
+!     free.
 !    
       data a27 /
-     &  2.37499d0,-0.12559d0,-0.16145d0,-0.12559d0,-0.12559d0,
-     & -0.16145d0, 0.11575d0,
-     & -0.16145d0, 0.32628d0, 0.11111d0, 0.11111d0, 0.32628d0,
-     &  0.11111d0,-0.10405d0,
-     & -0.10405d0, 0.11111d0, 0.32628d0, 0.11111d0,-0.10405d0,
-     &  0.11111d0,-0.31246d0,
-     & -0.31246d0, 0.31481d0, 0.31481d0, 0.31481d0, 0.31481d0,
-     & -0.16902d0,-0.16902d0,
-     &  1.28439d0,-0.27072d0,-0.19444d0,-0.27072d0,-0.19444d0,
-     &  0.15961d0,-0.00661d0,
-     &  0.15961d0,-0.27072d0,-0.27072d0, 0.15961d0, 0.15961d0,
-     & -0.12559d0, 2.37499d0,
-     & -0.12559d0,-0.16145d0,-0.16145d0,-0.12559d0,-0.16145d0,
-     &  0.11575d0, 0.32628d0,
-     &  0.32628d0, 0.11111d0, 0.11111d0, 0.11111d0, 0.11111d0,
-     & -0.10405d0,-0.10405d0,
-     &  0.11111d0, 0.32628d0, 0.11111d0,-0.10405d0,-0.31246d0,
-     &  0.31481d0, 0.31481d0,
-     & -0.31246d0, 0.31481d0,-0.16902d0,-0.16902d0, 0.31481d0,
-     & -0.27072d0,-0.19444d0,
-     & -0.27072d0, 1.28439d0, 0.15961d0,-0.00661d0, 0.15961d0,
-     & -0.19444d0,-0.27072d0,
-     &  0.15961d0, 0.15961d0,-0.27072d0,-0.48824d0,-0.48824d0,
-     & -0.48824d0,-0.48824d0,
-     &  0.22898d0, 0.22898d0, 0.22898d0, 0.22898d0, 0.05556d0,
-     &  0.05556d0, 0.05556d0,
-     &  0.05556d0, 0.05556d0, 0.05556d0, 0.05556d0, 0.05556d0,
-     & -0.22222d0,-0.22222d0,
-     & -0.22222d0,-0.22222d0, 0.31481d0,-0.31246d0,-0.31246d0,
-     &  0.31481d0,-0.16902d0,
-     &  0.31481d0, 0.31481d0,-0.16902d0,-0.27072d0, 1.28439d0,
-     & -0.27072d0,-0.19444d0,
-     &  0.15961d0,-0.19444d0, 0.15961d0,-0.00661d0, 0.15961d0,
-     & -0.27072d0,-0.27072d0,
-     &  0.15961d0,-0.12559d0,-0.16145d0,-0.12559d0, 2.37499d0,
-     & -0.16145d0, 0.11575d0,
-     & -0.16145d0,-0.12559d0, 0.11111d0, 0.11111d0, 0.32628d0,
-     &  0.32628d0,-0.10405d0,
-     & -0.10405d0, 0.11111d0, 0.11111d0, 0.11111d0,-0.10405d0,
-     &  0.11111d0, 0.32628d0,
-     &  0.31481d0, 0.31481d0,-0.31246d0,-0.31246d0,-0.16902d0,
-     & -0.16902d0, 0.31481d0,
-     &  0.31481d0,-0.19444d0,-0.27072d0, 1.28439d0,-0.27072d0,
-     & -0.00661d0, 0.15961d0,
-     & -0.19444d0, 0.15961d0, 0.15961d0, 0.15961d0,-0.27072d0,
-     & -0.27072d0,-0.16145d0,
-     & -0.12559d0, 2.37499d0,-0.12559d0, 0.11575d0,-0.16145d0,
-     & -0.12559d0,-0.16145d0,
-     &  0.11111d0, 0.32628d0, 0.32628d0, 0.11111d0,-0.10405d0,
-     &  0.11111d0, 0.11111d0,
-     & -0.10405d0,-0.10405d0, 0.11111d0, 0.32628d0, 0.11111d0,
-     & -0.31246d0, 0.31481d0,
-     & -0.16902d0, 0.31481d0,-0.31246d0, 0.31481d0,-0.16902d0,
-     &  0.31481d0,-0.27072d0,
-     &  0.15961d0, 0.15961d0,-0.27072d0,-0.27072d0, 0.15961d0,
-     &  0.15961d0,-0.27072d0,
-     &  1.28439d0,-0.19444d0,-0.00661d0,-0.19444d0,-0.48824d0,
-     & -0.48824d0, 0.22898d0,
-     &  0.22898d0,-0.48824d0,-0.48824d0, 0.22898d0, 0.22898d0,
-     &  0.05556d0,-0.22222d0,
-     &  0.05556d0,-0.22222d0, 0.05556d0,-0.22222d0, 0.05556d0,
-     & -0.22222d0, 0.05556d0,
-     &  0.05556d0, 0.05556d0, 0.05556d0, 0.31481d0,-0.31246d0,
-     &  0.31481d0,-0.16902d0,
-     &  0.31481d0,-0.31246d0, 0.31481d0,-0.16902d0,-0.27072d0,
-     & -0.27072d0, 0.15961d0,
-     &  0.15961d0,-0.27072d0,-0.27072d0, 0.15961d0, 0.15961d0,
-     & -0.19444d0, 1.28439d0,
-     & -0.19444d0,-0.00661d0,-0.48824d0, 0.22898d0, 0.22898d0,
-     & -0.48824d0,-0.48824d0,
-     &  0.22898d0, 0.22898d0,-0.48824d0,-0.22222d0, 0.05556d0,
-     & -0.22222d0, 0.05556d0,
-     & -0.22222d0, 0.05556d0,-0.22222d0, 0.05556d0, 0.05556d0,
-     &  0.05556d0, 0.05556d0,
-     &  0.05556d0,-0.29630d0,-0.29630d0,-0.29630d0,-0.29630d0,
-     & -0.29630d0,-0.29630d0,
-     & -0.29630d0,-0.29630d0,-0.11111d0,-0.11111d0,-0.11111d0,
-     & -0.11111d0,-0.11111d0,
-     & -0.11111d0,-0.11111d0,-0.11111d0,-0.11111d0,-0.11111d0,
-     & -0.11111d0,-0.11111d0,
-     &  0.22898d0,-0.48824d0,-0.48824d0, 0.22898d0, 0.22898d0,
-     & -0.48824d0,-0.48824d0,
-     &  0.22898d0,-0.22222d0, 0.05556d0,-0.22222d0, 0.05556d0,
-     & -0.22222d0, 0.05556d0,
-     & -0.22222d0, 0.05556d0, 0.05556d0, 0.05556d0, 0.05556d0,
-     &  0.05556d0, 0.31481d0,
-     & -0.16902d0, 0.31481d0,-0.31246d0, 0.31481d0,-0.16902d0,
-     &  0.31481d0,-0.31246d0,
-     &  0.15961d0, 0.15961d0,-0.27072d0,-0.27072d0, 0.15961d0,
-     &  0.15961d0,-0.27072d0,
-     & -0.27072d0,-0.19444d0,-0.00661d0,-0.19444d0, 1.28439d0,
-     &  0.22898d0, 0.22898d0,
-     & -0.48824d0,-0.48824d0, 0.22898d0, 0.22898d0,-0.48824d0,
-     & -0.48824d0, 0.05556d0,
-     & -0.22222d0, 0.05556d0,-0.22222d0, 0.05556d0,-0.22222d0,
-     &  0.05556d0,-0.22222d0,
-     &  0.05556d0, 0.05556d0, 0.05556d0, 0.05556d0,-0.16902d0,
-     &  0.31481d0,-0.31246d0,
-     &  0.31481d0,-0.16902d0, 0.31481d0,-0.31246d0, 0.31481d0,
-     &  0.15961d0,-0.27072d0,
-     & -0.27072d0, 0.15961d0, 0.15961d0,-0.27072d0,-0.27072d0
-     & , 0.15961d0,-0.00661d0,
-     & -0.19444d0, 1.28439d0,-0.19444d0,-0.12559d0,-0.16145d0,
-     &  0.11575d0,-0.16145d0,
-     &  2.37499d0,-0.12559d0,-0.16145d0,-0.12559d0, 0.11111d0,
-     & -0.10405d0,-0.10405d0,
-     &  0.11111d0, 0.32628d0, 0.11111d0, 0.11111d0, 0.32628d0,
-     &  0.32628d0, 0.11111d0,
-     & -0.10405d0, 0.11111d0, 0.31481d0, 0.31481d0,-0.16902d0,
-     & -0.16902d0,-0.31246d0,
-     & -0.31246d0, 0.31481d0, 0.31481d0,-0.19444d0, 0.15961d0,
-     & -0.00661d0, 0.15961d0,
-     &  1.28439d0,-0.27072d0,-0.19444d0,-0.27072d0,-0.27072d0,
-     & -0.27072d0, 0.15961d0,
-     &  0.15961d0,-0.16145d0,-0.12559d0,-0.16145d0, 0.11575d0,
-     & -0.12559d0, 2.37499d0,
-     & -0.12559d0,-0.16145d0, 0.11111d0, 0.11111d0,-0.10405d0,
-     & -0.10405d0, 0.32628d0,
-     &  0.32628d0, 0.11111d0, 0.11111d0, 0.11111d0, 0.32628d0,
-     &  0.11111d0,-0.10405d0,
-     &  0.31481d0,-0.16902d0,-0.16902d0, 0.31481d0,-0.31246d0,
-     &  0.31481d0, 0.31481d0,
-     & -0.31246d0, 0.15961d0,-0.00661d0, 0.15961d0,-0.19444d0,
-     & -0.27072d0,-0.19444d0,
-     & -0.27072d0, 1.28439d0,-0.27072d0, 0.15961d0, 0.15961d0,
-     & -0.27072d0, 0.22898d0,
-     &  0.22898d0, 0.22898d0, 0.22898d0,-0.48824d0,-0.48824d0,
-     & -0.48824d0,-0.48824d0,
-     &  0.05556d0, 0.05556d0, 0.05556d0, 0.05556d0, 0.05556d0,
-     &  0.05556d0, 0.05556d0,
-     &  0.05556d0,-0.22222d0,-0.22222d0,-0.22222d0,-0.22222d0,
-     & -0.16902d0, 0.31481d0,
-     &  0.31481d0,-0.16902d0, 0.31481d0,-0.31246d0,-0.31246d0,
-     &  0.31481d0, 0.15961d0,
-     & -0.19444d0, 0.15961d0,-0.00661d0,-0.27072d0, 1.28439d0,
-     & -0.27072d0,-0.19444d0,
-     &  0.15961d0,-0.27072d0,-0.27072d0, 0.15961d0,-0.16145d0,
-     &  0.11575d0,-0.16145d0,
-     & -0.12559d0,-0.12559d0,-0.16145d0,-0.12559d0, 2.37499d0,
-     & -0.10405d0,-0.10405d0,
-     &  0.11111d0, 0.11111d0, 0.11111d0, 0.11111d0, 0.32628d0,
-     &  0.32628d0, 0.11111d0,
-     & -0.10405d0, 0.11111d0, 0.32628d0,-0.16902d0,-0.16902d0,
-     &  0.31481d0, 0.31481d0,
-     &  0.31481d0, 0.31481d0,-0.31246d0,-0.31246d0,-0.00661d0,
-     &  0.15961d0,-0.19444d0,
-     &  0.15961d0,-0.19444d0,-0.27072d0, 1.28439d0,-0.27072d0,
-     &  0.15961d0, 0.15961d0,
-     & -0.27072d0,-0.27072d0, 0.11575d0,-0.16145d0,-0.12559d0,
-     & -0.16145d0,-0.16145d0,
-     & -0.12559d0, 2.37499d0,-0.12559d0,-0.10405d0, 0.11111d0,
-     &  0.11111d0,-0.10405d0,
-     &  0.11111d0, 0.32628d0, 0.32628d0, 0.11111d0,-0.10405d0,
-     &  0.11111d0, 0.32628d0,
-     &  0.11111d0/
+     & +2.3749905130142061615572402d+0, -1.2558804006385479445917497d-1,
+     & -1.6144899697318323972261567d-1, -1.2558804006385501650377989d-1,
+     & -1.2558804006385515528165797d-1, -1.6144899697318335074491813d-1,
+     & +1.1575022772654103286882332d-1, -1.6144899697318332298934251d-1,
+     & +3.2627685256707739114290234d-1, +1.1111111111111149352126404d-1,
+     & +1.1111111111111166005471773d-1, +3.2627685256707633643102895d-1,
+     & +1.1111111111111099392090296d-1, -1.0405463034485473716728166d-1,
+     & -1.0405463034485518125649151d-1, +1.1111111111111171556586896d-1,
+     & +3.2627685256707716909829742d-1, +1.1111111111111129923223473d-1,
+     & -1.0405463034485522288985493d-1, +1.1111111111111163229914212d-1,
+     & -3.1246265455939625788062131d-1, -3.1246265455939359334536221d-1,
+     & +3.1481481481481216055939854d-1, +3.1481481481481499162811133d-1,
+     & +3.1481481481481182749249115d-1, +3.1481481481481465856120394d-1,
+     & -1.6901882692208583525150800d-1, -1.6901882692208175518189250d-1,
+     & +1.2843861132567893079681198d+0, -2.7072129701152342073910972d-1,
+     & -1.9444444444444389241688498d-1, -2.7072129701152392033947081d-1,
+     & -1.9444444444444367037228005d-1, +1.5961018590041298192971908d-1,
+     & -6.6083354790133445183508698d-3, +1.5961018590041253784050923d-1,
+     & -2.7072129701152364278371465d-1, -2.7072129701152358727256342d-1,
+     & +1.5961018590041348153008016d-1, +1.5961018590041123332845530d-1,
+     & -1.2558804006385776430576584d-1, +2.3749905130142034970219811d+0,
+     & -1.2558804006385507201493112d-1, -1.6144899697318287890013266d-1,
+     & -1.6144899697318521036848438d-1, -1.2558804006385462792572127d-1,
+     & -1.6144899697318235154419597d-1, +1.1575022772654162961369906d-1,
+     & +3.2627685256707922301089297d-1, +3.2627685256707933403319544d-1,
+     & +1.1111111111111089677638830d-1, +1.1111111111111129923223473d-1,
+     & +1.1111111111111106330984200d-1, +1.1111111111111049432054187d-1,
+     & -1.0405463034485712414678460d-1, -1.0405463034485641637960640d-1,
+     & +1.1111111111111222904401785d-1, +3.2627685256707822381017081d-1,
+     & +1.1111111111111078575408584d-1, -1.0405463034485666617978694d-1,
+     & -3.1246265455939692401443608d-1, +3.1481481481481471407235517d-1,
+     & +3.1481481481481521367271625d-1, -3.1246265455939792321515824d-1,
+     & +3.1481481481481587980653103d-1, -1.6901882692208511360654200d-1,
+     & -1.6901882692208533565114692d-1, +3.1481481481481415896084286d-1,
+     & -2.7072129701152436442868066d-1, -1.9444444444444536346239261d-1,
+     & -2.7072129701152353176141219d-1, +1.2843861132567899741019346d+0,
+     & +1.5961018590041256559608485d-1, -6.6083354790130886466381632d-3,
+     & +1.5961018590041356479680701d-1, -1.9444444444444441977282167d-1,
+     & -2.7072129701152392033947081d-1, +1.5961018590041231579590431d-1,
+     & +1.5961018590041278764068977d-1, -2.7072129701152264358299249d-1,
+     & -4.8823919872290821864169175d-1, -4.8823919872290794108593559d-1,
+     & -4.8823919872291077215464838d-1, -4.8823919872290899579780898d-1,
+     & +2.2897993946364816886251958d-1, +2.2897993946364941786342229d-1,
+     & +2.2897993946364988970820775d-1, +2.2897993946365155504274469d-1,
+     & +5.5555555555554851643318415d-2, +5.5555555555555712066162499d-2,
+     & +5.5555555555555725943950307d-2, +5.5555555555556031255282079d-2,
+     & +5.5555555555555920232979616d-2, +5.5555555555555365121467304d-2,
+     & +5.5555555555554796132167183d-2, +5.5555555555555392877042920d-2,
+     & -2.2222222222222168253047414d-1, -2.2222222222222254295331823d-1,
+     & -2.2222222222222151599702045d-1, -2.2222222222222257070889384d-1,
+     & +3.1481481481481549122847241d-1, -3.1246265455939775668170455d-1,
+     & -3.1246265455939814525976317d-1, +3.1481481481481582429537980d-1,
+     & -1.6901882692208500258423953d-1, +3.1481481481481621287343842d-1,
+     & +3.1481481481481532469501872d-1, -1.6901882692208539116229815d-1,
+     & -2.7072129701152164438227032d-1, +1.2843861132567930827264036d+0,
+     & -2.7072129701152281011644618d-1, -1.9444444444444328179422143d-1,
+     & +1.5961018590041139986190899d-1, -1.9444444444444586306275369d-1,
+     & +1.5961018590041073372809421d-1, -6.6083354790144668844398268d-3,
+     & +1.5961018590041303744087031d-1, -2.7072129701152297664989987d-1,
+     & -2.7072129701152392033947081d-1, +1.5961018590041164966208953d-1,
+     & -1.2558804006385637652698506d-1, -1.6144899697318712550320186d-1,
+     & -1.2558804006384907681059815d-1, +2.3749905130142003883975121d+0,
+     & -1.6144899697317796616324870d-1, +1.1575022772654076919085497d-1,
+     & -1.6144899697318212949959104d-1, -1.2558804006385926310684908d-1,
+     & +1.1111111111111218741065443d-1, +1.1111111111111010574248326d-1,
+     & +3.2627685256707761318750727d-1, +3.2627685256707888994398559d-1,
+     & -1.0405463034485737394696514d-1, -1.0405463034485638862403079d-1,
+     & +1.1111111111111186822153485d-1, +1.1111111111111141025453719d-1,
+     & +1.1111111111111179883259581d-1, -1.0405463034485576412357943d-1,
+     & +1.1111111111110900939724644d-1, +3.2627685256707994465585898d-1,
+     & +3.1481481481481549122847241d-1, +3.1481481481481282669321331d-1,
+     & -3.1246265455939442601263067d-1, -3.1246265455939720157019224d-1,
+     & -1.6901882692208336500527821d-1, -1.6901882692208519687326884d-1,
+     & +3.1481481481481515816156502d-1, +3.1481481481481332629357439d-1,
+     & -1.9444444444444225483792366d-1, -2.7072129701152342073910972d-1,
+     & +1.2843861132567906402357494d+0, -2.7072129701152181091572402d-1,
+     & -6.6083354790147739304950747d-3, +1.5961018590041103903942599d-1,
+     & -1.9444444444444558550699753d-1, +1.5961018590041214926245061d-1,
+     & +1.5961018590041251008493361d-1, +1.5961018590041275988511416d-1,
+     & -2.7072129701152436442868066d-1, -2.7072129701152253256069002d-1,
+     & -1.6144899697318501607945507d-1, -1.2558804006385732021655599d-1,
+     & +2.3749905130142030529327712d+0, -1.2558804006385482221475058d-1,
+     & +1.1575022772654303127026765d-1, -1.6144899697318307318916197d-1,
+     & -1.2558804006385507201493112d-1, -1.6144899697318496056830384d-1,
+     & +1.1111111111111199312162512d-1, +3.2627685256707950056664913d-1,
+     & +3.2627685256707839034362451d-1, +1.1111111111111152127683965d-1,
+     & -1.0405463034485801232520430d-1, +1.1111111111111041105381503d-1,
+     & +1.1111111111111079963187365d-1, -1.0405463034485795681405307d-1,
+     & -1.0405463034485722129129925d-1, +1.1111111111111175719923239d-1,
+     & +3.2627685256707750216520481d-1, +1.1111111111111160454356650d-1,
+     & -3.1246265455939781219285578d-1, +3.1481481481481216055939854d-1,
+     & -1.6901882692208303193837082d-1, +3.1481481481481571327307734d-1,
+     & -3.1246265455939448152378191d-1, +3.1481481481481587980653103d-1,
+     & -1.6901882692208552994017623d-1, +3.1481481481481204953709607d-1,
+     & -2.7072129701152292113874864d-1, +1.5961018590041298192971908d-1,
+     & +1.5961018590041237130705554d-1, -2.7072129701152497505134420d-1,
+     & -2.7072129701152436442868066d-1, +1.5961018590041345377450455d-1,
+     & +1.5961018590041300968529470d-1, -2.7072129701152336522795849d-1,
+     & +1.2843861132567895300127248d+0, -1.9444444444444416997264113d-1,
+     & -6.6083354790131745154502241d-3, -1.9444444444444400343918744d-1,
+     & -4.8823919872290866273090160d-1, -4.8823919872291138277731193d-1,
+     & +2.2897993946365219342098385d-1, +2.2897993946365047257529568d-1,
+     & -4.8823919872290605370679373d-1, -4.8823919872290855170859913d-1,
+     & +2.2897993946364939010784667d-1, +2.2897993946364658679470949d-1,
+     & +5.5555555555556107583115022d-2, -2.2222222222222276499792315d-1,
+     & +5.5555555555554921032257454d-2, -2.2222222222222173804162537d-1,
+     & +5.5555555555554997360090397d-2, -2.2222222222222245968659138d-1,
+     & +5.5555555555556412894446794d-2, -2.2222222222222173804162537d-1,
+     & +5.5555555555554907154469646d-2, +5.5555555555556489222279737d-2,
+     & +5.5555555555554768376591568d-2, +5.5555555555555864721828385d-2,
+     & +3.1481481481481582429537980d-1, -3.1246265455939853383782179d-1,
+     & +3.1481481481481804474142905d-1, -1.6901882692208644587417155d-1,
+     & +3.1481481481481710105185812d-1, -3.1246265455939759014825086d-1,
+     & +3.1481481481481565776192610d-1, -1.6901882692208705649683509d-1,
+     & -2.7072129701152131131536294d-1, -2.7072129701152486402904174d-1,
+     & +1.5961018590041101128385037d-1, +1.5961018590041312070759716d-1,
+     & -2.7072129701152231051608510d-1, -2.7072129701152375380601711d-1,
+     & +1.5961018590041206599572376d-1, +1.5961018590041239906263115d-1,
+     & -1.9444444444444369812785567d-1, +1.2843861132567919725033789d+0,
+     & -1.9444444444444680675232462d-1, -6.6083354790126783845360947d-3,
+     & -4.8823919872290966193162376d-1, +2.2897993946365174933177400d-1,
+     & +2.2897993946364897377421244d-1, -4.8823919872290882926435529d-1,
+     & -4.8823919872291038357658977d-1, +2.2897993946364955664130036d-1,
+     & +2.2897993946364897377421244d-1, -4.8823919872290783006363313d-1,
+     & -2.2222222222222218213083522d-1, +5.5555555555556114522008926d-2,
+     & -2.2222222222222243193101576d-1, +5.5555555555555927171873520d-2,
+     & -2.2222222222222151599702045d-1, +5.5555555555555295732528265d-2,
+     & -2.2222222222222232090871330d-1, +5.5555555555554927971151358d-2,
+     & +5.5555555555555420632618535d-2, +5.5555555555554879398894030d-2,
+     & +5.5555555555556045133069887d-2, +5.5555555555555143076862379d-2,
+     & -2.9629629629629505860322070d-1, -2.9629629629629572473703547d-1,
+     & -2.9629629629629583575933793d-1, -2.9629629629629633535969901d-1,
+     & -2.9629629629629594678164040d-1, -2.9629629629629655740430394d-1,
+     & -2.9629629629629583575933793d-1, -2.9629629629629583575933793d-1,
+     & -1.1111111111111098004311515d-1, -1.1111111111111106330984200d-1,
+     & -1.1111111111111071636514680d-1, -1.1111111111111100779869076d-1,
+     & -1.1111111111111098004311515d-1, -1.1111111111111082738744926d-1,
+     & -1.1111111111111118820993227d-1, -1.1111111111111088289860049d-1,
+     & -1.1111111111111088289860049d-1, -1.1111111111111095228753953d-1,
+     & -1.1111111111111109106541761d-1, -1.1111111111111099392090296d-1,
+     & +2.2897993946365111095353484d-1, -4.8823919872291154931076562d-1,
+     & -4.8823919872290760801902820d-1, +2.2897993946365025053069076d-1,
+     & +2.2897993946365105544238361d-1, -4.8823919872290849619744790d-1,
+     & -4.8823919872290844068629667d-1, +2.2897993946364875172960751d-1,
+     & -2.2222222222222173804162537d-1, +5.5555555555555559410496613d-2,
+     & -2.2222222222222171028604976d-1, +5.5555555555555351243679496d-2,
+     & -2.2222222222222245968659138d-1, +5.5555555555555767577313730d-2,
+     & -2.2222222222222171028604976d-1, +5.5555555555555767577313730d-2,
+     & +5.5555555555555295732528265d-2, +5.5555555555555705127268595d-2,
+     & +5.5555555555555004298984301d-2, +5.5555555555555892477404001d-2,
+     & +3.1481481481481554673962364d-1, -1.6901882692208303193837082d-1,
+     & +3.1481481481481121686982760d-1, -3.1246265455939586930256269d-1,
+     & +3.1481481481481082829176898d-1, -1.6901882692208441971715160d-1,
+     & +3.1481481481481587980653103d-1, -3.1246265455939370436766467d-1,
+     & +1.5961018590041253784050923d-1, +1.5961018590041117781730406d-1,
+     & -2.7072129701152325420565603d-1, -2.7072129701152192193802648d-1,
+     & +1.5961018590041328724105085d-1, +1.5961018590041095577269914d-1,
+     & -2.7072129701152414238407573d-1, -2.7072129701152169989342156d-1,
+     & -1.9444444444444267117155789d-1, -6.6083354790150705682094667d-3,
+     & -1.9444444444444575204045123d-1, +1.2843861132567913063695642d+0,
+     & +2.2897993946365069461990061d-1, +2.2897993946365294282152547d-1,
+     & -4.8823919872291138277731193d-1, -4.8823919872290899579780898d-1,
+     & +2.2897993946364636475010457d-1, +2.2897993946364950113014913d-1,
+     & -4.8823919872290888477550652d-1, -4.8823919872290616472909619d-1,
+     & +5.5555555555555143076862379d-2, -2.2222222222222237641986453d-1,
+     & +5.5555555555556142277584541d-2, -2.2222222222222215437525961d-1,
+     & +5.5555555555556274116568716d-2, -2.2222222222222207110853276d-1,
+     & +5.5555555555554844704424511d-2, -2.2222222222222157150817168d-1,
+     & +5.5555555555556010438600367d-2, +5.5555555555555087565711148d-2,
+     & +5.5555555555556246360993100d-2, +5.5555555555554782254379376d-2,
+     & -1.6901882692208533565114692d-1, +3.1481481481481687900725319d-1,
+     & -3.1246265455939903343818287d-1, +3.1481481481481426998314532d-1,
+     & -1.6901882692208636260744470d-1, +3.1481481481481493611696010d-1,
+     & -3.1246265455939770117055332d-1, +3.1481481481481687900725319d-1,
+     & +1.5961018590041040066118683d-1, -2.7072129701152281011644618d-1,
+     & -2.7072129701152225500493387d-1, +1.5961018590041073372809421d-1,
+     & +1.5961018590041234355147992d-1, -2.7072129701152303216105111d-1,
+     & -2.7072129701152258807184126d-1, +1.5961018590041131659518214d-1,
+     & -6.6083354790149621479922182d-3, -1.9444444444444450303954852d-1,
+     & +1.2843861132567944149940331d+0, -1.9444444444444569652929999d-1,
+     & -1.2558804006385707041637545d-1, -1.6144899697318029763160041d-1,
+     & +1.1575022772653729974390302d-1, -1.6144899697318318421146444d-1,
+     & +2.3749905130141977238622530d+0, -1.2558804006385462792572127d-1,
+     & -1.6144899697318174092153242d-1, -1.2558804006385054785610578d-1,
+     & +1.1111111111111142413232500d-1, -1.0405463034485587514588190d-1,
+     & -1.0405463034485570861242820d-1, +1.1111111111111218741065443d-1,
+     & +3.2627685256707839034362451d-1, +1.1111111111111081350966145d-1,
+     & +1.1111111111110963389769779d-1, +3.2627685256707772420980973d-1,
+     & +3.2627685256707900096628805d-1, +1.1111111111111046656496626d-1,
+     & -1.0405463034485554207897451d-1, +1.1111111111111003635354422d-1,
+     & +3.1481481481481532469501872d-1, +3.1481481481481626838458965d-1,
+     & -1.6901882692208602954053731d-1, -1.6901882692208566871805431d-1,
+     & -3.1246265455939892241588041d-1, -3.1246265455939670196983116d-1,
+     & +3.1481481481481549122847241d-1, +3.1481481481481587980653103d-1,
+     & -1.9444444444444364261670444d-1, +1.5961018590041220477360184d-1,
+     & -6.6083354790139551410144136d-3, +1.5961018590041317621874839d-1,
+     & +1.2843861132567917504587740d+0, -2.7072129701152458647328558d-1,
+     & -1.9444444444444597408505615d-1, -2.7072129701152358727256342d-1,
+     & -2.7072129701152147784881663d-1, -2.7072129701152325420565603d-1,
+     & +1.5961018590041162190651391d-1, +1.5961018590041162190651391d-1,
+     & -1.6144899697318335074491813d-1, -1.2558804006385107521204247d-1,
+     & -1.6144899697318709774762624d-1, +1.1575022772654164349148687d-1,
+     & -1.2558804006385970719605893d-1, +2.3749905130142017206651417d+0,
+     & -1.2558804006385543283741413d-1, -1.6144899697317885434166840d-1,
+     & +1.1111111111111038329823941d-1, +1.1111111111111066085399557d-1,
+     & -1.0405463034485656903527229d-1, -1.0405463034485608331269901d-1,
+     & +3.2627685256708022221161514d-1, +3.2627685256707772420980973d-1,
+     & +1.1111111111110931470857821d-1, +1.1111111111111240945525935d-1,
+     & +1.1111111111111300620013509d-1, +3.2627685256707733563175111d-1,
+     & +1.1111111111111181271038362d-1, -1.0405463034485783191396280d-1,
+     & +3.1481481481481593531768226d-1, -1.6901882692208522462884446d-1,
+     & -1.6901882692208569647362992d-1, +3.1481481481481443651659902d-1,
+     & -3.1246265455939636890292377d-1, +3.1481481481481526918386749d-1,
+     & +3.1481481481481449202775025d-1, -3.1246265455939570276910899d-1,
+     & +1.5961018590041384235256317d-1, -6.6083354790135110518045636d-3,
+     & +1.5961018590041292641856785d-1, -1.9444444444444311526076774d-1,
+     & -2.7072129701152347625026096d-1, -1.9444444444444591857390492d-1,
+     & -2.7072129701152508607364666d-1, +1.2843861132567906402357494d+0,
+     & -2.7072129701152203296032894d-1, +1.5961018590041201048457253d-1,
+     & +1.5961018590041189946227007d-1, -2.7072129701152380931716834d-1,
+     & +2.2897993946365061135317376d-1, +2.2897993946364780804003658d-1,
+     & +2.2897993946365219342098385d-1, +2.2897993946364972317475406d-1,
+     & -4.8823919872290655330715481d-1, -4.8823919872290932886471637d-1,
+     & -4.8823919872290927335356514d-1, -4.8823919872291082766579962d-1,
+     & +5.5555555555555968805236944d-2, +5.5555555555555864721828385d-2,
+     & +5.5555555555555524716027094d-2, +5.5555555555555170832437994d-2,
+     & +5.5555555555554539393092739d-2, +5.5555555555555899416297905d-2,
+     & +5.5555555555555996560812559d-2, +5.5555555555555434510406343d-2,
+     & -2.2222222222222198784180591d-1, -2.2222222222222154375259606d-1,
+     & -2.2222222222222295928695246d-1, -2.2222222222222151599702045d-1,
+     & -1.6901882692208602954053731d-1, +3.1481481481481576878422857d-1,
+     & +3.1481481481481371487163301d-1, -1.6901882692208530789557130d-1,
+     & +3.1481481481481382589393547d-1, -3.1246265455939703503673854d-1,
+     & -3.1246265455939703503673854d-1, +3.1481481481481582429537980d-1,
+     & +1.5961018590041045617233806d-1, -1.9444444444444489161760714d-1,
+     & +1.5961018590041189946227007d-1, -6.6083354790152527141744443d-3,
+     & -2.7072129701152247704953879d-1, +1.2843861132567935268156134d+0,
+     & -2.7072129701152242153838756d-1, -1.9444444444444400343918744d-1,
+     & +1.5961018590041031739445998d-1, -2.7072129701152364278371465d-1,
+     & -2.7072129701152214398263141d-1, +1.5961018590041220477360184d-1,
+     & -1.6144899697318362830067429d-1, +1.1575022772653784097762752d-1,
+     & -1.6144899697318049192062972d-1, -1.2558804006385601570450206d-1,
+     & -1.2558804006385279605773064d-1, -1.6144899697318265685552774d-1,
+     & -1.2558804006385446139226758d-1, +2.3749905130141986120406727d+0,
+     & -1.0405463034485584739030628d-1, -1.0405463034485686046881625d-1,
+     & +1.1111111111111091065417611d-1, +1.1111111111111157678799088d-1,
+     & +1.1111111111111068860957118d-1, +1.1111111111111046656496626d-1,
+     & +3.2627685256707894545513682d-1, +3.2627685256707839034362451d-1,
+     & +1.1111111111111113269878103d-1, -1.0405463034485616657942586d-1,
+     & +1.1111111111110917593070013d-1, +3.2627685256707922301089297d-1,
+     & -1.6901882692208505809539076d-1, -1.6901882692208772263064986d-1,
+     & +3.1481481481481715656300935d-1, +3.1481481481481510265041379d-1,
+     & +3.1481481481481765616337043d-1, +3.1481481481481515816156502d-1,
+     & -3.1246265455939731259249470d-1, -3.1246265455939992161660257d-1,
+     & -6.6083354790145171914206301d-3, +1.5961018590041184395111884d-1,
+     & -1.9444444444444558550699753d-1, +1.5961018590041017861658190d-1,
+     & -1.9444444444444464181742660d-1, -2.7072129701152175540457279d-1,
+     & +1.2843861132567935268156134d+0, -2.7072129701152303216105111d-1,
+     & +1.5961018590041028963888436d-1, +1.5961018590041234355147992d-1,
+     & -2.7072129701152281011644618d-1, -2.7072129701152203296032894d-1,
+     & +1.1575022772654330882602380d-1, -1.6144899697318310094473759d-1,
+     & -1.2558804006385540508183851d-1, -1.6144899697318235154419597d-1,
+     & -1.6144899697318287890013266d-1, -1.2558804006385546059298974d-1,
+     & +2.3749905130141999443083023d+0, -1.2558804006385557161529221d-1,
+     & -1.0405463034485845641441415d-1, +1.1111111111111122984329569d-1,
+     & +1.1111111111111109106541761d-1, -1.0405463034485783191396280d-1,
+     & +1.1111111111111139637674938d-1, +3.2627685256708005567816144d-1,
+     & +3.2627685256708055527852252d-1, +1.1111111111111132698781034d-1,
+     & -1.0405463034485880335910934d-1, +1.1111111111111135474338596d-1,
+     & +3.2627685256708066630082499d-1, +1.1111111111111106330984200d-1/
 !
       data iflag /1/
 !
diff --git a/ccx/src/extrapolatefem.f b/ccx/src/extrapolatefem.f
index 4e98365..d5e564c 100644
--- a/ccx/src/extrapolatefem.f
+++ b/ccx/src/extrapolatefem.f
@@ -102,10 +102,20 @@
 !
       data a2 /  1.1455,-0.1455,1.1455,-0.1455,1.1455,-0.1455,
      &           -0.1455,1.1455,-0.1455,1.1455,-0.1455,1.1455/
-      data a4 /  1.92705, -0.30902, -0.30902, -0.30902,
-     &          -0.30902,  1.92705, -0.30902, -0.30902,
-     &          -0.30902, -0.30902,  1.92705, -0.30902,
-     &          -0.30902, -0.30902, -0.30902,  1.92705/
+
+!     Precision increased on these 8/2021 Pete Gustafson (after a course
+!     excercise demonstrating extrapolation). Current precision is
+!     ridiculous, but comes for free.
+      data a4 /
+     & +1.9270509831248423537886083d+0, -3.0901699437494745126286944d-1,
+     & -3.0901699437494745126286944d-1, -3.0901699437494745126286944d-1,
+     & -3.0901699437494745126286944d-1, +1.9270509831248423537886083d+0,
+     & -3.0901699437494745126286944d-1, -3.0901699437494745126286944d-1,
+     & -3.0901699437494745126286944d-1, -3.0901699437494745126286944d-1,
+     & +1.9270509831248423537886083d+0, -3.0901699437494745126286944d-1,
+     & -3.0901699437494745126286944d-1, -3.0901699437494745126286944d-1,
+     & -3.0901699437494745126286944d-1, +1.9270509831248423537886083d+0/
+
       data a9 / 1.63138,-0.32628,-0.32628,-0.52027, 0.10405, 0.10405,
      &         -0.32628, 1.63138,-0.32628, 0.10405,-0.52027, 0.10405,
      &         -0.32628,-0.32628, 1.63138, 0.10405, 0.10405,-0.52027,
@@ -117,19 +127,44 @@
      &          0.10405, 0.10405,-0.52027,-0.32628,-0.32628, 1.63138/
 !
 !     extrapolation from a 2x2x2=8 integration point scheme in a hex to
-!     the vertex nodes
-!    
-      data a8 /2.549,-.683,.183,-.683,-.683,.183,
-     &        -.04904,.183,-.683,2.549,-.683,.183,
-     &        .183,-.683,.183,-.04904,-.683,.183,
-     &        -.683,2.549,.183,-.04904,.183,-.683,
-     &        .183,-.683,2.549,-.683,-.04904,.183,
-     &        -.683,.183,-.683,.183,-.04904,.183,
-     &        2.549,-.683,.183,-.683,.183,-.683,
-     &        .183,-.04904,-.683,2.549,-.683,.183,
-     &        .183,-.04904,.183,-.683,-.683,.183,
-     &        -.683,2.549,-.04904,.183,-.683,.183,
-     &        .183,-.683,2.549,-.683/  
+!     the vertex nodes.  Precision increased on these 8/2021 Pete
+!     Gustafson (after a course excercise demonstrating
+!     extrapolation). Current precision is ridiculous, but comes for
+!     free.
+!     
+      data a8 /
+     & +2.5490381056766566736371260d+0, -6.8301270189221807704882394d-1,
+     & -6.8301270189221940931645349d-1, +1.8301270189221954809433157d-1,
+     & -6.8301270189221896522724364d-1, +1.8301270189221824358227764d-1,
+     & +1.8301270189221996442796581d-1, -4.9038105676658283460511711d-2,
+     & -6.8301270189221874318263872d-1, +2.5490381056766580059047556d+0,
+     & +1.8301270189221963136105842d-1, -6.8301270189221918727184857d-1,
+     & +1.8301270189221860440476064d-1, -6.8301270189221807704882394d-1,
+     & -4.9038105676659414500218048d-2, +1.8301270189221935380530226d-1,
+     & +1.8301270189222007545026827d-1, -6.8301270189221818807112641d-1,
+     & -6.8301270189221907624954611d-1, +2.5490381056766580059047556d+0,
+     & -4.9038105676658678977464234d-2, +1.8301270189221757744846286d-1,
+     & +1.8301270189221979789451211d-1, -6.8301270189221829909342887d-1,
+     & -6.8301270189221874318263872d-1, +1.8301270189221782724864340d-1,
+     & +2.5490381056766575618155457d+0, -6.8301270189222018647257073d-1,
+     & +1.8301270189222024198372196d-1, -4.9038105676655389941753782d-2,
+     & -6.8301270189222040851717566d-1, +1.8301270189221932604972665d-1,
+     & -6.8301270189221874318263872d-1, +1.8301270189221952033875596d-1,
+     & +1.8301270189221940931645349d-1, -4.9038105676658602649631291d-2,
+     & +2.5490381056766566736371260d+0, -6.8301270189221907624954611d-1,
+     & -6.8301270189221918727184857d-1, +1.8301270189221968687220965d-1,
+     & +1.8301270189221968687220965d-1, -6.8301270189221963136105842d-1,
+     & -4.9038105676657818554620150d-2, +1.8301270189221965911663403d-1,
+     & -6.8301270189221818807112641d-1, +2.5490381056766553413694965d+0,
+     & +1.8301270189222118567329289d-1, -6.8301270189221863216033626d-1,
+     & -4.9038105676657922638028708d-2, +1.8301270189221990891681457d-1,
+     & +1.8301270189221979789451211d-1, -6.8301270189222029749487319d-1,
+     & +1.8301270189221993667239019d-1, -6.8301270189221741091500917d-1,
+     & -6.8301270189222107465099043d-1, +2.5490381056766566736371260d+0,
+     & +1.8301270189221985340566334d-1, -4.9038105676658318154981231d-2,
+     & -6.8301270189221918727184857d-1, +1.8301270189222004769469265d-1,
+     & -6.8301270189221852113803379d-1, +1.8301270189221771622634094d-1,
+     & +2.5490381056766584499939654d+0, -6.8301270189221885420494118d-1/
 !
 !     extrapolation from a 14 integration point scheme in a hex to
 !     the vertex nodes
@@ -165,87 +200,282 @@
      &  0.2069E+00,0.2069E+00,0.2069E+00,0.2069E+00/
 !
 !     extrapolation from a 3x3x3=27 integration point scheme in a hex to
-!     the all nodes in a 20-node element
+!     the all nodes in a 20-node element. Precision increased on these
+!     8/2021 Pete Gustafson (after a course excercise demonstrating
+!     extrapolation). Current precision is ridiculous, but comes for
+!     free.
 !    
       data a27 /
-     &  2.37499,-0.12559,-0.16145,-0.12559,-0.12559,-0.16145, 0.11575,
-     & -0.16145, 0.32628, 0.11111, 0.11111, 0.32628, 0.11111,-0.10405,
-     & -0.10405, 0.11111, 0.32628, 0.11111,-0.10405, 0.11111,-0.31246,
-     & -0.31246, 0.31481, 0.31481, 0.31481, 0.31481,-0.16902,-0.16902,
-     &  1.28439,-0.27072,-0.19444,-0.27072,-0.19444, 0.15961,-0.00661,
-     &  0.15961,-0.27072,-0.27072, 0.15961, 0.15961,-0.12559, 2.37499,
-     & -0.12559,-0.16145,-0.16145,-0.12559,-0.16145, 0.11575, 0.32628,
-     &  0.32628, 0.11111, 0.11111, 0.11111, 0.11111,-0.10405,-0.10405,
-     &  0.11111, 0.32628, 0.11111,-0.10405,-0.31246, 0.31481, 0.31481,
-     & -0.31246, 0.31481,-0.16902,-0.16902, 0.31481,-0.27072,-0.19444,
-     & -0.27072, 1.28439, 0.15961,-0.00661, 0.15961,-0.19444,-0.27072,
-     &  0.15961, 0.15961,-0.27072,-0.48824,-0.48824,-0.48824,-0.48824,
-     &  0.22898, 0.22898, 0.22898, 0.22898, 0.05556, 0.05556, 0.05556,
-     &  0.05556, 0.05556, 0.05556, 0.05556, 0.05556,-0.22222,-0.22222,
-     & -0.22222,-0.22222, 0.31481,-0.31246,-0.31246, 0.31481,-0.16902,
-     &  0.31481, 0.31481,-0.16902,-0.27072, 1.28439,-0.27072,-0.19444,
-     &  0.15961,-0.19444, 0.15961,-0.00661, 0.15961,-0.27072,-0.27072,
-     &  0.15961,-0.12559,-0.16145,-0.12559, 2.37499,-0.16145, 0.11575,
-     & -0.16145,-0.12559, 0.11111, 0.11111, 0.32628, 0.32628,-0.10405,
-     & -0.10405, 0.11111, 0.11111, 0.11111,-0.10405, 0.11111, 0.32628,
-     &  0.31481, 0.31481,-0.31246,-0.31246,-0.16902,-0.16902, 0.31481,
-     &  0.31481,-0.19444,-0.27072, 1.28439,-0.27072,-0.00661, 0.15961,
-     & -0.19444, 0.15961, 0.15961, 0.15961,-0.27072,-0.27072,-0.16145,
-     & -0.12559, 2.37499,-0.12559, 0.11575,-0.16145,-0.12559,-0.16145,
-     &  0.11111, 0.32628, 0.32628, 0.11111,-0.10405, 0.11111, 0.11111,
-     & -0.10405,-0.10405, 0.11111, 0.32628, 0.11111,-0.31246, 0.31481,
-     & -0.16902, 0.31481,-0.31246, 0.31481,-0.16902, 0.31481,-0.27072,
-     &  0.15961, 0.15961,-0.27072,-0.27072, 0.15961, 0.15961,-0.27072,
-     &  1.28439,-0.19444,-0.00661,-0.19444,-0.48824,-0.48824, 0.22898,
-     &  0.22898,-0.48824,-0.48824, 0.22898, 0.22898, 0.05556,-0.22222,
-     &  0.05556,-0.22222, 0.05556,-0.22222, 0.05556,-0.22222, 0.05556,
-     &  0.05556, 0.05556, 0.05556, 0.31481,-0.31246, 0.31481,-0.16902,
-     &  0.31481,-0.31246, 0.31481,-0.16902,-0.27072,-0.27072, 0.15961,
-     &  0.15961,-0.27072,-0.27072, 0.15961, 0.15961,-0.19444, 1.28439,
-     & -0.19444,-0.00661,-0.48824, 0.22898, 0.22898,-0.48824,-0.48824,
-     &  0.22898, 0.22898,-0.48824,-0.22222, 0.05556,-0.22222, 0.05556,
-     & -0.22222, 0.05556,-0.22222, 0.05556, 0.05556, 0.05556, 0.05556,
-     &  0.05556,-0.29630,-0.29630,-0.29630,-0.29630,-0.29630,-0.29630,
-     & -0.29630,-0.29630,-0.11111,-0.11111,-0.11111,-0.11111,-0.11111,
-     & -0.11111,-0.11111,-0.11111,-0.11111,-0.11111,-0.11111,-0.11111,
-     &  0.22898,-0.48824,-0.48824, 0.22898, 0.22898,-0.48824,-0.48824,
-     &  0.22898,-0.22222, 0.05556,-0.22222, 0.05556,-0.22222, 0.05556,
-     & -0.22222, 0.05556, 0.05556, 0.05556, 0.05556, 0.05556, 0.31481,
-     & -0.16902, 0.31481,-0.31246, 0.31481,-0.16902, 0.31481,-0.31246,
-     &  0.15961, 0.15961,-0.27072,-0.27072, 0.15961, 0.15961,-0.27072,
-     & -0.27072,-0.19444,-0.00661,-0.19444, 1.28439, 0.22898, 0.22898,
-     & -0.48824,-0.48824, 0.22898, 0.22898,-0.48824,-0.48824, 0.05556,
-     & -0.22222, 0.05556,-0.22222, 0.05556,-0.22222, 0.05556,-0.22222,
-     &  0.05556, 0.05556, 0.05556, 0.05556,-0.16902, 0.31481,-0.31246,
-     &  0.31481,-0.16902, 0.31481,-0.31246, 0.31481, 0.15961,-0.27072,
-     & -0.27072, 0.15961, 0.15961,-0.27072,-0.27072, 0.15961,-0.00661,
-     & -0.19444, 1.28439,-0.19444,-0.12559,-0.16145, 0.11575,-0.16145,
-     &  2.37499,-0.12559,-0.16145,-0.12559, 0.11111,-0.10405,-0.10405,
-     &  0.11111, 0.32628, 0.11111, 0.11111, 0.32628, 0.32628, 0.11111,
-     & -0.10405, 0.11111, 0.31481, 0.31481,-0.16902,-0.16902,-0.31246,
-     & -0.31246, 0.31481, 0.31481,-0.19444, 0.15961,-0.00661, 0.15961,
-     &  1.28439,-0.27072,-0.19444,-0.27072,-0.27072,-0.27072, 0.15961,
-     &  0.15961,-0.16145,-0.12559,-0.16145, 0.11575,-0.12559, 2.37499,
-     & -0.12559,-0.16145, 0.11111, 0.11111,-0.10405,-0.10405, 0.32628,
-     &  0.32628, 0.11111, 0.11111, 0.11111, 0.32628, 0.11111,-0.10405,
-     &  0.31481,-0.16902,-0.16902, 0.31481,-0.31246, 0.31481, 0.31481,
-     & -0.31246, 0.15961,-0.00661, 0.15961,-0.19444,-0.27072,-0.19444,
-     & -0.27072, 1.28439,-0.27072, 0.15961, 0.15961,-0.27072, 0.22898,
-     &  0.22898, 0.22898, 0.22898,-0.48824,-0.48824,-0.48824,-0.48824,
-     &  0.05556, 0.05556, 0.05556, 0.05556, 0.05556, 0.05556, 0.05556,
-     &  0.05556,-0.22222,-0.22222,-0.22222,-0.22222,-0.16902, 0.31481,
-     &  0.31481,-0.16902, 0.31481,-0.31246,-0.31246, 0.31481, 0.15961,
-     & -0.19444, 0.15961,-0.00661,-0.27072, 1.28439,-0.27072,-0.19444,
-     &  0.15961,-0.27072,-0.27072, 0.15961,-0.16145, 0.11575,-0.16145,
-     & -0.12559,-0.12559,-0.16145,-0.12559, 2.37499,-0.10405,-0.10405,
-     &  0.11111, 0.11111, 0.11111, 0.11111, 0.32628, 0.32628, 0.11111,
-     & -0.10405, 0.11111, 0.32628,-0.16902,-0.16902, 0.31481, 0.31481,
-     &  0.31481, 0.31481,-0.31246,-0.31246,-0.00661, 0.15961,-0.19444,
-     &  0.15961,-0.19444,-0.27072, 1.28439,-0.27072, 0.15961, 0.15961,
-     & -0.27072,-0.27072, 0.11575,-0.16145,-0.12559,-0.16145,-0.16145,
-     & -0.12559, 2.37499,-0.12559,-0.10405, 0.11111, 0.11111,-0.10405,
-     &  0.11111, 0.32628, 0.32628, 0.11111,-0.10405, 0.11111, 0.32628,
-     &  0.11111/
+     & +2.3749905130142061615572402d+0, -1.2558804006385479445917497d-1,
+     & -1.6144899697318323972261567d-1, -1.2558804006385501650377989d-1,
+     & -1.2558804006385515528165797d-1, -1.6144899697318335074491813d-1,
+     & +1.1575022772654103286882332d-1, -1.6144899697318332298934251d-1,
+     & +3.2627685256707739114290234d-1, +1.1111111111111149352126404d-1,
+     & +1.1111111111111166005471773d-1, +3.2627685256707633643102895d-1,
+     & +1.1111111111111099392090296d-1, -1.0405463034485473716728166d-1,
+     & -1.0405463034485518125649151d-1, +1.1111111111111171556586896d-1,
+     & +3.2627685256707716909829742d-1, +1.1111111111111129923223473d-1,
+     & -1.0405463034485522288985493d-1, +1.1111111111111163229914212d-1,
+     & -3.1246265455939625788062131d-1, -3.1246265455939359334536221d-1,
+     & +3.1481481481481216055939854d-1, +3.1481481481481499162811133d-1,
+     & +3.1481481481481182749249115d-1, +3.1481481481481465856120394d-1,
+     & -1.6901882692208583525150800d-1, -1.6901882692208175518189250d-1,
+     & +1.2843861132567893079681198d+0, -2.7072129701152342073910972d-1,
+     & -1.9444444444444389241688498d-1, -2.7072129701152392033947081d-1,
+     & -1.9444444444444367037228005d-1, +1.5961018590041298192971908d-1,
+     & -6.6083354790133445183508698d-3, +1.5961018590041253784050923d-1,
+     & -2.7072129701152364278371465d-1, -2.7072129701152358727256342d-1,
+     & +1.5961018590041348153008016d-1, +1.5961018590041123332845530d-1,
+     & -1.2558804006385776430576584d-1, +2.3749905130142034970219811d+0,
+     & -1.2558804006385507201493112d-1, -1.6144899697318287890013266d-1,
+     & -1.6144899697318521036848438d-1, -1.2558804006385462792572127d-1,
+     & -1.6144899697318235154419597d-1, +1.1575022772654162961369906d-1,
+     & +3.2627685256707922301089297d-1, +3.2627685256707933403319544d-1,
+     & +1.1111111111111089677638830d-1, +1.1111111111111129923223473d-1,
+     & +1.1111111111111106330984200d-1, +1.1111111111111049432054187d-1,
+     & -1.0405463034485712414678460d-1, -1.0405463034485641637960640d-1,
+     & +1.1111111111111222904401785d-1, +3.2627685256707822381017081d-1,
+     & +1.1111111111111078575408584d-1, -1.0405463034485666617978694d-1,
+     & -3.1246265455939692401443608d-1, +3.1481481481481471407235517d-1,
+     & +3.1481481481481521367271625d-1, -3.1246265455939792321515824d-1,
+     & +3.1481481481481587980653103d-1, -1.6901882692208511360654200d-1,
+     & -1.6901882692208533565114692d-1, +3.1481481481481415896084286d-1,
+     & -2.7072129701152436442868066d-1, -1.9444444444444536346239261d-1,
+     & -2.7072129701152353176141219d-1, +1.2843861132567899741019346d+0,
+     & +1.5961018590041256559608485d-1, -6.6083354790130886466381632d-3,
+     & +1.5961018590041356479680701d-1, -1.9444444444444441977282167d-1,
+     & -2.7072129701152392033947081d-1, +1.5961018590041231579590431d-1,
+     & +1.5961018590041278764068977d-1, -2.7072129701152264358299249d-1,
+     & -4.8823919872290821864169175d-1, -4.8823919872290794108593559d-1,
+     & -4.8823919872291077215464838d-1, -4.8823919872290899579780898d-1,
+     & +2.2897993946364816886251958d-1, +2.2897993946364941786342229d-1,
+     & +2.2897993946364988970820775d-1, +2.2897993946365155504274469d-1,
+     & +5.5555555555554851643318415d-2, +5.5555555555555712066162499d-2,
+     & +5.5555555555555725943950307d-2, +5.5555555555556031255282079d-2,
+     & +5.5555555555555920232979616d-2, +5.5555555555555365121467304d-2,
+     & +5.5555555555554796132167183d-2, +5.5555555555555392877042920d-2,
+     & -2.2222222222222168253047414d-1, -2.2222222222222254295331823d-1,
+     & -2.2222222222222151599702045d-1, -2.2222222222222257070889384d-1,
+     & +3.1481481481481549122847241d-1, -3.1246265455939775668170455d-1,
+     & -3.1246265455939814525976317d-1, +3.1481481481481582429537980d-1,
+     & -1.6901882692208500258423953d-1, +3.1481481481481621287343842d-1,
+     & +3.1481481481481532469501872d-1, -1.6901882692208539116229815d-1,
+     & -2.7072129701152164438227032d-1, +1.2843861132567930827264036d+0,
+     & -2.7072129701152281011644618d-1, -1.9444444444444328179422143d-1,
+     & +1.5961018590041139986190899d-1, -1.9444444444444586306275369d-1,
+     & +1.5961018590041073372809421d-1, -6.6083354790144668844398268d-3,
+     & +1.5961018590041303744087031d-1, -2.7072129701152297664989987d-1,
+     & -2.7072129701152392033947081d-1, +1.5961018590041164966208953d-1,
+     & -1.2558804006385637652698506d-1, -1.6144899697318712550320186d-1,
+     & -1.2558804006384907681059815d-1, +2.3749905130142003883975121d+0,
+     & -1.6144899697317796616324870d-1, +1.1575022772654076919085497d-1,
+     & -1.6144899697318212949959104d-1, -1.2558804006385926310684908d-1,
+     & +1.1111111111111218741065443d-1, +1.1111111111111010574248326d-1,
+     & +3.2627685256707761318750727d-1, +3.2627685256707888994398559d-1,
+     & -1.0405463034485737394696514d-1, -1.0405463034485638862403079d-1,
+     & +1.1111111111111186822153485d-1, +1.1111111111111141025453719d-1,
+     & +1.1111111111111179883259581d-1, -1.0405463034485576412357943d-1,
+     & +1.1111111111110900939724644d-1, +3.2627685256707994465585898d-1,
+     & +3.1481481481481549122847241d-1, +3.1481481481481282669321331d-1,
+     & -3.1246265455939442601263067d-1, -3.1246265455939720157019224d-1,
+     & -1.6901882692208336500527821d-1, -1.6901882692208519687326884d-1,
+     & +3.1481481481481515816156502d-1, +3.1481481481481332629357439d-1,
+     & -1.9444444444444225483792366d-1, -2.7072129701152342073910972d-1,
+     & +1.2843861132567906402357494d+0, -2.7072129701152181091572402d-1,
+     & -6.6083354790147739304950747d-3, +1.5961018590041103903942599d-1,
+     & -1.9444444444444558550699753d-1, +1.5961018590041214926245061d-1,
+     & +1.5961018590041251008493361d-1, +1.5961018590041275988511416d-1,
+     & -2.7072129701152436442868066d-1, -2.7072129701152253256069002d-1,
+     & -1.6144899697318501607945507d-1, -1.2558804006385732021655599d-1,
+     & +2.3749905130142030529327712d+0, -1.2558804006385482221475058d-1,
+     & +1.1575022772654303127026765d-1, -1.6144899697318307318916197d-1,
+     & -1.2558804006385507201493112d-1, -1.6144899697318496056830384d-1,
+     & +1.1111111111111199312162512d-1, +3.2627685256707950056664913d-1,
+     & +3.2627685256707839034362451d-1, +1.1111111111111152127683965d-1,
+     & -1.0405463034485801232520430d-1, +1.1111111111111041105381503d-1,
+     & +1.1111111111111079963187365d-1, -1.0405463034485795681405307d-1,
+     & -1.0405463034485722129129925d-1, +1.1111111111111175719923239d-1,
+     & +3.2627685256707750216520481d-1, +1.1111111111111160454356650d-1,
+     & -3.1246265455939781219285578d-1, +3.1481481481481216055939854d-1,
+     & -1.6901882692208303193837082d-1, +3.1481481481481571327307734d-1,
+     & -3.1246265455939448152378191d-1, +3.1481481481481587980653103d-1,
+     & -1.6901882692208552994017623d-1, +3.1481481481481204953709607d-1,
+     & -2.7072129701152292113874864d-1, +1.5961018590041298192971908d-1,
+     & +1.5961018590041237130705554d-1, -2.7072129701152497505134420d-1,
+     & -2.7072129701152436442868066d-1, +1.5961018590041345377450455d-1,
+     & +1.5961018590041300968529470d-1, -2.7072129701152336522795849d-1,
+     & +1.2843861132567895300127248d+0, -1.9444444444444416997264113d-1,
+     & -6.6083354790131745154502241d-3, -1.9444444444444400343918744d-1,
+     & -4.8823919872290866273090160d-1, -4.8823919872291138277731193d-1,
+     & +2.2897993946365219342098385d-1, +2.2897993946365047257529568d-1,
+     & -4.8823919872290605370679373d-1, -4.8823919872290855170859913d-1,
+     & +2.2897993946364939010784667d-1, +2.2897993946364658679470949d-1,
+     & +5.5555555555556107583115022d-2, -2.2222222222222276499792315d-1,
+     & +5.5555555555554921032257454d-2, -2.2222222222222173804162537d-1,
+     & +5.5555555555554997360090397d-2, -2.2222222222222245968659138d-1,
+     & +5.5555555555556412894446794d-2, -2.2222222222222173804162537d-1,
+     & +5.5555555555554907154469646d-2, +5.5555555555556489222279737d-2,
+     & +5.5555555555554768376591568d-2, +5.5555555555555864721828385d-2,
+     & +3.1481481481481582429537980d-1, -3.1246265455939853383782179d-1,
+     & +3.1481481481481804474142905d-1, -1.6901882692208644587417155d-1,
+     & +3.1481481481481710105185812d-1, -3.1246265455939759014825086d-1,
+     & +3.1481481481481565776192610d-1, -1.6901882692208705649683509d-1,
+     & -2.7072129701152131131536294d-1, -2.7072129701152486402904174d-1,
+     & +1.5961018590041101128385037d-1, +1.5961018590041312070759716d-1,
+     & -2.7072129701152231051608510d-1, -2.7072129701152375380601711d-1,
+     & +1.5961018590041206599572376d-1, +1.5961018590041239906263115d-1,
+     & -1.9444444444444369812785567d-1, +1.2843861132567919725033789d+0,
+     & -1.9444444444444680675232462d-1, -6.6083354790126783845360947d-3,
+     & -4.8823919872290966193162376d-1, +2.2897993946365174933177400d-1,
+     & +2.2897993946364897377421244d-1, -4.8823919872290882926435529d-1,
+     & -4.8823919872291038357658977d-1, +2.2897993946364955664130036d-1,
+     & +2.2897993946364897377421244d-1, -4.8823919872290783006363313d-1,
+     & -2.2222222222222218213083522d-1, +5.5555555555556114522008926d-2,
+     & -2.2222222222222243193101576d-1, +5.5555555555555927171873520d-2,
+     & -2.2222222222222151599702045d-1, +5.5555555555555295732528265d-2,
+     & -2.2222222222222232090871330d-1, +5.5555555555554927971151358d-2,
+     & +5.5555555555555420632618535d-2, +5.5555555555554879398894030d-2,
+     & +5.5555555555556045133069887d-2, +5.5555555555555143076862379d-2,
+     & -2.9629629629629505860322070d-1, -2.9629629629629572473703547d-1,
+     & -2.9629629629629583575933793d-1, -2.9629629629629633535969901d-1,
+     & -2.9629629629629594678164040d-1, -2.9629629629629655740430394d-1,
+     & -2.9629629629629583575933793d-1, -2.9629629629629583575933793d-1,
+     & -1.1111111111111098004311515d-1, -1.1111111111111106330984200d-1,
+     & -1.1111111111111071636514680d-1, -1.1111111111111100779869076d-1,
+     & -1.1111111111111098004311515d-1, -1.1111111111111082738744926d-1,
+     & -1.1111111111111118820993227d-1, -1.1111111111111088289860049d-1,
+     & -1.1111111111111088289860049d-1, -1.1111111111111095228753953d-1,
+     & -1.1111111111111109106541761d-1, -1.1111111111111099392090296d-1,
+     & +2.2897993946365111095353484d-1, -4.8823919872291154931076562d-1,
+     & -4.8823919872290760801902820d-1, +2.2897993946365025053069076d-1,
+     & +2.2897993946365105544238361d-1, -4.8823919872290849619744790d-1,
+     & -4.8823919872290844068629667d-1, +2.2897993946364875172960751d-1,
+     & -2.2222222222222173804162537d-1, +5.5555555555555559410496613d-2,
+     & -2.2222222222222171028604976d-1, +5.5555555555555351243679496d-2,
+     & -2.2222222222222245968659138d-1, +5.5555555555555767577313730d-2,
+     & -2.2222222222222171028604976d-1, +5.5555555555555767577313730d-2,
+     & +5.5555555555555295732528265d-2, +5.5555555555555705127268595d-2,
+     & +5.5555555555555004298984301d-2, +5.5555555555555892477404001d-2,
+     & +3.1481481481481554673962364d-1, -1.6901882692208303193837082d-1,
+     & +3.1481481481481121686982760d-1, -3.1246265455939586930256269d-1,
+     & +3.1481481481481082829176898d-1, -1.6901882692208441971715160d-1,
+     & +3.1481481481481587980653103d-1, -3.1246265455939370436766467d-1,
+     & +1.5961018590041253784050923d-1, +1.5961018590041117781730406d-1,
+     & -2.7072129701152325420565603d-1, -2.7072129701152192193802648d-1,
+     & +1.5961018590041328724105085d-1, +1.5961018590041095577269914d-1,
+     & -2.7072129701152414238407573d-1, -2.7072129701152169989342156d-1,
+     & -1.9444444444444267117155789d-1, -6.6083354790150705682094667d-3,
+     & -1.9444444444444575204045123d-1, +1.2843861132567913063695642d+0,
+     & +2.2897993946365069461990061d-1, +2.2897993946365294282152547d-1,
+     & -4.8823919872291138277731193d-1, -4.8823919872290899579780898d-1,
+     & +2.2897993946364636475010457d-1, +2.2897993946364950113014913d-1,
+     & -4.8823919872290888477550652d-1, -4.8823919872290616472909619d-1,
+     & +5.5555555555555143076862379d-2, -2.2222222222222237641986453d-1,
+     & +5.5555555555556142277584541d-2, -2.2222222222222215437525961d-1,
+     & +5.5555555555556274116568716d-2, -2.2222222222222207110853276d-1,
+     & +5.5555555555554844704424511d-2, -2.2222222222222157150817168d-1,
+     & +5.5555555555556010438600367d-2, +5.5555555555555087565711148d-2,
+     & +5.5555555555556246360993100d-2, +5.5555555555554782254379376d-2,
+     & -1.6901882692208533565114692d-1, +3.1481481481481687900725319d-1,
+     & -3.1246265455939903343818287d-1, +3.1481481481481426998314532d-1,
+     & -1.6901882692208636260744470d-1, +3.1481481481481493611696010d-1,
+     & -3.1246265455939770117055332d-1, +3.1481481481481687900725319d-1,
+     & +1.5961018590041040066118683d-1, -2.7072129701152281011644618d-1,
+     & -2.7072129701152225500493387d-1, +1.5961018590041073372809421d-1,
+     & +1.5961018590041234355147992d-1, -2.7072129701152303216105111d-1,
+     & -2.7072129701152258807184126d-1, +1.5961018590041131659518214d-1,
+     & -6.6083354790149621479922182d-3, -1.9444444444444450303954852d-1,
+     & +1.2843861132567944149940331d+0, -1.9444444444444569652929999d-1,
+     & -1.2558804006385707041637545d-1, -1.6144899697318029763160041d-1,
+     & +1.1575022772653729974390302d-1, -1.6144899697318318421146444d-1,
+     & +2.3749905130141977238622530d+0, -1.2558804006385462792572127d-1,
+     & -1.6144899697318174092153242d-1, -1.2558804006385054785610578d-1,
+     & +1.1111111111111142413232500d-1, -1.0405463034485587514588190d-1,
+     & -1.0405463034485570861242820d-1, +1.1111111111111218741065443d-1,
+     & +3.2627685256707839034362451d-1, +1.1111111111111081350966145d-1,
+     & +1.1111111111110963389769779d-1, +3.2627685256707772420980973d-1,
+     & +3.2627685256707900096628805d-1, +1.1111111111111046656496626d-1,
+     & -1.0405463034485554207897451d-1, +1.1111111111111003635354422d-1,
+     & +3.1481481481481532469501872d-1, +3.1481481481481626838458965d-1,
+     & -1.6901882692208602954053731d-1, -1.6901882692208566871805431d-1,
+     & -3.1246265455939892241588041d-1, -3.1246265455939670196983116d-1,
+     & +3.1481481481481549122847241d-1, +3.1481481481481587980653103d-1,
+     & -1.9444444444444364261670444d-1, +1.5961018590041220477360184d-1,
+     & -6.6083354790139551410144136d-3, +1.5961018590041317621874839d-1,
+     & +1.2843861132567917504587740d+0, -2.7072129701152458647328558d-1,
+     & -1.9444444444444597408505615d-1, -2.7072129701152358727256342d-1,
+     & -2.7072129701152147784881663d-1, -2.7072129701152325420565603d-1,
+     & +1.5961018590041162190651391d-1, +1.5961018590041162190651391d-1,
+     & -1.6144899697318335074491813d-1, -1.2558804006385107521204247d-1,
+     & -1.6144899697318709774762624d-1, +1.1575022772654164349148687d-1,
+     & -1.2558804006385970719605893d-1, +2.3749905130142017206651417d+0,
+     & -1.2558804006385543283741413d-1, -1.6144899697317885434166840d-1,
+     & +1.1111111111111038329823941d-1, +1.1111111111111066085399557d-1,
+     & -1.0405463034485656903527229d-1, -1.0405463034485608331269901d-1,
+     & +3.2627685256708022221161514d-1, +3.2627685256707772420980973d-1,
+     & +1.1111111111110931470857821d-1, +1.1111111111111240945525935d-1,
+     & +1.1111111111111300620013509d-1, +3.2627685256707733563175111d-1,
+     & +1.1111111111111181271038362d-1, -1.0405463034485783191396280d-1,
+     & +3.1481481481481593531768226d-1, -1.6901882692208522462884446d-1,
+     & -1.6901882692208569647362992d-1, +3.1481481481481443651659902d-1,
+     & -3.1246265455939636890292377d-1, +3.1481481481481526918386749d-1,
+     & +3.1481481481481449202775025d-1, -3.1246265455939570276910899d-1,
+     & +1.5961018590041384235256317d-1, -6.6083354790135110518045636d-3,
+     & +1.5961018590041292641856785d-1, -1.9444444444444311526076774d-1,
+     & -2.7072129701152347625026096d-1, -1.9444444444444591857390492d-1,
+     & -2.7072129701152508607364666d-1, +1.2843861132567906402357494d+0,
+     & -2.7072129701152203296032894d-1, +1.5961018590041201048457253d-1,
+     & +1.5961018590041189946227007d-1, -2.7072129701152380931716834d-1,
+     & +2.2897993946365061135317376d-1, +2.2897993946364780804003658d-1,
+     & +2.2897993946365219342098385d-1, +2.2897993946364972317475406d-1,
+     & -4.8823919872290655330715481d-1, -4.8823919872290932886471637d-1,
+     & -4.8823919872290927335356514d-1, -4.8823919872291082766579962d-1,
+     & +5.5555555555555968805236944d-2, +5.5555555555555864721828385d-2,
+     & +5.5555555555555524716027094d-2, +5.5555555555555170832437994d-2,
+     & +5.5555555555554539393092739d-2, +5.5555555555555899416297905d-2,
+     & +5.5555555555555996560812559d-2, +5.5555555555555434510406343d-2,
+     & -2.2222222222222198784180591d-1, -2.2222222222222154375259606d-1,
+     & -2.2222222222222295928695246d-1, -2.2222222222222151599702045d-1,
+     & -1.6901882692208602954053731d-1, +3.1481481481481576878422857d-1,
+     & +3.1481481481481371487163301d-1, -1.6901882692208530789557130d-1,
+     & +3.1481481481481382589393547d-1, -3.1246265455939703503673854d-1,
+     & -3.1246265455939703503673854d-1, +3.1481481481481582429537980d-1,
+     & +1.5961018590041045617233806d-1, -1.9444444444444489161760714d-1,
+     & +1.5961018590041189946227007d-1, -6.6083354790152527141744443d-3,
+     & -2.7072129701152247704953879d-1, +1.2843861132567935268156134d+0,
+     & -2.7072129701152242153838756d-1, -1.9444444444444400343918744d-1,
+     & +1.5961018590041031739445998d-1, -2.7072129701152364278371465d-1,
+     & -2.7072129701152214398263141d-1, +1.5961018590041220477360184d-1,
+     & -1.6144899697318362830067429d-1, +1.1575022772653784097762752d-1,
+     & -1.6144899697318049192062972d-1, -1.2558804006385601570450206d-1,
+     & -1.2558804006385279605773064d-1, -1.6144899697318265685552774d-1,
+     & -1.2558804006385446139226758d-1, +2.3749905130141986120406727d+0,
+     & -1.0405463034485584739030628d-1, -1.0405463034485686046881625d-1,
+     & +1.1111111111111091065417611d-1, +1.1111111111111157678799088d-1,
+     & +1.1111111111111068860957118d-1, +1.1111111111111046656496626d-1,
+     & +3.2627685256707894545513682d-1, +3.2627685256707839034362451d-1,
+     & +1.1111111111111113269878103d-1, -1.0405463034485616657942586d-1,
+     & +1.1111111111110917593070013d-1, +3.2627685256707922301089297d-1,
+     & -1.6901882692208505809539076d-1, -1.6901882692208772263064986d-1,
+     & +3.1481481481481715656300935d-1, +3.1481481481481510265041379d-1,
+     & +3.1481481481481765616337043d-1, +3.1481481481481515816156502d-1,
+     & -3.1246265455939731259249470d-1, -3.1246265455939992161660257d-1,
+     & -6.6083354790145171914206301d-3, +1.5961018590041184395111884d-1,
+     & -1.9444444444444558550699753d-1, +1.5961018590041017861658190d-1,
+     & -1.9444444444444464181742660d-1, -2.7072129701152175540457279d-1,
+     & +1.2843861132567935268156134d+0, -2.7072129701152303216105111d-1,
+     & +1.5961018590041028963888436d-1, +1.5961018590041234355147992d-1,
+     & -2.7072129701152281011644618d-1, -2.7072129701152203296032894d-1,
+     & +1.1575022772654330882602380d-1, -1.6144899697318310094473759d-1,
+     & -1.2558804006385540508183851d-1, -1.6144899697318235154419597d-1,
+     & -1.6144899697318287890013266d-1, -1.2558804006385546059298974d-1,
+     & +2.3749905130141999443083023d+0, -1.2558804006385557161529221d-1,
+     & -1.0405463034485845641441415d-1, +1.1111111111111122984329569d-1,
+     & +1.1111111111111109106541761d-1, -1.0405463034485783191396280d-1,
+     & +1.1111111111111139637674938d-1, +3.2627685256708005567816144d-1,
+     & +3.2627685256708055527852252d-1, +1.1111111111111132698781034d-1,
+     & -1.0405463034485880335910934d-1, +1.1111111111111135474338596d-1,
+     & +3.2627685256708066630082499d-1, +1.1111111111111106330984200d-1/
 !
       data iflag /1/
 !
diff --git a/ccx/src/frd.c b/ccx/src/frd.c
index 97b0520..bef1122 100644
--- a/ccx/src/frd.c
+++ b/ccx/src/frd.c
@@ -1,5 +1,7 @@
 /*     CalculiX - A 3-dimensional finite element program                 */
-/*              Copyright (C) 1998-2020 Guido Dhondt                          */
+/*              Copyright (C) 1998-2020 Guido Dhondt                     */
+/*     Modifications of this subroutine                                  */
+/*              Copyright (C) 2013-2021 Peter A. Gustafson               */
 
 /*     This program is free software; you can redistribute it and/or     */
 /*     modify it under the terms of the GNU General Public License as    */
@@ -20,6 +22,9 @@
 #include <stdio.h>
 #include <string.h>
 #include "CalculiX.h"
+#ifdef EXODUSII
+#include "exo.h"
+#endif
 
 #define min(a,b) ((a) <= (b) ? (a) : (b))
 #define max(a,b) ((a) >= (b) ? (a) : (b))
@@ -88,6 +93,20 @@ void frd(double *co,ITG *nk,ITG *kon,ITG *ipkon,char *lakon,ITG *ne0,
 
   double pi,oner,*errn=NULL;
 
+#ifdef EXODUSII
+  if(strcmp1(output,"exo")==0){
+    exo(co,nk,kon,ipkon,lakon,ne0,v,stn,inum,nmethod,kode,
+	filab,een,t1,fn,time,epn,ielmat,matname,enern,
+	xstaten,nstate_,istep,iinc,ithermal,qfn,mode,noddiam,
+	trab,inotr,ntrans,orab,ielorien,norien,description,
+	ipneigh,neigh,mi,stx,vr,vi,stnr,stni,vmax,stnmax,
+	ngraph,veold,ener,ne,cs,set,nset,istartset,iendset,
+	ialset,eenmax,fnr,fni,emn,thicke,jobnamec,output,qfx,
+	cdn,mortar,cdnr,cdni,nmat,ielprop,prop);
+    return;
+  }
+#endif
+
   strcpy(fneig,jobnamec);
   strcat(fneig,".frd");
 
diff --git a/ccx/src/frdselect.c b/ccx/src/frdselect.c
index 6992a30..003e33a 100644
--- a/ccx/src/frdselect.c
+++ b/ccx/src/frdselect.c
@@ -50,7 +50,8 @@ void frdselect(double *field1,double *field2,ITG *iset,ITG *nkcoords,ITG *inum,
   
   if(*iset==0){
     for(i=0;i<*nkcoords;i++){
-
+      /* i+1 is the node ID */
+      
       /* check whether output is requested for solid nodes or
          network nodes */
 
@@ -75,11 +76,11 @@ void frdselect(double *field1,double *field2,ITG *iset,ITG *nkcoords,ITG *inum,
 	for(n=1;n<=(ITG)((*ncomp+5)/6);n++){
 	  if(n==1){
 	    if(strcmp1(output,"asc")==0){
-	      fprintf(f1,"%3s%10" ITGFORMAT "",m1,i+1);
+	      fprintf(f1,"%3s%10" ITGFORMAT "",m1,i+1); /* i+1 is the node number??? */
 	    }else{
 	      iw=(int)(i+1);fwrite(&iw,sizeof(int),1,f1);
 	    }
-	    for(j=0;j<min(6,*ncomp);j++){
+	    for(j=0;j<min(6,*ncomp);j++){ /* j is the number of components */
 	      if(ifield[j]==1){
 		if(strcmp1(output,"asc")==0){
 		    fprintf(f1,"%12.5E",(float)field1[i*nfield[0]+icomp[j]]);
diff --git a/ccx/src/ini_cal.c b/ccx/src/ini_cal.c
index e49689f..6f71ad2 100644
--- a/ccx/src/ini_cal.c
+++ b/ccx/src/ini_cal.c
@@ -188,7 +188,9 @@ void ini_cal(char *jobnamec,char *output,char *fneig,char *kind1,char *kind2,
   *nprint=0;
   *kode=0;
   for(i=0;i<2;i++) jout[i]=1;
-  strcpy1(output,"asc ",4);
+
+  /* Deactived by Pete to allow exodus output */
+  /* strcpy1(output,"asc ",4); */
 
   /* info on the kind of nonlinearity of the step */
   
