/[MITgcm]/MITgcm_contrib/gael/profilesMatlabProcessing/profiles_devel/profiles_subgrid_dist.m
ViewVC logotype

Annotation of /MITgcm_contrib/gael/profilesMatlabProcessing/profiles_devel/profiles_subgrid_dist.m

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


Revision 1.2 - (hide annotations) (download)
Tue Apr 21 20:15:40 2015 UTC (11 years, 4 months ago) by gforget
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +0 -0 lines
FILE REMOVED
- remove old codes.

1 gforget 1.1 function [myStat]=profiles_subgrid_stats(choiceProf,choiceLevel,choiceVar,choiceGrid,subGridN);
2    
3     %note : in wod09 ctd, we will use kk=2 to 33
4     % subGridN=8;%on the cs24 grid
5     % subGridN=18;%on the llc grid
6     %next things I need to do :
7     % IO/memory store
8     % loop over subGridN
9     % proper treatment of choiceLevel
10    
11    
12     %global variables:
13     gcmfaces_global;
14     global lon lat obs point;
15     global choiceGridOld; if isempty(choiceGridOld); choiceGridOld='x'; end;
16     global choiceProfOld; if isempty(choiceProfOld); choiceProfOld='x'; end;
17     global choiceVarOld; if isempty(choiceVarOld); choiceVarOld='x'; end;
18     global choiceLevelOld; if isempty(choiceLevelOld); choiceLevelOld=0; end;
19    
20     %choice of time/depth ranges
21     RC=squeeze(rdmds('/net/weddell/raid3/gforget/ecco_v4/GRID/RC'));
22     RF=squeeze(rdmds('/net/weddell/raid3/gforget/ecco_v4/GRID/RF'));
23     depth0=-RF(1:end-1); depth0(2:end-1)=depth0(1:end-2);
24     depth1=-RF(2:end); depth1(2:end-1)=depth1(3:end);
25     if choiceLevel>=1;
26     depth0=depth0(choiceLevel); depth1=depth1(choiceLevel);%choice of depth range
27     else;
28     error('not implemented');
29     end;
30     %
31     date0=datenum(1950,1,1); date1=datenum(2049,12,31);%choice of time range
32    
33     %get the grid:
34     if strcmp(choiceGrid,'cs24')&~strcmp(choiceGrid,choiceGridOld);
35     dirGrid='/net/weddell/raid3/gforget/grids/gridCompleted/cube_FM/';
36     % dirGrid='/Users/gforget/mywork/projects_inprogress/2012mayInputs/insitu/processed/';
37     grid_load_native([dirGrid 'cube_24/'],6);
38     gcmfaces_bindata;
39     elseif strcmp(choiceGrid,'cs96')&~strcmp(choiceGrid,choiceGridOld);
40     dirGrid='/net/weddell/raid3/gforget/grids/gridCompleted/cube_FM/';
41     % dirGrid='/Users/gforget/mywork/projects_inprogress/2012mayInputs/insitu/processed/';
42     grid_load_native([dirGrid 'cube_96/'],6);
43     gcmfaces_bindata;
44     elseif strcmp(choiceGrid,'v4')&~strcmp(choiceGrid,choiceGridOld);
45     dirGrid='/net/nares/raid10/gforget/2012julyIters/GRID/';
46     %dirGrid='/net/weddell/raid3/gforget/ecco_v4/GRID/';
47     grid_load(dirGrid,5,'compact');
48     gcmfaces_bindata;
49     end;
50     %
51     choiceGridOld=choiceGrid;
52    
53     %get the data:
54     if strcmp(choiceProf,'argo');
55     dirData='/net/nares/raid11/ecco-shared/ecco-version-4/input/input_insitu/';
56     listData={'argo_june2012_1992_to_2007*','argo_june2012_2008_to_2010*','argo_june2012_2011_to_2012*'};
57     suffOut='argo';
58     elseif strcmp(choiceProf,'model');
59     dirData='/net/nares/raid10/gforget/2012julyIters/ecco_it0003_link/mat/profiles/output/';
60     listData={'argo_june2012_1992_to_2007_model*','argo_june2012_2008_to_2010_model*','argo_june2012_2011_to_2012_model*'};
61     suffOut='model';
62     else;
63     dirData='/net/nares/raid11/gforget/2012mayInputs/insitu/processed/';
64     listData=dir([dirData '*.nc']);
65     for ii=1:length(listData); listData(ii).name=[listData(ii).name(1:end-3) '*']; end;
66     listData={listData(:).name};
67     %for ~ backward compatibility:
68     %listRm={'itp_MITprof*','bobbers_MITprof*','CLIMODE_Talley_ctd*','WOD09_XBT*'};
69     %remove climode and xbts
70     listRm={'bobbers_MITprof*','CLIMODE_Talley_ctd*','WOD09_XBT*'};
71     for ii=1:length(listRm); tmp1=find(~strcmp(listData,listRm{ii})); listData={listData{tmp1}}; end;
72     suffOut='';
73     end;
74     %
75     test1=~strcmp(choiceProf,choiceProfOld)|...
76     ~strcmp(choiceVar,choiceVarOld)|...
77     choiceLevel~=choiceLevelOld;
78     if test1;
79     [MITprof]=MITprof_stats_load(dirData,listData,choiceVar,['prof_' choiceVar suffOut]);
80     %in DRHODR case, switch to log10 :
81     if strcmp(choiceVar,'DRHODR'); MITprof.prof=log10(MITprof.prof); end;
82     %mask out values that are not in year range:
83     ii=find(MITprof.prof_date<date0|MITprof.prof_date>date1);
84     MITprof.prof(ii,:)=NaN;
85     %restrict to depth range of interest:
86     kk=find(MITprof.prof_depth>=depth0&MITprof.prof_depth<=depth1);
87     lon=MITprof.prof_lon; lat=MITprof.prof_lat; obs=MITprof.prof(:,kk);
88     end;
89     %
90     choiceLevelOld=choiceLevel;
91     choiceVarOld=choiceVar;
92     choiceProfOld=choiceProf;
93    
94     %get indices in full grid: (not necessarily the ecco_v4 one)
95     point=gcmfaces_bindata(lon,lat);
96    
97     if strcmp(choiceVar,'RHOP');
98     binE=[1015:0.2:1030];
99     else;
100     error('need to specify binE');
101     end;
102     nB=length(binE);
103    
104     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
105    
106     %check nearest neighbor mapping:
107     % [tmp1,tmp2]=gcmfaces_bindata(lon,lat,lon);
108     % figureL; qwckplot(log10(tmp2));
109     % figureL; qwckplot(tmp1./tmp2);
110    
111     %map of full grid indices: (consistent with prof_point2)
112     indGrid=convert2array(mygrid.XC);
113     indGrid(:)=[1:length(indGrid(:))];
114     indGrid=convert2array(indGrid);
115    
116     %reduce grid and map reduced grid indices:
117     indBox=mygrid.XC;
118     if subGridN==0;
119     %global computation
120     indBox(:)=1;
121     else;
122     %regional computation
123     boxMax=0;
124     for iFace=1:mygrid.nFaces;
125     tmp1=indBox{iFace};
126     tmp3=ceil([1:size(tmp1,1)]'/subGridN);
127     tmp4=(ceil([1:size(tmp1,2)]/subGridN)-1)*max(tmp3);
128     tmp5=tmp3*ones(1,size(tmp1,2))+ones(size(tmp1,1),1)*tmp4;
129     indBox{iFace}=tmp5+boxMax;
130     boxMax=boxMax+max(tmp5(:));
131     end;
132     end;
133    
134     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
135    
136     %main computational loop: stats for each region in indBox
137     %
138     indBox=convert2array(indBox);
139     box=repmat(indBox(point),[1 size(obs,2)]);
140     box(isnan(obs))=NaN;
141     boxList=unique(indBox(:));
142     boxList=boxList(find(~isnan(boxList)));
143     %
144     myStat.nb=zeros(360*360,1);
145     myStat.dist=zeros(360*360,nB);
146     %
147     for ii=boxList';
148     if mod(ii,1000)==0;
149     [choiceLevel ii length(boxList)]
150     end;
151     jj=find(box==ii);
152     if length(jj)>=10;
153     n=histc(obs(jj),binE);
154     jj=find(indBox==ii);
155     myStat.nb(jj)=sum(n);
156     myStat.dist(jj,:)=ones(length(jj),1)*n'/sum(n);
157     end;
158     end;
159     %
160     myStat.nb=reshape(myStat.nb,[360 360]);
161     myStat.dist=reshape(myStat.dist,[360 360 nB]);
162     myStat.nb=convert2array(myStat.nb);
163     myStat.dist=convert2array(myStat.dist);
164    

  ViewVC Help
Powered by ViewVC 1.1.22