/[MITgcm]/MITgcm_contrib/high_res_cube/input/lookat_uv.m
ViewVC logotype

Annotation of /MITgcm_contrib/high_res_cube/input/lookat_uv.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.2 - (hide annotations) (download)
Sat Jan 10 00:50:25 2004 UTC (21 years, 6 months ago) by dimitri
Branch: MAIN
CVS Tags: hrcube4, hrcube5, hrcube_1, hrcube_2, hrcube_3, HEAD
Changes since 1.1: +6 -2 lines
checking in code for first multi-year 510x510x6 integration

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

  ViewVC Help
Powered by ViewVC 1.1.22