1 |
edhill |
1.1 |
function ncep2cube |
2 |
|
|
% Interpolates NCEP data onto cubed sphere |
3 |
|
|
% Edit this script with path of NCEP data |
4 |
|
|
|
5 |
|
|
% Location of NCEP monthly climatology files |
6 |
|
|
archivedir='/u/u2/cheisey/stephd/ncep'; |
7 |
|
|
disp(['Reading NCEP data from ',archivedir]); |
8 |
|
|
|
9 |
|
|
% data path for grids |
10 |
|
|
path='./'; |
11 |
|
|
|
12 |
|
|
|
13 |
|
|
% Load grids |
14 |
|
|
load FMT.mat |
15 |
|
|
load HGRID.mat nxc nyc lon_lo lon_hi lat_lo lat_hi xc yc zA GRID yg |
16 |
|
|
load MASKS |
17 |
|
|
if ndims(msk)==4 |
18 |
|
|
msk=msk(:,:,:,1); |
19 |
|
|
else |
20 |
|
|
msk=msk(:,:,1); |
21 |
|
|
end |
22 |
|
|
|
23 |
|
|
% file list - ncep data |
24 |
|
|
flist={'ncep_downsolar_monthly.cdf',... |
25 |
|
|
'ncep_precip_monthly.cdf',... |
26 |
|
|
'ncep_tair_monthly.cdf',... |
27 |
|
|
'ncep_downlw_monthly.cdf',... |
28 |
|
|
'ncep_qair_monthly.cdf',... |
29 |
|
|
'ncep_runoff_monthly.cdf'}; |
30 |
|
|
|
31 |
|
|
% Label for plots |
32 |
|
|
lbl={'solar - incoming solar radiation',... |
33 |
|
|
'rain - precipitation',... |
34 |
|
|
'tair - air temperature', ... |
35 |
|
|
'flw - downward longwave flux',... |
36 |
|
|
'qair - specific humidity',... |
37 |
|
|
'runoff'}; |
38 |
|
|
|
39 |
|
|
% Name of corresponding variable in pkg/therm_seaice |
40 |
|
|
nm={'solr', 'prate', 'temp', 'lwflx', 'qa','runoff'}; |
41 |
|
|
|
42 |
|
|
% Scale Factors |
43 |
|
|
% Note reverse sign convention for some fields, see pkg/bulk_forcing/BULKF.h |
44 |
|
|
scaleFac=[-1, -1.0E-3, 1., -1., 1.0, 1.0]; |
45 |
|
|
|
46 |
|
|
|
47 |
|
|
% delete file for plots |
48 |
|
|
delete ncep2cube.ps |
49 |
|
|
|
50 |
|
|
for ifile=1:length(flist) |
51 |
|
|
if ifile > 1 |
52 |
|
|
clear Qshift; |
53 |
|
|
clear xs; |
54 |
|
|
end |
55 |
|
|
filename=flist{ifile}; |
56 |
|
|
disp(filename); |
57 |
|
|
ncepFile=fullfile(archivedir,filename); |
58 |
|
|
ncload(ncepFile ); |
59 |
|
|
|
60 |
|
|
% Rotate data in x so that it goes from -180 to 180 |
61 |
|
|
% instead of 0 to 360 |
62 |
|
|
xs=X-X(97); |
63 |
|
|
xs(end+1)=(xs(end)-xs(end-1))+xs(end); |
64 |
|
|
|
65 |
|
|
ys=Y; |
66 |
|
|
Qshi = permute(sq(eval(nm{ifile})),[3 2 1]); |
67 |
|
|
Qshift(97:192,:,:)=Qshi(1:96, :,:); |
68 |
|
|
Qshift(1:96, :,:)=Qshi(97:192, :,:); |
69 |
|
|
|
70 |
|
|
|
71 |
|
|
%Wrap one bin around |
72 |
|
|
Qshift(end+1, :,:) = Qshift(1,:,:); |
73 |
|
|
|
74 |
|
|
|
75 |
|
|
% Interpolate data |
76 |
|
|
Q=redo(Qshift, xs,ys,yc,xc,msk); |
77 |
|
|
Q=Q*scaleFac(ifile); |
78 |
|
|
|
79 |
|
|
% Permute indcies for cubed sphere 32x32x6 becomes 32x6x32 |
80 |
|
|
Qpermute=permute(Q,[1 3 2 4 5]); |
81 |
|
|
|
82 |
|
|
% write output file |
83 |
|
|
fout=[filename(1:end-12),'_cubed.bin']; |
84 |
|
|
wrda(fout,Qpermute,1,fmt,Ieee); |
85 |
|
|
|
86 |
|
|
% display min and max |
87 |
|
|
mx=max(max(max(max(Qpermute)))); |
88 |
|
|
mn=min(min(min(min(Qpermute)))); |
89 |
|
|
disp([' min = ',num2str(mn),' max = ',num2str(mx)]) |
90 |
|
|
|
91 |
|
|
% plot data |
92 |
|
|
if 1==1 |
93 |
|
|
|
94 |
|
|
figure(ifile) |
95 |
|
|
imonth = 3; |
96 |
|
|
subplot(2,1,1) |
97 |
|
|
|
98 |
|
|
imagesc(xs,ys,sq(Qshift(:,:,imonth))'); |
99 |
|
|
set(gca,'yDir','normal'); |
100 |
|
|
title([strrep(filename,'_','\_'),' month=',num2str(imonth)]); |
101 |
|
|
colorbar('horiz') |
102 |
|
|
|
103 |
|
|
xcs=rdmds(fullfile(path,'XC')); |
104 |
|
|
ycs=rdmds(fullfile(path,'YC')); |
105 |
|
|
xcg=rdmds(fullfile(path,'XG')); |
106 |
|
|
ycg=rdmds(fullfile(path,'YG')); |
107 |
|
|
ny=32; |
108 |
|
|
for n=1:6, |
109 |
|
|
v2d([(n-1)*ny+1:n*ny],[1:ny])=Q([1:ny],[1:ny],n, imonth); |
110 |
|
|
end |
111 |
|
|
|
112 |
|
|
subplot(2,1,2) |
113 |
|
|
shift=-1; |
114 |
|
|
c1=0; |
115 |
|
|
c2=0; |
116 |
|
|
grph_CS(sq(v2d),xcs,ycs,xcg,ycg,c1,c2,shift) |
117 |
|
|
h=title([strrep(fout,'_','\_'),' ',lbl{ifile}]); |
118 |
|
|
print -dpsc2 -append ncep2cube.ps |
119 |
|
|
end |
120 |
|
|
|
121 |
|
|
|
122 |
|
|
end |
123 |
|
|
return |
124 |
|
|
|
125 |
|
|
function Q=redo(Qshi, xs,ys,yc,xc,msk) |
126 |
|
|
|
127 |
|
|
|
128 |
|
|
for jj=1:12, |
129 |
|
|
for t=1:size(xc,3); |
130 |
|
|
Q(:,:,t,jj)=interp2(ys,xs,sq(Qshi(:,:,jj)),yc(:,:,t),xc(:,:,t)) .*msk(:,:,t); |
131 |
|
|
end |
132 |
|
|
end |
133 |
|
|
Q( isnan(Q) )=0; |
134 |
|
|
|
135 |
|
|
return |
136 |
|
|
|
137 |
|
|
|
138 |
|
|
|