function [load, rc] = mfs_new_load(msg, nodes, elements, dofs, inload);

# usage: [load, rc] = mfs_new_load(msg, nodes, elements, dofs, inload);
#
# Input  msg         File handle of message file
#        nodes       Structure with nodal point data
#        elements    Structure with element data
#        dofs        Structure with degree of fredom data
#        inload      Structure with input load data
#
# Output load        Structure with loads
#        rc          Return code: 0 means no erros
#                                 1 means errors
#
# The function builds the load structure.
#
# Input data:
#
#   inload.point(:)   Nodal point forces or moments
#   inload.disp(:)    Prescribed displacements
#   inload.velo(:)    Prescribed velocities
#   inload.acce(:)    Prescribed accelerations
#   inload.inertia(:) Inertia loads
#
# Output data:
#
#   load             see Programmer's Manual
#
# -----------------------------------------------------------------------------

# Initialize

  ntypes   = 5;
  typelist = {"point", "disp", "velo", "acce", "inertia"};
  fields   = {"id", "data", "lc"};
  needf    = {[1, 2], [1, 2], [1, 2], [1, 2], 2};
  legalf   = {1 : 3, 1 : 3, 1 : 3, 1 : 3, 2 : 3}; 
  lbits    = 1 : ntypes;
  nloads   = zeros(ntypes, 1);
  uva      = struct("name", {"u", "v", "a"});

  rc       = 0;
  load     = struct();

# Check arguments

  if (nargin != 5 || nargout != 2)
     print_usage();
  end

# Check which loads are defined

  typesdef = fieldnames(inload);
  for k = 1 : length(typesdef)
      t = typesdef{k};
      ix = find(strcmp(t, typelist) == 1);
      if (length(ix) == 1)
         ix1 = ix(1);
         nloads(ix1) = 1; 
      else
         fprintf(msg, "*E* Unknown load type ""%s"" found\n", t);
         rc = 1;
      end
  end

# Get number of load cases and check definition of required fields

  nofldc = 1;
  rcf    = 0;

  for k = 1 : ntypes

      if (nloads(k))

         loadtype = getfield(inload, typelist{k});
         errmsg = sprintf("loads.%s", typelist{k});
         rcf = mfs_checkfields1(msg, loadtype, fields(needf{k}), errmsg);
         rcf = rcf || mfs_checkfields4(msg, loadtype, fields(legalf{k}), 
                                       errmsg);

         if (rcf)
            rc = 1;
         else
            nloads(k) = length(loadtype);
            if (! isfield(loadtype, "lc"))
               loadtype(1).lc = 1;
            end
            for n = 1 : nloads(k)
                needid = ismember(1, needf{k});
                if (needid && isempty(loadtype(n).id))
                   fprintf(msg, "*E* Nodal point identifier not defined ");
                   fprintf(msg, " for loads.%s(%2d)\n", typelist{k}, n);
                   rc = 1;
                end
                if (isempty(loadtype(n).data))
                   fprintf(msg, "*E* Data not defined ");
                   fprintf(msg, " for loads.%s(%2d)\n", typelist{k}, n);
                   rc = 1;
                end
                if (isempty(loadtype(n).lc))
                   loadtype(n).lc = 1;
                else
                   nofldc = max(nofldc, loadtype(n).lc);
                end
            end
            inload = setfield(inload, typelist{k}, loadtype);
         end

      end

  end

  enfmot = nloads(2) || nloads(3) || nloads(4);

  if (enfmot && ! isfield(dofs, "dofp"))
     fprintf(msg, "*E* Prescribed motion defined but");
     fprintf(msg, " prescribed degrees of freedom undefined\n")
     rc = 1;
  end

  if (rc) return; end

# Initialize load case information (bits correspond to loadtypes)

  inflc = zeros(nofldc, 1);

