-
Notifications
You must be signed in to change notification settings - Fork 0
/
qg2c_step.m
160 lines (133 loc) · 4.5 KB
/
qg2c_step.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
function []= qg2c_step(q1, q2, by, U1, U2, beta1, beta2, qforce)
global L W nx ny a11 a12 a21 a22 b11 b12 b21 b22 massfac Hfree del F1 F2 rd
dt=1/128; % time step
tmax=10; % max time
tpl=1/dt; % time plot in number of time steps
dx=L/nx; % grid spacing
dy=W/ny; % grid spacing y
[xh,yh]=meshgrid([0:nx]*dx,[0:ny]*dy); % extended grid with end points
t=0; % time
tc=0; % count
q1_p=q1;q2_p=q2;
zy=zeros(1,nx);
psimax=[];
ts=[];
ubs=[];
qps=[];
while t<=tmax+dt/2
% second order Adams-Bashforth
tmp=q1;q1=1.5*q1-0.5*q1_p;q1_p=tmp; % extrapolation
tmp=q2;q2=1.5*q2-0.5*q2_p;q2_p=tmp;
% p1 = psi1, p2 = psi2 total psi
[p1,p2]=invert(q1,q2,a11,a12,a21,a22,b11,b12,b21,b22,massfac,Hfree,del);
u1=-p(p1,dy); % differentiate
v1=p(p1',dx)';
u2=-p(p2,dy);
v2=p(p2',dx)';
if(rem(tc,tpl)==0) % plot t_0 time point
ff=[rs(q2_p);zy;zy;rs(q1_p)];
imagesc(ff);axis('xy');title(num2str(t));drawnow();
ts=[ts,t];
ubs=[ubs,mean(u1')'];
qp=q1_p-mean(q1_p')'*ones(1,nx);
qps=[qps,max(qp(:))];
end
qt1=q1+by; % add background PV
qt2=q2+by;
Fy=flux(q1,v1,dy,dt,0); % compute northward flux of total pv vq
Fy(1,:)=0; Fy(ny+1,:)=0;
Fx=flux(q1',u1'+U1',dx,dt,1)'; % compute eastward flux of total pv uq
dq1dt=-p(Fx',dx)'-p(Fy,dy)-beta1.*v1(2:end,:)-rd*q1-rd*F1*qforce; % PV equation including linear damping and forcing
Fy=flux(q2,v2,dy,dt,0);
Fy(1,:)=0; Fy(ny+1,:)=0;
Fx=flux(q2',u2'+U2',dx,dt,1)';
dq2dt=-p(Fx',dx)'-p(Fy,dy)-beta2.*v2(2:end,:)-rd.*q2+rd.*F2*qforce; % other layer
% keyboard;
q1=q1_p+dt*dq1dt; % give new PV
q2=q2_p+dt*dq2dt;
tc=tc+1; % time step
t=tc*dt;
end
plot(ts,max(ubs));
function df=p(f,dy) % compute derivative in y
df=(1/dy)*diff(f);
return
end
function qp=rs(q) % rescaling
del=max(max(q))-min(min(q));
if(del==0)
qp=q;
return;
end
qp=(q-min(min(q)))/del;
end
% compute psi, q sep. and compute correct advection, remove zonal mean from
% dq/dt
function fa=flux(f,v,dy,dt,bcf) % compute fluxes with flux correction scheme (limit magnitude)
if(bcf==0)
f=[-f(1,:);f;-f(end,:)];
else
f=[f(end,:);f;f(1,:)];
end
n=size(f,1);nm1=n-1;
fbar=0.5*(f(1:nm1,:)+f(2:n,:));
delf=f(2:n,:)-f(1:nm1,:);
absv=abs(v);
fup=v.*fbar-0.5*absv.*delf;
flw=v.*fbar-0.5*dt/dy*absv.^2.*delf;
% fbar=shift(f,1);
delfp=circshift(delf,1);
delfm=circshift(delf,-1);
r=((v>=0).*delfp+(v<0).*delfm)./delf;
if(sum(sum(isnan(r)))>0)
r(isnan(r))=0;
end
if(sum(sum(isinf(r)))>0)
r(isinf(r))=1e20*sign(r(isinf(r)));
end
% van leer
psi=(r+abs(r))./(1+abs(r));
% superbee
% psi=max(0,max(min(1,2*r),min(2,r)));
% if(sum(sum(isnan(psi)))>0)
% psi(isnan(psi))=0;
% endif
fa=fup+psi.*(flw-fup);
% fa=flw;
end
% compute PV inversion
function [p1,p2]=invert(q1,q2,a11,a12,a21,a22,b11,b12,b21,b22,massfac,Hfree,del)
nx=size(q1,2);
ny=size(q1,1);
z1b=mean(q1')'; % compute means
z2b=mean(q2')';
q1=q1-z1b*ones(1,nx);
q2=q2-z2b*ones(1,nx); % interpolates to corners of box
z1=(q1(1:ny-1,:)+q1(1:ny-1,[nx,1:nx-1])...
+q1(2:ny,:)+q1(2:ny,[nx,1:nx-1]))/4;
z2=(q2(1:ny-1,:)+q2(1:ny-1,[nx,1:nx-1])...
+q2(2:ny,:)+q2(2:ny,[nx,1:nx-1]))/4;
z1b=dct(z1b); % inversions for mean for cosine transform
z2b=dct(z2b);
p1b=idct(b11.*z1b+b12.*z2b);
p2b=idct(b21.*z1b+b22.*z2b); % average back because of cosine
p1b=[0.5*p1b(1)+0.5*p1b(2);0.5*(p1b(1:ny-1)+p1b(2:ny));0.5*p1b(ny-1)+0.5*p1b(ny)];
p2b=[0.5*p2b(1)+0.5*p2b(2);0.5*(p2b(1:ny-1)+p2b(2:ny));0.5*p2b(ny-1)+0.5*p2b(ny)];
delmass=massfac*(p1b-p2b); % conserve mass
p1b=p1b-delmass*Hfree;
p2b=p2b+del*delmass*Hfree;
z1=fft(dst(z1)')'; % Fourier transform for perturbations
z2=fft(dst(z2)')';
p1=a11.*z1+a12.*z2;
p2=a21.*z1+a22.*z2;
p1=real([zeros(1,nx);...
idst(ifft(p1')');...
zeros(1,nx)]);
p2=real([zeros(1,nx);...
idst(ifft(p2')');...
zeros(1,nx)]);
p1=[p1,p1(:,1)]+p1b*ones(1,nx+1);
p2=[p2,p2(:,1)]+p2b*ones(1,nx+1);
return
end
end