/[MITgcm]/MITgcm_contrib/gael/matlab_class/sample_analysis/example_MITprof.m
ViewVC logotype

Annotation of /MITgcm_contrib/gael/matlab_class/sample_analysis/example_MITprof.m

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


Revision 1.2 - (hide annotations) (download)
Tue Apr 4 15:43:27 2017 UTC (8 years, 3 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint66o, checkpoint66f, HEAD
Changes since 1.1: +4 -0 lines
- example_MITprof.m: also store msk2D and msk3D in myprofmyop

1 gforget 1.1 function []=example_MITprof();
2     % EXAMPLE_MITPROF illustrates the use of MITprof datasets by computing
3     % Argo stats over a regional domain that includes the ACC core
4    
5     gcmfaces_global;
6    
7     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
8     if myenv.verbose>0;
9     gcmfaces_msg('===============================================');
10     gcmfaces_msg(['*** entering example_MITprof: illustrate ' ...
11     'the use of MITprof datasets by computing Argo ACC stats'],'');
12     end;
13    
14     %%%%%%%%%%%%%%%%%
15     %load grid and setup test case:
16     %%%%%%%%%%%%%%%%%
17    
18     if isempty(mygrid);
19     grid_load;
20     end;
21    
22     myenv.nctilesdir=fullfile('release2_climatology',filesep,'nctiles_climatology',filesep);
23     myenv.mitprofdir=fullfile('release2_climatology',filesep,'profiles',filesep);
24    
25     if ~isdir(myenv.nctilesdir);
26     diags=[];
27     help gcmfaces_demo;
28     warning(['skipping example_MITprof (missing ' myenv.nctilesdir ')']);
29     return;
30     end;
31    
32     if ~isdir(myenv.mitprofdir);
33     diags=[];
34     help gcmfaces_demo;
35     warning(['skipping example_MITprof (missing ' myenv.mitprofdir ')']);
36     return;
37     end;
38    
39     %%%%%%%%%%%%%%%%%%%%%%%
40     %define regional mask
41    
42     if myenv.verbose>0;
43     gcmfaces_msg('* call calc_barostream and define regional mask accordingly');
44     end;
45    
46     %compute mean streamfunction to define mask:
47     fldU=mygrid.mskW.*mean(read_nctiles([myenv.nctilesdir 'UVELMASS/UVELMASS']),4);
48     fldV=mygrid.mskS.*mean(read_nctiles([myenv.nctilesdir 'VVELMASS/VVELMASS']),4);
49     [fldBAR]=calc_barostream(fldU,fldV);
50    
51     %define 2D mask:
52     msk2D=fldBAR;
53     msk2D(fldBAR>100)=0; msk2D(fldBAR<30)=0; msk2D(mygrid.YC>0)=0;
54     msk2D(msk2D>0)=1; msk2D=msk2D.*mygrid.mskC(:,:,1);
55    
56     %max climatological MLD could aternatively be used to define mask:
57     % maxMLD=mygrid.mskC(:,:,1).*max(read_nctiles('release2_climatology/nctiles_climatology/MXLDEPTH/MXLDEPTH'),[],3);
58    
59     %define 3D mask:
60     msk3D=repmat(msk2D,[1 1 50]).*mygrid.mskC;
61     msk3D(:,:,26:50)=0.*msk3D(:,:,26:50);
62    
63    
64     %%%%%%%%%%%%%%%%%%%%%%%
65     %read in MITprof data sets and matching climatology profiles
66    
67     if myenv.verbose>0;
68     gcmfaces_msg('* read-in and select Argo profiles');
69     end;
70    
71     listData={'MITprof_mar2016_argo9506*','MITprof_mar2016_argo0708*','MITprof_mar2016_argo0910*',...
72     'MITprof_mar2016_argo1112*','MITprof_mar2016_argo1314*','MITprof_mar2016_argo1515*'};
73     listVar={'prof_T'}; listV={'T'};
74    
75     for vv=1:length(listVar);
76     tmpprof=MITprof_stats_load(myenv.mitprofdir,listData,listV{vv},listVar{vv});
77     %
78     loc_tile=gcmfaces_loc_tile(90,90,tmpprof.prof_lon,tmpprof.prof_lat);
79     tmpprof.prof_msk=convert2array(msk2D);
80     tmpprof.prof_msk=tmpprof.prof_msk(loc_tile.point);
81     tmpprof=MITprof_subset(tmpprof,'msk',1);
82     %
83     if vv==1; myprof=tmpprof; end;
84     myprof.(listVar{vv})=tmpprof.prof;
85     %
86     clear tmpprof;
87     end;
88     myprof=rmfield(myprof,{'prof','prof_msk'});
89    
90     if myenv.verbose>0;
91     gcmfaces_msg('* subsample monthly climatology to profile locations');
92     end;
93    
94     fldIn.fil=fullfile(myenv.nctilesdir,filesep,'THETA');
95     fldIn.tim='monclim'; fldIn.name='prof_Tclim';
96     myprof=MITprof_resample(myprof,fldIn);
97    
98    
99     %%%%%%%%%%%%%%%%%%%%%%%
100     %compute monthly averages
101    
102     if myenv.verbose>0;
103     gcmfaces_msg('* call MITprof_wrapper to compute monthly climatology');
104     end;
105    
106     %specify operation to compute via MITprof_wrapper
107     myop=[];
108     %myop.op_name='mean';
109     myop.op_name='cycle'; myop.op_tim=[0:7:365];
110     myop.op_vars={'prof_T','prof_Tclim'};
111    
112     %compute cycle via MITprof_wrapper
113     [cy_T,cy_Tclim]=MITprof_wrapper(myprof,myop);
114    
115     if myenv.verbose>0;
116     gcmfaces_msg('* call MITprof_wrapper via bootstrp to compute std');
117     end;
118    
119     %store myprof and myop into global variable
120     global myprofmyop; myprofmyop=myprof;
121     for ii=fieldnames(myop)'; myprofmyop.(ii{1})=myop.(ii{1}); end;
122    
123 gforget 1.2 %store mask since it could be useful for plotting or other purposes
124     myprofmyop.msk2D=msk2D;
125     myprofmyop.msk3D=msk3D;
126    
127 gforget 1.1 %compute boostrap samples using MITprof_wrapper
128     K=[1:myprofmyop.np]; N=100;
129     myout = bootstrp(N, @MITprof_wrapper,K);
130    
131     %compute standard deviation
132     bo_cy_T=reshape(cell2mat(myout(:,1)'),[size(cy_T) N]);
133     bo_cy_Tclim=reshape(cell2mat(myout(:,2)'),[size(cy_T) N]);
134     std_cy_T=std(bo_cy_T,[],3);
135     std_cy_Tclim=std(bo_cy_Tclim,[],3);
136    
137     %%%%%%%%%%%%%%%%%%%%%%%
138     %display results
139    
140     figure;
141     subplot(2,2,1); imagescnan(cy_T'); colorbar; title('mean using Argo');
142     subplot(2,2,2); imagescnan(cy_Tclim'); colorbar; title('mean using climatology');
143     subplot(2,2,3); imagescnan(std_cy_T'); colorbar; title('std using Argo');
144     subplot(2,2,4); imagescnan(std_cy_Tclim'); colorbar; title('std using climatology');
145    
146     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
147     if myenv.verbose>0;
148     gcmfaces_msg('*** leaving example_MITprof');
149     gcmfaces_msg('===============================================');
150     end;
151    

  ViewVC Help
Powered by ViewVC 1.1.22