function OUT = writeSWMSIN(parRSink,parKs,mTime,presHead52,presHead1)

includeWaterFlow = true;
includeRain = true;
includeRoot = true;
includeSeep = true;

includeAtmos = includeRain|includeRoot;

verboseMode = false;
uniformMaterial = true;
writeToFile = true;

hydrusFileName = fullfile(pwd,'MESHTRIA.TXT');
swmsFileName = fullfile(pwd,filesep,'SWMS_3D.IN',filesep,'GRID.IN');
% First, read in data from the mesh file generated by HYDRUS based on the
% hydrusFileName  = 'coord ground and bedrock-inclined base.txt' file.
[NumNP,NumEl] = textread(hydrusFileName,'%d%d',1,'headerlines',5);
[nodeID,nodeX,nodeY,nodeZ]=textread(hydrusFileName,'%d%f%f%f',NumNP,'headerlines',7);
[elemID,elemI,elemJ,elemK,elemL,elemM,elemN,elemO,elemP,elemSub]=...
    textread(hydrusFileName,'%d%d%d%d%d%d%d%d%d%d',NumEl,'headerlines',NumNP+9);
xSet = 9:3:27;
ySet = 1:3:43;
IO = ismember(nodeX,xSet) &...
     ismember(nodeY,ySet);
tmp = [nodeX(IO),nodeY(IO),nodeZ(IO)];

nodeX = tmp(:,1);
nodeY = tmp(:,2);
nodeZ = tmp(:,3);

clear tmp

nodeID = [1:sum(IO)]';
xv = unique(nodeX)';
yv = unique(nodeY);
nRows = numel(yv);
nCols = numel(xv);
laySel = logical([0,0,0,0,0,1,1,1,1,1]);
nLays = sum(laySel);

X = repmat(xv,[nRows,1,nLays]);
Y = repmat(yv,[1,nCols,nLays]);
Z = repmat(NaN,[nRows,nCols,nLays]);
ID = repmat(NaN,[nRows,nCols,nLays]);
Layer = repmat(shiftdim(1:nLays,-1),[nRows,nCols,1]);

if uniformMaterial
    Material(1:nRows,1:nCols,1:nLays) = 1;
else
    Material(1:nRows,1:nCols,[nLays:-1:1]) = repmat(shiftdim([1,1,2,2,3],-1),[nRows,nCols,1]);
end

iID = 0;

for iRow = 1:nRows
    for iCol = 1:nCols        
        tmp = sort(nodeZ(nodeX==X(iRow,iCol)&nodeY==Y(iRow,iCol)),'ascend');
        Z(iRow,iCol,1:nLays) = tmp(laySel);
        soilDepthIncrease = 0.15;%(tmp(10)-tmp(6)-0.30)*-0.7;
        Z(iRow,iCol,1:nLays) = linspace(tmp(6),tmp(10)+soilDepthIncrease,nLays)';
        iID  = iID + 1;
        ID(iRow,iCol,nLays:-1:1) = (iID-1)*nLays+1:iID*nLays;
    end
end
clear tmp
soilDepth = Z(:,:,nLays)-Z(:,:,1);

% % save('Z.mat','Z')
% presHeadIni = repmat(Z(:,:,1),[1,1,5])-Z
% save('preshead-0-at-lower-bound.mat','presHeadIni')
% 
% figure
% hist(soilDepth(:),20)
% save('soildepth.mat','soilDepth')


% select which nodes should be part of the seepage face:
nodeType = zeros(nRows,nCols,nLays);
if includeSeep
    nodeType = nodeType+((nodeType==0)&...
                         (Y==min(Y(:))))*-2;
end

if includeRain
    nodeType = nodeType+(nodeType==0&...
               ismember(Y,[4:3:46])&...
                        Layer==nLays)*-4;
end

