function mfs_mergex(fnin1, fnin2, fnout)

# usage: mfs_mergex(fnin1, fnin2, fnout)
#
# Input  fnin1     Name of first input file
#        fnin2     Name of second input file
#        fnout     Name of output file
#
# The function merges two msh-files into one.
#
# The following datablocks are supported:
#
#   PhysicalNames
#   Nodes
#   Elements
#   NodeData
#
# -----------------------------------------------------------------------------

# Check arguments

  if (nargin != 3 || nargout != 0)
     print_usage();
  end

# Read the first input file

  [fid, vers, rc] = mfs_mshfileopen(fnin1);
  if (rc) return; end

  printf("    mfs_merge: reading file %s\n", fnin1);
  mfs_paths("add", "mfs_mergex.m", vers);
  data1 = mfs_merge_read(fid);
  mfs_paths("remove");

  fclose(fid);

# Read the second input file

  [fid, vers, rc] = mfs_mshfileopen(fnin2);
  if (rc) return; end

  printf("    mfs_merge: reading file %s\n", fnin2);
  mfs_paths("add", "mfs_mergex.m", vers);
  data2 = mfs_merge_read(fid);
  mfs_paths("remove");

  fclose(fid);

# Merge the data

  data = struct();

  if (isfield(data1, "nprops"))
     prpinc = data1.nprops;
     if (isfield(data2, "nprops"))
        data.nprops = data1.nprops + data2.nprops;        
        for n = 1 : data2.nprops
            data2.props{2, n} += prpinc;
        end
        data.props = [data1.props, data2.props];
     else
        data.nprops = data1.nprops;
        data.props  = data1.props;
     end
  else
     prpinc = 0;
     if (isfield(data2, "nprops"))
        data.nprops = data2.nprops;
        data.props  = data2.props;
     end
  end

  if (isfield(data1, "nnodes"))
     nodinc = data1.maxnod;
     if (isfield(data2, "nnodes"))
        data.nnodes = data1.nnodes + data2.nnodes;
        for n = 1 : data2.nnodes
            data2.nodes{1, n} += nodinc;
        end
        data.nodes = [data1.nodes, data2.nodes];
     else
        fprintf("*E* mfs_merge: $Nodes only in file %s\n", ...
                fnin1);
        return;
     end
  end

  if (isfield(data1, "nelem"))
     elminc = data1.maxelt;
     if (isfield(data2, "nelem"))
        data.nelem = data1.nelem + data2.nelem;
        for n = 1 : data2.nelem
            data2.elements{n}{1}(1) += elminc;
            data2.elements{n}{1}(3) += prpinc;
            data2.elements{n}{1}(4) += prpinc;
            nofk = length(data2.elements{n}{2});
            for k = 1 : nofk
                data2.elements{n}{2}(k) += nodinc;
            end
        end
        data.elements = [data1.elements, data2.elements];
     else
        fprintf("*E* mfs_merge: $Elements only in file %s\n", ...
                fnin1);
        return;
     end
  end

  if (data1.nblocks)
     if (data2.nblocks)
        if (data1.nblocks != data2.nblocks)
           printf("*E* mfs_merge: inconsistent number of $NodeData blocks\n");
           printf("               file %s contains %d blocks\n", ...
                  fnin1, data1.nblocks);
           printf("               file %s contains %d blocks\n", ...
                  fnin2, data2.nblocks);
           return;
        end
        nb     = data1.nblocks;
        nodinc = data1.nodedata{1}.id(end);
        n1     = data1.nodedata{1}.nofnod;
        n2     = data2.nodedata{1}.nofnod;
        data.nblocks = nb;
        for n = 1 : nb
            if (data1.nodedata{n}.tstpno != data2.nodedata{n}.tstpno)
               printf("*E* mfs_merge: inconsistent time step numbers");
               printf(" (%d/%d)\n", data1.nodedata{n}.tstpno,
                                    data2.nodedata{n}.tstpno);
               printf("               Data block titles = \n");
               printf("               %s\n", data1.nodedata{n}.title);
               printf("               %s\n", data2.nodedata{n}.title);
               return;
            end
            nodedata.title  = data2.nodedata{n}.title;
            nodedata.tstpno = data1.nodedata{n}.tstpno;
            nodedata.ndata  = data1.nodedata{n}.ndata;
            nodedata.nofnod = n1 + n2;
            for k = 1 : n2
                id(k) = data2.nodedata{n}.id(k) + nodinc;
            end
            nodedata.id = [data1.nodedata{n}.id, id];
            nodedata.data = [data1.nodedata{n}.data, ...
                             data2.nodedata{n}.data];
            data.nodedata{n} = nodedata;
        end
     else
        printf("*E* mfs_merge: $NodeData only in file %s\n", fnin1);
        return;
     end
  end

# Write the output file

  if ((fout = fopen(fnout, "wt")) == -1)
     printf("*E* mfs_merge: could not open file %s\n", fnout);
     return; 
  end

  printf("    mfs_merge: writing file %s\n", fnout);

  fprintf(fout, "$MeshFormat\n2.2 0 8\n$EndMeshFormat\n");

  if (isfield(data, "nprops"))
     fprintf(fout, "$PhysicalNames\n");
     fprintf(fout, "%d\n", data.nprops);
     for n = 1 : data.nprops
         fprintf(fout, "%d %d %s\n", data.props{:, n});
     end
     fprintf(fout, "$EndPhysicalNames\n");
  end

  if (isfield(data, "nnodes"))
     fprintf(fout, "$Nodes\n");
     fprintf(fout, "%d\n", data.nnodes);
     for n = 1 : data.nnodes
         fprintf(fout, "%d %f %f %f\n", data.nodes{:, n});
     end
     fprintf(fout ,"$EndNodes\n");
  end

  if (isfield(data, "nelem"))
     fprintf(fout, "$Elements\n");
     fprintf(fout, "%d\n", data.nelem);
     for n = 1 : data.nelem
         fprintf(fout, "%d %d 2 %d %d", data.elements{n}{1});
         nofk = length(data.elements{n}{2});
         for k = 1 : nofk
             fprintf(fout, " %d", data.elements{n}{2}(k));
         end
         fprintf(fout, "\n");
     end
     fprintf(fout, "$EndElements\n");
  end

  if (isfield(data, "nodedata"))
     for b = 1 : data.nblocks
         fprintf(fout, "$NodeData\n");
         fprintf(fout, "1\n%s\n", data.nodedata{b}.title);
         fprintf(fout, "1\n0.000\n3\n%d\n%d\n", 
                 data.nodedata{b}.tstpno, data.nodedata{b}.ndata); 
         nofnod = data.nodedata{b}.nofnod;
         fprintf(fout, "%d\n", nofnod);
         for k = 1 : nofnod
             fprintf(fout, "%d", data.nodedata{b}.id(k));
             for l = 1 : data.nodedata{b}.ndata
                 fprintf(fout, " %e", data.nodedata{b}.data{k}(l));
             end
             fprintf(fout, "\n");
         end
         fprintf(fout, "$EndNodeData\n");
     end
  end

  fclose(fout);

end
