function [mStatesKF,mStatesNOKF]=swmsRunFun(sodaPar,mConstants,mStatesKF,mStatesNOKF,mPars,mTime)

sodaUnpack

% presHead1Old = presHead1;
% 
% if ismember(34,get(0,'Children'))
%     close(34)
% end
% subplotScreen(1,1,1,'figureNumber',34)
% clf
% set(gcf,'inverthardcopy','off','paperpositionmode','auto')
% subplot(1,3,1)
% imagesc(presHead1Old)
% caxis([-1,1]*max(abs(presHead1(:))))
% axis image
% hColorbar=colorbar;
% set(get(hColorbar,'title'),'string',num2str(max(presHead1Old(:))))
% set(get(hColorbar,'xlabel'),'string',num2str(min(presHead1Old(:))))
% set(gca,'ydir','normal')
% colormap([autumn;flipud(winter)])
% title(['presHead1Old: Time = ',sprintf('%7.2f',mTime{1,2}(1)), ' ',swmsTimeUnits])
% 


% rSink = 0.01
% Ks = 0.35 % 10 times smaller than before
% disp('hardcoded parameters')

rSink = 10^LOG10rSink
Ks = 10^LOG10Ks

saveFigs = false;
visualizeTransects = false;
viz = false;

% state updates have been suggested by the EnKF. These have been applied to
% variable presHead1. First, determine the sign and magnitude of the
% update:
dPresHead1 = presHead1-presHead1Old;
% then, apply the update to the rest of the profile also:
presHead52New = presHead52+repmat(dPresHead1,[1,1,size(presHead52,3)]);


% delete all old files in ./SWMS_3D.IN and ./SWMS_3D.OUT:
% swmsCleanup('IMG','OUT','MAT')
OUT = writeSWMSIN(rSink,Ks,mTime,presHead52New,presHead1);
for iVar = 2:numel(OUT)
    eval([OUT{1}{iVar-1},' = OUT{iVar};'])
end
clear iVar
clear OUT

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% visualize the layout of the grid:
vizScriptBoundaries
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% visualize the materials properties:
vizScriptMaterials
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% visualize the Feddes root extraction/ bedrock sink parameters
vizScriptFeddes

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

% call the SWMS binary:
showSTDOUT = true;
[status,result] = swmsCallBat('12',showSTDOUT);
clear status
clear showSTDOUT
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

if isempty(strfind(result,'error')) &&...
        isempty(strfind(result,'*ERR*')) &&...
        isempty(strfind(result,'is exceeded'))
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
    vizScriptTransectsHeads
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
    vizScriptTransectsThetas
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
    vizScriptFlow
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
    vizScriptInOut
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
    vizScriptCondensedFlow
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
    vizScriptMesh
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
    vizScriptSoilDepth
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
    vizScriptGroundWater
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
    vizScriptRootExtractionPresHead
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
    vizScriptRootExtractionTheta
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
else
    disp('EXE exited with error --no output files read.')
    return
end % if isempty(findstr(result,'error'))



% save([datestr(now,30),'-results.mat'])


presHeads = readOutput(pwd,'h.out',false);
presHeads4D = conv2dto4d(presHeads.values,nRows,nCols,nLays,nodes);

[B,IX] = sort([nodes.ID(:)]);
presHeads3D = repmat(NaN,[nRows,nCols,nLays]);
presHeads3D(IX) = presHeads.values(:,end);

presHead52 = presHeads3D(:,:,2:5);
presHead1 = presHeads3D(:,:,1);

% store the prior state values in variable presHead1Old, to be able to
% calculate the sign and magnitude of the state update when the program
% returns to this function.
presHead1Old = presHead1;

[vInit,vEnd] = readBALANCEOUT();
[vPrec,vSink,vSeep] = readQOUT();

discharge = -vSeep/(mTime{1,2}(2)-mTime{1,2}(1));

sodaPack

