matlab_class/sample_analysis/example_transports.m

function [diags]=example_transports();
%EXAMPLE_TRANSPORTS computes transports (and zonal averages)
%
% Note: this demonstration routine includes its own call to grid_load
%
% Example:
%          addpath gcmfaces/; gcmfaces_global; myenv.verbose=1;
%          [diags]=example_transports();
%          example_transports_disp(diags);

gcmfaces_global;
if myenv.verbose>0;
    gcmfaces_msg('===============================================');
    gcmfaces_msg(['*** entering example_transports: ' ...
        'load flow field form nctiles file and compute transports'],'');
end;

%%%%%%%%%%%%%%%%%
%load grid:
%%%%%%%%%%%%%%%%%

if isempty(mygrid);
   grid_load;
end;

myenv.nctilesdir=fullfile('release1',filesep,'nctiles_climatology',filesep);

if ~isdir(myenv.nctilesdir);
  myenv.nctilesdir=fullfile('release2_climatology',filesep,'nctiles_climatology',filesep);
end;

if ~isdir(myenv.nctilesdir);
    diags=[];
    help gcmfaces_demo;
    warning(['skipping example_transports (missing ' myenv.nctilesdir ')']);
    return;
end;

if myenv.verbose>0;
    gcmfaces_msg('* call gcmfaces_lines_zonal : determine grid lines that closely follow lines of');
    gcmfaces_msg('constant latitudes and will be used in zonal mean and meridional transport computations','  ');
end;
gcmfaces_lines_zonal;

if myenv.verbose>0;
    gcmfaces_msg('* call gcmfaces_lines_transp : determine grid lines that closely follow');
    gcmfaces_msg('great circles and will be used to compute transports through select sections','  ');
end;
% warning('skipping gcmfaces_lines_transp\n');
[lonPairs,latPairs,names]=gcmfaces_lines_pairs;
gcmfaces_lines_transp(lonPairs,latPairs,names);

%%%%%%%%%%%%%%%%%
%do computations:
%%%%%%%%%%%%%%%%%

listDiags={};

% [listTimes]=diags_list_times;
diags.listTimes=1;

%part 1:

listVars={'UVELMASS','VVELMASS','GM_PsiX','GM_PsiY'};
missingVars={};
for vv=1:length(listVars);
  tmp1=[myenv.nctilesdir listVars{vv} filesep listVars{vv} '*nc'];
  if isempty(dir(tmp1)); missingVars={missingVars{:},tmp1}; end;
end;

