<p><b>mpetersen@lanl.gov</b> 2012-09-07 13:04:33 -0600 (Fri, 07 Sep 2012)</p><p>Adding Matlab code to create cross-sections by interpolation. This was used to create images of the Equatorial Undercurrent for the paper, and fixes the previous problem of jagged contours by following cells. The code is written generally for any cross-section, but there are lines in sub_plot_cross_section.m that are specific to the formatting of the EUC right now, like contour intervals and axis labels.<br>
</p><hr noshade><pre><font color="gray">Added: branches/tools/cross_section_interp/cross_section_interp.m
===================================================================
--- branches/tools/cross_section_interp/cross_section_interp.m (rev 0)
+++ branches/tools/cross_section_interp/cross_section_interp.m 2012-09-07 19:04:33 UTC (rev 2146)
@@ -0,0 +1,253 @@
+%function cross_section_interp
+
+% Plot cross-sections of MPAS fields.
+%
+% This is the main function, where the user can specify data files,
+% coordinates and text, then call functions to find sections, load
+% data, and plot cross-sections.
+%
+% The final product is a set of plots as jpg files, a latex file,
+% and a compiled pdf file of the plots, if desired.
+%
+% Mark Petersen, MPAS-Ocean Team, LANL, Sept 2012
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Specify data files
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% all plots are placed in the f directory. Comment out if not needed.
+unix('mkdir -p f docs');
+
+% 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 = '/';
+
+% 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_uReconstructMeridional,acc_uReconstructZonal, \
+% nAccumulate,latVertex,lonVertex,verticesOnEdge,edgesOnVertex,hZLevel,\
+% dvEdge,latCell,lonCell,cellsOnCell,nEdgesOnCell \
+% file_in.nc file_out.nc
+
+sim(1).dir = '/var/tmp/mpeterse/runs/todds_runs/P1/x5.NA.15km/';
+sim(1).netcdf_file = ['total_avg_o.x5.NA.15km.0020-01-01_00.00.00.nc'];
+
+sim(2).dir = '/var/tmp/mpeterse/runs/todds_runs/P1/x5.NA.7.5km/';
+sim(2).netcdf_file = ['total_avg_o.x5.NA.7.5km.0021-01-01_00.00.00.nc'];
+
+sim(3).dir = '/var/tmp/mpeterse/runs/todds_runs/P1/x1.15km/';
+sim(3).netcdf_file = ['total_avg_o.x1.15km.0020-01-01_00.00.00.nc'];
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Specify section coordinates and text
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% sectionText a cell array with text describing each section
+sectionText = {
+'N Atlantic 26N lat',...
+'N Atlantic 36N lat',...
+'N Atlantic 41N lat',...
+'N Atlantic 46N lat',...
+'N Atlantic 56N lat',...
+'N Atlantic 70W lon',...
+'N Atlantic 65W lon',...
+'N Atlantic 60W lon',...
+'N Atlantic 50W lon',...
+'N Atlantic 40W lon',...
+'N Atlantic 30W lon',...
+'Eq Pacific 140W lon',...
+'Eq Pacific 0 lat ',...
+'Drake Pass 65W lon',...
+'ACC Tasman 147E lon',...
+ };
+
+% coord(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.
+coord = [...
+ 26 -80 26 -15;... % N Atl Zonal
+ 36 -76 36 -10;... % N Atl Zonal
+ 41 -72 41 -10;... % N Atl Zonal
+ 46 -60 46 -10;... % N Atl Zonal
+ 56 -60 56 -10;... % N Atl Zonal
+ 20 -70 44 -70;... % N Atl Meridional
+ 19 -65 44 -65;... % N Atl Meridional
+ 8.5 -60 46 -60;... % N Atl Meridional
+ 1.8 -50 62 -50;... % N Atl Meridional
+ -3 -40 65 -40;... % N Atl Meridional
+ -5 -30 68.2 -30;... % N Atl Meridional
+ -8 -140 8 -140;... % Eq Pac Meridional
+ 0 140 0 -95;... % Eq Pac Zonal
+ -65 -65 -55 -65;... % Drake
+ -67 147 -43.5 147;... % S Oc Tas
+ ];
+
+
+% These are the cross-sections for the current plots:
+
+% sectionText a cell array with text describing each section
+sectionText = {
+'Eq Pacific 140W lon',...
+'Eq Pacific 0 lat ',...
+ };
+
+coord = [...
+ -8 -140 10 -140;... % Eq Pac Meridional
+ 0 143 0 -95;... % Eq Pac Zonal
+ ];
+
+% number of points to plot for each figure
+nPoints = 300;
+
+% plotDepth(nSections) depth to which to plot each section
+plotDepth = 5000*ones(1,size(coord,1));
+plotDepth = 400*ones(1,size(coord,1));
+
+% For plotting, only four plots are allowed per row.
+% Choose sections above for each page.
+% page.name name of this group of sections
+% sectionID section numbers for each row of this page
+page(1).name = 'EUC';
+page(1).sectionID = [1:2];
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Specify variables to view
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% 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_uReconstructZonal',...
+'acc_uReconstructMeridional',...
+'ke_acc_uReconstruct'};
+
+var_name = {...
+'acc_uReconstructZonal',...
+};
+
+var_conv_factor = [100 100 1]; % convert m/s to cm/s for velocities
+
+var_lims = [-10 10 1.0; -10 10 1.0; 0 20 2.5];
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Specify latex command
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% This matlab script will invoke a latex compiler in order to
+% produce a pdf file. Specify the unix command-line latex
+% executable, or 'none' to not compile the latex document.
+
+latex_command = 'latex';
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Specify actions to be taken
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+find_cell_weights_flag = true ;
+plot_section_locations_flag = true ;
+load_large_variables_flag = true ;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Begin main code. Normally this does not need to change.
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% change the coordinate range to be 0 to 360.
+coord(:,2) = mod(coord(:,2),360);
+coord(:,4) = mod(coord(:,4),360);
+
+for iSim = 1:length(sim)
+
+ fprintf(['**** simulation: ' sim(iSim).dir '</font>
<font color="blue">'])
+ fid_latex = fopen('temp.tex','w');
+ fprintf(fid_latex,['%% file created by plot_mpas_cross_sections, ' date '</font>
<font color="black"></font>
<font color="blue">']);
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %
+ % Find cells that connect beginning and end points of section
+ %
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ if find_cell_weights_flag
+ [sim(iSim).cellsOnVertexSection, sim(iSim).cellWeightsSection, ...
+ sim(iSim).latSection,sim(iSim).lonSection, ...
+ depth, botDepth, sim(iSim).maxLevelCellSection] ...
+ = find_cell_weights(wd,sim(iSim).dir,sim(iSim).netcdf_file, ...
+ sectionText,coord,nPoints);
+ end
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %
+ % Plot cell section locations on world map
+ %
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ if plot_section_locations_flag
+ sub_plot_section_locations(sim(iSim).dir,coord, ...
+ sim(iSim).latSection,sim(iSim).lonSection,fid_latex)
+ end
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %
+ % Load large variables from netcdf file
+ %
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ if load_large_variables_flag
+ [sim(iSim).sectionData] = load_large_variables ...
+ (wd,sim(iSim).dir,sim(iSim).netcdf_file, var_name, var_conv_factor, ...
+ sim(iSim).cellsOnVertexSection, sim(iSim).cellWeightsSection);
+ end
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %
+ % Plot data on cross-sections
+ %
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ for iPage = 1:length(page)
+
+ sub_plot_cross_sections(sim(iSim).netcdf_file,sectionText, ...
+ page(iPage).name, page(iPage).sectionID,sim(iSim).sectionData,depth,botDepth,...
+ sim(iSim).latSection,sim(iSim).lonSection,sim(iSim).maxLevelCellSection,coord,plotDepth,...
+ var_name,var_lims,fid_latex)
+
+ end % iPage
+ fprintf(fid_latex,['</font>
<font color="black">\\end{document}</font>
<font color="black"></font>
<font color="blue">']);
+ fclose(fid_latex);
+
+% doc_dir = ['docs/' regexprep(sim(iSim).dir,'/','_') '_' ...
+% sim(iSim).netcdf_file '_dir' ];
+% unix(['mkdir -p ' doc_dir '/f']);
+% unix(['mv f/*jpg ' doc_dir '/f']);
+
+% filename = [ regexprep(sim(iSim).dir,'/','_') '_' sim(iSim).netcdf_file '.tex'];
+% unix(['cat mpas_sections.head.tex temp.tex > ' doc_dir '/' filename ]);
+
+% if not(strcmp(latex_command,'none'))
+% fprintf('*** Compiling latex document </font>
<font color="gray">')
+% cd(doc_dir);
+% unix([latex_command ' ' filename]);
+% cd('../..');
+% end
+
+end % iSim
+
+
Added: branches/tools/cross_section_interp/find_cell_weights.m
===================================================================
--- branches/tools/cross_section_interp/find_cell_weights.m (rev 0)
+++ branches/tools/cross_section_interp/find_cell_weights.m 2012-09-07 19:04:33 UTC (rev 2146)
@@ -0,0 +1,159 @@
+function [cellsOnVertexSection, cellWeightsSection, latSection,lonSection, ...
+ depth, botDepth, maxLevelCellSection] = find_cell_weights ...
+ (wd,dir,netcdf_file,sectionText,coord,nPoints)
+
+% This function reads grid data from an MPAS-Ocean grid or restart
+% netCDF file, and finds a path of cells that connect the endpoints
+% specified in coord. The path is forced to travel through cells
+% that are closest to the line connecting the beginning and end
+% cells.
+%
+% 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
+% coord(nSections,4) endpoints of sections, with one section per row as
+% [startlat startlon endlat endlon]
+%
+%%%%%%%%%% output arguments %%%%%%%%%
+% cellsOnVertexSection(vertexDegree,nPoints,nSections) cells neighboring nearest vertex
+% cellWeightsSection(vertexDegree,nPoints,nSections) weights for each cell
+% latSection(nPoints,nSections) lat coordinates of each section
+% lonSection(nPoints,nSections) lon coordinates of each section
+% depth(nVertLevels) depth of center of each layer, for plotting
+% latCellDeg(nCells) lat arrays for all cells
+% lonCellDeg(nCells) lon arrays for all cells
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Read data from file
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+fprintf(['** find_cell_sections, simulation: ' dir '</font>
<font color="blue">'])
+
+filename = [wd '/' dir '/' netcdf_file ];
+ncid = netcdf.open(filename,'nc_nowrite');
+
+xCell = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'xCell'));
+yCell = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'yCell'));
+zCell = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'zCell'));
+latVertex = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'latVertex'));
+lonVertex = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'lonVertex'));
+xVertex = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'xVertex'));
+yVertex = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'yVertex'));
+zVertex = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'zVertex'));
+cellsOnVertex = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'cellsOnVertex'));
+hZLevel = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'hZLevel'));
+maxLevelCell = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'maxLevelCell'));
+sphere_radius = netcdf.getAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'sphere_radius');
+[dimname,nVertices]= netcdf.inqDim(ncid,netcdf.inqDimID(ncid,'nVertices'));
+[dimname,vertexDegree]= netcdf.inqDim(ncid,netcdf.inqDimID(ncid,'vertexDegree'));
+[dimname,nVertLevels]= netcdf.inqDim(ncid,netcdf.inqDimID(ncid,'nVertLevels'));
+netcdf.close(ncid)
+
+nSections = size(coord,1);
+
+% Compute depth of center of each layer, for plotting
+depth(1) = hZLevel(1)/2;
+botDepth(1) = hZLevel(1);
+for i=2:nVertLevels
+ depth(i) = depth(i-1) + 0.5*(hZLevel(i) + hZLevel(i-1));
+ botDepth(i) = botDepth(i-1) + hZLevel(i);
+end
+depth(1)=0; % make top layer plot to surface
+
+latSection = zeros(nPoints,nSections);
+lonSection = zeros(nPoints,nSections);
+latVertexSection = zeros(nPoints,nSections);
+lonVertexSection = zeros(nPoints,nSections);
+maxLevelCellSection = zeros(nPoints,nSections);
+nearestVertexSection = zeros(nPoints,nSections);
+cellsOnVertexSection = zeros(vertexDegree,nPoints,nSections);
+cellWeightsSection = zeros(vertexDegree,nPoints,nSections);
+margin=.5;
+
+for iSection=1:nSections
+ fprintf('Finding nearest vertex for Section %g </font>
<font color="blue">',iSection)
+ latSection(:,iSection) = linspace(coord(iSection,1),coord(iSection,3),nPoints);
+ lonSection(:,iSection) = linspace(coord(iSection,2),coord(iSection,4),nPoints);
+
+ maxLat = (max(latSection(:,iSection))+margin)*pi/180;
+ minLat = (min(latSection(:,iSection))-margin)*pi/180;
+ maxLon = (max(lonSection(:,iSection))+margin)*pi/180;
+ minLon = (min(lonSection(:,iSection))-margin)*pi/180;
+
+ vInd = find(latVertex>minLat&latVertex<maxLat ...
+ &lonVertex>minLon&lonVertex<maxLon);
+ distToVertex = zeros(length(vInd),1);
+
+ for iPt=1:nPoints
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %
+ % Find x,y,z coordinates of each section point
+ %
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ latPt = latSection(iPt,iSection)*pi/180; % lat of point, in radians
+ lonPt = lonSection(iPt,iSection)*pi/180; % lon of point, in radians
+ zPt = sphere_radius*sin(latPt);
+ r = sphere_radius*cos(latPt);
+ xPt = r*cos(lonPt);
+ yPt = r*sin(lonPt);
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %
+ % Find nearest vertex to each section point
+ %
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ for i=1:length(vInd)
+ distToVertex(i) = sqrt(...
+ (xPt - xVertex(vInd(i)))^2 ...
+ + (yPt - yVertex(vInd(i)))^2 ...
+ + (zPt - zVertex(vInd(i)))^2 );
+ end
+
+ [mindist,i] = min(distToVertex);
+ iVertex = vInd(i);
+ nearestVertexSection(iPt,iSection) = iVertex;
+ cellsOnVertexSection(:,iPt,iSection) = cellsOnVertex(:,iVertex);
+ latVertexSection(iPt,iSection) = latVertex(iVertex)*180/pi;
+ lonVertexSection(iPt,iSection) = lonVertex(iVertex)*180/pi;
+
+ maxLevelCellSection(iPt,iSection) = min(...
+ maxLevelCell(cellsOnVertexSection(:,iPt,iSection)));
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %
+ % Compute area weights for each cell
+ %
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ vPt = [xPt,yPt,zPt];
+ for iCell=1:3
+ vCell(:,iCell) = [...
+ xCell(cellsOnVertexSection(iCell,iPt,iSection)) ...
+ yCell(cellsOnVertexSection(iCell,iPt,iSection)) ...
+ zCell(cellsOnVertexSection(iCell,iPt,iSection))];
+ end
+
+ for iCell=1:3
+ cellWeightsSection(iCell,iPt,iSection) = triArea(vPt,...
+ vCell(:,mod(iCell,3)+1),vCell(:,mod(iCell+1,3)+1),sphere_radius);
+ end
+ cellWeightsSection(:,iPt,iSection) = cellWeightsSection(:,iPt,iSection)/sum(cellWeightsSection(:,iPt,iSection));
+
+ end
+
+ % print out if desired:
+ fprintf('lat: desired, vertex lon: desired, vertex</font>
<font color="blue">')
+ %fprintf('%10.2f %10.2f %10.2f %10.2f</font>
<font color="blue">', ...
+ %[latSection(:,iSection) latVertexSection(:,iSection) ...
+ % lonSection(:,iSection) lonVertexSection(:,iSection) ]')
+
+ %'cell weights'
+ %cellWeightsSection(:,:,iSection)
+end % iSections
+
+
+fprintf(['</font>
<font color="gray">'])
Added: branches/tools/cross_section_interp/load_large_variables.m
===================================================================
--- branches/tools/cross_section_interp/load_large_variables.m (rev 0)
+++ branches/tools/cross_section_interp/load_large_variables.m 2012-09-07 19:04:33 UTC (rev 2146)
@@ -0,0 +1,82 @@
+function [sectionData] = load_large_variables ...
+ (wd,dir,netcdf_file, var_name, var_conv_factor, ...
+ cellsOnVertexSection, cellWeightsSection)
+
+% 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.
+% cellsOnVertexSection(vertexDegree,nPoints,nSections) cells neighboring nearest vertex
+% cellWeightsSection(vertexDegree,nPoints,nSections) weights for each cell
+%
+%%%%%%%%%% output arguments %%%%%%%%%
+% sectionData(nVertLevels,nPoints,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,nCells]= netcdf.inqDim(ncid,netcdf.inqDimID(ncid,'nCells'));
+
+% use only when averaging in matlab:
+%[dimname,nTimeSlices]= netcdf.inqDim(ncid,netcdf.inqDimID(ncid,'Time'));
+
+nVars = length(var_name);
+vertexDegree = size(cellsOnVertexSection,1);
+nPoints = size(cellsOnVertexSection,2);
+nSections = size(cellsOnVertexSection,3);
+
+sectionData = zeros(nVertLevels,nPoints,nSections,nVars);
+
+for iVar=1:nVars
+ temptext = char(var_name(iVar));
+ fprintf(['loading: ' temptext '</font>
<font color="blue">'])
+ if (temptext(1:6)=='ke_acc')
+ % assume zonal and meridional velocity are in index 1 and 2.
+ sectionData(:,:,:,iVar) = sqrt(sectionData(:,:,:,1).^2 + sectionData(:,:,:,2).^2)/2;
+ else
+
+ mean_var = netcdf.getVar(ncid,netcdf.inqVarID(ncid,char(var_name(iVar))))*var_conv_factor(iVar);
+
+ % use only when averaging in matlab:
+ %acc_var = netcdf.getVar(ncid,netcdf.inqVarID(ncid,char(var_name(iVar))));
+ %mean_var = zeros(nVertLevels, nCells);
+ %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 iPt = 1:nPoints
+ for k=1:nVertLevels
+ sectionData(k,iPt,iSection,iVar) = ...
+ mean_var(k,cellsOnVertexSection(1,iPt,iSection))*cellWeightsSection(1,iPt,iSection)...
+ + mean_var(k,cellsOnVertexSection(2,iPt,iSection))*cellWeightsSection(2,iPt,iSection)...
+ + mean_var(k,cellsOnVertexSection(3,iPt,iSection))*cellWeightsSection(3,iPt,iSection);
+ end
+ end
+ end
+ end
+
+end
+netcdf.close(ncid)
+
+fprintf('</font>
<font color="gray">')
+
Added: branches/tools/cross_section_interp/sub_plot_cross_sections.m
===================================================================
--- branches/tools/cross_section_interp/sub_plot_cross_sections.m (rev 0)
+++ branches/tools/cross_section_interp/sub_plot_cross_sections.m 2012-09-07 19:04:33 UTC (rev 2146)
@@ -0,0 +1,210 @@
+function sub_plot_cross_sections(dir,sectionText,pageName,sectionID,sectionData,depth,botDepth,...
+ latSection,lonSection, maxLevelCellSection,coord, plotDepth,...
+ var_name,var_lims,fid_latex)
+
+% Plot cross-sections of MPAS fields.
+%
+% Mark Petersen, MPAS-Ocean Team, LANL, May 2012
+%
+%%%%%%%%%% input arguments %%%%%%%%%
+% dir text string, name of simulation
+% sectionText a cell array with text describing each section
+% sectionID section numbers for each row of this page
+% pageName name of this group of sections
+% sectionData(nVertLevels,nPoints,nSections,nVars)
+% data in each cross-section for each variable
+% depth(nVertLevels) depth of center of each layer, for plotting
+% latSection(nPoints,nSections) lat coordinates of each section
+% lonSection(nPoints,nSections) lon coordinates of each section
+% coord(nSections,4) endpoints of sections, with one section per row as
+% [startlat startlon endlat endlon]
+% plotDepth(nSections) depth to which to plot each section
+% var_lims(nVars,3) contour line definition: min, max, interval
+% var_name(nVars) a cell array with text for each variable to
+% load or compute.
+% fid_latex file ID of latex file
+
+fprintf(['** sub_plot_cross_sections simulation: ' dir ...
+ ' plotting page: ' pageName '</font>
<font color="blue">'])
+
+px = [.28 .78];
+py=linspace(.84,.13,4); % Midpoint position of plots
+pw = [.4]; ph=[.18]; % width and height of plots
+
+nPoints = size(sectionData,2);
+nSections = size(sectionData,3);
+nVars = size(sectionData,4);
+
+for iVar=1:nVars
+ %figure(iVar+1); clf
+ %set(gcf,'Position',[100+(iVar*100) 600-iVar*100 550 400])
+ %temptext2=char(var_name(iVar));
+
+ for iRow = 1:length(sectionID)
+ figure(iRow+1); clf
+ set(gcf,'Position',[100+(iRow*100) 600-iRow*100 550 400])
+ iSection = sectionID(iRow);
+ if coord(iSection,1)==coord(iSection,3) % meridional section
+ xtext = 'longitude';
+ xaxis = lonSection(:,iSection);
+ else % zonal section
+ xtext = 'latitude';
+ xaxis = latSection(:,iSection);
+ end
+
+ % left column
+ % ha=subplot('position',[px(1)-pw/2 py(iRow)-ph/2 pw ph]);
+ % temptext = char(sectionText(iSection));
+ % if temptext2(1:6)=='ke_acc'
+ % h=surf(xaxis, depth,log10(sectionData(:,1:nCellsInSection(iSection),iSection,iVar)));
+ % set(gca,'CLim',[-1 1.2])
+ % else
+ % h=surf(xaxis, depth,sectionData(:,1:nCellsInSection(iSection),iSection,iVar));
+ % end
+
+ % set(h,'EdgeColor','none')
+ % view(0,-90)
+ % title([temptext ', cm/s'])
+ % ylabel('depth, m')
+ % xlabel(xtext)
+ % axis tight
+ % set(gca,'YLim',[0 plotDepth(iSection)])
+ % h=colorbar ;
+ % if temptext2(1:6)=='ke_acc'
+% set(h,'YTick',[-1:1:1.2],'YTickLabel',[0.1 1 10])
+ % end
+
+ % right column
+
+px = [.54];
+py=.53;
+pw = [.87]; ph=[.83]; % width and height of plots
+ ha=subplot('position',[px(1)-pw/2 py-ph/2 pw ph]);
+ temptext = char(sectionText(iSection));
+ hold on
+% contour(xaxis, depth,sectionData(:,1:nCellsInSection(iSection),iSection,iVar), ...
+% [var_lims(iVar,1):var_lims(iVar,3):var_lims(iVar,2)]);
+% set(gca,'CLim',var_lims(iVar,1:2))
+
+ %%%%%% special commands for DWBC mrp
+ % imitating colorbar at http://www.agu.org/journals/jc/jc1203/2011JC007586/figures.shtml#fig10
+
+ %xaxis = xaxis - 360*ones(size(xaxis));
+ % xaxis is now in longitude. Convert to Distance (km)
+ % along 26.5N east of 77W
+ %xaxis = (xaxis+77)*99; % for DWBC only
+ %contour_lims = [-25 -20 -15 -10 -5 -2 -1 1 2 5 10 15 20 25];
+ %contour_lims = [-20 -15 -10 -5 -2 0 2 5 10 15 20 25 30];
+ contour_lims = [-110:10:110];
+ [cout,h]=contourf(xaxis, depth,sectionData(:,:,iSection,iVar), ...
+ contour_lims);
+ set(gca,'CLim',[min(contour_lims) max(contour_lims)])
+ set(h,'LineColor',[.5 .5 .5])
+ cbfrac=0;
+
+ % Text labels on countours
+ [cout,h]=contour(xaxis, depth,sectionData(:,:,iSection,iVar),...
+ [-100 -50 -30:10:30 50 100]);
+ ls=[100 300];
+ clabel(cout,h,'fontsize',10,'color','k','rotation',0,'LabelSpacing',ls(iRow));
+ set(h,'LineColor',[.5 .5 .5])
+
+ % Black lines
+ [cout,h]=contour(xaxis, depth,sectionData(:,:,iSection,iVar),[-100:50:100]);
+ set(h,'LineColor',[0 0 0],'LineWidth',1)
+
+
+ % stretched colorbar using contour_lims:
+ cmin=min(contour_lims);
+ cmax=max(contour_lims);
+% cvalue = cmin-.5*dc:dc:cmax+.5*dc;
+ nc_orig = 256;
+ nc = length(contour_lims);
+ cmap_orig = ColdHot(nc_orig);
+ cmap_orig_short = zeros(nc-1,3);
+ ind=(.5:1:nc-.5);
+ for j=1:nc-1
+ cmap_orig_short(j,:) = cmap_orig( floor((j-.5)/(nc-1)*nc_orig),:);
+ end
+
+ cvalue = linspace(cmin,cmax,256);
+ nc_inc = length(cvalue);
+
+ cmapnew = zeros(nc_inc,3);
+ for jnew=2:nc_inc
+ jold = max(min(min(find(contour_lims>=cvalue(jnew))),nc)-1,1);
+ cmapnew(jnew-1,:) = cmap_orig_short(jold,:);
+ end
+ cmapnew(nc_inc,:) = cmap_orig_short(nc-1,:);
+
+ colormap(cmapnew)
+ %colorbarf_spec(cout,h,'vert',contour_lims);
+ %xlabel('Distance (km) along 26.5N east of 77W') % for DWBC only
+
+ axis tight
+ %set(gca,'YLim',[0 plotDepth(iSection)],'XLim',[0 175]) % for DWBC only
+ %set(gca,'YTick',[0:1000:5000],'XTick',[0:25:175])
+ set(gca,'YLim',[0 plotDepth(iRow)])
+ set(gca,'YTick',[0:100:400])
+ xlabel(xtext)
+ % set(gca,'XTick',-1*[80:.5:70])
+ %%%%%% special commands for DWBC mrp end
+
+ %%%%%% special commands for EUC mrp end
+ if iRow==2
+ set(gca,'XTick',[143 156 165 180 190 205 220 235 250 265])
+ set(gca,'XTickLabel',{'143' '156' '165E' '180' '170W' '155' '140' '125' '110' '95'})
+ end
+
+ %%%%%% special commands for EUC mrp end
+
+ set(gca,'YDir','reverse')
+ %title([temptext ', ' num2str(var_lims(iVar,3)) 'cm/s int'])
+ ylabel('depth, m')
+ grid on
+ set(gca,'layer','top');
+ h=colorbar;
+
+ % mrp draw bottom based on zero contour
+ hold on
+ n = nPoints;
+ % old way: maxLevelCell=zeros(1,n);
+ x(2:n) = (xaxis(1:n-1)+xaxis(2:n))/2;
+ x(1) = xaxis(1) - (xaxis(2)-xaxis(1))/2;
+ x(n+1) = xaxis(n) + (xaxis(n)-xaxis(n-1))/2;
+ b = max(botDepth);
+ for j=1:n
+ % old way: maxLevelCell(j)=max(min(find(sectionData(:,j,iSection,iVar)==0.0))-1,1);
+ depthline(j) = botDepth(maxLevelCellSection(j,iSection));
+ h=patch([x(j) x(j+1) x(j+1) x(j) x(j)],...
+ [b b depthline(j) depthline(j) b], [.5 .5 .5]);
+ set(h,'LineStyle','none')
+ end
+
+ % mrp draw bottom based on zero contour end
+
+ set(gcf,'PaperPositionMode','auto','color',[.8 1 .8], ...
+ 'PaperPosition',[0.25 0.25 5.5 3.2])
+ subplot('position',[0 .95 1 .05]); axis off
+ title_txt = [regexprep(char(var_name(iVar)),'_','\\_') ', ' regexprep(dir,'_','\\_')];
+% h=text(.55,.4,title_txt);
+% set(h,'HorizontalAlignment','center','FontWeight','bold','FontSize',14)
+% text(.005,.7,[ date ]);
+
+ dir_name1 = regexprep(dir,'\.','_');
+ dir_name2 = regexprep(dir_name1,'/','_');
+ %filename=['f/' dir_name2 '_' pageName '_var' num2str(iVar)];
+ filename=['f/' dir_name2 '_' pageName num2str(iRow) '_var' num2str(iVar)];
+ print('-djpeg',[filename '.jpg']);
+ print('-depsc2',[filename '.eps']);
+ unix(['convert ' filename '.eps ' filename '.pdf']);
+ 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">']);
+
+ % print('-depsc2',[filename '.eps'])
+
+ end
+
+
+end
+
Added: branches/tools/cross_section_interp/sub_plot_section_locations.m
===================================================================
--- branches/tools/cross_section_interp/sub_plot_section_locations.m (rev 0)
+++ branches/tools/cross_section_interp/sub_plot_section_locations.m 2012-09-07 19:04:33 UTC (rev 2146)
@@ -0,0 +1,98 @@
+function sub_plot_section_locations(dir,coord, ...
+ latSection,lonSection,fid_latex)
+
+% Plot section locations on world map
+
+% Mark Petersen, MPAS-Ocean Team, LANL, Sep 2012
+
+%%%%%%%%%% input arguments %%%%%%%%%
+% dir text string, name of simulation
+% coord(nSections,4) endpoints of sections, with one section per row as
+% [startlat startlon endlat endlon]
+% latSection(nPoints,nSections) lat coordinates of each section
+% lonSection(nPoints,nSections) lon coordinates of each section
+% fid_latex file ID of latex file
+
+fprintf(['** sub_plot_cell_sections, on figure 1.</font>
<font color="blue">'])
+
+nSections = size(coord,1);
+
+figure(1); clf
+
+ minLon = 0.0;
+ latTrans = 360;
+
+ % 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([-360+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])
+
+ hold on
+ grid on
+
+ for iSection=1:nSections
+ h=plot(lonSection(:,iSection),latSection(:,iSection),'r-');
+ h=text(lonSection(1,iSection),latSection(1,iSection), ...
+ num2str(iSection))
+ get(h)
+
+ set(h,'Color',[1 1 1],'FontWeight','bold')
+ %h=plot(lonSection(:,iSection),latSection(:,iSection),'y.');
+ %set(h,'Color','y','LineWidth',1)
+ end
+
+ ylabel('latitude')
+ xlabel('longitude')
+ title(['Domain: ' regexprep(dir,'_','\\_') ' Cells in cross 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 '_cell_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/cross_section_interp/triArea.m
===================================================================
--- branches/tools/cross_section_interp/triArea.m (rev 0)
+++ branches/tools/cross_section_interp/triArea.m 2012-09-07 19:04:33 UTC (rev 2146)
@@ -0,0 +1,24 @@
+function area=triArea(A,B,C,R)
+% - This function calculates the area of the triangle A,B,C on the
+% surface of a sphere.
+%
+% Input: A, B, C
+% A: vertex 1 of triangle
+% B: vertex 2 of triangle
+% C: vertex 3 of triangle
+% R: radius of sphere
+% Output: (returned value area)
+% area: surface area of triangle on sphere.
+
+R2inv = 1/R/R;
+
+a = acos(dot(B,C)*R2inv);
+b = acos(dot(C,A)*R2inv);
+c = acos(dot(A,B)*R2inv);
+
+s = 0.5*(a+b+c);
+
+tanqe = sqrt(tan(0.5*s)*tan(0.5*(s-a))*tan(0.5*(s-b))*tan(0.5*(s-c)));
+
+area = abs(4.0*atan(tanqe));
+
</font>
</pre>