% AOSC 614 - Intro to Data Assimilation

% MATLAB Examples by Steven Greybush

% to complement lectures by Eugenia Kalnay



% =========================================================================

% Toy Example 1: 2 thermometers 

clear all



Tt = 273 % truth [K]



T1 = 271 % thermometer 1 (forecast) [K]

T2 = 274 % thermometer 2 (observation) [K]

sigma1 = 2^0.5 % error standard deviation of thermometer 1 (fcst) [K]

sigma2 = 1^0.5 % error standard deviation of thermometer 2 (obs) [K]



a1 = sigma2^2 / (sigma1^2 + sigma2^2) % first weight

a2 = 1 - a1 % second weight

Ta = a1*T1+a2*T2 % temperature analysis [K]

sigma_a = (sigma1^2*sigma2^2)/(sigma1^2+sigma2^2) % analysis uncertainty [K]



% =========================================================================

% Toy Example 2: stone in space



sb = 5.67e-8; % W m^-2 K^-4; Stefan-Boltzman constant

h = inline('5.67e-8*x^4'); % observation operator

hlin = inline('4*5.67e-8*x^3'); % dh/dT



Tb = 271 % background temperature [K]

yo = 300 % radiance observation [W/m^2]

H = hlin(Tb) % linearize the observation operator around Tb [(W/m^2)/K]

sigmab = 2 ^ 0.5 % background error [K]

sigmao = 10 ^ 0.5 % observation error [W/m^2]



obsincr = yo - h(Tb) % observation increment [W/m^2]

w = sigmab^2*H*(sigmao^2 + H * sigmab^2 * H)^-1 % weight [K/(W/m^2)]

Ta = Tb + w * obsincr % temperature analysis [K]



clear('H');

% =========================================================================

% Toy Example 3: model with 3 grid points

xb = [265 272 270]' % 3 model grid points

yo = [269 271]' % 2 observations: between grid points 2 & 3, 1 & 2

obserr = 1; % 1 K

R = [obserr^2 0; 0 obserr^2]

B = 2*[1 -0.5 -0.25; -0.5 1 0.5; -0.25 0.5 1]

H = [0 0.5 0.5; 0.5 0.5 0];

K = B*H'*inv(H*B*H' + R)

xa = xb + B*H'*inv(H*B*H' + R)*(yo-H*xb)



plot(xb,'r','Linewidth',2);

hold on

plot(xa,'b','Linewidth',2);

scatter(2.5, yo(1), 9,'g','filled');

scatter(1.5, yo(2), 9,'g','filled');

legend('Background','Analysis','Observation')

xlabel('Model Grid Point')

ylabel('Temperature [K]')

title('Data Assimilation Toy Example 3')

ylim([260 280])



