<p><b>mpetersen@lanl.gov</b> 2012-05-24 14:06:44 -0600 (Thu, 24 May 2012)</p><p>Added C Shell script averaging_script to take averages of variables in netcdf files. See notes at top of averaging_script. Updated matlab code to work with the output file created by averaging_script.<br>
</p><hr noshade><pre><font color="gray">Added: branches/tools/plot_mpas_cross_sections/averaging_script
===================================================================
--- branches/tools/plot_mpas_cross_sections/averaging_script (rev 0)
+++ branches/tools/plot_mpas_cross_sections/averaging_script 2012-05-24 20:06:44 UTC (rev 1937)
@@ -0,0 +1,123 @@
+#!/bin/tcsh
+
+# This C Shell script uses nco commands to do the following:
+# 1. Create a small netcdf file that only contains relevant grid variables.
+# 2. Average timeslices within each output file.
+# 3. Average among output files created in step 2.
+# In step 2, output files already averaged will not be recomputed.
+# Step 3 does not use nAccumulate to weight the average. This is not a
+# problem if all output files use the same number of timeslices
+# (like 12 months in an annual file).
+#
+# The final product of this script is a single netcdf file named
+# total_avg_${lastFileName}
+# (for example, total_avg_o.x5.NA.75km_15km.0029-02-01_00.00.00.nc)
+# which contains grid information and all averaged variables, and
+# where ${lastFileName} is the final file used in averaging.
+
+# LANL notes:
+# This script may be run from an interactive (llogin) node or using the
+# supplied msub scripts.
+#
+# For the 15km (1.8M cells), I had to run on a mustang compute node
+# (64GB, versus 24 or 32 on the others), and the averaging step (step 2)
+# is limited to 16 timeslices in a file (18 fails). Processing a
+# 12-timeslice average took 8 minutes, so five years of monthly data
+# would require about 40 minutes per variable.
+#
+# For all other runs, including the x5.NA.37.5km_7.5km (1.0M cells), I
+# could run on conejo or lobo with no trouble, and it took 2 minutes
+# to average each 10-slice file for one variable.
+
+# Mark Petersen, MPAS-Ocean Team, LANL, May 2012
+
+########################################################################
+#
+# Specify variables and output files
+#
+########################################################################
+
+# variables
+set varList = (acc_uReconstructZonal acc_uReconstructMeridional)
+# To compute transport through sections, include acc_u as follows:
+#set varList = (acc_uReconstructZonal acc_uReconstructMeridional acc_u)
+
+# This is the list of output files that averages should be computed from.
+set outputFileList = o*00.nc
+
+########################################################################
+#
+# Machine specific paths
+#
+########################################################################
+
+# On LANL machines like lobo and mapache, loading pgi library will
+# load the correct libpgc.so.
+module purge
+module load pgi
+
+# Specify path of nco operators.
+# The following paths work on LANL turquoise machines.
+alias ncks '/usr/projects/cesm/software/conejo/nco/bin/ncks'
+alias ncra '/usr/projects/cesm/software/conejo/nco/bin/ncra'
+alias ncwa '/usr/projects/cesm/software/conejo/nco/bin/ncwa'
+alias ncea '/usr/projects/cesm/software/conejo/nco/bin/ncea'
+
+########################################################################
+#
+# Begin main code. Normally this does not need to change.
+#
+########################################################################
+
+# Find the last output file
+set lastFileName = `ls $outputFileList | tail -n 1`
+
+# Step 1. Create small grid file. This only needs to be done with one output file.'
+echo 'Create file with grid variables only:'
+set newFile = total_avg_${lastFileName}
+if ( -f $newFile ) then
+ echo " File $newFile exists. Skipping grid file generation."
+else
+ echo " File $newFile does not exists. Generating grid file..."
+
+ ncks -v latVertex,lonVertex,verticesOnEdge,edgesOnVertex,hZLevel,dvEdge,latCell,lonCell,cellsOnCell,nEdgesOnCell \
+ $lastFileName $newFile
+endif
+echo
+
+foreach varName ($varList)
+ echo 'Processing variable '$varName':'
+
+ foreach ncfile ($outputFileList)
+
+ # Step 2. Average timeslices within a single file
+ set newFile = ${ncfile}_${varName}_internal_avg.nc
+ if ( -f $newFile ) then
+ echo " File $newFile exists. Skipping internal averaging."
+ else
+ echo " File $newFile does not exists. Computing internal average..."
+
+ ncwa -y avg -v $varName -w nAccumulate -a Time $ncfile $newFile
+ ncwa -y ttl -v nAccumulate -w nAccumulate -a Time -A $ncfile $newFile
+ endif
+
+ end
+
+ # Step 3. Average variable among grid files.
+ # Note that this is not a weighted average. If there are a different number
+ # of timeslices per file, this will not weight the average correctly.
+ # I would like to change this to a weighted average, but don't know how at
+ # this point.
+
+ set finalFile = total_avg_${lastFileName}
+ echo " Computing total average among files. Appending to "$finalFile" (expect warnings for repeat dims)"
+ ncea -y avg -v $varName -A *_${varName}_internal_avg.nc $finalFile
+ ncea -y ttl -v nAccumulate -A *_${varName}_internal_avg.nc $finalFile
+ echo
+
+end
+
+echo "averaging_script complete!"
+
+
+
Property changes on: branches/tools/plot_mpas_cross_sections/averaging_script
___________________________________________________________________
Added: svn:executable
+ *
Modified: branches/tools/plot_mpas_cross_sections/load_large_variables.m
===================================================================
--- branches/tools/plot_mpas_cross_sections/load_large_variables.m 2012-05-24 19:05:33 UTC (rev 1936)
+++ branches/tools/plot_mpas_cross_sections/load_large_variables.m 2012-05-24 20:06:44 UTC (rev 1937)
@@ -34,7 +34,10 @@
[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'));
+
+% use only when averaging in matlab:
+%[dimname,nTimeSlices]= netcdf.inqDim(ncid,netcdf.inqDimID(ncid,'Time'));
+
nVars = length(var_name);
nSections = length(nCellsInSection);
@@ -48,13 +51,16 @@
% 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);
+
+ 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 i=1:nCellsInSection(iSection)
Added: branches/tools/plot_mpas_cross_sections/load_large_variables_avg.m
===================================================================
--- branches/tools/plot_mpas_cross_sections/load_large_variables_avg.m (rev 0)
+++ branches/tools/plot_mpas_cross_sections/load_large_variables_avg.m 2012-05-24 20:06:44 UTC (rev 1937)
@@ -0,0 +1,73 @@
+function [sectionData] = load_large_variables_avg ...
+ (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/msub_script
===================================================================
--- branches/tools/plot_mpas_cross_sections/msub_script (rev 0)
+++ branches/tools/plot_mpas_cross_sections/msub_script 2012-05-24 20:06:44 UTC (rev 1937)
@@ -0,0 +1,6 @@
+#! /bin/tcsh
+#MSUB -N nco_post_processing
+#MSUB -l walltime=2:00:00
+#MSUB -l nodes=1:ppn=1
+
+averaging_script
Modified: branches/tools/plot_mpas_cross_sections/plot_mpas_cross_sections.m
===================================================================
--- branches/tools/plot_mpas_cross_sections/plot_mpas_cross_sections.m 2012-05-24 19:05:33 UTC (rev 1936)
+++ branches/tools/plot_mpas_cross_sections/plot_mpas_cross_sections.m 2012-05-24 20:06:44 UTC (rev 1937)
@@ -34,21 +34,17 @@
% 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(1).netcdf_file = ['total_avg_o.x5.NA.75km_15km.0029-02-01_00.00.00.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(2).netcdf_file = ['total_avg_o.x5.NA.75km_15km.0027-01-01_00.00.00.nc'];
sim(3).dir = 'x1.15km';
-sim(3).netcdf_file = ['o.x1.15km.0018-11-03_00.00.00_acc_uReconstruct.nc'];
+sim(3).netcdf_file = ['total_avg_o.x1.15km.0018-11-03_00.00.00.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'];
+sim(4).netcdf_file = ['total_avg_o.x5.NA.37.5km_7.5km.0007-01-01_00.00.00.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
</font>
</pre>