leakageMap = false;
D = zeros(nRows+1,nCols+1);
D(2:end-1,2:end-1) = 0.25*ones(nRows-1,nCols-1);
relativeResolution = D(1:end-1,1:end-1)+D(2:end,2:end)+D(1:end-1,2:end)+D(2:end,1:end-1);
clear D
if includeRoot
    leakageMap = ismember(Y(:,:,1),[4:3:46]);
    applyRootExtrHere = (nodeType==0)&...
                (Layer==1)&...
                repmat(leakageMap,[1,1,nLays]);

    % make extraction rates in certains nodes stronger:
%     extraLeakageNodes = [4,7,1;...
%                          6,4,1;...
%                          7,1,1;...
%                         11,3,1;...
%                         14,6,1]
    extraLeakageNodes = []
    
    TMP = repmat(relativeResolution,[1,1,nLays]).*applyRootExtrHere;
    betaTmp = TMP./sum(TMP(:));
    
    nodeType = nodeType + applyRootExtrHere.*0;
    nodeBeta = applyRootExtrHere.*betaTmp+...
               ~applyRootExtrHere*0.0;
           
    for iRecord = 1:size(extraLeakageNodes,1)
        iRow = extraLeakageNodes(iRecord,1);
        iCol = extraLeakageNodes(iRecord,2);
        iLay = extraLeakageNodes(iRecord,3);
        nodeBeta(iRow,iCol,iLay) = nodeBeta(iRow,iCol,iLay)*30;
    end           
    clear TMP
else
    nodeBeta = zeros(nRows,nCols,nLays);
end

% presHeadIni = repmat(Z(:,:,1),[1,1,nLays])-Z;
% save('preshead-0-at-lower-bound.mat','presHeadIni')

presHeadIni = cat(3,presHead1,presHead52);

M = sortrows([ID(:),X(:),Y(:),Z(:),nodeType(:),presHeadIni(:),Material(:),Layer(:),nodeBeta(:)],1);

clear elemID
clear elemSub

dRow = [0,0,1,1,0,0,1,1];
dCol = [0,1,0,1,0,1,0,1];
dLay = [1,1,1,1,0,0,0,0];

iElem = 0;
for iCol = 1:nCols-1
    for iRow = 1:nRows-1
        for iLay = 1:nLays-1
            iElem = iElem+1;
            elemID(iElem,1)=iElem;
            for iCorner = 1:8
                r = iRow+dRow(iCorner);
                c = iCol+dCol(iCorner);
                l = iLay+dLay(iCorner);

                hexas(iElem,iCorner)=ID(r,c,l);
            end
            test = (isodd(iCol)&isodd(iRow))|(~isodd(iCol)&~isodd(iRow));
            if test
                elemSub(iElem,1) = 1;
            else
                elemSub(iElem,1) = 2;
            end
        end
    end
end

nHexas=size(hexas,1);
tetras = repmat(NaN,[nHexas*6,4]);
for iHexa=1:nHexas
    tetras((iHexa-1)*6+1:iHexa*6,1:4) = hexas2tetras(hexas(iHexa,:),elemSub(iHexa,1),'tetraMode1');
end

NumBP = sum(nodeType(:)~=0);

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % % % % % % %             BLOCK H           % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

% beta
blockHStr = swmsPrintBlockHStr(numel(ID),size(tetras,1),NumBP,M(:,1),M(:,2),M(:,3),M(:,4),M(:,5),M(:,6),M(:,7),M(:,8),M(:,9));

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % % % % % % %             BLOCK I           % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

blockIStr = swmsPrintBlockIStr(tetras);

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % % % % % % %             BLOCK J           % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

totalArea = (max(X(:))-min(X(:)))*(max(Y(:))-min(Y(:)));
if includeRoot
    transLength = sum(sum(applyRootExtrHere(:,:,1).*...
        (relativeResolution./sum(relativeResolution(:)))*totalArea));
else
    transLength  = 0;
end

