| 1 |
function []=example_budget(); |
| 2 |
%EXAMPLE_TRANSPORTS illustrates ECCO v4 budgets |
| 3 |
% |
| 4 |
%stand-alone call: addpath gcmfaces/sample_analysis/; example_budget; |
| 5 |
% |
| 6 |
%needed input files: |
| 7 |
% wget --recursive ftp://mit.ecco-group.org/gforget/nctiles_budget_2d |
| 8 |
% mv mit.ecco-group.org/gforget/nctiles_budget_2d sample_input/. |
| 9 |
% rm -rf mit.ecco-group.org |
| 10 |
% |
| 11 |
%note: in the nctiles_budget_2d files from 03-Feb-2015 the sign convention |
| 12 |
% for trWtop and trWbot was positive downward; this sign convention was |
| 13 |
% reversed in later versions for consistency with MITgcm convention. |
| 14 |
|
| 15 |
gcmfaces_global; |
| 16 |
|
| 17 |
if myenv.verbose>0; |
| 18 |
gcmfaces_msg('==============================================='); |
| 19 |
gcmfaces_msg(['*** entering example_budget ' ... |
| 20 |
'load budget terms from nctiles files ' ... |
| 21 |
'and compute hemispheric budgets'],''); |
| 22 |
end; |
| 23 |
|
| 24 |
%%%%%%%%%%%%%%%%% |
| 25 |
%load grid: |
| 26 |
%%%%%%%%%%%%%%%%% |
| 27 |
|
| 28 |
%expected location: |
| 29 |
myenv.nctilesdir=fullfile('sample_input',filesep,'nctiles_budget_2d',filesep); |
| 30 |
%if nctiles_budget_2d is not found then try old location: |
| 31 |
if ~isdir(myenv.nctilesdir); |
| 32 |
%if not found then try old location: |
| 33 |
tmpdir=fullfile(myenv.gcmfaces_dir,'/sample_input/nctiles_budget_2d/'); |
| 34 |
if isdir(tmpdir); myenv.nctilesdir=tmpdir; end; |
| 35 |
end; |
| 36 |
%if nctiles_budget_2d is still not found then issue warning and skip example_budget |
| 37 |
if ~isdir(myenv.nctilesdir); |
| 38 |
diags=[]; |
| 39 |
help example_budget; |
| 40 |
warning(['skipping example_budget (missing ' myenv.nctilesdir ')']); |
| 41 |
return; |
| 42 |
end; |
| 43 |
|
| 44 |
dirName=fullfile(myenv.nctilesdir,filesep,'budgMo',filesep); |
| 45 |
|
| 46 |
if ~isdir(dirName); |
| 47 |
help gcmfaces_demo; |
| 48 |
warning(['skipping example_budget (missing ' dirName ')']); |
| 49 |
return; |
| 50 |
end; |
| 51 |
|
| 52 |
if isempty(which('ncload')); |
| 53 |
warning(['skipping example_budget (missing ncload that is part of MITprof)']); |
| 54 |
return; |
| 55 |
end; |
| 56 |
|
| 57 |
%%%%%%%%%%%%%%%%%% |
| 58 |
%main computation: |
| 59 |
%%%%%%%%%%%%%%%%%% |
| 60 |
|
| 61 |
%load grid: |
| 62 |
if isempty(mygrid); |
| 63 |
grid_load; |
| 64 |
end; |
| 65 |
|
| 66 |
%select budget of interest |
| 67 |
nameBudg='budgMo'; |
| 68 |
|
| 69 |
%load budget terms |
| 70 |
fileName=[myenv.nctilesdir nameBudg filesep nameBudg(1:6)]; |
| 71 |
tend=read_nctiles(fileName,'tend'); |
| 72 |
trU=read_nctiles(fileName,'trU'); |
| 73 |
trV=read_nctiles(fileName,'trV'); |
| 74 |
trWtop=read_nctiles(fileName,'trWtop'); |
| 75 |
trWbot=read_nctiles(fileName,'trWbot'); |
| 76 |
|
| 77 |
%load dt (time increments) vector (not a gcmfaces object) |
| 78 |
ncload([fileName '.0001.nc'],'dt'); |
| 79 |
|
| 80 |
%get budget descitption and units |
| 81 |
ncid = netcdf.open([fileName '.0001.nc'],'NC_NOWRITE'); |
| 82 |
descr = netcdf.getAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'description'); |
| 83 |
varid = netcdf.inqVarID(ncid,'tend'); |
| 84 |
units = netcdf.getAtt(ncid,varid,'units'); |
| 85 |
netcdf.close(ncid); |
| 86 |
|
| 87 |
if myenv.verbose>0; gcmfaces_msg(['* display Northern Hemisphere budget for ' nameBudg ' (' descr ')']);end; |
| 88 |
|
| 89 |
%define northern hemisphere as domain of integration |
| 90 |
nameMask='Northern Hemisphere'; |
| 91 |
mask=mygrid.mskC(:,:,1).*(mygrid.YC>0); |
| 92 |
areaMask=mygrid.RAC.*mask; |
| 93 |
|
| 94 |
%edit plot title accordingly |
| 95 |
tmp1=strfind(descr,'-- ECCO v4'); |
| 96 |
descr=[descr(1:tmp1-1) 'for: ' nameMask]; |
| 97 |
|
| 98 |
%compute northern hemisphere integrals |
| 99 |
budg.tend=NaN*dt; |
| 100 |
budg.hconv=NaN*dt; |
| 101 |
budg.zconv=NaN*dt; |
| 102 |
for tt=1:length(dt); |
| 103 |
%compute flux convergence |
| 104 |
hconv=calc_UV_conv(trU(:,:,tt),trV(:,:,tt)); |
| 105 |
zconv=trWtop(:,:,tt)-trWbot(:,:,tt); |
| 106 |
%compute sum over domain |
| 107 |
budg.tend(tt)=nansum(tend(:,:,tt).*mask)/nansum(areaMask); |
| 108 |
budg.hconv(tt)=nansum(hconv.*mask)/nansum(areaMask); |
| 109 |
budg.zconv(tt)=nansum(zconv.*mask)/nansum(areaMask); |
| 110 |
end; |
| 111 |
|
| 112 |
%display result |
| 113 |
figureL; |
| 114 |
subplot(3,1,1); set(gca,'FontSize',12); |
| 115 |
plot(cumsum(dt.*budg.tend)); |
| 116 |
grid on; xlabel('month'); ylabel([units ' x s']); |
| 117 |
legend('content anomaly'); |
| 118 |
subplot(3,1,2); set(gca,'FontSize',12); |
| 119 |
plot(cumsum(dt.*budg.tend)); hold on; |
| 120 |
plot(cumsum(dt.*budg.hconv),'r'); |
| 121 |
plot(cumsum(dt.*budg.zconv),'g'); |
| 122 |
grid on; xlabel('month'); ylabel([units ' x s']); |
| 123 |
legend('content anomaly','horizontal convergence','vertical convergence'); |
| 124 |
subplot(3,1,3); set(gca,'FontSize',12); |
| 125 |
plot(cumsum(dt.*(budg.tend-budg.hconv-budg.zconv))); |
| 126 |
grid on; xlabel('month'); ylabel([units ' x s']); |
| 127 |
legend('budget residual'); |
| 128 |
|
| 129 |
if myenv.verbose>0; |
| 130 |
gcmfaces_msg('*** leaving example_budget ',''); |
| 131 |
gcmfaces_msg('==============================================='); |
| 132 |
end; |
| 133 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph