<p><b>mpetersen@lanl.gov</b> 2012-05-22 10:18:06 -0600 (Tue, 22 May 2012)</p><p>Adding branches/tools/transport_sections. This provides the matlab code to find connecting edges along a section, and to compute the transport through a section. See README file.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/tools/mpas_draw/Makefile
===================================================================
--- branches/tools/mpas_draw/Makefile 2012-05-21 22:12:55 UTC (rev 1926)
+++ branches/tools/mpas_draw/Makefile 2012-05-22 16:18:06 UTC (rev 1927)
@@ -1,7 +1,5 @@
CXX=g++
-#CXXFLAGS= -O3
CXXFLAGS= -O3 -m64
-NETCDF=/opt/netcdf/3.6.3/64
DBGFLAGS= -g -m64
@@ -11,11 +9,11 @@
PLATFORM=_LINUX
ifeq ($(PLATFORM),_LINUX)
- CXXLIBS = -I$(NETCDF)/include -L$(NETCDF)/lib -lGL -lglut -lnetcdf_c++ -lnetcdf -lGLU -lstdc++
+ CXXLIBS = -I$(NETCDF)/include -I${GLUT}/include -L${GLUT}/lib -L$(NETCDF)/lib -lGL -lglut -lnetcdf_c++ -lnetcdf -lGLU -lstdc++
endif
ifeq ($(PLATFORM),_MACOS)
- CXXLIBS = -L$(NETCDF)/lib -lnetcdf -lnetcdf_c++ -I$(NETCDF)/include -lstdc++ -framework OpenGL -framework GLUT
+ CXXLIBS = -L$(NETCDF)/lib -lnetcdf -lnetcdf_c++ -I$(NETCDF)/include -lstdc++ -framework OpenGL -framework GLUT
endif
all:
Added: branches/tools/transport_sections/README
===================================================================
--- branches/tools/transport_sections/README (rev 0)
+++ branches/tools/transport_sections/README 2012-05-22 16:18:06 UTC (rev 1927)
@@ -0,0 +1,32 @@
+transport_sections README
+
+This is a matlab tool to find sections that connect two points on the
+globe. These sections are a sequence of connected edges, and the
+edges and other variables are output as both a netcdf and text files.
+The transport can then be measured using this matlab code using output
+files, or in MPAS-Ocean during runtime.
+
+To begin, change the parameters at the top of transport_sections.m.
+You will need to change the text strings wd, sim(i).dir, and
+sim(i).netcdf_file so that the text string
+
+[wd '/' sim(i).dir '/' sim(i).netcdf_file ] is the file path,
+
+where wd is the working directory and dir is the run directory.
+Details of the section coordinates and variables may be specified in
+transport_sections.m.
+
+The data files only need to contain a small number of variables.
+You may need to reduce the file size before copying to a local
+machine using:
+
+ncks -v acc_u, \
+nAccumulate,latVertex,lonVertex,verticesOnEdge,edgesOnVertex,hZLevel,\
+dvEdge \
+file_in.nc file_out.nc
+
+This plotting tool has been tested for ocean output files, but should
+nearly work for other cores as well. A few lines will need to be
+changed.
+
+Mark Petersen, MPAS-Ocean Team, LANL, May 2012
Added: branches/tools/transport_sections/compute_transport.m
===================================================================
--- branches/tools/transport_sections/compute_transport.m (rev 0)
+++ branches/tools/transport_sections/compute_transport.m 2012-05-22 16:18:06 UTC (rev 1927)
@@ -0,0 +1,86 @@
+function compute_transport ...
+ (wd,dir,netcdf_file, var_name, ...
+ sectionEdgeIndex, sectionEdgeSign, ...
+ nEdgesInSection, sectionData,sectionText,sectionAbbreviation)
+
+% Load large variables from netcdf file
+
+% Mark Petersen, MPAS-Ocean Team, LANL, May 2012
+
+%%%%%%%%%% input arguments %%%%%%%%%
+% The text string [wd '/' dir '/' netcdf_file ] is the file path,
+% where wd is the working directory and dir is the run directory.
+% var_name(nVars) a cell array with text for each variable to
+% load or compute.
+% var_conv_factor multiply each variable by this unit conversion.
+% sectionEdgeIndex(maxEdges,nSections) cell index of each section
+% nEdgesInSection(nSections) number of cells in each section
+% sectionData(nVertLevels,max(nEdgesInSection),nSections,nVars)
+% data in each cross-section for each variable
+% sectionText a cell array with text describing each section
+% sectionAbbreviation an 8-character title for each section
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Compute transport through each section
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+fprintf('</font>
<font color="blue">')
+fprintf(['** Compute transport: ' dir '</font>
<font color="blue">'])
+
+filename = [wd '/' dir '/' netcdf_file ];
+ncid = netcdf.open(filename,'nc_nowrite');
+
+hZLevel = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'hZLevel'));
+dvEdge = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'dvEdge'));
+[dimname,nVertLevels]= netcdf.inqDim(ncid,netcdf.inqDimID(ncid,'nVertLevels'));
+netcdf.close(ncid)
+
+nSections = length(nEdgesInSection);
+maxNEdgesInSection = max(nEdgesInSection);
+
+m3ps_to_Sv = 1e-6; % m^3/sec flux to Sverdrups
+
+% the volume transport
+tr = zeros(nVertLevels,maxNEdgesInSection,nSections);
+tr_total = zeros(nSections,1);
+header=' ';
+data_str=' ';
+
+for iSection = 1:nSections
+ for i=1:nEdgesInSection(iSection)
+ iEdge = sectionEdgeIndex(i,iSection);
+ for k=1:nVertLevels
+ % Compute transport.
+ % I am assuming here that sectionData(:,:,:,1) contains acc_u
+ tr(k,i,iSection) = sectionEdgeSign(i,iSection)...
+ *sectionData(k,i,iSection,1)*dvEdge(iEdge)* ...
+ hZLevel(k)*m3ps_to_Sv;
+ tr_total(iSection) = tr_total(iSection) + tr(k,i,iSection);
+
+ % This is edge velocity
+ tr(k,i,iSection) = sectionEdgeSign(i,iSection)*sectionData(k,i,iSection,1);
+ end
+ end
+ %figure(iSection+1)
+ %imagesc(log(abs(tr(:,1:nEdgesInSection(iSection),iSection))))
+ %imagesc(tr(:,1:nEdgesInSection(iSection),iSection))
+ %colorbar
+
+ temptext = char(sectionText(iSection));
+ fprintf(['Section %3i, ' temptext(1:20) 'observed flow:' ...
+ temptext(63:75) ' mpas flow: %20.15f Sv</font>
<font color="blue">'],iSection,tr_total(iSection))
+
+% fprintf(['Section %3i, ' temptext(1:20) 'observed flow:' ...
+% temptext(63:75) ' mpas flow: %6.1f Sv</font>
<font color="blue">'],iSection,tr_total(iSection))
+
+ header = [header sectionAbbreviation(iSection,:) ' '];
+ data_str = [data_str num2str_fixed(tr_total(iSection),'%4.1f',8)...
+ ' '];
+end
+
+fprintf(['</font>
<font color="black"> Summary, in Sv: </font>
<font color="black">' header ' Simulation: </font>
<font color="black">' data_str ' ' dir '</font>
<font color="blue">'])
+
+
+fprintf('</font>
<font color="gray">')
+
Added: branches/tools/transport_sections/find_edge_sections.m
===================================================================
--- branches/tools/transport_sections/find_edge_sections.m (rev 0)
+++ branches/tools/transport_sections/find_edge_sections.m 2012-05-22 16:18:06 UTC (rev 1927)
@@ -0,0 +1,219 @@
+function [sectionEdgeIndex, sectionEdgeSign, nEdgesInSection, ...
+ latSectionVertex,lonSectionVertex, ...
+ latVertexDeg,lonVertexDeg] = find_edge_sections ...
+ (wd,dir,netcdf_file,sectionText,sectionCoord)
+
+% This function reads grid data from an MPAS-Ocean grid or restart
+% netCDF file, and finds a path of edges that connect the endpoints
+% specified in sectionCoord. The path is forced to travel through edges
+% that are closest to the line connecting the beginning and end
+% edges.
+%
+% Mark Petersen, MPAS-Ocean Team, LANL, May 2012
+%
+%%%%%%%%%% input arguments %%%%%%%%%
+% The text string [wd '/' dir '/' netcdf_file ] is the file path,
+% where wd is the working directory and dir is the run directory.
+% sectionText a cell array with text describing each section
+% sectionCoord(nSections,4) endpoints of sections, with one section per row as
+% [startlat startlon endlat endlon]
+%
+%%%%%%%%%% output arguments %%%%%%%%%
+% sectionEdgeIndex(maxNEdgesInSection,nSections) edge index of each section
+% sectionEdgeSign(maxNEdgesInSection,nSections) sign of each
+% section, positive is to right of path direction.
+% nEdgesInSection(nSections) number of edges in each section
+% latSectionVertex(maxNEdgesInSection,nSections) lat coordinates of each section
+% lonSectionVertex(maxNEdgesInSection,nSections) lon coordinates of each section
+% latVertexDeg(nEdges) lat arrays for all edges
+% lonVertexDeg(nEdges) lon arrays for all edges
+
+%%%%%%%%%% parameters internal to this function %%%%%%%%%%
+
+% maxEdges specifies the maximum number of Edges attempted along
+% the path to the end-edge before stopping with a warning.
+maxEdges = 1500;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Read edge and edge data from grid file
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+fprintf(['** find_edge_sections, simulation: ' dir '</font>
<font color="blue">'])
+
+filename = [wd '/' dir '/' netcdf_file ];
+ncid = netcdf.open(filename,'nc_nowrite');
+
+latVertex = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'latVertex'));
+lonVertex = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'lonVertex'));
+verticesOnEdge = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'verticesOnEdge'));
+edgesOnVertex = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'edgesOnVertex'));
+[dimname,nEdges]= netcdf.inqDim(ncid,netcdf.inqDimID(ncid,'nEdges'));
+[dimname,nVertices]= netcdf.inqDim(ncid,netcdf.inqDimID(ncid,'nVertices'));
+[dimname,vertexDegree]= netcdf.inqDim(ncid,netcdf.inqDimID(ncid,'vertexDegree'));
+netcdf.close(ncid)
+
+nSections = size(sectionCoord,1);
+
+% Grid variables should be:
+% lat varies from -pi/2:pi/2
+% lon varies from 0:2*pi
+if (min(lonVertex)<-1e-8)
+ lonVertex = mod(lonVertex,2*pi);
+end
+% convert to degrees for plotting:
+latVertexDeg = latVertex*180/pi;
+lonVertexDeg = lonVertex*180/pi;
+
+sectionVertexIndex = zeros(maxEdges,nSections);
+sectionEdgeIndex = zeros(maxEdges,nSections);
+sectionEdgeSign = zeros(maxEdges,nSections);
+nEdgesInSection = zeros(1,nSections);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Find edges that connect beginning and ending points
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+for iSection=1:nSections
+ latCoord = [sectionCoord(iSection,1) sectionCoord(iSection,3)]/180*pi;
+ lonCoord = [sectionCoord(iSection,2) sectionCoord(iSection,4)]/180*pi;
+
+ % Find vertex closest to start and end coordinates.
+ % The seed vertex array simply stores start and end index.
+ minDist = 1e10*ones(1,2);
+ seedVertexIndex = zeros(1,2);
+ for iVertex = 1:nVertices
+ for i=1:2
+ dist = sqrt( ...
+ (lonCoord(i) - lonVertex(iVertex))^2 ...
+ + (latCoord(i) - latVertex(iVertex))^2);
+ if (dist<minDist(i))
+ minDist(i) = dist;
+ seedVertexIndex(i) = iVertex;
+% fprintf([' iVertex %i, i %i, dist %g lonVertex %g latVertex' ...
+% ' %g </font>
<font color="blue">'],iVertex,i,dist,latVertex(iVertex),lonVertex(iVertex))
+ end
+ end
+ end
+
+ % Rename first vertex on route:
+ lonBeginVertex = lonVertex(seedVertexIndex(1));
+ latBeginVertex = latVertex(seedVertexIndex(1));
+ beginVertexIndex = seedVertexIndex(1);
+
+ % Assign first vertex on route to vertex index array:
+ sectionVertexIndex(1,iSection) = beginVertexIndex;
+
+ % Rename last vertex on route:
+ lonEndVertex = lonVertex(seedVertexIndex(2));
+ latEndVertex = latVertex(seedVertexIndex(2));
+ endVertexIndex = seedVertexIndex(2);
+
+ % Assign distance from beginVertex to endVertex. This is the
+ % distance to 'beat' by candidate neighbor vertexs.
+ distLastVertex = sqrt( ...
+ (lonBeginVertex - lonEndVertex)^2 ...
+ + (latBeginVertex - latEndVertex)^2 );
+
+ % loop through edges in search for section
+ for i=1:maxEdges
+ additionalEdgeFound = 0;
+ minDistToLine = 1e10;
+
+ for j = 1:vertexDegree
+ iEdge = edgesOnVertex(j,sectionVertexIndex(i,iSection));
+ if (iEdge>0)
+ % Find the vertex on the other side of iEdge
+ if (verticesOnEdge(1,iEdge)==sectionVertexIndex(i,iSection))
+ iVertex = verticesOnEdge(2,iEdge);
+ % Going from vertex 1 to vertex 2. Leave positive.
+ edgeSign = 1;
+ else
+ iVertex = verticesOnEdge(1,iEdge);
+ % Going from vertex 2 to vertex 1. Make negative.
+ edgeSign = -1;
+ end
+
+ % I am using lat/lon Cartesian distance.
+ % This is distance to the final vertex location.
+ dist = sqrt( ...
+ (lonVertex(iVertex) - lonVertex(endVertexIndex))^2 ...
+ + (latVertex(iVertex) - latVertex(endVertexIndex))^2 );
+
+%fprintf('%6i %6i %8.4f %8.4f h1=plot(%g,%g); h2=plot(%g,%g); </font>
<font color="blue">',...
+%i,j,dist,distLastVertex,...
+% lonVertex(iVertex)*180/pi,latVertex(iVertex)*180/pi,...
+% lonVertex(endVertexIndex)*180/pi,latVertex(endVertexIndex)*180/pi)
+ % check if this vertex is closer to the end vertex than the
+ % last vertex. If so, it is a candidate, and we can continue.
+ if (dist<distLastVertex)
+
+ % distance formula from http://mathworld.wolfram.com/Point-LineDistance2-Dimensional.html
+ % distToLine = abs((x2-x1)*(y1-y0) - (x1-x0)*(y2-y1)) ...
+ % / sqrt((x2-x1)^2 + (y2-y1)^2);
+
+ distToLine = abs(...
+ (lonEndVertex -lonBeginVertex )*(latBeginVertex-latVertex(iVertex)) ...
+ - (lonBeginVertex-lonVertex(iVertex))*(latEndVertex -latBeginVertex) ) ...
+ / sqrt((lonEndVertex-lonBeginVertex)^2 + (latEndVertex-latBeginVertex)^2);
+
+ % Find the vertex that is closest to the line connecting
+ % beginning and ending vertexs.
+ if (distToLine<minDistToLine)
+ additionalEdgeFound = 1;
+ distChosenVertex = dist;
+ minDistToLine = distToLine;
+ minDistVertexIndex = iVertex;
+ minDistEdgeIndex = iEdge;
+ minDistEdgeSign = edgeSign;
+ end
+
+ end
+ end
+
+ end % nEdgesOnEdge
+
+ distLastVertex = distChosenVertex;
+
+ sectionVertexIndex(i+1,iSection) = minDistVertexIndex;
+ sectionEdgeIndex (i,iSection) = minDistEdgeIndex;
+ sectionEdgeSign (i,iSection) = minDistEdgeSign ;
+ nEdgesInSection(iSection) = i;
+
+ if (minDistVertexIndex==endVertexIndex)
+ break
+ end
+
+ end % maxEdges
+
+ temptext = char(sectionText(iSection));
+ fprintf(['Section %3i, ' temptext(1:49) ' '],iSection)
+ if (minDistVertexIndex==endVertexIndex)
+ fprintf(' complete with %g edges.</font>
<font color="blue">',nEdgesInSection(iSection))
+ else
+ fprintf(['WARNING: Did not complete path to end' ...
+ ' vertex with %g maxEdges. </font>
<font color="blue">'],maxEdges)
+ end
+
+end % iSections
+
+maxNEdgesInSection = max(nEdgesInSection);
+
+% reduce size of sectionEdgeIndex array
+sectionEdgeIndex = sectionEdgeIndex(1:maxNEdgesInSection,:);
+sectionEdgeSign = sectionEdgeSign(1:maxNEdgesInSection,:);
+
+% assign lat/lon on section
+latSection = zeros(maxNEdgesInSection+1,nSections);
+lonSection = zeros(maxNEdgesInSection+1,nSections);
+for iSection = 1:nSections
+ for i=1:nEdgesInSection(iSection)+1
+ iVertex = sectionVertexIndex(i,iSection);
+ latSectionVertex(i,iSection) = latVertexDeg(iVertex);
+ lonSectionVertex(i,iSection) = lonVertexDeg(iVertex);
+ end
+end
+
+fprintf(['</font>
<font color="gray">'])
Added: branches/tools/transport_sections/load_large_variables_edge.m
===================================================================
--- branches/tools/transport_sections/load_large_variables_edge.m (rev 0)
+++ branches/tools/transport_sections/load_large_variables_edge.m 2012-05-22 16:18:06 UTC (rev 1927)
@@ -0,0 +1,68 @@
+function [sectionData] = load_large_variables_edge ...
+ (wd,dir,netcdf_file, var_name, var_conv_factor, ...
+ sectionEdgeIndex, nEdgesInSection)
+
+% Load large variables from netcdf file
+%
+% Mark Petersen, MPAS-Ocean Team, LANL, May 2012
+%
+%%%%%%%%%% input arguments %%%%%%%%%
+% The text string [wd '/' dir '/' netcdf_file ] is the file path,
+% where wd is the working directory and dir is the run directory.
+% var_name(nVars) a cell array with text for each variable to
+% load or compute.
+% var_conv_factor multiply each variable by this unit conversion.
+% sectionEdgeIndex(maxEdges,nSections) cell index of each section
+% nEdgesInSection(nSections) number of cells in each section
+%
+%%%%%%%%%% output arguments %%%%%%%%%
+% sectionData(nVertLevels,max(nEdgesInSection),nSections,nVars)
+% data in each cross-section for each variable
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Load large variables
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+fprintf('</font>
<font color="blue">')
+fprintf(['** load_large_variables simulation: ' dir '</font>
<font color="blue">'])
+
+filename = [wd '/' dir '/' netcdf_file ];
+ncid = netcdf.open(filename,'nc_nowrite');
+
+nAccumulate = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'nAccumulate'));
+
+[dimname,nVertLevels]= netcdf.inqDim(ncid,netcdf.inqDimID(ncid,'nVertLevels'));
+[dimname,nEdges]= netcdf.inqDim(ncid,netcdf.inqDimID(ncid,'nEdges'));
+[dimname,nTimeSlices]= netcdf.inqDim(ncid,netcdf.inqDimID(ncid,'Time'));
+nVars = length(var_name);
+nSections = length(nEdgesInSection);
+
+maxNumberEdges = max(nEdgesInSection);
+sectionData = zeros(nVertLevels,maxNumberEdges,nSections,nVars);
+
+for iVar=1:nVars
+ temptext = char(var_name(iVar));
+ fprintf(['loading: ' temptext '</font>
<font color="blue">'])
+
+ acc_var = netcdf.getVar(ncid,netcdf.inqVarID(ncid,char(var_name(iVar))));
+ mean_var = zeros(nVertLevels, nEdges);
+ for i=1:nTimeSlices
+ mean_var = mean_var + nAccumulate(i)*squeeze(acc_var(:,:,i));
+ end
+ mean_var = mean_var/sum(nAccumulate)*var_conv_factor(iVar);
+
+ for iSection = 1:nSections
+ for i=1:nEdgesInSection(iSection)
+ iEdge = sectionEdgeIndex(i,iSection);
+ for k=1:nVertLevels
+ sectionData(k,i,iSection,iVar) = mean_var(k,iEdge);
+ end
+ end
+ end
+
+end
+netcdf.close(ncid)
+
+fprintf('</font>
<font color="gray">')
+
Added: branches/tools/transport_sections/sub_plot_edge_sections.m
===================================================================
--- branches/tools/transport_sections/sub_plot_edge_sections.m (rev 0)
+++ branches/tools/transport_sections/sub_plot_edge_sections.m 2012-05-22 16:18:06 UTC (rev 1927)
@@ -0,0 +1,116 @@
+function sub_plot_edge_sections(dir,sectionCoord, ...
+ latSectionVertex,lonSectionVertex, ...
+ latVertexDeg,lonVertexDeg, ...
+ sectionEdgeIndex, nEdgesInSection,...
+ fid_latex)
+
+% Plot cell section locations on world map
+
+% Mark Petersen, MPAS-Ocean Team, LANL, May 2012
+
+%%%%%%%%%% input arguments %%%%%%%%%
+% dir text string, name of simulation
+% sectionCoord(nSections,4) endpoints of sections, with one section per row as
+% [startlat startlon endlat endlon]
+% latCellDeg(nCells) lat arrays for all cells
+% lonCellDeg(nCells) lon arrays for all cells
+% sectionEdgeIndex(maxEdges,nSections) cell index of each section
+% nEdgesInSection(nSections) number of cells in each section
+% fid_latex file ID of latex file
+
+fprintf(['** sub_plot_cell_sections, on figure 1.</font>
<font color="blue">'])
+
+nSections = size(sectionCoord,1);
+
+figure(1); clf
+
+if (min(lonVertexDeg)>-1e-8)
+ minLon = 0.0;
+ latTrans = 360;
+else
+ minLon = -180.0;
+ latTrans = 0.0;
+end
+
+ % plot topo data of the earth. This is just low-rez one deg
+ % data for visual reference.
+ load('topo.mat','topo','topomap1');
+ if minLon==-180
+ topoNew(:,1:180) = topo(:,181:360);
+ topoNew(:,181:360) = topo(:,1:180);
+ image([-180 180],[-90 90],topoNew,'CDataMapping', 'scaled');
+ else
+ image([0 360],[-90 90],topo,'CDataMapping', 'scaled');
+ end
+
+ colormap(topomap1);
+ set(gca,'YDir','normal')
+
+ hold on
+
+ % world
+ axis tight
+ set(gca,'XTick',30*[-10:12])
+ set(gca,'YTick',15*[-20:20])
+
+ % half world
+% axis([-240+latTrans 0+latTrans -80 70])
+% set(gca,'XTick',20*[-10:20])
+% set(gca,'YTick',10*[-20:20])
+
+ % N Atlantic
+% axis([-90+latTrans -5+latTrans -5 70])
+% set(gca,'XTick',[-100:5:360])
+% set(gca,'YTick',[-90:5:90])
+
+ % Drake passage
+% axis([-90+latTrans,-50+latTrans,-75,-50])
+% set(gca,'XTick',[-100:2:360])
+ % set(gca,'YTick',[-200:2:200])
+
+ % Pacific
+% axis([130 260 -10 10])
+% set(gca,'XTick',[0:1:300])
+% set(gca,'YTick',[-20:.1:20])
+
+
+ % plot vertexs. This is just done for debugging.
+ %plot(lonVertexDeg,latVertexDeg,'.y')
+
+ grid on
+
+ for iSection=1:nSections
+ latCoordDeg = [sectionCoord(iSection,1) sectionCoord(iSection,3)];
+ lonCoordDeg = [sectionCoord(iSection,2) sectionCoord(iSection,4)];
+
+ h=plot([mod(lonCoordDeg,360)],[latCoordDeg],'*-');
+ set(h,'Color','y','LineWidth',1)
+
+ for i=1:nEdgesInSection(iSection)
+ h = line([lonSectionVertex(i,iSection) lonSectionVertex(i+1,iSection)],...
+ [latSectionVertex(i,iSection) latSectionVertex(i+1,iSection)]);
+ set(h,'Color','r','LineWidth',2)
+ %plot([lonVertexDeg(sectionVertexIndex(i+1,iSection))], ...
+ % [latVertexDeg(sectionVertexIndex(i+1,iSection))],'sk')
+ end
+ end
+
+ ylabel('latitude')
+ xlabel('longitude')
+ title(['Domain: ' regexprep(dir,'_','\\_') ' Edges of transport sections. '])
+
+ set(gcf,'PaperPositionMode','auto','color',[.8 1 .8], ...
+ 'PaperPosition',[0.25 0.25 8 8])
+
+ subplot('position',[0 .95 1 .05]); axis off
+ text(.005,.7,[ date ]);
+
+ dir_name1 = regexprep(dir,'\.','_');
+ dir_name2 = regexprep(dir_name1,'/','_');
+ filename=['f/' dir_name2 '_vertex_map' ];
+ print('-djpeg',[filename '.jpg'])
+
+ % put printing text in a latex file
+ fprintf(fid_latex,...
+ ['\\begin{figure}[btp] \\center </font>
<font color="blue"> \\includegraphics[width=7.5in]{'...
+ filename '.jpg} </font>
<font color="black">\\end{figure} </font>
<font color="gray">']);
Added: branches/tools/transport_sections/transport_sections.m
===================================================================
--- branches/tools/transport_sections/transport_sections.m (rev 0)
+++ branches/tools/transport_sections/transport_sections.m 2012-05-22 16:18:06 UTC (rev 1927)
@@ -0,0 +1,227 @@
+function transport_sections
+
+% Specify data files, coordinates and text, then call functions
+% to find edge sections, load data, and compute transport through
+% each section.
+%
+% This script produces new netcdf files in the subdirectory
+% netcdf_files which can then be merged with grid.nc or restart.nc
+% files to collect transport data in MPAS-Ocean
+%
+% To merge the new *_section_edge_data.nc with an existing grid or
+% restart file, use:
+% ncks -A -v sectionEdgeIndex,sectionEdgeSign,nEdgesInSection,\
+% sectionText,sectionAbbreviation,sectionCoord \
+% your_file_section_edge_data.nc your_restart_file.nc
+%
+% Mark Petersen, MPAS-Ocean Team, LANL, May 2012
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Specify data files
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% all plots are placed in the f directory. Comment out if not needed.
+unix('mkdir -p f netcdf_files docs text_files');
+
+% The text string [wd '/' sim(i).dir '/' sim(i).netcdf_file ] is the file path,
+% where wd is the working directory and dir is the run directory.
+
+wd = '/var/tmp/mpeterse/runs/todds_runs/P1';
+
+% These files only need to contain a small number of variables.
+% You may need to reduce the file size before copying to a local
+% machine using:
+% ncks -v acc_u, \
+% nAccumulate,latVertex,lonVertex,verticesOnEdge,edgesOnVertex,hZLevel,\
+% dvEdge \
+% file_in.nc file_out.nc
+%
+% and acc_u is only required if the flux is computed by this matlab script.
+
+sim(1).dir = 'x5.NA.75km_15km';
+sim(1).netcdf_file = ['o.x5.NA.75km_15km.0027-01-01_00.00.00_acc_u.nc'];
+
+sim(2).dir = 'x5.NA.75km_15km/LEITH';
+sim(2).netcdf_file = ['o.x5.NA.75km_15km.0027-01-01_00.00.00_acc_u.nc'];
+
+sim(3).dir = 'x1.15km';
+sim(3).netcdf_file = ['o.x1.15km.0018-11-03_00.00.00_acc_u.nc'];
+
+sim(4).dir = 'x5.NA.37.5km_7.5km';
+sim(4).netcdf_file = ['o.x5.NA.37.5km_7.5km.0007-01-01_00.00.00_acc_u.nc'];
+
+%clear sim
+%sim(1).dir='x1.120km';
+%sim(1).netcdf_file = 'output120km.0001-10-22_12:30:00.nc';
+
+%clear sim
+%sim(1).dir='p91l';
+%sim(1).netcdf_file = 'output120km.0001-10-22_00.00.00.nc';
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Specify section coordinates and text
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% sectionText a cell array with text describing each section
+sectionText = {
+'Drake Passage, -56 to -63 lat, 68W lon, section A21, 140+/- 6 Sv in Ganachaud_Wunsch00n and Ganachaud99thesis ',...
+'S Ocean, Tasmania-Ant, -44 to -66 lat, 140E lon, section P12, 157+/-10 Sv in Ganachaud_Wunsch00n and Ganachaud99thesis ',...
+'S Ocean, Africa-Ant, -31.3to -70 lat, 30E lon, section I6, Sv in Ganachaud99thesis ',...
+'Florida Current, 27 lat, -80 to -78.8lon, 31.5+/-1.5Sv in Johns_ea02dsr, 32.3+/-3.2Sv Larsen92rslpt '...
+'Indonesian Throughflow, -9 to -18 lat, 116E lon, section J89, 16+/- 5 Sv in Ganachaud_Wunsch00n and Ganachaud99thesis ',...
+'Mozambique Channel, -25 lat, 35 to 44E lon, section I4 , 14+/- 6 Sv in Ganachaud_Wunsch00n and Ganachaud99thesis ',...
+ };
+
+% sectionAbbreviation an 8-character title for each section
+sectionAbbreviation = [...
+ 'Drake Pa';...
+ 'Tasm-Ant';...
+ 'Afri-Ant';...
+ 'Flor Cur';...
+ 'Ind Thru';...
+ 'Mozam Ch';...
+];
+
+% sectionCoord(nSections,4) endpoints of sections, with one section per row as
+% [startlat startlon endlat endlon]
+% Traverse from south to north, and from east to west.
+% Then positive velocities are eastward and northward.
+sectionCoord = [...
+ -64.5 -64 -55 -65.3;... % Drake
+ -67 140 -43.5 147;... % Tasm-Ant
+ -31.3 30 -70.0 30;... % Afri-Ant
+ 26.52 -78.78 26.7 -80;... % Flor Cur
+ -21 116 -8.8 116;... % Ind Thru
+ -25 44 -25 34;... % Mozam Ch
+ ];
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Specify variables
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% var_name(nVars) a cell array with text for each variable to
+% load or compute.
+% var_conv_factor multiply each variable by this unit conversion.
+% var_lims(nVars,3) contour line definition: min, max, interval
+
+var_name = {...
+'acc_u',...
+}
+
+var_conv_factor = [1 1 1]; % No conversion here.
+
+var_lims = [-10 10 2.5; -10 10 2.5; 0 20 2.5];
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Specify actions to be taken
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+find_edge_sections_flag = true ;
+write_edge_sections_text_flag = false;
+write_edge_sections_netcdf_flag = false;
+plot_edge_sections_flag = false;
+compute_transport_flag = true ;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Begin main code. Normally this does not need to change.
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%close all
+
+% change the coordinate range to be 0 to 360.
+sectionCoord(:,2) = mod(sectionCoord(:,2),360);
+sectionCoord(:,4) = mod(sectionCoord(:,4),360);
+
+for iSim = 1:length(sim)
+
+ fprintf(['**** simulation: ' sim(iSim).dir '</font>
<font color="blue">'])
+ unix(['mkdir -p docs/' sim(iSim).netcdf_file '_dir/f']);
+ fid_latex = fopen('temp.tex','w');
+ fprintf(fid_latex,['%% file created by plot_mpas_cross_sections, ' date '</font>
<font color="black"></font>
<font color="gray">']);
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %
+ % Find edges that connect beginning and end points of section
+ %
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ if find_edge_sections_flag
+ [sim(iSim).sectionEdgeIndex, sim(iSim).sectionEdgeSign, sim(iSim).nEdgesInSection, ...
+ sim(iSim).latSectionVertex,sim(iSim).lonSectionVertex, ...
+ sim(iSim).latVertexDeg,sim(iSim).lonVertexDeg] ...
+ = find_edge_sections(wd,sim(iSim).dir,sim(iSim).netcdf_file, ...
+ sectionText,sectionCoord);
+ end
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %
+ % Write section edge information to a netcdf file
+ %
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ if write_edge_sections_text_flag
+ write_edge_sections_text...
+ (sim(iSim).dir, sim(iSim).sectionEdgeIndex, ...
+ sim(iSim).sectionEdgeSign, sim(iSim).nEdgesInSection, ...
+ sectionText,sectionAbbreviation,sectionCoord)
+ end
+
+ if write_edge_sections_netcdf_flag
+ write_edge_sections_netcdf...
+ (sim(iSim).dir, sim(iSim).sectionEdgeIndex, ...
+ sim(iSim).sectionEdgeSign, sim(iSim).nEdgesInSection, ...
+ sectionText,sectionAbbreviation,sectionCoord)
+ end
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %
+ % Plot edge section locations on world map
+ %
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ if plot_edge_sections_flag
+ sub_plot_edge_sections(sim(iSim).dir,sectionCoord, ...
+ sim(iSim).latSectionVertex,sim(iSim).lonSectionVertex, ...
+ sim(iSim).latVertexDeg,sim(iSim).lonVertexDeg, ...
+ sim(iSim).sectionEdgeIndex, sim(iSim).nEdgesInSection,...
+ fid_latex);
+ end
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %
+ % Load large variables from netcdf file
+ %
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ if compute_transport_flag
+ [sim(iSim).sectionData] = load_large_variables_edge ...
+ (wd,sim(iSim).dir,sim(iSim).netcdf_file, var_name, var_conv_factor, ...
+ sim(iSim).sectionEdgeIndex, sim(iSim).nEdgesInSection);
+ end
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %
+ % Compute transport through each section
+ %
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ if compute_transport_flag
+ compute_transport ...
+ (wd,sim(iSim).dir,sim(iSim).netcdf_file, var_name, ...
+ sim(iSim).sectionEdgeIndex, sim(iSim).sectionEdgeSign, sim(iSim).nEdgesInSection, ...
+ sim(iSim).sectionData,sectionText,sectionAbbreviation)
+ end
+
+ close(fid_latex)
+
+end % iSim
+
+
Added: branches/tools/transport_sections/write_edge_sections_netcdf.m
===================================================================
--- branches/tools/transport_sections/write_edge_sections_netcdf.m (rev 0)
+++ branches/tools/transport_sections/write_edge_sections_netcdf.m 2012-05-22 16:18:06 UTC (rev 1927)
@@ -0,0 +1,76 @@
+function write_edge_sections_netcdf ...
+ (dir, ...
+ sectionEdgeIndex, sectionEdgeSign, nEdgesInSection,...
+ sectionText,sectionAbbreviation,sectionCoord)
+
+% Write section edge information to the netcdf file
+% netcdf_files/your_dir_transport_section_edges.nc
+%
+% Mark Petersen, MPAS-Ocean Team, LANL, May 2012
+%
+%%%%%%%%%% input arguments %%%%%%%%%
+% dir string with run directory name
+% sectionEdgeIndex(maxNEdgesInSection,nSections) edge index of each section
+% sectionEdgeSign(maxNEdgesInSection,nSections) sign of each
+% section, positive is to right of path direction.
+% nEdgesInSection(nSections) number of cells in each section
+% sectionText a cell array with text describing each section
+% sectionAbbreviation an 8-character title for each section
+% sectionCoord(nSections,4) endpoints of sections, with one section per row as
+% [startlat startlon endlat endlon]
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Write section edge information to a netcdf file
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+fprintf('</font>
<font color="blue">')
+fprintf(['** Write edge information to file: ' dir '</font>
<font color="blue">'])
+
+nSections = length(nEdgesInSection);
+maxNEdgesInSection = max(nEdgesInSection);
+
+dir_name1 = regexprep(dir,'/','_');
+filename = ['netcdf_files/' dir_name1 '_transport_section_edges.nc'];
+ncid = netcdf.create(filename,'nc_clobber');
+
+% Define the dimensions of the variable.
+dimid_nSections = netcdf.defDim(ncid,'nSections',nSections);
+dimid_maxNEdgesInSection = netcdf.defDim(ncid,'maxNEdgesInSection',maxNEdgesInSection);
+dimid_latLonPairs = netcdf.defDim(ncid,'latLonPairs',4);
+dimid_CharLength8 = netcdf.defDim(ncid,'CharLength8',8);
+dimid_CharLength120 = netcdf.defDim(ncid,'CharLength120',120);
+
+% Define a new variable in the file.
+sectionEdgeIndex_varID = netcdf.defVar(ncid,'sectionEdgeIndex',...
+ 'int',[dimid_maxNEdgesInSection dimid_nSections]);
+
+sectionEdgeSign_varID = netcdf.defVar(ncid,'sectionEdgeSign',...
+ 'int',[dimid_maxNEdgesInSection dimid_nSections]);
+
+nEdgesInSection_varID = netcdf.defVar(ncid,'nEdgesInSection',...
+ 'int', [dimid_nSections]);
+
+sectionText_varID = netcdf.defVar(ncid,'sectionText',...
+ 'char',[dimid_CharLength120 dimid_nSections]);
+sectionAbbreviation_varID = netcdf.defVar(ncid,'sectionAbbreviation',...
+ 'char',[dimid_CharLength8 dimid_nSections]);
+sectionCoord_varID = netcdf.defVar(ncid,'sectionCoord',...
+ 'double',[dimid_latLonPairs dimid_nSections]);
+
+
+% Leave define mode and enter data mode to write data.
+netcdf.endDef(ncid)
+
+% Write data to variable.
+netcdf.putVar(ncid,sectionEdgeIndex_varID,sectionEdgeIndex);
+netcdf.putVar(ncid,sectionEdgeSign_varID,sectionEdgeSign);
+netcdf.putVar(ncid,nEdgesInSection_varID,nEdgesInSection);
+netcdf.putVar(ncid,sectionText_varID,char(sectionText)');
+netcdf.putVar(ncid,sectionAbbreviation_varID,sectionAbbreviation');
+netcdf.putVar(ncid,sectionCoord_varID,sectionCoord);
+
+netcdf.close(ncid)
+
+fprintf('</font>
<font color="gray">')
+
Added: branches/tools/transport_sections/write_edge_sections_text.m
===================================================================
--- branches/tools/transport_sections/write_edge_sections_text.m (rev 0)
+++ branches/tools/transport_sections/write_edge_sections_text.m 2012-05-22 16:18:06 UTC (rev 1927)
@@ -0,0 +1,104 @@
+function write_edge_sections_text ...
+ (dir, ...
+ sectionEdgeIndex, sectionEdgeSign, nEdgesInSection,...
+ sectionText,sectionAbbreviation,sectionCoord)
+
+% Write section edge information to the text file
+% text_files/your_dir_transport_section_edges.nc
+%
+% Mark Petersen, MPAS-Ocean Team, LANL, May 2012
+%
+%%%%%%%%%% input arguments %%%%%%%%%
+% dir string with run directory name
+% sectionEdgeIndex(maxNEdgesInSection,nSections) edge index of each section
+% sectionEdgeSign(maxNEdgesInSection,nSections) sign of each
+% section, positive is to right of path direction.
+% nEdgesInSection(nSections) number of cells in each section
+% sectionText a cell array with text describing each section
+% sectionAbbreviation an 8-character title for each section
+% sectionCoord(nSections,4) endpoints of sections, with one section per row as
+% [startlat startlon endlat endlon]
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Write section edge information to a text file
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+fprintf('</font>
<font color="blue">')
+fprintf(['** Write edge information to file: ' dir '</font>
<font color="blue">'])
+
+nSections = length(nEdgesInSection);
+maxNEdgesInSection = max(nEdgesInSection);
+
+dir_name1 = regexprep(dir,'/','_');
+unix(['mkdir text_files/' dir_name1 ]);
+
+% sectionEdgeIndex
+filename = ['text_files/' dir_name1 '/sectionEdgeIndex.txt'];
+fid = fopen(filename,'w');
+for j=1:nSections
+ fprintf(fid,' %10i',sectionEdgeIndex(:,j));
+ fprintf(fid,' </font>
<font color="blue">');
+end
+fclose(fid);
+
+% sectionEdgeIndex
+dir_name1 = regexprep(dir,'/','_');
+filename = ['text_files/' dir_name1 '/sectionEdgeIndex.txt'];
+fid = fopen(filename,'w');
+for j=1:nSections
+ fprintf(fid,' %10i',sectionEdgeIndex(:,j));
+ fprintf(fid,' </font>
<font color="blue">');
+end
+fclose(fid);
+
+% sectionEdgeSign
+dir_name1 = regexprep(dir,'/','_');
+filename = ['text_files/' dir_name1 '/sectionEdgeSign.txt'];
+fid = fopen(filename,'w');
+for j=1:nSections
+ fprintf(fid,' %10i',sectionEdgeSign(:,j));
+ fprintf(fid,' </font>
<font color="blue">');
+end
+fclose(fid);
+
+% nEdgesInSection
+dir_name1 = regexprep(dir,'/','_');
+filename = ['text_files/' dir_name1 '/nEdgesInSection.txt'];
+fid = fopen(filename,'w');
+for j=1:nSections
+ fprintf(fid,' %10i',nEdgesInSection(j));
+ fprintf(fid,' </font>
<font color="blue">');
+end
+fclose(fid);
+
+% sectionText
+dir_name1 = regexprep(dir,'/','_');
+filename = ['text_files/' dir_name1 '/sectionText.txt'];
+fid = fopen(filename,'w');
+for j=1:nSections
+ fprintf(fid,' %s',char(sectionText(j)));
+ fprintf(fid,' </font>
<font color="blue">');
+end
+fclose(fid);
+
+% sectionAbbreviation
+dir_name1 = regexprep(dir,'/','_');
+filename = ['text_files/' dir_name1 '/sectionAbbreviation.txt'];
+fid = fopen(filename,'w');
+for j=1:nSections
+ fprintf(fid,' %s',sectionAbbreviation(j,:));
+ fprintf(fid,' </font>
<font color="blue">');
+end
+fclose(fid);
+
+% sectionCoord
+dir_name1 = regexprep(dir,'/','_');
+filename = ['text_files/' dir_name1 '/sectionCoord.txt'];
+fid = fopen(filename,'w');
+for j=1:nSections
+ fprintf(fid,' %10.3f',sectionCoord(j,:));
+ fprintf(fid,' </font>
<font color="blue">');
+end
+fclose(fid);
+
</font>
</pre>