<p><b>mpetersen@lanl.gov</b> 2012-05-11 09:55:34 -0600 (Fri, 11 May 2012)</p><p>Adding new tool:  branches/tools/plot_mpas_cross_sections<br>
<br>
Adding a new tool to view cross-sections of mpas output data in<br>
matlab, by finding a path of cells that most closely follows a line in<br>
lat/lon coordinates.  This plotting tool has been tested for ocean<br>
output files, but should nearly work for other cores as well.  See<br>
README file.<br>
</p><hr noshade><pre><font color="gray">Added: branches/tools/plot_mpas_cross_sections/README
===================================================================
--- branches/tools/plot_mpas_cross_sections/README                                (rev 0)
+++ branches/tools/plot_mpas_cross_sections/README        2012-05-11 15:55:34 UTC (rev 1898)
@@ -0,0 +1,33 @@
+plot_mpas_cross_sections README
+
+This is a matlab post-processing tool to view cross sections of
+MPAS-Ocean data.  To begin, change the parameters at the top of
+plot_mpas_cross_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
+plot_mpas_cross_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_uReconstructMeridional,acc_uReconstructZonal, \
+nAccumulate,latVertex,lonVertex,verticesOnEdge,edgesOnVertex,hZLevel,\
+dvEdge,latCell,lonCell,cellsOnCell,nEdgesOnCell \
+file_in.nc file_out.nc
+
+After plot_mpas_cross_sections is executed in matlab, plots are placed
+in the docs/*.nc_dir/f subdirectory, and a latex file is created and
+compiled to produce a docs/*.nc_dir/*.nc.pdf file that includes all
+the figures.
+
+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.  We compute depth from the hZLevel variable, and plot cross
+sections with depth increasing downward.
+
+Mark Petersen, MPAS-Ocean Team, LANL, May 2012

Added: branches/tools/plot_mpas_cross_sections/find_cell_sections.m
===================================================================
--- branches/tools/plot_mpas_cross_sections/find_cell_sections.m                                (rev 0)
+++ branches/tools/plot_mpas_cross_sections/find_cell_sections.m        2012-05-11 15:55:34 UTC (rev 1898)
@@ -0,0 +1,208 @@
+function [sectionCellIndex, nCellsInSection, latSection,lonSection, ...
+          depth, latCellDeg,lonCellDeg] = find_cell_sections ...
+   (wd,dir,netcdf_file,sectionText,coord)
+
+% 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 %%%%%%%%%
+% sectionCellIndex(maxCells,nSections)       cell index of each section
+% nCellsInSection(nSections)                 number of cells in each section
+% latSection(max(nCellsInSection),nSections) lat coordinates of each section
+% lonSection(max(nCellsInSection),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  
+
+%%%%%%%%%% parameters internal to this function %%%%%%%%%%
+
+% maxCells specifies the maximum number of Cells attempted along
+% the path to the end-cell before stopping with a warning.
+maxCells = 1500;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  Read cell and edge data from grid file
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+fprintf(['** find_cell_sections, simulation: ' dir '</font>
<font color="blue">'])
+ 
+filename = [wd '/' dir '/' netcdf_file ];
+ncid = netcdf.open(filename,'nc_nowrite');
+
+latCell = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'latCell'));
+lonCell = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'lonCell'));
+cellsOnCell = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'cellsOnCell'));
+nEdgesOnCell = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'nEdgesOnCell'));
+hZLevel = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'hZLevel'));
+[dimname,nCells]= netcdf.inqDim(ncid,netcdf.inqDimID(ncid,'nCells'));
+[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;
+for i=2:nVertLevels
+  depth(i) = depth(i-1) + 0.5*(hZLevel(i) + hZLevel(i-1));
+end
+
+if (min(lonCell)&gt;-1e-8)
+  for iSection=1:nSections
+    if (coord(iSection,2)&lt;0)
+      coord(iSection,2) = coord(iSection,2) + 360;
+    end
+    if (coord(iSection,4)&lt;0)
+      coord(iSection,4) = coord(iSection,4) + 360;
+    end
+  end
+end
+
+% The native grid variables use:
+% lat varies from -pi/2:pi/2
+% lon varies from -pi:pi
+% convert to degrees for plotting:
+latCellDeg = latCell*180/pi;
+lonCellDeg = lonCell*180/pi;
+
+sectionCellIndex = zeros(maxCells,nSections);
+nCellsInSection = zeros(1,nSections);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  Find cells that connect beginning and ending points
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+for iSection=1:nSections
+   latCoord = [coord(iSection,1) coord(iSection,3)]/180*pi;
+   lonCoord = [coord(iSection,2) coord(iSection,4)]/180*pi;
+
+   % Find cell closest to start and end coordinates.
+   % The seed cell array simply stores start and end index.
+   minDist = 1e10*ones(1,2);
+   seedCellIndex = zeros(2,nSections);
+   for iCell = 1:nCells
+     for i=1:2
+       dist = sqrt( ...
+          (lonCoord(i) - lonCell(iCell))^2 ...
+        + (latCoord(i) - latCell(iCell))^2);
+       if (dist&lt;minDist(i))
+         minDist(i) = dist;
+         seedCellIndex(i,iSection) = iCell;
+       end
+     end
+   end
+
+   % Rename first cell on route:
+   lonBeginCell = lonCell(seedCellIndex(1,iSection));
+   latBeginCell = latCell(seedCellIndex(1,iSection));
+   beginCellIndex = seedCellIndex(1,iSection);
+
+   % Assign first cell on route to cell index array:
+   sectionCellIndex(1,iSection) = beginCellIndex;
+
+   % Rename last cell on route:
+   lonEndCell = lonCell(seedCellIndex(2,iSection));
+   latEndCell = latCell(seedCellIndex(2,iSection));
+   endCellIndex = seedCellIndex(2,iSection);
+   
+   % Assign distance from beginCell to endCell.  This is the
+   % distance to 'beat' by candidate neighbor cells.
+   distLastCell = sqrt( ...
+       (lonBeginCell - lonEndCell)^2 ...
+     + (latBeginCell - latEndCell)^2 );
+
+   % loop through cells in search for section
+   for i=1:maxCells
+     additionalCellFound = 0;
+     minDistToLine = 1e10;
+
+     for j = 1:nEdgesOnCell(sectionCellIndex(i,iSection))
+       iCell = cellsOnCell(j,sectionCellIndex(i,iSection));
+
+       if (iCell&lt;nCells+1 &amp;&amp; iCell&gt;0)
+          % I am using lat/lon Cartesian distance.
+          % This is distance to the final cell location.
+          dist = sqrt( ...
+             (lonCell(iCell) - lonEndCell)^2 ...
+           + (latCell(iCell) - latEndCell)^2 );
+
+          % check if this cell is closer to the end cell than the
+          % last cell.  If so, it is a candidate, and we can continue.
+          if (dist&lt;distLastCell)
+            
+            % 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(...
+                 (lonEndCell  -lonBeginCell  )*(latBeginCell-latCell(iCell)) ...
+               - (lonBeginCell-lonCell(iCell))*(latEndCell  -latBeginCell) ) ...
+               / sqrt((lonEndCell-lonBeginCell)^2 + (latEndCell-latBeginCell)^2);
+
+            % Find the cell that is closest to the line connecting
+            % beginning and ending cells.
+            if (distToLine&lt;minDistToLine)
+              additionalCellFound = 1;
+              distChosenCell = dist;
+              minDistToLine = distToLine;
+              minDistCellIndex = iCell;
+            end
+            
+          end
+       end
+       
+     end  % nEdgesOnCell
+
+     distLastCell = distChosenCell;
+
+     sectionCellIndex(i+1,iSection) = minDistCellIndex;
+     nCellsInSection(iSection) = i+1;
+
+     if (minDistCellIndex==seedCellIndex(2,iSection))
+        break
+     end
+
+   end % maxCells
+
+   temptext = char(sectionText(iSection));
+   fprintf(['Section %3i, ' temptext ' '],iSection)
+   if (minDistCellIndex==seedCellIndex(2,iSection))
+     fprintf(' complete with %g cells.</font>
<font color="blue">',nCellsInSection(iSection))
+   else
+     fprintf(['WARNING: Did not complete path to end' ...
+              ' cell with %g maxCells. </font>
<font color="blue">'],maxCells)
+   end
+
+end % iSections
+
+maxNumberCells = max(nCellsInSection);
+
+% reduce size of sectionCellIndex array
+sectionCellIndex = sectionCellIndex(1:maxNumberCells,:);
+
+% assign lat/lon on section
+latSection = zeros(maxNumberCells,nSections);
+lonSection = zeros(maxNumberCells,nSections);
+for iSection = 1:nSections
+  for i=1:nCellsInSection(iSection)
+    iCell = sectionCellIndex(i,iSection);
+    latSection(i,iSection) = latCellDeg(iCell);
+    lonSection(i,iSection) = lonCellDeg(iCell);
+  end
+end
+
+fprintf(['</font>
<font color="gray">'])

Added: branches/tools/plot_mpas_cross_sections/load_large_variables.m
===================================================================
--- branches/tools/plot_mpas_cross_sections/load_large_variables.m                                (rev 0)
+++ branches/tools/plot_mpas_cross_sections/load_large_variables.m        2012-05-11 15:55:34 UTC (rev 1898)
@@ -0,0 +1,73 @@
+function [sectionData] = load_large_variables ...
+   (wd,dir,netcdf_file, var_name, var_conv_factor, ...
+    sectionCellIndex, nCellsInSection)
+
+% 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.
+% sectionCellIndex(maxCells,nSections)       cell index of each section
+% nCellsInSection(nSections)                 number of cells in each section
+
+%%%%%%%%%% output arguments %%%%%%%%%
+% sectionData(nVertLevels,max(nCellsInSection),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'));
+[dimname,nTimeSlices]= netcdf.inqDim(ncid,netcdf.inqDimID(ncid,'Time'));
+nVars = length(var_name);
+nSections = length(nCellsInSection);
+
+maxNumberCells = max(nCellsInSection);
+sectionData = zeros(nVertLevels,maxNumberCells,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
+  
+    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 i=1:nCellsInSection(iSection)
+        iCell = sectionCellIndex(i,iSection);
+        for k=1:nVertLevels
+          sectionData(k,i,iSection,iVar) = mean_var(k,iCell);
+        end
+      end
+    end
+  end
+
+end
+netcdf.close(ncid)
+
+fprintf('</font>
<font color="gray">')
+

Added: branches/tools/plot_mpas_cross_sections/mpas_sections.head.tex
===================================================================
--- branches/tools/plot_mpas_cross_sections/mpas_sections.head.tex                                (rev 0)
+++ branches/tools/plot_mpas_cross_sections/mpas_sections.head.tex        2012-05-11 15:55:34 UTC (rev 1898)
@@ -0,0 +1,16 @@
+\documentclass[11pt]{article}
+\usepackage{color,amsmath,amsfonts,amssymb,amsthm}
+ 
+\usepackage[pdftex]{graphicx}
+\DeclareGraphicsExtensions{.pdf, .jpg, .tif, .png}
+
+\textwidth=7.5in
+\textheight=9.8in
+\topmargin=-0.5in
+\headsep=0.0in
+\headheight=0.0in
+\oddsidemargin=-0.25in
+\evensidemargin=-0.25in
+
+\begin{document}
+

Added: branches/tools/plot_mpas_cross_sections/plot_mpas_cross_sections.m
===================================================================
--- branches/tools/plot_mpas_cross_sections/plot_mpas_cross_sections.m                                (rev 0)
+++ branches/tools/plot_mpas_cross_sections/plot_mpas_cross_sections.m        2012-05-11 15:55:34 UTC (rev 1898)
@@ -0,0 +1,205 @@
+% script plot_mpas_cross_sections
+
+% Specify data files, coordinates and text, then call functions
+% to find sections, load data, and plot cross-sections.
+
+% 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 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 = '/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_uReconstructMeridional,acc_uReconstructZonal, \
+% nAccumulate,latVertex,lonVertex,verticesOnEdge,edgesOnVertex,hZLevel,\
+% dvEdge,latCell,lonCell,cellsOnCell,nEdgesOnCell \
+% file_in.nc file_out.nc
+
+sim(1).dir = 'x5.NA.75km_15km';
+sim(1).netcdf_file = ['o.x5.NA.75km_15km.0027-01-01_00.00.00_acc_uReconstruct.nc'];
+
+sim(2).dir = 'x5.NA.75km_15km.leith';
+sim(2).netcdf_file = ['o.x5.NA.75km_15km.leith.0027-01-01_00.00.00_acc_uReconstruct.nc'];
+
+sim(3).dir = 'x1.15km';
+sim(3).netcdf_file = ['o.x1.15km.0018-11-03_00.00.00_acc_uReconstruct.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_uReconstruct.nc'];
+
+%clear sim
+%sim(1).dir='x1.120km';
+%sim(1).netcdf_file = 'output120km.0001-10-22_12:30: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
+  ];
+
+% plotDepth(nSections) depth to which to plot each section
+plotDepth = 4000*ones(1,size(coord,1));
+plotDepth(12:13) = 400;
+
+% 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 = 'nAtlantic_Zon_p1';
+page(1).sectionID = [1:4];
+page(2).name = 'nAtlantic_Zon_p2';
+page(2).sectionID = [5];
+
+page(3).name = 'nAtlantic_Mer_p1';
+page(3).sectionID = [6:9];
+page(4).name = 'nAtlantic_Mer_p2';
+page(4).sectionID = [10:11];
+                    
+page(5).name = 'eq_Pac_ACC';
+page(5).sectionID = [12:15];
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  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_conv_factor = [100 100 1]; % convert m/s to cm/s for velocities
+
+var_lims = [-10 10 2.5; -10 10 2.5; 0 20 2.5];
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  Begin main code.  Normally this does not need to change.
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%close all
+
+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="blue">']);
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  %
+  %  Find cells that connect beginning and end points of section
+  %
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+  [sim(iSim).sectionCellIndex, sim(iSim).nCellsInSection, sim(iSim).latSection,sim(iSim).lonSection, ...
+   depth, sim(iSim).latCellDeg,sim(iSim).lonCellDeg] ...
+     = find_cell_sections(wd,sim(iSim).dir,sim(iSim).netcdf_file, ...
+     sectionText,coord);
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  %
+  %  Plot cell section locations on world map
+  %
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+  sub_plot_cell_sections(sim(iSim).dir,coord, ...
+     sim(iSim).latCellDeg,sim(iSim).lonCellDeg, ...
+     sim(iSim).sectionCellIndex, sim(iSim).nCellsInSection,fid_latex)
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  %
+  %  Load large variables from netcdf file
+  %
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+  [sim(iSim).sectionData] = load_large_variables ...
+     (wd,sim(iSim).dir,sim(iSim).netcdf_file, var_name, var_conv_factor, ...
+      sim(iSim).sectionCellIndex, sim(iSim).nCellsInSection);
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  %
+  %  Plot data on cross-sections
+  %
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+  for iPage = 1:length(page)
+
+    sub_plot_cross_sections(sim(iSim).dir,sectionText, ...
+       page(iPage).name, page(iPage).sectionID,sim(iSim).sectionData,depth,...
+       sim(iSim).sectionCellIndex, sim(iSim).nCellsInSection, ...
+       sim(iSim).latSection,sim(iSim).lonSection,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="gray">']);
+  fclose(fid_latex);
+  dir = ['docs/' sim(iSim).netcdf_file '_dir/'];
+  filename = [sim(iSim).netcdf_file '.tex'];
+  unix(['cat mpas_sections.head.tex temp.tex &gt; ' dir filename ]);
+  unix(['mv f/*jpg ' dir 'f']);
+  cd(dir);
+  unix(['latex ' filename]);
+  cd('../..');
+end % iSim
+
+

Added: branches/tools/plot_mpas_cross_sections/sub_plot_cell_sections.m
===================================================================
--- branches/tools/plot_mpas_cross_sections/sub_plot_cell_sections.m                                (rev 0)
+++ branches/tools/plot_mpas_cross_sections/sub_plot_cell_sections.m        2012-05-11 15:55:34 UTC (rev 1898)
@@ -0,0 +1,110 @@
+function sub_plot_cell_sections(dir,coord, ...
+   latCellDeg,lonCellDeg, ...
+   sectionCellIndex, nCellsInSection,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
+% coord(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  
+% sectionCellIndex(maxCells,nSections)       cell index of each section
+% nCellsInSection(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(coord,1);
+
+figure(1); clf
+
+if (min(lonCellDeg)&gt;-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 cells.  This is just done for debugging.
+   %plot(lonCellDeg,latCellDeg,'.y')
+
+   grid on
+
+   for iSection=1:nSections
+     latCoordDeg = [coord(iSection,1) coord(iSection,3)];
+     lonCoordDeg = [coord(iSection,2) coord(iSection,4)];
+
+     h=plot([mod(lonCoordDeg,360)],[latCoordDeg],'*-');
+     set(h,'Color','y','LineWidth',1)
+
+     for i=1:nCellsInSection(iSection)-1
+        h = line([lonCellDeg(sectionCellIndex(i,iSection)) lonCellDeg(sectionCellIndex(i+1,iSection))],...
+                 [latCellDeg(sectionCellIndex(i,iSection)) latCellDeg(sectionCellIndex(i+1,iSection))]);
+        set(h,'Color',[1 0 0],'LineWidth',2)
+      end
+   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 ]);
+
+   filename=['f/' regexprep(dir,'\.','_') '_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/plot_mpas_cross_sections/sub_plot_cross_sections.m
===================================================================
--- branches/tools/plot_mpas_cross_sections/sub_plot_cross_sections.m                                (rev 0)
+++ branches/tools/plot_mpas_cross_sections/sub_plot_cross_sections.m        2012-05-11 15:55:34 UTC (rev 1898)
@@ -0,0 +1,109 @@
+function sub_plot_cross_sections(dir,sectionText,pageName,sectionID,sectionData,depth,...
+   sectionCellIndex, nCellsInSection, latSection,lonSection, coord, plotDepth,...
+   var_name,var_lims,fid_latex)
+
+% Plot data on cross-sections
+
+% 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,max(nCellsInSection),nSections,nVars)
+%   data in each cross-section for each variable
+% depth(nVertLevels)                         depth of center of each layer, for plotting
+% sectionCellIndex(maxCells,nSections)       cell index of each section
+% nCellsInSection(nSections)                 number of cells in each section
+% latSection(max(nCellsInSection),nSections) lat coordinates of each section
+% lonSection(max(nCellsInSection),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
+nSections = length(nCellsInSection);
+nVars = length(var_name);
+
+for iVar=1:nVars
+   figure(iVar+1); clf
+   set(gcf,'Position',[100+(iVar*100) 600-iVar*100   715   975])
+   temptext2=char(var_name(iVar));
+
+   for iRow = 1:length(sectionID)
+      iSection = sectionID(iRow);
+      if coord(iSection,1)==coord(iSection,3) % meridional section
+        xtext = 'longitude';
+        xaxis = lonSection(1:nCellsInSection(iSection),iSection);
+      else % zonal section
+        xtext = 'latitude';
+        xaxis = latSection(1:nCellsInSection(iSection),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
+        
+      ha=subplot('position',[px(2)-pw/2 py(iRow)-ph/2 pw ph]);
+      temptext = char(sectionText(iSection));
+      hold on
+        contour(xaxis, depth,sectionData(:,1:nCellsInSection(iSection),iSection,iVar), ...
+               [-30:var_lims(iVar,3):30]);
+        set(gca,'CLim',var_lims(iVar,1:2))
+        set(gca,'YDir','reverse')
+        title([temptext ', ' num2str(var_lims(iVar,3)) 'cm/s int'])
+      ylabel('depth, m')
+      xlabel(xtext)
+      axis tight
+      set(gca,'YLim',[0 plotDepth(iSection)])
+      grid on
+      h=colorbar;
+   
+   end
+
+   set(gcf,'PaperPositionMode','auto','color',[.8 1 .8], ...
+     'PaperPosition',[0.25 0.25 8 10])
+   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 ]);
+
+      filename=['f/' regexprep(dir,'\.','_') '_' pageName '_var' num2str(iVar)];
+      print('-djpeg',[filename '.jpg'])
+      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="blue">']);
+   %   print('-depsc2',[filename '.eps'])
+
+end
+

</font>
</pre>