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

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

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


Revision 1.1 - (hide annotations) (download)
Tue Jan 20 13:31:44 2015 UTC (10 years, 6 months ago) by gforget
Branch: MAIN
- illustrates budget analysis

1 gforget 1.1 function []=example_budget(choiceV3orV4);
2     %object: illustrates budget analysis
3     %inputs: choiceV3orV4 ('v3' or 'v4') selects the sample GRID
4    
5     gcmfaces_global;
6     if myenv.verbose>0;
7     gcmfaces_msg('===============================================');
8     gcmfaces_msg(['*** entering example_budget ' ...
9     'that will load a set of budget terms from file ' ...
10     'and compute a series of derived diagnostics '],'');
11     end;
12    
13     %%%%%%%%%%%%%%%%%
14     %load parameters:
15     %%%%%%%%%%%%%%%%%
16    
17     if myenv.verbose>0;
18     gcmfaces_msg('* set grid files path and format, and number of faces for this grid');
19     end;
20     myenv.nctiles=1;
21     if strcmp(choiceV3orV4,'v4');
22     nF=5;
23     fileFormat='compact';
24     dir0=fullfile([myenv.gcmfaces_dir '../']);
25     dirGrid=fullfile(dir0,'GRID/');
26     myenv.nctilesdir=fullfile(dir0,'release1/nctiles_budget/');
27     else;
28     error('files are missing');
29     nF=1;
30     dir0=fullfile(myenv.gcmfaces_dir,'sample_input/');
31     fileFormat='straight';
32     dirGrid=fullfile(dir0,'grid_occa/');
33     myenv.nctilesdir=fullfile(dir0,'nctiles_occa/');
34     end;
35     if myenv.verbose>0;
36     gcmfaces_msg('* call grid_load : load grid to memory (mygrid) according to');
37     gcmfaces_msg(['dirGrid = ' dirGrid],' ');
38     gcmfaces_msg(['nFaces = ' num2str(nF)],' ');
39     gcmfaces_msg(['fileFormat = ' fileFormat],' ');
40     % fprintf([' > dirGrid = ' dirGrid '\n']);
41     end;
42    
43     %%%%%%%%%%%%%%%%%%
44     %main computation:
45     %%%%%%%%%%%%%%%%%%
46    
47     %select budget of interest
48     nameBudg='budgMo';
49    
50     %load budget terms
51     fileName=[myenv.nctilesdir nameBudg filesep nameBudg(1:6)];
52     tend=read_nctiles(fileName,'tend');
53     trU=read_nctiles(fileName,'trU');
54     trV=read_nctiles(fileName,'trV');
55     trWtop=read_nctiles(fileName,'trWtop');
56     trWbot=read_nctiles(fileName,'trWbot');
57    
58     %load dt (time increments) vector (not a gcmfaces object)
59     ncload([fileName '.0001.nc'],'dt');
60    
61     %get budget descitption and units
62     ncid = netcdf.open([fileName '.0001.nc'],'NC_NOWRITE');
63     descr = netcdf.getAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'description');
64     varid = netcdf.inqVarID(ncid,'tend');
65     units = netcdf.getAtt(ncid,varid,'units');
66     netcdf.close(ncid);
67    
68     %define northern hemisphere as domain of integration
69     nameMask='Northern Hemisphere';
70     mask=mygrid.mskC(:,:,1).*(mygrid.YC>0);
71     areaMask=mygrid.RAC.*mask;
72    
73     %edit plot title accordingly
74     tmp1=strfind(descr,'-- ECCO v4');
75     descr=[descr(1:tmp1-1) 'for: ' nameMask];
76    
77     %compute northern hemisphere integrals
78     budg.tend=NaN*dt;
79     budg.hconv=NaN*dt;
80     budg.zconv=NaN*dt;
81     for tt=1:length(dt);
82     %compute flux convergence
83     hconv=calc_UV_conv(trU(:,:,tt),trV(:,:,tt));
84     zconv=trWtop(:,:,tt)-trWbot(:,:,tt);
85     %compute sum over domain
86     budg.tend(tt)=nansum(tend(:,:,tt).*mask)/nansum(areaMask);
87     budg.hconv(tt)=nansum(hconv.*mask)/nansum(areaMask);
88     budg.zconv(tt)=nansum(zconv.*mask)/nansum(areaMask);
89     end;
90    
91     %display result
92     figureL;
93     subplot(3,1,1); set(gca,'FontSize',12);
94     plot(cumsum(dt.*budg.tend));
95     grid on; xlabel('month'); ylabel([units '.s']);
96     legend('content anomaly'); title(descr);
97     subplot(3,1,2); set(gca,'FontSize',12);
98     plot(cumsum(dt.*budg.tend)); hold on;
99     plot(cumsum(dt.*budg.hconv),'r');
100     plot(cumsum(dt.*budg.zconv),'g');
101     grid on; xlabel('month'); ylabel([units '.s']);
102     legend('content anomaly','horizontal convergence','vertical convergence');
103     subplot(3,1,3); set(gca,'FontSize',12);
104     plot(cumsum(dt.*(budg.tend-budg.hconv-budg.zconv)));
105     grid on; xlabel('month'); ylabel([units '.s']);
106     legend('budget residual');
107    

  ViewVC Help
Powered by ViewVC 1.1.22