# Process point loads

  if (nloads(1))

     [idp, ixp]   = sort([inload.point.id]);
     inload.point = inload.point(ixp);

     nodno = lookup(nodes.ids, idp, "m");
     badno = find(nodno == 0);

     if (! isempty(badno))
        rc = 1;
        for ixb = badno
            fprintf(msg, "*E* loads.%s(%2d) ", typelist{1}, ixp(ixb));
            fprintf(msg, "references undefined node %5d\n", idp(ixb));
        end
     end

     if (! rc)

        [point(1 : nofldc).type]   = deal(0);
        [point(1 : nofldc).nofnod] = deal(0);

        f = sparse(dofs.ndofg, nofldc);

        for n = 1 : nloads(1)

            lc        = inload.point(n).lc;
            inflc(lc) = bitset(inflc(lc), lbits(1));

            np = ++point(lc).nofnod;
            point(lc).id(np) = nodno(n);

            point(lc).d(np, :) = zeros(1, dofs.mxdofpnt);
            fin  = inload.point(n).data;
            lnf  = min(length(fin), dofs.mxdofpnt);
            fvec = fin(1 : lnf);
            point(lc).d(np, 1 : lnf) += fvec;

            dof1 = (nodno(n) - 1) * dofs.mxdofpnt + 1;
            dof2 = dof1 + lnf - 1;
            f(dof1 : dof2, lc) += fvec';

         end 

         point = mfs_resnodedata(point, nodes);

     end

  end

# Process enforced motion

  if (enfmot)
     [motion(1 : nofldc).type]   = deal(0);
     [motion(1 : nofldc).nofnod] = deal(0);
  end

  for kk = 2 : 4

      k = kk - 1;

      if (nloads(kk))

         inp = getfield(inload, typelist{kk});

         [idp, ixp] = sort([inp.id]);
         inp        = inp(ixp);

         nodno = lookup(nodes.ids, idp, "m");
         badno = find(nodno == 0);

         if (! isempty(badno))
            rc = 1;
            for ixb = badno
                fprintf(msg, "*E* loads.%s(%2d) ", typelist{kk}, ixp(ixb));
                fprintf(msg, "references undefined node %5d\n", idp(ixb));
            end
         end

         if (! rc)

            u = sparse(dofs.ndofg, nofldc);

            for n = 1 : nloads(kk)

               lc        = inp(n).lc;
               inflc(lc) = bitset(inflc(lc), lbits(kk));

               if (motion(lc).type)
                  if (motion(lc).type != k)
                     fprintf(msg, "*E* Loadcase %2d contains mixed ", lc)
                     fprintf(msg, "types of enforced motion, ");
                     fprintf(msg, " node = %5d\n", inp(n).id);
                     rc = 1;
                  end
               else
                  motion(lc).type = k;
               end

               np = ++motion(lc).nofnod;
               motion(lc).id(np) = nodno(n);

               motion(lc).d(np, :) = zeros(1, dofs.mxdofpnt);
               uin  = inp(n).data;
               lnf  = min(length(uin), dofs.mxdofpnt);
               uvec = uin(1 : lnf);
               motion(lc).d(np, 1 : lnf) += uvec;

               dof1 = (nodno(n) - 1) * dofs.mxdofpnt + 1;
               dof2 = dof1 + lnf - 1;
               u(dof1 : dof2, lc) += uvec';

               ixp = lookup(dofs.dofp, [dof1 : dof2], "m");
               ixz = find(ixp == 0);
               if ((ln = length(ixz)))
                   for j = ixz
                       if (uvec(j))
                          fprintf(msg, "*W* Node %5d, dof %1d: ",
                                  inp(n).id, j);
                          fprintf(msg, "Enforced motion ignored\n");
                      end
                  end
               end

            end

            uva(k).mat = u;

         end

      end

  end

# Process inertia loads

  if (nloads(5))
     [inertia(1 : nofldc).type] = deal(4);
     ar = zeros(nodes.ncoor, nofldc);
     for n = 1 : nloads(5)
         lc        = inload.inertia(n).lc;
         inflc(lc) = bitset(inflc(lc), lbits(5));
         inertia(lc).d = zeros(1, nodes.ncoor);
         ain  = inload.inertia(n).data;
         lna  = min(length(ain), nodes.ncoor);
         avec = ain(1 : lna);
         inertia(lc).d(1 : lna) += avec;
         ar(:, lc) += avec';
     end
  end

# Build the load structure

  if (! rc)
     load = struct("nofldc", nofldc, "inflc", inflc);
     if (nloads(1))
        load = setfield(load, "point", point);
        load = setfield(load, "f", f);
     end
     if (enfmot)
        load = setfield(load, "motion", motion);
        for k = 1 : 3
           if (nloads(k + 1))
              load = setfield(load, uva(k).name, uva(k).mat);
           end
        end
     end
     if (nloads(5))
        load = setfield(load, "inertia", inertia);
        load = setfield(load, "ar", ar);
     end
  end

end
