1 |
% |
2 |
% function [XT,YT,TADV] = calc_hadv(T,U,V,DX,DY,Ymin) |
3 |
% |
4 |
% Computes horizontal advection of scalar T |
5 |
% on the XT, YT (physics) grid at a given vertical level. |
6 |
% NB: (T,U,V) are 2-D fields |
7 |
% |
8 |
% (U,V) on a C-grid (XU,YU) (XV,YV) |
9 |
% NB: The way it is computed is consistent with |
10 |
% the flux form used by MIT-GCM (i.e. if |
11 |
% T DIV is added one recovers the flux form) |
12 |
% |
13 |
% Ymin is the southern latitude (negative, in degree) |
14 |
% DX is the longitudinal resolution (in degree) |
15 |
% DY is the latitudinal resolution (in degree) |
16 |
% |
17 |
% (c) acz, Jul. 2003 |
18 |
|
19 |
|
20 |
function [XT,YT,TADV] = calc_hadv(T,U,V,DX,DY,Ymin) |
21 |
|
22 |
% C-grid |
23 |
% |
24 |
[NX NY] = size(U); %or V |
25 |
XU = [0:DX:(DX*NX-DX)]; |
26 |
XV = XU + DX/2; |
27 |
YU = [(Ymin+DY/2):DY:(-Ymin-DY/2)]; |
28 |
YV = [Ymin:DY:-Ymin-DY]; |
29 |
XT = XV; YT = YU; |
30 |
|
31 |
% Constants |
32 |
RADIUS = 6371 * 1000; |
33 |
DYG = RADIUS * DY * pi/180; |
34 |
DXG = RADIUS * DX * pi/180; |
35 |
|
36 |
% Calculate zonal advection on U-grid |
37 |
advu = zeros(NX,NY); |
38 |
for i = 1:NX-1 |
39 |
advu(i+1,:) = U(i+1,:) .* (T(i+1,:)-T(i,:)); |
40 |
end |
41 |
advu(1,:) = U(1,:) .* (T(1,:)-T(NX,:)); |
42 |
|
43 |
% Average advu on T-grid |
44 |
advuTG = zeros(NX,NY); |
45 |
AG = cos(YT*pi/180) * DYG * DXG; |
46 |
for i = 1:NX-1 |
47 |
advuTG(i,:) = DYG*( advu(i,:)+advu(i+1,:) ) ./ (2*AG); |
48 |
end |
49 |
advuTG(NX,:) = DYG*( advu(NX,:)+advu(1,:) ) ./ (2*AG); |
50 |
|
51 |
% Calculate meridional advection on V-grid |
52 |
advv = zeros(NX,NY); %note advv(:,1) = 0 because v(:,1)=0 |
53 |
for j = 2:NY |
54 |
advv(:,j) = V(:,j) .* (T(:,j)-T(:,j-1)); |
55 |
end |
56 |
|
57 |
% Average advv on T-grid |
58 |
advvTG = zeros(NX,NY); |
59 |
for j = 1:NY-1 |
60 |
AG(j) = DYG * DXG * ( cos(YV(j+1)*pi/180)+cos(YV(j)*pi/180) )/2; |
61 |
end |
62 |
AG(NY) = DYG * DXG * ( cos((YV(NY)+DY)*pi/180)+cos(YV(NY)*pi/180) )/2; |
63 |
DDXG = DXG * cos((YV+DY)*pi/180); |
64 |
for j = 1:NY-1 |
65 |
advvTG(:,j) = DDXG(j)*( advv(:,j)+advv(:,j+1) ) ./ (2*AG(j)); |
66 |
end |
67 |
advvTG(:,NY) = DDXG(NY)*( advv(:,NY)+0 ) ./ (2*AG(NY)); |
68 |
|
69 |
% Horizontal advection |
70 |
TADV = advuTG + advvTG; |
71 |
|