diff --git a/5. Data Assimilation/Utilities/Ensrf/jointENSRF.m b/5. Data Assimilation/Utilities/Ensrf/jointENSRF.m index 95ba6830..c833dd99 100644 --- a/5. Data Assimilation/Utilities/Ensrf/jointENSRF.m +++ b/5. Data Assimilation/Utilities/Ensrf/jointENSRF.m @@ -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 + [K, Kdenom] = jointKalman( 'K', Knum(:,obs), Ydev(obs,:), yloc(obs,obs), Rt(obs) ); %#ok % Update the mean - Amean(:,t) = Mmean + K * ( D(obs,t) - Ymean(obs) ); %#ok + Amean(:,t) = Mmean + K * ( Dt(obs) - Ymean(obs) ); %#ok K = []; %#ok (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 \ No newline at end of file diff --git a/5. Data Assimilation/Utilities/Ensrf/serialENSRF.m b/5. Data Assimilation/Utilities/Ensrf/serialENSRF.m index c6dbb9f6..954d0f68 100644 --- a/5. Data Assimilation/Utilities/Ensrf/serialENSRF.m +++ b/5. Data Assimilation/Utilities/Ensrf/serialENSRF.m @@ -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. @@ -62,7 +62,7 @@ Mpsm = Am(F{d}.H) + Ad(F{d}.H,:); %#ok % 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 + [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 @@ -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 + [K, a] = serialKalman( Ad, Ydev, w(:,d), 1, R(d,t) ); %#ok % Update Am = Am + K*( D(d,t) - Ymean );