1 |
dimitri |
1.1 |
% generate rotated u/v vectors and compare to model output |
2 |
|
|
% and to those used by global_ocean.cs32x15 |
3 |
|
|
|
4 |
|
|
clear all, close all |
5 |
|
|
|
6 |
|
|
% load model grid |
7 |
|
|
xc=readbin('XC.data',[32,6,32]);xc=permute(xc,[1 3 2]); |
8 |
|
|
xg=readbin('XG.data',[32,6,32]);xg=permute(xg,[1 3 2]); |
9 |
|
|
yc=readbin('YC.data',[32,6,32]);yc=permute(yc,[1 3 2]); |
10 |
|
|
yg=readbin('YG.data',[32,6,32]);yg=permute(yg,[1 3 2]); |
11 |
|
|
%tmp=yg; cx=[min(tmp(:)) max(tmp(:))]; plot_cube |
12 |
|
|
|
13 |
|
|
% load input velocity on Gaussian NCEP grid |
14 |
dimitri |
1.2 |
u=readbin('uwind_192_94_12.bin',[192 94]); |
15 |
|
|
v=readbin('vwind_192_94_12.bin',[192 94]); |
16 |
dimitri |
1.1 |
lon=0:1.875:358.125; |
17 |
|
|
lat=[-88.5420, -86.6532, -84.7532, -82.8508, -80.9474, -79.0435, ... |
18 |
|
|
-77.1393, -75.2351, -73.3307, -71.4262, -69.5217, -67.6171, ... |
19 |
|
|
-65.7125, -63.8079, -61.9033, -59.9986, -58.0940, -56.1893, ... |
20 |
|
|
-54.2846, -52.3799, -50.4752, -48.5705, -46.6658, -44.7611, ... |
21 |
|
|
-42.8564, -40.9517, -39.0470, -37.1422, -35.2375, -33.3328, ... |
22 |
|
|
-31.4281, -29.5234, -27.6186, -25.7139, -23.8092, -21.9044, ... |
23 |
|
|
-19.9997, -18.0950, -16.1902, -14.2855, -12.3808, -10.4760, ... |
24 |
|
|
- 8.5713, - 6.6666, - 4.7618, - 2.8571, - 0.9524, ... |
25 |
|
|
0.9524, 2.8571, 4.7618, 6.6666, 8.5713, 10.4760, 12.3808, ... |
26 |
|
|
14.2855, 16.1902, 18.0950, 19.9997, 21.9044, 23.8092, ... |
27 |
|
|
25.7139, 27.6186, 29.5234, 31.4281, 33.3328, 35.2375, ... |
28 |
|
|
37.1422, 39.0470, 40.9517, 42.8564, 44.7611, 46.6658, ... |
29 |
|
|
48.5705, 50.4752, 52.3799, 54.2846, 56.1893, 58.0940, ... |
30 |
|
|
59.9986, 61.9033, 63.8079, 65.7125, 67.6171, 69.5217, ... |
31 |
|
|
71.4262, 73.3307, 75.2351, 77.1393, 79.0435, 80.9474, ... |
32 |
|
|
82.8508, 84.7532, 86.6532, 88.5420]; |
33 |
|
|
|
34 |
|
|
% interpolate u and v on xc, yc |
35 |
|
|
u1=[u(192,:); u; u(1,:)]; |
36 |
|
|
u1=[u1(:,1) u1 u1(:,94)]; |
37 |
|
|
v1=[v(192,:); v; v(1,:)]; |
38 |
|
|
v1=[v1(:,1) v1 v1(:,94)]; |
39 |
|
|
lon1=-1.875:1.875:360; |
40 |
|
|
lat1=[-90 lat 90]; |
41 |
|
|
u2=interp2(lon1,lat1,u1',xc(:),yc(:)); |
42 |
|
|
v2=interp2(lon1,lat1,v1',xc(:),yc(:)); |
43 |
|
|
u2=reshape(u2,size(xc)); v2=reshape(v2,size(xc)); |
44 |
|
|
|
45 |
|
|
% compute adjacent grid coordinates for dx |
46 |
|
|
nx=1:31; ny=1:31; |
47 |
|
|
x1=xg(nx,ny,:); |
48 |
|
|
x2=xg(nx+1,ny,:); |
49 |
|
|
x3=xg(nx,ny+1,:); |
50 |
|
|
x4=xg(nx+1,ny+1,:); |
51 |
|
|
ix=find(x2-x1>180); x2(ix)=x2(ix)-360; |
52 |
|
|
ix=find(x1-x2>180); x2(ix)=x2(ix)+360; |
53 |
|
|
ix=find(x3-x1>180); x3(ix)=x3(ix)-360; |
54 |
|
|
ix=find(x1-x3>180); x2(ix)=x3(ix)+360; |
55 |
|
|
ix=find(x4-x1>180); x4(ix)=x4(ix)-360; |
56 |
|
|
ix=find(x1-x4>180); x4(ix)=x4(ix)+360; |
57 |
|
|
|
58 |
|
|
% compute adjacent grid coordinates for dy |
59 |
|
|
y1=yg(nx,ny,:); |
60 |
|
|
y2=yg(nx+1,ny,:); |
61 |
|
|
y3=yg(nx,ny+1,:); |
62 |
|
|
y4=yg(nx+1,ny+1,:); |
63 |
|
|
|
64 |
|
|
% compute cube-sphere u |
65 |
|
|
dx=0.5*(x2+x4-x1-x3).*cos(yc(nx,ny,:)*pi/180); minmax(dx) |
66 |
|
|
dy=0.5*(y2+y4-y1-y3); minmax(dy) |
67 |
|
|
denom=1./sqrt(dx.*dx+dy.*dy); |
68 |
|
|
u3=(u2(nx,ny,:).*dx+v2(nx,ny,:).*dy).*denom; |
69 |
|
|
|
70 |
|
|
% compute cube-sphere v |
71 |
|
|
dx=0.5*(x3+x4-x1-x2).*cos(yc(nx,ny,:)*pi/180); minmax(dx) |
72 |
|
|
dy=0.5*(y3+y4-y1-y2); minmax(dy) |
73 |
|
|
denom=1./sqrt(dx.*dx+dy.*dy); |
74 |
|
|
v3=(u2(nx,ny,:).*dx+v2(nx,ny,:).*dy).*denom; |
75 |
|
|
|
76 |
|
|
cx=[-10 10]; |
77 |
|
|
figure(1), tmp=u3; plot_cube, title('u from matlab') |
78 |
|
|
figure(2), tmp=v3; plot_cube, title('v from matlab') |
79 |
|
|
|
80 |
|
|
% load global_ocean.cs32x15 input files |
81 |
|
|
un=permute(readbin('../inp_thsice/ncep_uwind_cs.bin',[32,6,32],1,'real*8'),[1 3 2]); |
82 |
|
|
figure(3), tmp=un; cx=[min(tmp(:)) max(tmp(:))]; plot_cube |
83 |
|
|
title('u from ocean.cs32x15') |
84 |
dimitri |
1.2 |
print -djpeg u_cs32 |
85 |
dimitri |
1.1 |
vn=permute(readbin('../inp_thsice/ncep_vwind_cs.bin',[32,6,32],1,'real*8'),[1 3 2]); |
86 |
|
|
figure(4), tmp=vn; cx=[min(tmp(:)) max(tmp(:))]; plot_cube |
87 |
|
|
title('v from ocean.cs32x15') |
88 |
dimitri |
1.2 |
print -djpeg v_cs32 |
89 |
dimitri |
1.1 |
|
90 |
|
|
% load model output files |
91 |
|
|
um=permute(readbin('UWIND.0000000020.data',[32,6,32]),[1 3 2]); |
92 |
|
|
figure(5), tmp=um; plot_cube, title('u from model output') |
93 |
dimitri |
1.2 |
print -djpeg u_exf |
94 |
dimitri |
1.1 |
vm=permute(readbin('VWIND.0000000020.data',[32,6,32]),[1 3 2]); |
95 |
|
|
figure(6), tmp=vm; plot_cube, title('v from model output') |
96 |
dimitri |
1.2 |
print -djpeg v_exf |