function npn_array_calc(i,i7)
% CALCULATE NPN ARRAY with information about unique nodes
global join
global lnp
global npn nplc xlc ylc
global lnr
global ndiv
global inp
global meltyp

jbgn=1;
jend=i7; % number of nodes per loop

npn=zeros(1,i7);

for j=1:4

    if(join(i,j,1)~=0)  % if this side of the loop connects to another loop

        j1=join(i,j,1); % get the subdomain
        j2=join(i,j,2); % get the side of the other SD

        kbgn=1;

        switch(meltyp)
            case 3
                kend = ndiv(i,j)+1; % kend = number of nodes on this boundary
            case 4
                kend = ndiv(i,j)+1;
            case 6
                kend = 2*ndiv(i,j)+1;
            case 9
                kend = 2*ndiv(i,j)+1;
        end

        k2=kend+1;

        for k=kbgn:kend
            k1=lnr(i,j,k); % k1 is the local node number
            k2=k2-1;
            npn(k1)=lnp(j1,j2,k2);
        end

    end % end if joined
end % end loop over sides

% check to see if loop corners have appeared before. if so, use the already
% existing global node number. if not, assign a new global node number.
corn2npn=[1,ndiv(i,1)+1,i7,i7-ndiv(i,1)];
for jlp=1:(i-1)
    for corn=1:4
        for jcorn=1:4
            if (xlc(i,corn)==xlc(jlp,jcorn))&(ylc(i,corn)==ylc(jlp,jcorn))
                npn(corn2npn(corn))=nplc(jlp,jcorn);
                break
            end
        end
    end
end

jbgn=1;
jend=i7;

for j=jbgn:jend

    if(npn(j)==0)
        inp=inp+1; % inp is zeroed before first SD, then it increments (here) 
                   % for each node that isn't joined to another subdomain
        npn(j)=inp;
    end
    % by the way, for the first loop, npn is just a vector from 1 to the number of
    % nodes in the loop
end
% update nplc with correct global node numbers
nplc(i,:)=npn(corn2npn);