function [OUT,sodaPar]=sodaEnKF(sodaPar,parVec)

mConstants = evalin('base','mConstants');
mStatesExcluded = evalin('base','mStatesExcluded');

for k=1:numel(parVec)
    mPars{k,1} = sodaPar.parNames{k};
    mPars{k,2} = parVec(k);
end
clear k

measOperator = 1;

nIter = numel(sodaPar.iterator{1,2});
nStates = sodaPar.nStates;
nMembers = sodaPar.nEnsembleMembers;
nOutputs = sodaPar.nOutputs;
nOptPars = sodaPar.nOptPars;


TMP=strcmp(sodaPar.stochForcePars,'optimized');
stochForceParsFlag = any(TMP(:,2));

if stochForceParsFlag
    Ix = find(TMP(:,2));
    sodaPar.stochForcePars{Ix,2} = parVec(sodaPar.stochForceParsIx);
end
clear TMP

mOutputs = repmat(NaN,[nIter,nOutputs,nMembers]);

defaultNaNArray = repmat(NaN,[nIter,nStates,nMembers]);

mStatesPrior = defaultNaNArray;
mStatesPert = defaultNaNArray;
mStatesPost = defaultNaNArray;

perturbedObs = defaultNaNArray;
stateInn = defaultNaNArray;


iIter = 1;
stochForce = defaultNaNArray;
eval(sodaPar.errModelCallStr)


kalmanGain = repmat(NaN,[nIter,nStates,nStates]);

switch sodaPar.initState
    case 'unifrand'
        mStatesPost(1,1:nStates,1:nMembers) = sodaUnifRandDraw(sodaPar,'stateSpace');
    case 'reference'
        for iMember=1:nMembers
            mStatesPost(1,1:nStates,iMember) =...
            sodaPar.initStateRef + randn(1,nStates).*sodaPar.initStateErr;
        
            statesTooLow = sodaPar.stateSpaceLoBound > mStatesPost(1,1:nStates,iMember);
            mStatesPost(1,statesTooLow,iMember) = sodaPar.stateSpaceLoBound(statesTooLow);
            
            statesTooHigh = sodaPar.stateSpaceHiBound < mStatesPost(1,1:nStates,iMember);
            mStatesPost(1,statesTooHigh,iMember) = sodaPar.stateSpaceHiBound(statesTooHigh);
        end
end

for iIter = 2:nIter
    
    % Propagate all ensemble members one time step...
    % further using the specified model: 
    [mStatesPrior,mStatesExcluded,mOutputs]=sodaPropStates(sodaPar,iIter,...
       mConstants,mStatesPrior,mStatesPost,mPars,mStatesExcluded,mOutputs);


    % determine the stochastic forcing value:
    eval(sodaPar.errModelCallStr)
       
    % Add a model error to the simulated state values:
    mStatesPert = sodaPerturbModel(sodaPar,mStatesPrior,stochForce,mStatesPert,iIter);
    
    % Compute the ensemble covariance matrix, Pe (See...
    % Equation 47 of Evensen et al. 2002):
    ensembleCov = sodaCalcEnsembleCov(mStatesPert,iIter);

    % For all ensemble members, perturb the current observation...
    % using a randn distribution:
    [perturbedObs,sodaPar] = sodaPerturbObs(sodaPar,mConstants,perturbedObs,iIter);
    
    % Compute the measurement error covariance matrix, Re (See...
    % Equation 51 of Evensen et al. 2002):
    measErrCov = sodaCalcMeasErrCov(sodaPar,perturbedObs,iIter);


    % Caclulate the classical Kalman Gain (See...
    % Equation 52 of Evensen et al. 2002):
    kalmanGain = sodaCalcKalmanGain(ensembleCov,measOperator,measErrCov,kalmanGain,iIter);
    

    % Calculate the state innovations:
    stateInn = sodaCalcInnovations(sodaPar,kalmanGain,perturbedObs,...
        measOperator,mStatesPert,stateInn,iIter);
    
    % Add innovation to ensemble of states (See Equation 54 of ...
    % Evensen et al. 2002):
    mStatesPost = sodaUpdateStates(sodaPar,mStatesPert,...
        mStatesPost,stateInn,iIter);
    
    
end

meanEnPredPrior = measOperator*mean(mStatesPrior,3);
meanEnPredPost = measOperator*mean(mStatesPost,3);
OUT = measOperator*mean(mStatesPert,3);


if sodaPar.plotEnsemble

    trueObs = mConstants{sodaPar.isMeas,2};

    figure(873)
    clf
    set(gcf,'numbertitle','off','name','ensemble results')

    for iStateTMP = 1:nStates+1
        
        if iStateTMP>nStates
            iState = iStateTMP-1;
        else
            iState = iStateTMP;
        end

        if nStates>1
            subplot(ceil(sqrt(nStates+1)),ceil(sqrt(nStates+1)),iStateTMP)
        else
            subplot(1,2,iStateTMP)
        end
        
        for iMember=1:nMembers
            xOffset = [1:nIter]-0.125+(((iMember-1)/(nMembers-1))*0.25);
            plotHandles(:,iMember) = plot(xOffset,mStatesPrior(1:nIter,iState,iMember),'o',...
                xOffset,mStatesPert(1:nIter,iState,iMember),'+',...
                xOffset,mStatesPost(1:nIter,iState,iMember),'s',...
                xOffset,perturbedObs(1:nIter,iState,iMember),'d');
            hold on
        end
        plotHandles2 =plot(1:nIter,trueObs(:,iState),'k*',...
             1:nIter,meanEnPredPrior(:,iState),...
             1:nIter,meanEnPredPost(:,iState));

        set(plotHandles,'markersize',5)
        colorArray={'b';'r';[0,0.5,0];'m'};
        for k=1:4
            set(plotHandles(k,:),'markerfacecolor',colorArray{k},'markeredgecolor',colorArray{k})
        end


        set(plotHandles2(2),'color',[0,0,1])
        set(plotHandles2(3),'color',[0,0.5,0])
        set(gca,'ylim',[sodaPar.stateSpaceLoBound(iState),...
            sodaPar.stateSpaceHiBound(iState)],'xlim',[1,nIter])

        xlabel(sodaPar.iterator{1})
        ylabel(sodaPar.stateNames{iState})
        
        if iStateTMP==nStates+1
            set(gca,'xlim',min(get(gca,'xlim'))+[-2,-1])
%             hLegend = legend([plotHandles(1:4,1);plotHandles2],...
%                 '\eta(\Psi_e^a,X,\theta)',...
%                 '\Psi_e^f',...
%                 '\Psi_e^a',...
%                 'H(\Psi_t)+\epsilon_t',...
%                 'true obs',...
%                 'mean(\Psi_e^f)',...
%                 'mean(\Psi_e^a)');
%             set(hLegend,'interpreter','tex');
            hLegend = legend([plotHandles(1:4,1);plotHandles2],...
                'mStatesPrior',...
                'mStatesPert',...
                'mStatesPost',...
                'perturbedObs',...
                'trueObs',...
                'mean(mStatesPrior)',...
                'mean(mStatesPost)');
            set(hLegend,'interpreter','none');

            axis off
            
        end
    end

    drawnow
end

% (perturbedObs<mStatesPost & mStatesPost<mStatesPert) |...
% (mStatesPert<mStatesPost & mStatesPost<perturbedObs)


if stochForceParsFlag
    
    sodaPar.stochForcePars{Ix,2}='optimized';
    
end