blockJStr = swmsPrintBlockJStr(ID,X,Y,Z,nodeType,transLength);
eoGridInStr = '***** End of file Grid.In ***********************************************************************************';

% monitor
monitor.x = 15;%X(1,1,1);
monitor.y = 19;%Y(1,1,1);
monitor.layer = 1;
monitor.z = Z(X==monitor.x & Y==monitor.y & Layer==monitor.layer);


% transect
ind = find(X==monitor.x);
[a,Ix,b] = ind2sub(size(X),ind(1));
transect.x = X(1,Ix,1);
transect.y = Y(:,Ix,1);
transect.zBottom = Z(:,Ix,1);
transect.zSurf = Z(:,Ix,nLays);
clear a b Ix ind


nodes.n = numel(ID);
nodes.ID = ID;
nodes.X = X;
nodes.Y = Y;
nodes.Z = Z;
nodes.Layer = Layer;
nodes.type = nodeType;
nodes.beta = nodeBeta;
nodes.presHeadIni = presHeadIni;
nodes.material = Material;
nodes.description.n = 'number of nodes in the grid';
nodes.description.ID = 'unique identifier for each node in the grid';
nodes.description.X = ['distance in SWMS/HYDRUS length units from the left ',char(10),...
                       'side of the grid; first row is lowest part of the hillslope'];
nodes.description.Y = ['distance in SWMS/HYDRUS length units from the bottom',char(10),...
                       'of the grid; first row is lowest part of the hillslope'];
nodes.description.Z = 'elevation in SWMS/HYDRUS length units, relative to a horizontal plane';
nodes.description.type = 'type of node/type of boundary condition.';
nodes.description.material = 'identifier for the material (soil physical properties). See  block B in SELECTOR.IN.';


elems.n = size(tetras,1);
elems.tetras = tetras;
elems.description.n = ['number of elements in the grid. Each element is a',char(10),...
                        'tetrahedron consisting of 4 points.'];
elems.description.tetras = ['Nx4 array of node identifiers; each row in ',char(10),...
                            'the array contains the 4 corners of a tetrahedron.'];
                        

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % %                 END OF GRID.IN PART                       % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 


% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % %                 START OF SELECTOR.IN PART                 % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 


% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % % % % % % %             BLOCK A           % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
headingLine = 'test with Luisas data --automated writing of block A,B,C,D,E,H,I,J,K';
lengthUnits = 'm';
timeUnits = 'hr';
maxIterations = 31;
tolTheta = 0.000001;
tolHead = 0.000001;

strA = swmsPrintBlockAStr(headingLine,lengthUnits,timeUnits,...
    maxIterations,tolTheta,tolHead,includeWaterFlow,includeAtmos,includeSeep);

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % % % % % % %             BLOCK B           % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

%  excerpt from luisa's data
%   thr    ths   Alfa     n         Ks      l
%   0.28  0.475      4      2        3.5    0.5 
%   0.28   0.46      4      2        1.5    0.5 
%  0.325   0.45      4      2       0.65    0.5 
%    0.3   0.45   3.25   1.75      0.006    0.5 
%   0.28    0.4      3    1.5     0.0006    0.5 
% 

if uniformMaterial
    materials.captions = {'thr','ths','Alfa','n','Ks'};
    materials.values = [0.28,0.475,4,2,parKs];
else
    materials.captions = {'thr','ths','Alfa','n','Ks'};
    materials.values = [0.28,0.475,4,2,3.5;...
                        0.28,0.46,4,2,1.5;...
                        0.325,0.45,4,2,0.65];
%                     ;...
%                         0.3,0.45,3.25,1.75,0.006;...
%                         0.28,0.4,3,1.5,0.0006];
end
nMaterials = size(materials.values,1);
strB = swmsPrintBlockBStr(materials);


% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % % % % % % %             BLOCK C           % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% integration
timeDeltaMin = 1/60;
timeDelta = 6/60;
timeDeltaMax = 20/60;
multiplierHigher = 1.2;
multiplierLower = 0.7;

