Skip to content

Commit

Permalink
Corrected sliced variables
Browse files Browse the repository at this point in the history
  • Loading branch information
JonKing93 committed Apr 30, 2019
1 parent 0b3d429 commit 9d46624
Show file tree
Hide file tree
Showing 2 changed files with 12 additions and 8 deletions.
14 changes: 9 additions & 5 deletions 5. Data Assimilation/Utilities/Ensrf/jointENSRF.m
Original file line number Diff line number Diff line change
Expand Up @@ -42,26 +42,30 @@
% For each time step
parfor t = 1:nTime

% Slice variables
Rt = R(:,t);
Dt = D(:,t);

% Update using obs that are not NaN
update(:,t) = update(:,t) & ~isnan( D(:,t) );
update(:,t) = update(:,t) & ~isnan( Dt );
obs = update(:,t);

% Get the full kalman gain and kalman denominator
[K, Kdenom] = jointKalman( 'K', Knum(:,obs), Ydev(obs,:), yloc(obs,obs), R(obs,t) ); %#ok<PFBNS>
[K, Kdenom] = jointKalman( 'K', Knum(:,obs), Ydev(obs,:), yloc(obs,obs), Rt(obs) ); %#ok<PFBNS>

% Update the mean
Amean(:,t) = Mmean + K * ( D(obs,t) - Ymean(obs) ); %#ok<PFBNS>
Amean(:,t) = Mmean + K * ( Dt(obs) - Ymean(obs) ); %#ok<PFBNS>
K = []; %#ok<NASGU> (Free up memory)

% If returning variance
if ~meanOnly

% Get the adjusted kalman gain
Ka = jointKalman( 'Ka', Knum(:,obs), Kdenom, R(obs,t) );
Ka = jointKalman( 'Ka', Knum(:,obs), Kdenom, Rt(obs) );

% Update the deviations and get the variance
Avar(:,t) = var( Mdev - Ka * Ydev(obs,:), 0, 2 );
end
end

end
end
6 changes: 3 additions & 3 deletions 5. Data Assimilation/Utilities/Ensrf/serialENSRF.m
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
function[Amean, Avar, Ye, update] = serialENSRF( M, D, R, F, w )
function[Amean, Avar, Ye, R, update] = serialENSRF( M, D, R, F, w )
%% Implements data assimilation using an ensemble square root filter with serial
% updates for individual observations. Time steps are assumed independent
% and processed in parallel.
Expand Down Expand Up @@ -62,7 +62,7 @@
Mpsm = Am(F{d}.H) + Ad(F{d}.H,:); %#ok<PFBNS>

% Run the PSM. Generate R. Error check.
[Ye(d,:,t), R(d,t), update(d,t)] = getPSMOutput( F{d}, Mpsm, R(d,t), t, d ); %#ok<PFOUS>
[Ye(d,:,t), R(d,t), update(d,t)] = getPSMOutput( F{d}, Mpsm, R(d,t), t, d );

% If no errors occured in the PSM, and the R value is valid,
% update the analysis
Expand All @@ -72,7 +72,7 @@
[Ymean, Ydev] = decomposeEnsemble( Ye(d,:,t) );

% Get the Kalman gain and alpha scaling factor
[K, a] = serialKalman( Ad, Ydev, w(:,d), 1, R(d,t) ); %#ok<PFBNS>
[K, a] = serialKalman( Ad, Ydev, w(:,d), 1, R(d,t) ); %#ok<PFBNS>

% Update
Am = Am + K*( D(d,t) - Ymean );
Expand Down

0 comments on commit 9d46624

Please sign in to comment.