function data = mfs_merge_read(fid)

# usage: data = mfs_merge_read(fid)
#
# Input  fname  File handle of file to read
#
# Output data   Structure with contents of input file
#
# The function reads a file in Gmsh msh format 4.1 The following
# datablocks are supported:
#
#   PhysicalNames
#   Entities
#   Nodes
#   Elements
#   NodeData
#
# --------------------------------------------------------------------

  nphysTagPos = [5, 8, 8, 8]; % position of number of physical tags

# Check arguments

  if (nargin != 1 || nargout != 1)
     print_usage();
  end

# Initialize

  data         = struct();
  data.nblocks = 0;         % Number of $NodeData blocks

# Process the contents

  while ((line = fgetl(fid)) != -1)

     key = sscanf(line, "%s");

     switch key

     case "$PhysicalNames"
        printf("               found datablock $PhyscialNames\n");
        data.nprops = fscanf(fid, "%d\n", 1);
        for n = 1 : data.nprops
            [props{1, n}, props{2, n}, props{3, n}] = ...
            fscanf(fid, "%d %d %s\n", "C");
        end
        data.props = props;
        line = fgetl(fid);

     case "$Entities"
        printf("               found datablock $Entities\n");
        num = fscanf(fid, "%d %d %d %d\n", 4);
        for n = 1 : 4
            np  = nphysTagPos(n);
            ids = zeros(num(n), 2);
            entities(n) = struct("num", num(n), "ids", ids);
            for m = 1 : num(n)
                line  = fgetl(fid);
                words = strsplit(line);
                ids(m, 1) = sscanf(words{1}, "%d"); 
                ntags = sscanf(words{np}, "%d");
                if (ntags)
                   ids(m, 2) = sscanf(words{np+1}, "%d");
                end
            end
            entities(n).ids = sortrows(ids);
        end
        line = fgetl(fid);

     case "$Nodes"
        printf("               found datablock $Nodes\n");
        line  = fgetl(fid);
        words = strsplit(line);
        numBlocks   = sscanf(words{1}, "%d");
        data.nnodes = sscanf(words{2}, "%d");
        data.maxnod = sscanf(words{4}, "%d");
        nbase = 0;
        for nb = 1 : numBlocks
            line  = fgetl(fid);
            words = strsplit(line);
            np    = sscanf(words{4}, "%d");
            for n = 1 : np
                line = fgetl(fid);
                ids(nbase + n) = sscanf(line, "%d");
            end
            for n = 1 : np
                line = fgetl(fid);
                coor(nbase + n, 1 : 3) = sscanf(line, "%f", 3);
            end
            nbase += np;
        end
        [ids, ix] = sort(ids);
        coor      = coor(ix, :);
        data.nodes = [num2cell(ids); num2cell(coor')];
        line = fgetl(fid);

     case "$Elements"
        printf("               found datablock $Elements\n");
        line  = fgetl(fid);
        words = strsplit(line);
        numBlocks   = sscanf(words{1}, "%d");
        data.nelem  = sscanf(words{2}, "%d");
        data.maxelt = sscanf(words{4}, "%d");
        data.elements = cell(1, data.nelem);
        nbase = 0;
        for nb = 1 : numBlocks
            [entityDim, entityTag, elemType, ne] = ...
            fscanf(fid, "%d %d %d %d\n", "C");
            ixe = entityDim + 1;
            ixg = lookup(entities(ixe).ids(:, 1), entityTag);
            grpid = entities(ixe).ids(ixg, 2);
            [ids, nofnod, nodes] = mfs_readelt(fid, ne);
            elem = ids';
            elem(:, 2) = elemType; elem(:, 3 : 4) = grpid;
            for n = 1 : ne
                element = {elem(n, :), nodes(n, :)};
                data.elements{nbase + n} = element;
            end
            nbase += ne;
        end
        line = fgetl(fid);

     case "$NodeData"
        printf("               found datablock $NodeData\n");
        nb = ++data.nblocks;
        line = fgetl(fid);
        nodedata.title = fgetl(fid);
        line = fgetl(fid); line = fgetl(fid); line = fgetl(fid);
        nodedata.tstpno = fscanf(fid, "%d/n", 1);
        nodedata.ndata  = fscanf(fid, "%d\n", 1);
        nodedata.nofnod = fscanf(fid, "%d\n", 1);
        for n = 1 : nodedata.nofnod
            nodedata.id(n) = fscanf(fid, "%d", 1);
            val = fscanf(fid, "%f\n", nodedata.ndata);
            nodedata.data{n} = val;
        end
        data.nodedata{nb} = nodedata;
        line = fgetl(fid);

     otherwise

        mfs_skip(fid);

     end

  end

end
