/[MITgcm]/MITgcm_contrib/high_res_cube/eddy_flux/ncep2cube.m
ViewVC logotype

Annotation of /MITgcm_contrib/high_res_cube/eddy_flux/ncep2cube.m

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


Revision 1.2 - (hide annotations) (download)
Tue Aug 3 15:54:26 2004 UTC (20 years, 11 months ago) by edhill
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +0 -0 lines
FILE REMOVED
 o sanitizing and checking in parts of the eddy flux diagnostics

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    

  ViewVC Help
Powered by ViewVC 1.1.22