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