if ~isempty(missingVars);
    fprintf('\n example_transports could not find the following files ---> skipping related computation!\n');
    disp(missingVars');
else;

if myenv.verbose>0; gcmfaces_msg('* call read_nctiles : load velocity GM streamfunction fields');end;

for vvv=1:length(listVars);
    vv=listVars{vvv};
    tmp1=read_nctiles([myenv.nctilesdir vv '/' vv],vv);
    tmp1=mean(tmp1,4);
    tmp1(mygrid.mskC==0)=NaN;
    eval([vv '=tmp1;']);
end;

%apply NaN-masks:
UVELMASS=UVELMASS.*mygrid.mskW;
VVELMASS=VVELMASS.*mygrid.mskS;

if myenv.verbose>0; gcmfaces_msg('* call calc_bolus :  compute bolus component from GM streamfunction');end;

%add bolus component from GM scheme:
[UVELbol,VVELbol,fldWbolus]=calc_bolus(GM_PsiX,GM_PsiY);
UVELbol=UVELbol.*mygrid.mskW; VVELbol=VVELbol.*mygrid.mskS;
UVELtot=UVELMASS+UVELbol; VVELtot=VVELMASS+VVELbol;

%compute transports:
listDiags={listDiags{:},'fldBAR','gloOV','fldTRANSPORTS','gloMT_FW'};
if myenv.verbose>0; gcmfaces_msg('* call calc_barostream : comp. barotropic stream function');end;
[fldBAR]=calc_barostream(UVELMASS,VVELMASS);
if myenv.verbose>0; gcmfaces_msg('* call calc_overturn : comp. residual overturning stream function');end;
[gloOV]=calc_overturn(UVELtot,VVELtot);
if myenv.verbose>0; gcmfaces_msg('* call calc_transports : comp. transects transports');end;
[fldTRANSPORTS]=1e-6*calc_transports(UVELMASS,VVELMASS,mygrid.LINES_MASKS,{'dh','dz'});
if myenv.verbose>0; gcmfaces_msg('* call calc_MeridionalTransport : comp. meridional seawater transport');end;
[gloMT_FW]=1e-6*calc_MeridionalTransport(UVELMASS,VVELMASS,1);

end;%if ~isempty(missingVars);

%part 2:

listVars={'THETA','SALT','ADVx_TH','ADVy_TH','ADVx_SLT','ADVy_SLT'};
listVars={listVars{:},'DFxE_TH','DFyE_TH','DFxE_SLT','DFyE_SLT'};

missingVars={};
for vv=1:length(listVars);
  tmp1=[myenv.nctilesdir listVars{vv} filesep listVars{vv} '*nc'];
  if isempty(dir(tmp1)); missingVars={missingVars{:},tmp1}; end;
end;

if ~isempty(missingVars);
    fprintf('\n example_transports could not find the following files ---> skipping related computation!\n');
    disp(missingVars');
else;

if myenv.verbose>0; gcmfaces_msg('* load tracer and transports fields');end;
listDiags={listDiags{:},'fldTzonmean','fldSzonmean','gloMT_H','gloMT_SLT'};
for vvv=1:length(listVars);
    vv=listVars{vvv};
    tmp1=read_nctiles([myenv.nctilesdir vv '/' vv],vv);
    tmp1=mean(tmp1,4);
    tmp1(mygrid.mskC==0)=NaN;
    eval([vv '=tmp1;']);
end;

if myenv.verbose>0; gcmfaces_msg('* call calc_zonmean_T : comp. zonal mean temperature');end;
[fldTzonmean]=calc_zonmean_T(THETA);
if myenv.verbose>0; gcmfaces_msg('* call calc_zonmean_T : comp. zonal mean salinity');end;
[fldSzonmean]=calc_zonmean_T(SALT);

if myenv.verbose>0; gcmfaces_msg('* call calc_MeridionalTransport : comp. meridional heat transport');end;
tmpU=(ADVx_TH+DFxE_TH); tmpV=(ADVy_TH+DFyE_TH);
[gloMT_H]=1e-15*4e6*calc_MeridionalTransport(tmpU,tmpV,0);
if myenv.verbose>0; gcmfaces_msg('* call calc_MeridionalTransport : comp. meridional salt transport');end;
tmpU=(ADVx_SLT+DFxE_SLT); tmpV=(ADVy_SLT+DFyE_SLT);
[gloMT_SLT]=1e-6*calc_MeridionalTransport(tmpU,tmpV,0);

end;%if ~isempty(missingVars);

%part 3: format output

for ddd=1:length(listDiags);
    dd=listDiags{ddd};
    eval(['diags.' dd '=' dd ';']);
end;

if myenv.verbose>0;
    gcmfaces_msg('*** leaving example_transports');
    gcmfaces_msg('===============================================','');
end;

1	function [diags]=example_transports();
2	%EXAMPLE_TRANSPORTS computes transports (and zonal averages)
3	%
4	% Note: this demonstration routine includes its own call to grid_load
5	%
6	% Example:
7	% addpath gcmfaces/; gcmfaces_global; myenv.verbose=1;
8	% [diags]=example_transports();
9	% example_transports_disp(diags);
10
11	gcmfaces_global;
12	if myenv.verbose>0;
13	gcmfaces_msg('===============================================');
14	gcmfaces_msg(['*** entering example_transports: ' ...
15	'load flow field form nctiles file and compute transports'],'');
16	end;
17
18	%%%%%%%%%%%%%%%%%
19	%load grid:
20	%%%%%%%%%%%%%%%%%
21
22	if isempty(mygrid);
23	grid_load;
24	end;
25
26	myenv.nctilesdir=fullfile('release1',filesep,'nctiles_climatology',filesep);
27
28	if ~isdir(myenv.nctilesdir);
29	myenv.nctilesdir=fullfile('release2_climatology',filesep,'nctiles_climatology',filesep);
30	end;
31
32	if ~isdir(myenv.nctilesdir);
33	diags=[];
34	help gcmfaces_demo;
35	warning(['skipping example_transports (missing ' myenv.nctilesdir ')']);
36	return;
37	end;
38
39	if myenv.verbose>0;
40	gcmfaces_msg('* call gcmfaces_lines_zonal : determine grid lines that closely follow lines of');
41	gcmfaces_msg('constant latitudes and will be used in zonal mean and meridional transport computations',' ');
42	end;
43	gcmfaces_lines_zonal;
44
45	if myenv.verbose>0;
46	gcmfaces_msg('* call gcmfaces_lines_transp : determine grid lines that closely follow');
47	gcmfaces_msg('great circles and will be used to compute transports through select sections',' ');
48	end;
49	% warning('skipping gcmfaces_lines_transp\n');
50	[lonPairs,latPairs,names]=gcmfaces_lines_pairs;
51	gcmfaces_lines_transp(lonPairs,latPairs,names);
52
53	%%%%%%%%%%%%%%%%%
54	%do computations:
55	%%%%%%%%%%%%%%%%%
56
57	listDiags={};
58
59	% [listTimes]=diags_list_times;
60	diags.listTimes=1;
61
62	%part 1:
63
64	listVars={'UVELMASS','VVELMASS','GM_PsiX','GM_PsiY'};
65	missingVars={};
66	for vv=1:length(listVars);
67	tmp1=[myenv.nctilesdir listVars{vv} filesep listVars{vv} '*nc'];
68	if isempty(dir(tmp1)); missingVars={missingVars{:},tmp1}; end;
69	end;
70
71	if ~isempty(missingVars);
72	fprintf('\n example_transports could not find the following files ---> skipping related computation!\n');
73	disp(missingVars');
74	else;
75
76	if myenv.verbose>0; gcmfaces_msg('* call read_nctiles : load velocity GM streamfunction fields');end;
77
78	for vvv=1:length(listVars);
79	vv=listVars{vvv};
80	tmp1=read_nctiles([myenv.nctilesdir vv '/' vv],vv);
81	tmp1=mean(tmp1,4);
82	tmp1(mygrid.mskC==0)=NaN;
83	eval([vv '=tmp1;']);
84	end;
85
86	%apply NaN-masks:
87	UVELMASS=UVELMASS.*mygrid.mskW;
88	VVELMASS=VVELMASS.*mygrid.mskS;
89
90	if myenv.verbose>0; gcmfaces_msg('* call calc_bolus : compute bolus component from GM streamfunction');end;
91
92	%add bolus component from GM scheme:
93	[UVELbol,VVELbol,fldWbolus]=calc_bolus(GM_PsiX,GM_PsiY);
94	UVELbol=UVELbol.mygrid.mskW; VVELbol=VVELbol.mygrid.mskS;
95	UVELtot=UVELMASS+UVELbol; VVELtot=VVELMASS+VVELbol;
96
97	%compute transports:
98	listDiags={listDiags{:},'fldBAR','gloOV','fldTRANSPORTS','gloMT_FW'};
99	if myenv.verbose>0; gcmfaces_msg('* call calc_barostream : comp. barotropic stream function');end;
100	[fldBAR]=calc_barostream(UVELMASS,VVELMASS);
101	if myenv.verbose>0; gcmfaces_msg('* call calc_overturn : comp. residual overturning stream function');end;
102	[gloOV]=calc_overturn(UVELtot,VVELtot);
103	if myenv.verbose>0; gcmfaces_msg('* call calc_transports : comp. transects transports');end;
104	[fldTRANSPORTS]=1e-6*calc_transports(UVELMASS,VVELMASS,mygrid.LINES_MASKS,{'dh','dz'});
105	if myenv.verbose>0; gcmfaces_msg('* call calc_MeridionalTransport : comp. meridional seawater transport');end;
106	[gloMT_FW]=1e-6*calc_MeridionalTransport(UVELMASS,VVELMASS,1);
107
108	end;%if ~isempty(missingVars);
109
110	%part 2:
111
112	listVars={'THETA','SALT','ADVx_TH','ADVy_TH','ADVx_SLT','ADVy_SLT'};
113	listVars={listVars{:},'DFxE_TH','DFyE_TH','DFxE_SLT','DFyE_SLT'};
114
115	missingVars={};
116	for vv=1:length(listVars);
117	tmp1=[myenv.nctilesdir listVars{vv} filesep listVars{vv} '*nc'];
118	if isempty(dir(tmp1)); missingVars={missingVars{:},tmp1}; end;
119	end;
120
121	if ~isempty(missingVars);
122	fprintf('\n example_transports could not find the following files ---> skipping related computation!\n');
123	disp(missingVars');
124	else;
125
126	if myenv.verbose>0; gcmfaces_msg('* load tracer and transports fields');end;
127	listDiags={listDiags{:},'fldTzonmean','fldSzonmean','gloMT_H','gloMT_SLT'};
128	for vvv=1:length(listVars);
129	vv=listVars{vvv};
130	tmp1=read_nctiles([myenv.nctilesdir vv '/' vv],vv);
131	tmp1=mean(tmp1,4);
132	tmp1(mygrid.mskC==0)=NaN;
133	eval([vv '=tmp1;']);
134	end;
135
136	if myenv.verbose>0; gcmfaces_msg('* call calc_zonmean_T : comp. zonal mean temperature');end;
137	[fldTzonmean]=calc_zonmean_T(THETA);
138	if myenv.verbose>0; gcmfaces_msg('* call calc_zonmean_T : comp. zonal mean salinity');end;
139	[fldSzonmean]=calc_zonmean_T(SALT);
140
141	if myenv.verbose>0; gcmfaces_msg('* call calc_MeridionalTransport : comp. meridional heat transport');end;
142	tmpU=(ADVx_TH+DFxE_TH); tmpV=(ADVy_TH+DFyE_TH);
143	[gloMT_H]=1e-154e6calc_MeridionalTransport(tmpU,tmpV,0);
144	if myenv.verbose>0; gcmfaces_msg('* call calc_MeridionalTransport : comp. meridional salt transport');end;
145	tmpU=(ADVx_SLT+DFxE_SLT); tmpV=(ADVy_SLT+DFyE_SLT);
146	[gloMT_SLT]=1e-6*calc_MeridionalTransport(tmpU,tmpV,0);
147
148	end;%if ~isempty(missingVars);
149
150	%part 3: format output
151
152	for ddd=1:length(listDiags);
153	dd=listDiags{ddd};
154	eval(['diags.' dd '=' dd ';']);
155	end;
156
157	if myenv.verbose>0;
158	gcmfaces_msg('*** leaving example_transports');
159	gcmfaces_msg('===============================================','');
160	end;
161