%hydrusPrec = readHydrusRain();
% rainTable format =
% rainTable = [  a,b;...
%                c,d;...
%                e,f];   "From a through a+c, use intensity d"

rainTable = [  0,0.000;...
              96,0.000;...
               6,0.006;...
             500,0.000];
         
rainTable = [rainTable(:,1),cumsum(rainTable(:,1)),rainTable(:,2)];

tmp.t = [rainTable(1,2):0.25:rainTable(end,2)]';
tmp.precRate = repmat(0,[numel(tmp.t),1]);
for iRecord = 2:size(rainTable,1)
    rainNow = rainTable(iRecord-1,2)<tmp.t & tmp.t<=rainTable(iRecord,2);
    tmp.precRate(rainNow) = rainTable(iRecord,3);
end

% figure(36)
% plot(tmp.t,tmp.precRate)
% ylabel('prec [m/h]')

% [C,IA,IB] = intersect(tmp.t,mTime{1,2}(2));
% prec.time = tmp.time(IA);
% prec.rate = tmp.rate(IA);
IO = (mTime{1,2}(1)<tmp.t)&(tmp.t<=mTime{1,2}(end));
prec.time = tmp.t(IO);
prec.rate = tmp.precRate(IO);
clear tmp

printTimes = mTime{1,2}(2:end);

strC = swmsPrintBlockCStr(timeDelta,timeDeltaMin,timeDeltaMax,...
    multiplierHigher,multiplierLower,printTimes);

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % % % % % % %             BLOCK D           % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

if includeRoot
    % feddes
    P0 = 11;
    P2H = 1.00;%0.50;
    P2L = P2H;
    P3 = 0.0;
    POptm = ones(1,nMaterials)*P0;
    strD = swmsPrintBlockDStr(P3,P2L,P2H,POptm,P0);
    
    feddes.P0 = P0;
    feddes.POptm = POptm;
    feddes.P2H = P2H;
    feddes.P2L = P2L;
    feddes.P3 = P3;
    clear P0 POptm P2L P2H P3
else
    feddes = [];
    strD = '';
end


% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % % % % % % %             BLOCK E           % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

if includeSeep
    strE = swmsPrintBlockEStr(nodes);
else
    strE = '';
end

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % %                 SELECTOR.IN END-OF-FILE                   % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

eoSELECTORINStr = '*** END OF INPUT FILE ''SELECTOR.IN'' ************************************';

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % %                 END OF SELECTOR.IN PART                   % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 


% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % %                 START OF ATMOSPH.IN PART                  % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
if includeRain|includeRoot
    templK = readTemplate('k');
    if includeRoot
        sinkF = 't';
    else
        sinkF = 'f';
    end
    qGWLF = 'f';

    if strcmp(qGWLF,'f')
        Aqh = 0;
        Bqh = 0;
    end

    hCritS = 0;
    
    rPrec = prec.rate; % swmsLengthUnits/swmsTimeUnits
    % Prec = zeros(numel(printTimes),1);
   
    rSink = zeros(numel(prec.time),1);
    
    if includeRoot
        % rSink(:) = 0.0065888889;
        rSink(:) = parRSink;
    end

    tInit = mTime{1,2}(1);

    MaxAl = numel(rPrec);
    tAtm = prec.time;

    dataStr = sprintf('    %f       %.9f       0       0    %f  100     0     0     0     0 \n',[tAtm,rPrec,rSink]');
    
    blockKStr = sprintf(templK,sinkF,qGWLF,Aqh,Bqh,tInit,MaxAl,hCritS,dataStr);
    
    eoATMOSPHStr = '*** END OF INPUT FILE ''ATMOSPH.IN'' *************************************';
    
else
    blockKStr = '';
    tAtm = [];
    rSink = [];
    rPrec = [];
end




% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % %                 END OF ATMOSPH.IN PART                    % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 





% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % %                     WRITE TO FILE                       % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
if writeToFile
    fid=fopen(fullfile(pwd,'SWMS_3D.IN','GRID.IN'),'wt');
    fprintf(fid,'%s%s%s%s',blockHStr,blockIStr,blockJStr,eoGridInStr);
    fclose(fid);

    fid = fopen(fullfile(pwd,'SWMS_3D.IN','SELECTOR.IN'),'wt');
    fprintf(fid,'%s%s%s%s%s',strA,strB,strC,strD,strE,eoSELECTORINStr);
    fclose(fid);

    if includeRain|includeRoot
        fid = fopen(fullfile(pwd,'SWMS_3D.IN','ATMOSPH.IN'),'wt');
        fprintf(fid,'%s%s',blockKStr,eoATMOSPHStr);
        fclose(fid);
    end
    if verboseMode
        disp([char(10),'Writing of input files completed.',char(10)])
    end
else
    if verboseMode
        disp([char(10),'Writing of input files skipped.',char(10)])
    end
end




OUT{1} = {'nodes','elems','leakageMap','soilDepth','relativeResolution',...
                'atmosphTable','nRows','nCols','nLays','totalArea','transLength',...
                'materials','feddes','includeRain','includeSeep','includeRoot',...
                'transect','monitor'};
OUT{2} = nodes;
OUT{3} = elems;
OUT{4} = leakageMap;
OUT{5} = soilDepth;
OUT{6} = relativeResolution;
OUT{7} = [tAtm,rPrec,rSink];
OUT{8} = nRows;
OUT{9} = nCols;
OUT{10} = nLays;
OUT{11} = totalArea;
OUT{12} = transLength;
OUT{13} = materials;
OUT{14} = feddes;
OUT{15} = includeRain;
OUT{16} = includeSeep;
OUT{17} = includeRoot;
OUT{18} = transect;
OUT{19} = monitor;



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


function tetras = hexas2tetras(hexas,elemSub,biasMode)

switch biasMode
    case 'tetraMode1'
        nElems = size(hexas,1);
        for iElem=1:nElems
            clear tmp
            if elemSub(iElem,1)==1
                nodeOrder = [1,4,7,3;...
                             1,7,4,8;...
                             1,5,7,8;...
                             1,2,6,4;...
                             1,6,5,8;...
                             1,8,4,6]; 
                for k=1:size(nodeOrder,1)
                    tetras(k,1:4) = hexas(iElem,nodeOrder(k,:));
                end
            elseif elemSub(iElem,1)==2
                nodeOrder = [3,8,4,5;...
                             3,5,7,8;...
                             1,3,4,5;...
                             1,4,2,5;...
                             5,6,8,2;...
                             5,2,8,4];
                for k=1:size(nodeOrder,1)
                    tetras(k,1:4) = hexas(iElem,nodeOrder(k,:));
                end
            else
            end
        end
    case 'tetraMode2'
        % striped bias, works only on grids with one layer.
        nElems = size(hexas,1);
        for iElem=1:nElems
            clear tmp
            if elemSub(iElem,1)==1
                nodeOrder = [7,2,4,3;...
                             7,6,8,4;...
                             6,2,4,7;... %3a
                             6,7,5,2;...
                             2,7,5,1;...
                             2,7,1,3];   %3c
                for k=1:size(nodeOrder,1)
                    tetras(k,1:4) = hexas(iElem,nodeOrder(k,:));
                end
            elseif elemSub(iElem,1)==2
                nodeOrder = [8,3,1,4;...
                             5,7,1,8;...
                             5,7,3,8;...
                             6,2,8,1;...
                             1,6,5,8;...
                             1,2,8,4]; % manually mirrored version of geometry 
                                       % created under "elemSub(iElem,1)==1"
                for k=1:size(nodeOrder,1)
                    tetras(k,1:4) = hexas(iElem,nodeOrder(k,:));
                end
            else
            end
        end        
        
        
end

