%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% %%%%%%%%%%%%%%%
%%%%%%%%%%% ENIAC %%%%%%%%%%%%%%%
%%%%%%%%%%% %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% %%%
%%% Solve the Barotropic Vorticity Equation %%%
%%% %%%
%%% d/dt (Del^2(z)) = J(h*Del^2(z)+f,z) (1) %%%
%%% %%%
%%% where h = g*m^2/f with map factor m. %%%
%%% %%%
%%% The formulation is the same as used by %%%
%%% Charney, Fjortoft and von Neumann (1950) in %%%
%%% the ENIAC integrations (Tellus, 2, 237-254). %%%
%%% %%%
%%% An alternative formulation of the BVE using %%%
%%% the stream function is used if the indicator %%%
%%% StreamFunction = 1 %%%
%%% It is formally identical to (1) above, but now %%%
%%% z stands for the streamfunction and h has a %%%
%%% different definition: h = m^2. %%%
%%% %%%
%%% The initial data are as used in tho original %%%
%%% integrations in 1950. The boundary conditions %%%
%%% are formulated in the same way. %%%
%%% %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% %%%
%%% The Barotropic Vorticity Equation is approximated %%%
%%% by centered finite differences. The domain has %%%
%%% M * N points. %%%
%%% %%%
%%% The boundary conditions are as described %%%
%%% in the Tellus paper: z is given on the boundary %%%
%%% and vorticity (xi) is given at inflow points. %%%
%%% %%%
%%% The time scheme is the leapfrog method. %%%
%%% (the first step is a forward step). %%%
%%% %%%
%%% The quantity to be stepped forward is %%%
%%% xi = Del^2(z) %%%
%%% Thus, after each time-step, it is necessary %%%
%%% to solve a Poisson equation to get z. %%%
%%% This is done by the Fourier Transform method, %%%
%%% as in the Tellus paper. (For a description of %%%
%%% the method, see Numerical Recipes Sec. 19.4.) %%%
%%% %%%
%%% A variety of initial conditions may be specified %%%
%%% determined by the parameter Ncase (see below). %%%
%%% %%%
%%% The data were downloaded from the NCEP/NCAR %%%
%%% Reanalysis site. See, for example, %%%
%%% http://nomad3.ncep.noaa.gov/ %%%
%%% cgi-bin/ftp2u_6p_r1.sh %%%
%%% %%%
%%% They were converted from GRIB to ASCII and %%%
%%% saved in files with the name YYYYMMDDHH.asc %%%
%%% where YYYYMMDDHH is the date and time of the %%%
%%% analysis. Thus, for example, for Ncase=1, %%%
%%% the file is 1949010503.asc. %%%
%%% %%%
%%% The re-analysis is on a 2.5 x 2.5 degree grid. %%%
%%% The data were recovered on a domain covering %%%
%%% the northern hemisphere. The fields were %%%
%%% interpolated to the 19 x 16 polar stereographic %%%
%%% grid used by Charney et al, and the results %%%
%%% saved in files with the names YYYYMMDDHH.z00 %%%
%%% %%%
%%% The four cases in question are (all in 1949): %%%
%%% %%%
%%% (1) 0300 Z, January 5th, 1949 %%%
%%% with verification analysis Jan 6 %%%
%%% (2) 0300 Z, January 30th, 1949 %%%
%%% with verification analysis Jan 31 %%%
%%% (3) 0300 Z, January 31st, 1949 %%%
%%% with verification analysis Feb 1 %%%
%%% (4) 0300 Z, February 13th, 1949 %%%
%%% with verification analysis Feb 14 %%%
%%% %%%
%%% Note: Charney et al. used the analyses valid %%%
%%% at 0300 Z. These were not immediately available %%%
%%% from the re-analysis site. %%%
%%% Thanks to Chi-Fan Shih for assistance in %%%
%%% acquiring these 0300Z fields. %%%
%%% %%%
%%% Ncase Analysis file Verification file %%%
%%% 1 Case1-1949010503.z00 Case1-1949010603.z00 %%%
%%% 2 Case2-1949013003.z00 Case2-1949013103.z00 %%%
%%% 3 Case2-1949013103.z00 Case2-1949020103.z00 %%%
%%% 4 Case2-1949002133.z00 Case2-1949021403.z00 %%%
%%% %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% %%%
%%% Date: February, 2007 %%%
%%% %%%
%%% Author: Peter Lynch, UCD Met & Climate Centre %%%
%%% University College Dublin %%%
%%% %%%
%%% Email: Peter.Lynch@ucd.ie %%%
%%% %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% diary DIARY % SAVE OUTPUT IN DIARY IF REQUIRED.
clear % Clear the memory.
clf % clear the display.
%%% Define the graphic display colours.
whitebg([1 1 0.50]); % Beige background.
colormap('jet');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Basic parameters that may be varied
Ncase = 1; % Forecast case number.
DAYLEN = 1; % Forecast length in days.
DtHours = 1; % Timestep in hours.
%% Indicator for psi-form of equation
StreamFunction = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Details of polar stereographic grid
% M: (= 19) points in x-direction
% N: (= 16) points in y-direction
% Ds: (= 736 km) Grid step at North Pole
% Xp: x-coordinate of North Pole (varies)
% Yp: y-coordinate of North Pole (varies)
% Centerangle: Angle between Date-line and
% positive y axis (varies).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set the PS grid parameters for 4 cases.
NUMx = [ 19 19 19 19 ]; % x dimension
NUMy = [ 16 16 16 16 ]; % y dimension
Xpol = [ 10 10 10 10 ]; % North Pole
Ypol = [ 14 12 14 14 ]; % North Pole
DELs = [ 736E+3 736E+3 736E+3 736E+3];
Cang = [ -90 -70 -35 -85 ];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set the parameters for the forecast (Ncase)
M = NUMx(Ncase); % Points in x direction
N = NUMy(Ncase); % Points in y direction
Xp = Xpol(Ncase); % Coords of North Pole
Yp = Ypol(Ncase); % Coords of North Pole
Ds = DELs(Ncase); % Grid step at NP (metres)
Centerangle = Cang(Ncase); % Central angle of map.
MN = M*N; % Total number of grid points.
%% Define the spatial grid (not used in this form).
%% x = (0:M-1) * Ds; % x-coordinates of grid points
%% y = (0:N-1) * Ds; % y-coordinates of grid points
% Define the (X,Y) grid (for plotting)
[X, Y ] = meshgrid(1:M,1:N);
X = X'; Y=Y'; % Transpose for standard plot.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Define the time variable
daylen = DAYLEN; % Integration time (in days)
seclen = daylen*24*60*60; % Integration time (in seconds)
Dt = DtHours*60*60; % Timestep in seconds
nt = seclen/Dt; % Total number of time-steps.
t = (0:nt)*Dt; % Time variable
time = t/(60*60); % Time in hours (for plotting)
nspd = (24*60*60)/Dt; % Time steps per day
%% Print out information on space and time grid.
fprintf('Ncase = %i \n', Ncase)
fprintf('Grid size, M=%i N=%i \n',M,N)
fprintf('Total grid points, M*N=%i \n',MN)
fprintf('Timesteps per day, nspd=%i \n',nspd)
fprintf('Stream function, %i \n',StreamFunction)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Define some geophysical constants (SI units)
a = (4*10^7)/(2*pi); % Radius of the Earth
grav = 9.80665; % Gravitational acceleration
Omega = 2*pi/(24*60*60); % Angular velocity of Earth.
f0 = 2*Omega*sin(pi/4); % Mean Coriolis parameter.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Compute latitude and longitude
%% on the polar stereographic grid.
r2d = 180/pi; % Conversion factor: radians to degrees.
d2r = pi/180; % Conversion factor: degrees to radians.
for ny=1:N
for nx=1:M
xx = (nx-Xp)*Ds;
yy = (ny-Yp)*Ds;
rr = sqrt(xx^2+yy^2);
theta = atan2(yy,xx);
lambda = theta + d2r*(90+Centerangle);
if(lambda>pi) lambda = lambda - 2*pi; end
phi = 2*((pi/4)-atan(rr/(2*a)));
LON(nx,ny) = lambda; % Longitude (radians)
LAT(nx,ny) = phi; % Latitude (radians)
%% Angles in degrees (for plotting)
lamd = r2d*lambda; phid = r2d*phi;
LONDEG(nx,ny) = lamd; % Longitude (degrees)
LATDEG(nx,ny) = phid; % Latitude (degrees)
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Record the latitude and longitude of the
%% corners of the outer domain and of the
%% inner domain (for checking purposes).
Corners = 0;
if(Corners == 1 )
latSW = LATDEG(1,1); lonSW = LONDEG(1,1);
latSE = LATDEG(M,1); lonSE = LONDEG(M,1);
latNE = LATDEG(M,N); lonNE = LONDEG(M,N);
latNW = LATDEG(1,N); lonNW = LONDEG(1,N);
fprintf(' Corner latitudes and longitudes, OUTER DOMAIN \n');
fprintf(' Lat and long of SW corner %g %g \n',latSW,lonSW);
fprintf(' Lat and long of SE corner %g %g \n',latSE,lonSE);
fprintf(' Lat and long of NE corner %g %g \n',latNE,lonNE);
fprintf(' Lat and long of NW corner %g %g \n',latNW,lonNW);
latSW = LATDEG(3 ,2 ); lonSW = LONDEG(3 ,2 );
latSE = LATDEG(M-2,2 ); lonSE = LONDEG(M-2,2 );
latNE = LATDEG(M-2,N-2); lonNE = LONDEG(M-2,N-2);
latNW = LATDEG(3 ,N-2); lonNW = LONDEG(3 ,N-2);
fprintf(' Corner latitudes and longitudes, INNER DOMAIN \n');
fprintf(' Lat and long of SW inner %g %g \n',latSW,lonSW);
fprintf(' Lat and long of SE inner %g %g \n',latSE,lonSE);
fprintf(' Lat and long of NE inner %g %g \n',latNE,lonNE);
fprintf(' Lat and long of NW inner %g %g \n',latNW,lonNW);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Compute Coriolis Parameter and Map Factor
%% and parameter h = g*m^2/f used in the BVE
%% (For psi-equation, h = m^2 is used).
for ny=1:N
for nx=1:M
lambda = LON(nx,ny);
phi = LAT(nx,ny);
map = 2 / (1+sin(phi));
f = 2*Omega*sin(phi);
MAP(nx,ny) = map;
FCOR(nx,ny) = f;
h(nx,ny) = grav * map^2 / f;
if (StreamFunction == 1) h(nx,ny) = map^2; end
end
end
clf; % Clear the plot window
PlotAll = 0;
if(PlotAll==1)
subplot(2,2,1)
contourf(X,Y,MAP); title('MAP FACTOR');
subplot(2,2,2)
contourf(X,Y,FCOR); title('f CORIOLIS');
subplot(2,2,3)
contourf(X,Y,h); title('FACTOR h');
pause
clf
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute Sine Matrices for the Poisson solver
%% Coefficients for x-transformation
for m1=1:M-2
for m2=1:M-2
SM(m1,m2) = sin(pi*m1*m2/(M-1));
end
end
%% Coefficients for y-transformation
for n1=1:N-2
for n2=1:N-2
SN(n1,n2) = sin(pi*n1*n2/(N-1));
end
end
%% Eigenvalues of Laplacian operator
for mm=1:M-2
for nn=1:N-2
eigen = (sin(pi*mm/(2*(M-1))))^2 +(sin(pi*nn/(2*(N-1))))^2;
EIGEN(mm,nn) = (-4/Ds^2) * eigen;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Read and plot the initial and
%%% verification height data
%%% Define the initial field sizes.
z0 = zeros(M,N); % Initial height
z24 = zeros(M,N); % Verifying Analysis
xi0 = zeros(M,N); % Initial Laplacian of height
eta0 = zeros(M,N); % Initial absolute vorticity
%% File name tag.
Case(1,1:6) = 'Case1-';
Case(2,1:6) = 'Case2-';
Case(3,1:6) = 'Case3-';
Case(4,1:6) = 'Case4-';
%% Initial analysis
YMDH1(1,1:10) = '1949010503';
YMDH1(2,1:10) = '1949013003';
YMDH1(3,1:10) = '1949013103';
YMDH1(4,1:10) = '1949021303';
%% Verifying analysis
YMDH2(1,1:10) = '1949010603';
YMDH2(2,1:10) = '1949013103';
YMDH2(3,1:10) = '1949020103';
YMDH2(4,1:10) = '1949021403';
%% Initial and verification analysis on PS grid
File1 = [Case(Ncase,:) YMDH1(Ncase,:) '.z00'];
File2 = [Case(Ncase,:) YMDH2(Ncase,:) '.z00'];
%% Read in the analyses on the PS grid
z0 = load(File1); % Initial height analysis
z24 = load(File2); % Verifying height analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Plot the analysis fields
zcontours = 4500:50:6000;
contourf(X,Y,z0,zcontours); title('Initial Data')
pause
contourf(X,Y,z24,zcontours); title('Final Analysis')
pause
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Plot inner region: omit one row at bottom
%% and two rows on the other three sides
%% (as plotted by Charney et al.)
rx = 3:M-2; ry = 2:N-2;
subplot(2,2,1)
contour(X(rx,ry),Y(rx,ry),z0(rx,ry),zcontours);
title('Initial Analysis')
subplot(2,2,2)
contour(X(rx,ry),Y(rx,ry),z24(rx,ry),zcontours);
title('Verifying Analysis')
subplot(2,2,3)
contour(X(rx,ry),Y(rx,ry),z24(rx,ry)-z0(rx,ry),21);
title('Analysed Change')
pause
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Plot the analysis on the inner region
%% omit one row at bottom and two rows
%% on the other three sides.
clf; % Clear the plot window.
rx = 3:M-2; ry = 2:N-2;
contour(X(rx,ry),Y(rx,ry),z0(rx,ry),zcontours);
title('INITIAL HEIGHT (inner area)')
hold on
rx = 3:M-2; ry = 2:N-4;
loncontours = linspace(-180,+180,37);
latcontours = linspace(0,90,10);
contour(X(rx,ry),Y(rx,ry),LONDEG(rx,ry),loncontours);
contour(X(rx,ry),Y(rx,ry),LATDEG(rx,ry),latcontours);
hold off
pause
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define working arrays to have correct size
ddxz = zeros(M,N); % x derivative of z
ddyz = zeros(M,N); % y derivative of z
gradzsq = zeros(M,N); % squared gradiant of z
d2dx2z = zeros(M,N); % second x derivative of z
d2dy2z = zeros(M,N); % second y derivative of z
xi = zeros(M,N); % Laplacian of z
eta = zeros(M,N); % absolute vorticity
ddxeta = zeros(M,N); % x derivative of eta
ddyeta = zeros(M,N); % y derivative of eta
Jacobi = zeros(M,N); % Jacobian J(eta,z)
ddtz = zeros(M,N); % Tendency of z
ddtxi = zeros(M,N); % Tendency of xi
zdot = zeros(M-2,N-2); % Interior values of ddtz;
xidot = zeros(M-2,N-2); % Interior values of ddtxi;
ZDOT = zeros(M-2,N-2); % Fourier transform of zdot;
XIDOT = zeros(M-2,N-2); % Fourier transform of xidot;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Install the initial conditions for the forecast
if(StreamFunction==1)
% z0 = grav * z0 ./ FCOR; % Initial streamfunction
z0 = grav * z0 / f0; % Initial streamfunction
end
z = z0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Compute the Laplacian of the height/psi
% Second x-derivative of z
d2dx2z(2:M-1,:) = (z(3:M,:)+z(1:M-2,:)-2*z(2:M-1,:))/(Ds^2);
% Second y-derivative of z
d2dy2z(:,2:N-1) = (z(:,3:N)+z(:,1:N-2)-2*z(:,2:N-1))/(Ds^2);
%% Laplacian of height (or vorticity)
xi0(2:M-1,2:N-1) = d2dx2z(2:M-1,2:N-1)+d2dy2z(2:M-1,2:N-1);
%% Extend xi0 to boundaries (for Jacobian).
xi0(1,:) = 2*xi0(2 ,:)-xi0(3 ,:); % West
xi0(M,:) = 2*xi0(M-1,:)-xi0(M-2,:); % East
xi0(:,1) = 2*xi0(:,2 )-xi0(:,3 ); % South
xi0(:,N) = 2*xi0(:,N-1)-xi0(:,N-2); % North
%% Absolute vorticity
eta0 = h.*xi0 + FCOR;
%% Compute statistics of initial fields
zmin = min(min(z0));
zave = mean(mean(z0));
zmax = max(max(z0));
fprintf(' z0: min, ave, max %g %g %g \n', zmin, zave, zmax)
ximin = min(min(xi0));
xiave = mean(mean(xi0));
ximax = max(max(xi0));
fprintf(' xi0: min, ave, max %g %g %g \n', ximin, xiave, ximax)
etamin = min(min(eta0));
etaave = mean(mean(eta0));
etamax = max(max(eta0));
fprintf('eta0: min, ave, max %g %g %g \n',etamin,etaave,etamax)
% Plot the initial height, Laplacian and Vorticity
subplot(2,2,1)
contourf(X,Y,z0,21);
title('INITIAL HEIGHT/STREAMFUNCTION');
subplot(2,2,2)
contourf(X,Y,xi0,21);
title('INITIAL LAPLACIAN OF HEIGHT/PSI');
subplot(2,2,3)
contourf(X,Y,eta0,21);
title('INITIAL ABSOLUTE VORTICITY');
subplot(2,2,4)
contourf(X,Y,h);
title('FACTOR h');
fprintf('Press RETURN to continue \n');
pause;
clf; % Clear the plot window.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% MAIN LOOP %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t1elapsed = clock; % Starting time
t1cpu = cputime; % For CPU time.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% %%%
%%% Integrate the BVE in time %%%
%%% %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%
%%% Time-stepping is by leapfrog method
%%% (first step is forward)
%%%
%%% Define xi = Del^2(z). The BVE is
%%%
%%% (d/dt)xi = J( h*xi+f, z)
%%%
%%% We approximate the time derivative by
%%% ( xi(n+1)-xi(n-1) ) / (2*Dt)
%%%
%%% The Jacobian term J(eta,z) is approximated
%%% by centered space differences
%%%
%%% Then the value of xi at the new time (n+1)*Dt is:
%%% xi(n+1) = xi(n-1) + 2*Dt * J(n)
%%%
%%% When we have ddtxi, we have to solve a Poisson
%%% equation to get ddtz. Then both xi and z are stepped
%%% forward to (n+1)*Dt and the cycle is repeated.
%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(' Total Number of Steps: %i \n', nt)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% Start of the time-stepping loop %%%%%%%%%
for n=1:nt
% Print out tracking information.
PrintStep = 0;
if ( PrintStep == 1 )
fprintf(' Step number, time(h), time(d) %i %g %g \n', ...
n,n*Dt/3600,n*Dt/(3600*24));
end
%% First time through the loop:
if(n==1)
dt = Dt/2; % First step is forward
znm1 = z0; % Copy initial height field
xinm1 = xi0; % Copy initial vorticity field
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Compute the derivatives, Laplacian and Jacobian.
% x-derivative of z
ddxz(2:M-1,:) = (z(3:M,:)-z(1:M-2,:))/(2*Ds);
% y-derivative of z
ddyz(:,2:N-1) = (z(:,3:N)-z(:,1:N-2))/(2*Ds);
% Square of the gradient of z
gradzsq = ddxz.^2+ddyz.^2;
% Second x-derivative of z
d2dx2z(2:M-1,:) = (z(3:M,:)+z(1:M-2,:)-2*z(2:M-1,:))/(Ds^2);
% Second y-derivative of z
d2dy2z(:,2:N-1) = (z(:,3:N)+z(:,1:N-2)-2*z(:,2:N-1))/(Ds^2);
%% Laplacian of height
%% xi = d2dx2z + d2dy2z;
xi(2:M-1,2:N-1) = d2dx2z(2:M-1,2:N-1)+d2dy2z(2:M-1,2:N-1);
%% Extend xi to boundaries (to compute Jacobian).
%% (First time step only; xi from BCs after that)
if ( n == 1 )
xi(1,:) = 2*xi(2 ,:)-xi(3 ,:); % West
xi(M,:) = 2*xi(M-1,:)-xi(M-2,:); % East
xi(:,1) = 2*xi(:,2 )-xi(:,3 ); % South
xi(:,N) = 2*xi(:,N-1)-xi(:,N-2); % North
end
%% Absolute vorticity
eta = h.*xi + FCOR;
% x-derivative of eta
ddxeta(2:M-1,:) = (eta(3:M,:)-eta(1:M-2,:))/(2*Ds);
% y-derivative of eta
ddyeta(:,2:N-1) = (eta(:,3:N)-eta(:,1:N-2))/(2*Ds);
%% Compute the Jacobian J(eta,z)
Jacobi = ddxeta .* ddyz - ddyeta .* ddxz;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the energy and enstrophy integrals
E(n) = 0.5 * sum(sum(gradzsq));
S(n) = 0.5 * sum(sum(xi.^2));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Solve the Poisson Equation Del^2(ddtz) = ddtxi
%% with homogeneous boundary conditions: z is
%% constant on the boundaries, so ddtz vanishes.
%% Note: Fourier transform of xidot denoted XIDOT
%% Fourier transform of zdot denoted ZDOT.
%% Forward Fourier transform:
%% XIDOT = SM*xidot*SN;
%% Inverse transform:
%% zdot = (4/((M-1)*(N-1)))*SM*ZDOT*SN;
%%
% Tendency values in interior.
xidot = Jacobi(2:M-1,2:N-1);
% Compute the transform of the solution
XIDOT = SM*xidot*SN;
% Convert transform of d(xi)/dt to transform of d(z)/dt
ZDOT = XIDOT ./ EIGEN;
% Compute inverse transform to get the height tendency.
zdot = (4/((M-1)*(N-1))) * SM*ZDOT*SN;
ddtz (2:M-1,2:N-1) = zdot; % Insert inner values
ddtxi(2:M-1,2:N-1) = xidot; % Insert inner values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Compute ddtxi on the boundaries.
%% If fluid is entering through the boundary,
%% we set ddtxi to zero. If fluid is leaving
%% the region, we extrapolate ddtxi linearly
%% from the interior. Charney, et al., Eqn (21).
% Western boundary
sigW = max(0,+sign(ddyz(1,:)));
ddtxi(1,:) = sigW.*(2*ddtxi(2 ,:)-ddtxi(3 ,:));
% Eastern boundary
sigE = max(0,-sign(ddyz(M,:)));
ddtxi(M,:) = sigE.*(2*ddtxi(M-1,:)-ddtxi(M-2,:));
% Southern boundary
sigS = max(0,-sign(ddxz(:,1)));
ddtxi(:,1) = sigS.*(2*ddtxi(:,2 )-ddtxi(:,3 ));
% Northern boundary
sigN = max(0,+sign(ddxz(:,N)));
ddtxi(:,N) = sigN.*(2*ddtxi(:,N-1)-ddtxi(:,N-2));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step forward one time-step (leapfrog scheme).
xinp1 = xinm1 + (2*dt) * ddtxi;
znp1 = znm1 + (2*dt) * ddtz;
% Move the old values into arrays znm1 and xinm1
% Move the new values into arrays z and xi
znm1 = z;
xinm1 = xi;
z = znp1;
xi = xinp1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Save the fields at quarterpoints of the integration
if(n==1*nt/4) zq1 = z; end
if(n==2*nt/4) zq2 = z; end
if(n==3*nt/4) zq3 = z; end
if(n==4*nt/4) zq4 = z; end
% Plot the height field at each time step (if required)
plotsteps = 0; % Turn plotting off for timing.
if(plotsteps==1)
surf(X,Y,z); shading('interp')
axis('off'); view(0,90);
drawnow;
end
% Restore the timestep (after the first step)
dt = Dt;
end %%%%% End of the time-stepping loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Get the elapsed ime the integration
t2elapsed = clock; % Starting time
ElapsedTime = etime(t2elapsed,t1elapsed); % Elapsed time
t2cpu = cputime;
CPUtime = t2cpu - t1cpu; % For CPU time.
fprintf(' Elapsed time %g seconds \n', ElapsedTime)
fprintf(' CPU time %g seconds \n', CPUtime)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% END MAIN LOOP %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the energy and enstrophy
ddxz(2:M-1,:) = (z(3:M,:)-z(1:M-2,:))/(2*Ds);
ddyz(:,2:N-1) = (z(:,3:N)-z(:,1:N-2))/(2*Ds);
gradzsq = ddxz.^2+ddyz.^2;
E(nt+1) = 0.5 * sum(sum(gradzsq));
S(nt+1) = 0.5 * sum(sum(xi.^2));
% Plot the Energy and Enstrophy Integrals.
plotnn=88; % Indicator for plotting
if(plotnn==88)
clf; % Clear the plot window.
plot(time,E); title('Total Energy');
axis([0,time(nt),0,2*E(1)]);
fprintf('Press RETURN to continue \n'); pause
plot(time,S); title('Total Enstrophy');
axis([0,time(nt),0,2*S(1)]);
fprintf('Press RETURN to continue \n'); pause
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Convert back from psi to z if necessary
if(StreamFunction==1)
z0 = f0 * z0 / grav; % Initial streamfunction
z = f0 * z / grav; % Initial streamfunction
xi = f0 * xi / grav; % Initial Laplacian
zq1 = f0 * zq1 / grav; % Solution at quarter-point
zq2 = f0 * zq2 / grav; % Solution at quarter-point
zq3 = f0 * zq3 / grav; % Solution at quarter-point
zq4 = f0 * zq4 / grav; % Solution at quarter-point
end
% Plot the final height field
contourf(X,Y,z,zcontours)
title('FORECAST HEIGHT FIELD');
fprintf('Press RETURN to continue \n');
pause
hold off
% Plot the final height, laplacian and vorticity
subplot(2,2,1)
contourf(X,Y,z,zcontours);
title('FORECAST GEOPOTENTIAL HEIGHT');
subplot(2,2,2)
contourf(X,Y,xi,21);
title('FORECAST LAPLACIAN OF HEIGHT');
subplot(2,2,3)
contourf(X,Y,eta,21);
title('FORECAST ABSOLUTE VORTICITY');
fprintf('Press RETURN to continue \n');
pause; clf; % Clear the plot window.
%% Plot the heights at check-points.
rx = 1:M ; ry = 1:N ; % Full area
subplot(2,2,1)
contour(X(rx,ry),Y(rx,ry),zq1(rx,ry),zcontours);
title('Z (6 hours)')
subplot(2,2,2)
contour(X(rx,ry),Y(rx,ry),zq2(rx,ry),zcontours);
title('Z (12 hours)')
subplot(2,2,3)
contour(X(rx,ry),Y(rx,ry),zq3(rx,ry),zcontours);
title('Z (18 hours)')
subplot(2,2,4)
contour(X(rx,ry),Y(rx,ry),zq4(rx,ry),zcontours);
title('Z (24 hours)')
pause
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Plot inner region: omit one row at bottom
%% and two rows on the other three sides.
%% These figures are in the same form as
%% Figs 2 to 5 in Charney et al. (1950)
rx = 3:M-2; ry = 2:N-2;
subplot(2,2,1)
contour(X(rx,ry),Y(rx,ry),z0(rx,ry),zcontours,'b');
title('(A) INITIAL ANALYSIS')
axis off; hold on;
plot([3 M-2 M-2 3 3],[2 2 N-2 N-2 2],'k');
subplot(2,2,2)
contour(X(rx,ry),Y(rx,ry),z24(rx,ry),zcontours,'b');
title('(B) VERIFYING ANALYSIS')
axis off; hold on;
plot([3 M-2 M-2 3 3],[2 2 N-2 N-2 2],'k');
subplot(2,2,3)
contour(X(rx,ry),Y(rx,ry),z24(rx,ry)-z0(rx,ry),11,'b-');
hold on
contour(X(rx,ry),Y(rx,ry),z(rx,ry)-z0(rx,ry),11,'r--');
title('(C) ANALYSED & FORECAST CHANGES')
axis off; hold on;
plot([3 M-2 M-2 3 3],[2 2 N-2 N-2 2],'k');
subplot(2,2,4)
contour(X(rx,ry),Y(rx,ry),z(rx,ry),zcontours,'b');
title('(D) FORECAST HEIGHT')
axis off; hold on;
plot([3 M-2 M-2 3 3],[2 2 N-2 N-2 2],'k');
pause
%% Save the four-panel image.
PrintFiles = 1;
if ( PrintFiles == 1 )
print -depsc 'NEW-Case-N.eps'
print -djpeg 'NEW-Case-N.jpg'
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calculate the Forecast & Persistence BIAS
sum0 = 0; sum1 = 0; sum2 = 0; sum3 = 0;
rx = 3:M-2; ry = 2:N-2;
for nx=rx
for ny=ry
sum0 = sum0 + 1;
sum1 = sum1 + ( z(nx,ny)-z24(nx,ny));
sum2 = sum2 + (z0(nx,ny)-z24(nx,ny));
sum3 = sum3 + ( z(nx,ny)- z0(nx,ny));
end
end
rms1 = (sum1/sum0);
rms2 = (sum2/sum0);
rms3 = (sum3/sum0);
fprintf(' Bias of 24h forecast: %g \n', rms1);
fprintf(' Bias of persistence: %g \n', rms2);
fprintf(' Mean 24h forecast change: %g \n', rms3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calculate the Forecast & Persistence RMS Errors
sum0 = 0; sum1 = 0; sum2 = 0; sum3 = 0;
rx = 3:M-2; ry = 2:N-2;
for nx=rx
for ny=ry
sum0 = sum0 + 1;
sum1 = sum1 + ( z(nx,ny)-z24(nx,ny))^2;
sum2 = sum2 + (z0(nx,ny)-z24(nx,ny))^2;
sum3 = sum3 + ( z(nx,ny)- z0(nx,ny))^2;
end
end
rms1 = sqrt(sum1/sum0);
rms2 = sqrt(sum2/sum0);
rms3 = sqrt(sum3/sum0);
fprintf(' RMS error of 24h forecast: %g \n', rms1);
fprintf(' RMS error of persistence: %g \n', rms2);
fprintf(' RMS 24h forecast change: %g \n', rms3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Compute the S1 Scores for the Forecast and
%% for a persistence forecast
ddxz(2:M-1,:) = (z(3:M,:)-z(1:M-2,:))/(2*Ds);
ddyz(:,2:N-1) = (z(:,3:N)-z(:,1:N-2))/(2*Ds);
gradzsq = ddxz.^2+ddyz.^2;
sgradz = sum(sum(sqrt(gradzsq)));
ddxz(2:M-1,:) = (z0(3:M,:)-z0(1:M-2,:))/(2*Ds);
ddyz(:,2:N-1) = (z0(:,3:N)-z0(:,1:N-2))/(2*Ds);
gradzsq = ddxz.^2+ddyz.^2;
sgradz0 = sum(sum(sqrt(gradzsq)));
ddxz(2:M-1,:) = (z24(3:M,:)-z24(1:M-2,:))/(2*Ds);
ddyz(:,2:N-1) = (z24(:,3:N)-z24(:,1:N-2))/(2*Ds);
gradzsq = ddxz.^2+ddyz.^2;
sgradz24 = sum(sum(sqrt(gradzsq)));
ddxz(2:M-1,:) = ( z(3:M,:)- z(1:M-2,:))/(2*Ds) ...
-(z24(3:M,:)-z24(1:M-2,:))/(2*Ds);
ddyz(:,2:N-1) = ( z(:,3:N)- z(:,1:N-2))/(2*Ds) ...
-(z24(:,3:N)-z24(:,1:N-2))/(2*Ds);
gradzsq = ddxz.^2+ddyz.^2;
sgradzmz24 = sum(sum(sqrt(gradzsq)));
S1fcst = (sgradzmz24/max(sgradz,sgradz24)) * 100;
ddxz(2:M-1,:) = ( z0(3:M,:)- z0(1:M-2,:))/(2*Ds) ...
-(z24(3:M,:)-z24(1:M-2,:))/(2*Ds);
ddyz(:,2:N-1) = ( z0(:,3:N)- z0(:,1:N-2))/(2*Ds) ...
-(z24(:,3:N)-z24(:,1:N-2))/(2*Ds);
gradzsq = ddxz.^2+ddyz.^2;
sgradz0mz24 = sum(sum(sqrt(gradzsq)));
S1Pers = (sgradz0mz24/max(sgradz0,sgradz24)) * 100;
fprintf(' S1 score for 24h forecast: %g \n', S1fcst);
fprintf(' S1 score for persistence: %g \n', S1Pers);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%