function strJ = swmsPrintBlockJStr(nodeID,nodeX,nodeY,nodeZ,nodeType,transLength)


templJ = readTemplate('j');

% calculate node area of seepage nodes:
[boundSeepNodeID,boundSeepNodeArea] = calcNodeArea(nodeID,nodeX,nodeY,nodeZ,nodeType==-2,'xz');

% calculate node area of precipitation nodes:
[boundAtmNodeID,boundAtmNodeArea] = calcNodeArea(nodeID,nodeX,nodeY,nodeZ,nodeType==-4,'xy');

TMP = [boundSeepNodeID,boundSeepNodeArea;...
       boundAtmNodeID,boundAtmNodeArea];
   
if ~isempty(TMP)
    boundAll = sortrows(TMP,1);
    boundAllNodeIDStr = sprintf('%d %d %d %d %d %d %d %d %d %d\n',boundAll(:,1));
    boundAllNodeAreaStr = sprintf('%f %f %f %f %f %f %f %f %f %f\n',boundAll(:,2));
    
    if mod(size(boundAll,1),10)~=0
        boundAllNodeIDStr = [boundAllNodeIDStr,char(10)];
        boundAllNodeAreaStr = [boundAllNodeAreaStr,char(10)];
    end
else
    boundAllNodeIDStr = char(10);
    boundAllNodeAreaStr = char(10);
end



strJ = sprintf(templJ,boundAllNodeIDStr,boundAllNodeAreaStr,transLength);




% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % %                 BEGIN LOCAL FUNCTIONS                     % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %





function [ID,nodeArea] = calcNodeArea(ID,X,Y,Z,T,plane)


[nRows,nCols,nLays] = size(T);

ind = find(T);
[r,c,p] = ind2sub([nRows,nCols,nLays],ind);

if ~isempty(ind)

    switch plane
        case 'xz'
            ur = unique(r);
            if numel(ur)>1
                error('multiple faces')
            end
            if ~ismember(ur,[1,nRows])
                error('face not on edge of domain')
            end
            clear nRows
            clear nCols
            clear nLays

            ID = shiftdim(ID(ur,:,:),1)';
            X = shiftdim(X(ur,:,:),1)';
            Y = shiftdim(Y(ur,:,:),1)';
            Z = shiftdim(Z(ur,:,:),1)';
            T = shiftdim(T(ur,:,:),1)';
            nodeArea = zeros(size(ID));

            faceArea = calcFaceArea(X,Z);

        case 'xy'
            up = unique(p);
            if numel(up)>1
                error('multiple faces')
            end
            if ~ismember(up,[1,nLays])
                error('face not on edge of domain')
            end
            clear nRows
            clear nCols
            clear nLays

            ID = ID(:,:,up);
            X = X(:,:,up);
            Y = Y(:,:,up);
            Z = Z(:,:,up);
            T = T(:,:,up);
            nodeArea = zeros(size(ID));

            faceArea = calcFaceArea(X,Y);

        otherwise 
            error('not implemented')
    end

    % now redistribute the face area to the nodes:

    [nFaceRows,nFaceCols] = size(faceArea);

    for iFaceRow = 1:nFaceRows
        for iFaceCol = 1:nFaceCols

            if T(iFaceRow,iFaceCol)
                nodeArea(iFaceRow,iFaceCol) = nodeArea(iFaceRow,iFaceCol) + 0.25*faceArea(iFaceRow,iFaceCol);
            end
            if T(iFaceRow+1,iFaceCol)
                nodeArea(iFaceRow+1,iFaceCol) = nodeArea(iFaceRow+1,iFaceCol) + 0.25*faceArea(iFaceRow,iFaceCol);
            end
            if T(iFaceRow,iFaceCol+1)
                nodeArea(iFaceRow,iFaceCol+1) = nodeArea(iFaceRow,iFaceCol+1) + 0.25*faceArea(iFaceRow,iFaceCol);
            end
            if T(iFaceRow+1,iFaceCol+1)
                nodeArea(iFaceRow+1,iFaceCol+1) = nodeArea(iFaceRow+1,iFaceCol+1) + 0.25*faceArea(iFaceRow,iFaceCol);
            end

        end
    end
    
    % delete nodes which are not part of the boundary condition    
    
    IxDelete = find(nodeArea==0);
    ID = ID(:);
    ID(IxDelete) = [];

    nodeArea = nodeArea(:);
    nodeArea(IxDelete) = [];

    TMP = sortrows([ID,nodeArea],1);
    ID = TMP(:,1);
    nodeArea = TMP(:,2);
else
    disp([char(39),plane,char(39),' : no face defined'])
    ID = [];
    nodeArea = [];
end



% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %




function faceArea  = calcFaceArea(X,Y)


nFaceRows = size(X,1)-1;
nFaceCols = size(X,2)-1;
faceArea = repmat(NaN,[nFaceRows,nFaceCols]);

for iFaceRow = 1:nFaceRows
    for iFaceCol = 1:nFaceCols
        
        yUL = Y(iFaceRow,iFaceCol);
        yUR = Y(iFaceRow,iFaceCol+1);
        yBR = Y(iFaceRow+1,iFaceCol+1);
        yBL = Y(iFaceRow+1,iFaceCol);

        xUL = X(iFaceRow,iFaceCol);
        xUR = X(iFaceRow,iFaceCol+1);
        xBR = X(iFaceRow+1,iFaceCol+1);
        xBL = X(iFaceRow+1,iFaceCol);
        
        faceArea(iFaceRow,iFaceCol) = parallArea([xUL,yUL;xUR,yUR;xBR,yBR;xBL,yBL]);
        
    end
end



function A = parallArea(M)
% calculate the area of a parallelogram
% coordinates are [xUL,zUL;xUR,zUR;xBR,zBR;xBL,zBL]

% hFig = figure;
% plot(M(:,1),M(:,2))

K = sortrows(M,2);
dz = K(4,2)-K(1,2);
dx = M(2,1)-M(1,1);

a = abs(M(2,2)-M(1,2)); %zUR-zUL
b = abs(M(3,2)-M(4,2)); %zBR-zBL

A = dx*dz-...
    0.5*(a*dx)-...
    0.5*(b*dx);

% close(hFig)






