% QCOM.m

% All array subscripts start from 1

clear all; 

global fv jt jv kt kv v dt ; 

 jt = 20+1;
 kt = 20+1;
 jv = jt;
 kv = kt;
 v=zeros(jv+1,kv+1);
 fv=zeros(jv,kv, 2);

%     Assign parameters for run

ittmax = 100;
dt = 1.;

%     Initialize all values of arrays for v, w, theta, pi, [fv, fw, ftheta, fpi,]
%     including values for boundary points for v, w, theta, pi.

% initialize all variables
init;

% itt is time step index
itt = 1;

%     USE FORWARD SCHEME TO do first step

a = 1.;
b = 0.;
n1 = rem( itt    , 2 ) + 1;
n2 = rem( itt - 1, 2 ) + 1;

% do first time step
[ n1, n2, a, b ]=step( n1, n2, a, b );

%     ADAMS - BASHFORTH TWO - LEVEL SCHEME

a =   3. ./ 2.;
b = - 1. ./ 2.;

ittnow = itt  + 1;

for itt = ittnow: ittmax;

n1 = rem( itt    , 2 ) + 1;
n2 = rem( itt - 1, 2 ) + 1;

% do subsequent time steps
[ n1, n2, a, b ]=step( n1, n2, a, b );

end
itt = ittmax+1

%     END-OF-RUN OUTPUT ROUTINES GO HERE


function [n1, n2, a, b]=step( n1, n2, a, b ,varargin);

%     This is the entire subroutine.

% calculate forcing terms from variables at current time

[ n2 ]=rcalc( n2 );
% update variables using a time scheme
[ n1, n2, a, b ]=ab( n1, n2, a, b );
% apply boundary conditions to variables
bound;
%
end %subroutine step

function [n2]=rcalc(  n2 ,varargin);
global fv jt kt ;
%     CALCULATES FORCING TERMS FOR V(J,K), ETC.; STORES THEM  IN FV(J,K,N2), ETC.

for k = 1: kt;
for j = 1: jt;
fv(j,k,n2)  = 1.;
end
end

%     ETC (forcing for w, theta, and pi)

end %subroutine rcalc

function [n1, n2, a, b]=ab( n1, n2, a, b ,varargin);
global fv jt kt v dt ;
%     THE FOLLOWING LOOP UPDATES V USING EITHER THE FORWARD OR THE ADAMS-BASHFORTH
%     SCHEME DEPENDING ON THE VALUES OF A, B.
%     SUBSCRIPT N2 OF FV ALWAYS REFERS TO THE MOST RECENTLY CALCULATED VALUES FOR FV.

for k = 1: kt;
for j = 1: jt;
v(j,k)  = v(j,k)  + dt .*( a .* fv(j,k,n2)  + b .* fv(j,k,n1) );
end
end

%     ETC (update w, theta, and pi)

end %subroutine ab

function bound(varargin)
global jt kt v ;
%     apply boundary conditions for v

% Matlab does not like 0 subscripts

for j = 1: jt;
% free slip b.c.
v(j,1) = v(j,2);
% free slip b.c.
v(j,kt+1) = v(j,kt);
end

for k = 1: kt+1;
% periodic b.c.
v(1,k) = v(jt,k);
% periodic b.c.
v(jt+1,k) = v(2,k);
end

%     ETC (apply b.c. for w, theta, and pi)

end %subroutine bound

function init(varargin)
global jt kt v ;
%     initialize all variables


for k = 1: kt;
for j = 1: jt;
v(j,k)  = 0.;
end; j = jt+1;
end; k = kt+1;

bound;

end %subroutine init