Commit 3ac27449 authored by Pierre Cazenave's avatar Pierre Cazenave

Prepare for a new release.

parent d6a0c75b
Pierre Cazenave <pica@pml.ac.uk> fvcom-toolbox ChangeLog
20190220
--------
Minor update.
Thanks to Ricardo Torres and Simon Waldman for their contributions.
README.md
Updated with the current release details.
doc:
Updated for the latest release.
fvcom_postproc:
do_surface_plotMatlabMap
* Cosmetic changes to the plot.
show_max_CFL
* New function to compute the CFL number for a given grid.
fvcom_prepro:
add_obc_node_list
* Fix checks for invalid open boundaries.
add_sponge_nodes_list
* Code cleanup.
add_weights_FVCOM_nested_forcing
* Fix shape of the weights for the nodes.
change_shallow_bathy
* New function to adjust bathymetry in the grid.
find_nesting_region
* Code cleanup.
hybrid_coordinate
* Fix output format for the sigma file.
interp_POLCOMS2FVCOM
* Merge the v1 version of this and remove the duplicate.
make_default_nml
* Tweak the default number format for bed roughness.
read_fabm_variables
* Code cleanup.
read_sigma
* Add complete support for generalised and hyperbolic tangent sigma distributions.
read_sms_mesh
* Allow a given depth file to override values in the .2dm.
set_elevtide_tmd
* Add argument for the location of the TMD data.
write_FVCOM_river_nml
* Add argument to supply a vertical distribution as fractional depths (0-1).
write_FVCOM_tsobcERSEM
* Code cleanup.
utilities:
ComputeMatricRx1_nodes
* New function to calculate the hydrostatic consistency condition.
Times2Datetime
* New function to convert FVCOM 'Times' to MATLAB datetime objects.
concat_by_struct
* New function to concatenate data contained within a structure.
mjul2str
* Make the input always be a double.
plot_fvcom_field
* Add quiver plotting option.
read_FVCOM_river_file
* New function to read in FVCOM river files.
read_netCDF_FVCOM
* Fix loading data with only a single time stamp.
restrict_spatial_indices
* New function to eliminate nodes and elements for some given mask.
sigma_gen
* Add comment about potential bug (not yet fully investigated).
sigma_tanh
* Add option to return the sigma data into the supplied mesh object.
tubine_area_sigma
* Add option to plot into subplots.
20180201
--------
......
......@@ -22,6 +22,7 @@ Notes:
The PML version of the toolbox includes tagged releases, which can be downloaded as standalone (and thus relatively stable) versions. See the PML_ChangeLog.txt for details. Links to the direct downloads are:
- v20190220: https://github.com/pwcazenave/fvcom-toolbox/releases/tag/20190120
- v20180201: https://github.com/pwcazenave/fvcom-toolbox/releases/tag/20180201
- v20160811: https://github.com/pwcazenave/fvcom-toolbox/releases/tag/20160811
- v20160218: https://github.com/pwcazenave/fvcom-toolbox/releases/tag/20160218
......
......@@ -59,6 +59,6 @@ This function is called by:
0007 <span class="comment">% Revision history</span>
0008 <span class="comment">%</span>
0009 <span class="comment">%==============================================================================</span></pre></div>
<hr><address>Generated on Thu 01-Feb-2018 09:49:00 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
<hr><address>Generated on Wed 20-Feb-2019 16:06:01 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>
\ No newline at end of file
......@@ -114,6 +114,6 @@ This function is called by:
0061 <span class="comment">% Calculate direction and magnitude.</span>
0062 rDir=atan2(uDiff,vDiff)*(180/pi); <span class="comment">% in degrees.</span>
0063 rMag=sqrt(uDiff.^2+vDiff.^2)/tideDuration; <span class="comment">% in units/s.</span></pre></div>
<hr><address>Generated on Thu 01-Feb-2018 09:49:00 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
<hr><address>Generated on Wed 20-Feb-2019 16:06:01 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>
\ No newline at end of file
......@@ -121,6 +121,6 @@ This function is called by:
0076 <span class="keyword">end</span>
0077 <span class="keyword">end</span>
0078</pre></div>
<hr><address>Generated on Thu 01-Feb-2018 09:49:00 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
<hr><address>Generated on Wed 20-Feb-2019 16:06:01 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>
\ No newline at end of file
......@@ -168,6 +168,6 @@ This function is called by:
0094
0095 <span class="keyword">return</span>
0096</pre></div>
<hr><address>Generated on Thu 01-Feb-2018 09:49:00 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
<hr><address>Generated on Wed 20-Feb-2019 16:06:01 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>
\ No newline at end of file
......@@ -150,90 +150,89 @@ This function is called by:
0076 <span class="keyword">if</span> isfield (plotOPTS,<span class="string">'PlotoutS'</span>) &amp;&amp; ~isempty(plotOPTS.PlotoutS(plotOPTS.figure).handles)
0077 <span class="keyword">else</span>
0078 <span class="comment">% coast=load(plotOPTS.coastline_file);</span>
0079 geoshow(plotOPTS.coastline_file,<span class="string">'Color'</span>,<span class="string">'black'</span>)
0080
0081 <span class="keyword">end</span>
0082
0083 <span class="keyword">end</span>
0084 <span class="comment">%------------------------------------------------------------------------------</span>
0085 <span class="comment">% Generate maps at a give level</span>
0086 <span class="comment">%------------------------------------------------------------------------------</span>
0087 mstruct = gcm;
0079 geoshow(plotOPTS.coastline_file,<span class="string">'Facecolor'</span>,[0.7 .7 .7])
0080 <span class="keyword">end</span>
0081
0082 <span class="keyword">end</span>
0083 <span class="comment">%------------------------------------------------------------------------------</span>
0084 <span class="comment">% Generate maps at a give level</span>
0085 <span class="comment">%------------------------------------------------------------------------------</span>
0086 mstruct = gcm;
0087
0088
0089
0090 [X, Y] = mfwdtran(mstruct,plotOPTS.mesh.lat, plotOPTS.mesh.lon);
0091
0092 aa=plotOPTS.Time_record;
0093 fprintf(<span class="string">'Time step %i of %i\n'</span>,aa,length(FVCOM.Time_record));
0094 <span class="comment">%plot map</span>
0095 hold on
0096 ND=ndims(FVCOM.(plotOPTS.var_plot))
0097 <span class="keyword">switch</span> ND
0098 <span class="keyword">case</span> 2
0099 Plots(plotOPTS.figure).handles=patch(<span class="string">'Vertices'</span>,[X,Y],<span class="string">'Faces'</span>,plotOPTS.mesh.tri,<span class="keyword">...</span>
0100 <span class="string">'Cdata'</span>,squeeze(FVCOM.(plotOPTS.var_plot)(:,aa)),<span class="string">'edgecolor'</span>,<span class="string">'interp'</span>,<span class="string">'facecolor'</span>,<span class="string">'interp'</span>);
0101
0102 <span class="keyword">case</span> 3
0103 Plots(plotOPTS.figure).handles=patch(<span class="string">'Vertices'</span>,[X,Y],<span class="string">'Faces'</span>,plotOPTS.mesh.tri,<span class="keyword">...</span>
0104 <span class="string">'Cdata'</span>,squeeze(FVCOM.(plotOPTS.var_plot)(:,plotOPTS.nz_plot,aa)),<span class="string">'edgecolor'</span>,<span class="string">'interp'</span>,<span class="string">'facecolor'</span>,<span class="string">'interp'</span>);
0105 <span class="keyword">end</span>
0106 caxis(plotOPTS.clims)
0107 colorbar h
0108 <span class="comment">% check if mesh elements are required</span>
0109 <span class="keyword">if</span> plotOPTS.do_mesh
0110 <span class="comment">%plot vertices</span>
0111 Plots(plotOPTS.figure).mesh=patch(<span class="string">'Vertices'</span>,[X,Y],<span class="string">'Faces'</span>,plotOPTS.mesh.tri,<span class="keyword">...</span>
0112 <span class="string">'EdgeColor'</span>,[0.6 0.6 0.6],<span class="string">'FaceColor'</span>,<span class="string">'none'</span>);hold on
0113
0114 <span class="keyword">end</span>
0115 <span class="comment">% setm(gca,'MapLatLimit',plotOPTS.range_lat,'MapLonLimit',[plotOPTS.range_lon])</span>
0116 tightmap
0117 <span class="comment">%-----------------------------------------------------------------------</span>
0118 <span class="comment">% Only in my case it needs adding 6 because proj automatically determines a</span>
0119 <span class="comment">% reference long in strides of 6deg while m_map doesn't</span>
0120 <span class="comment">%------------------------------------------------------------------------</span>
0121 pause(plotOPTS.pause)
0122 <span class="comment">% if aa~=length(plotOPTS.Time_record)</span>
0123 <span class="comment">% delete(Plots(plotOPTS.figure).handles(:))</span>
0124 <span class="comment">% end</span>
0125
0126 <span class="comment">% end</span>
0127
0089 [X, Y] = mfwdtran(mstruct,plotOPTS.mesh.lat, plotOPTS.mesh.lon);
0090
0091 aa=plotOPTS.Time_record;
0092 fprintf(<span class="string">'Time step %i of %i\n'</span>,aa,length(FVCOM.Time_record));
0093 <span class="comment">%plot map</span>
0094 hold on
0095 ND=ndims(FVCOM.(plotOPTS.var_plot))
0096 <span class="keyword">switch</span> ND
0097 <span class="keyword">case</span> 2
0098 Plots(plotOPTS.figure).handles=patch(<span class="string">'Vertices'</span>,[X,Y],<span class="string">'Faces'</span>,plotOPTS.mesh.tri,<span class="keyword">...</span>
0099 <span class="string">'Cdata'</span>,squeeze(FVCOM.(plotOPTS.var_plot)(:,aa)),<span class="string">'edgecolor'</span>,<span class="string">'interp'</span>,<span class="string">'facecolor'</span>,<span class="string">'interp'</span>);
0100
0101 <span class="keyword">case</span> 3
0102 Plots(plotOPTS.figure).handles=patch(<span class="string">'Vertices'</span>,[X,Y],<span class="string">'Faces'</span>,plotOPTS.mesh.tri,<span class="keyword">...</span>
0103 <span class="string">'Cdata'</span>,squeeze(FVCOM.(plotOPTS.var_plot)(:,plotOPTS.nz_plot,aa)),<span class="string">'edgecolor'</span>,<span class="string">'interp'</span>,<span class="string">'facecolor'</span>,<span class="string">'interp'</span>);
0104 <span class="keyword">end</span>
0105 caxis(plotOPTS.clims)
0106 colorbar(<span class="string">'peer'</span>,gca,<span class="string">'EastOutside'</span>)
0107 <span class="comment">% check if mesh elements are required</span>
0108 <span class="keyword">if</span> plotOPTS.do_mesh
0109 <span class="comment">%plot vertices</span>
0110 Plots(plotOPTS.figure).mesh=patch(<span class="string">'Vertices'</span>,[X,Y],<span class="string">'Faces'</span>,plotOPTS.mesh.tri,<span class="keyword">...</span>
0111 <span class="string">'EdgeColor'</span>,[0.6 0.6 0.6],<span class="string">'FaceColor'</span>,<span class="string">'none'</span>);hold on
0112
0113 <span class="keyword">end</span>
0114 <span class="comment">% setm(gca,'MapLatLimit',plotOPTS.range_lat,'MapLonLimit',[plotOPTS.range_lon])</span>
0115 tightmap
0116 <span class="comment">%-----------------------------------------------------------------------</span>
0117 <span class="comment">% Only in my case it needs adding 6 because proj automatically determines a</span>
0118 <span class="comment">% reference long in strides of 6deg while m_map doesn't</span>
0119 <span class="comment">%------------------------------------------------------------------------</span>
0120 pause(plotOPTS.pause)
0121 <span class="comment">% if aa~=length(plotOPTS.Time_record)</span>
0122 <span class="comment">% delete(Plots(plotOPTS.figure).handles(:))</span>
0123 <span class="comment">% end</span>
0124
0125 <span class="comment">% end</span>
0126
0127
0128
0129
0130
0131
0132
0133
0134 <span class="comment">% fprintf('Time step %i of %i\n',aa,length(plotOPTS.Time_record));</span>
0135 <span class="comment">% %plot map</span>
0136 <span class="comment">% hold on</span>
0137 <span class="comment">% Plots(plotOPTS.figure).handles=patch('Vertices',[lat,lon],'Faces',plotOPTS.mesh.tri,...</span>
0138 <span class="comment">% 'Cdata',squeeze(FVCOM.(plotOPTS.var_plot)(:,plotOPTS.nz_plot,aa)),'edgecolor','interp','facecolor','interp');</span>
0139 <span class="comment">%</span>
0140 <span class="comment">% caxis(plotOPTS.clims)</span>
0141 <span class="comment">% colorbar</span>
0142 <span class="comment">% % check if mesh elements are required</span>
0143 <span class="comment">% if plotOPTS.do_mesh</span>
0144 <span class="comment">% %plot vertices</span>
0145 <span class="comment">% Plots(plotOPTS.figure).mesh=patch('Vertices',[lon,lat],'Faces',plotOPTS.mesh.tri,...</span>
0146 <span class="comment">% 'EdgeColor',[0.6 0.6 0.6],'FaceColor','none');hold on</span>
0147 <span class="comment">%</span>
0148 <span class="comment">% end</span>
0149 <span class="comment">% setm(gca,'MapLatLimit',plotOPTS.range_lat,'MapLonLimit',[plotOPTS.range_lon])</span>
0150 <span class="comment">% %-----------------------------------------------------------------------</span>
0151 <span class="comment">% % Only in my case it needs adding 6 because proj automatically determines a</span>
0152 <span class="comment">% % reference long in strides of 6deg while m_map doesn't</span>
0153 <span class="comment">% %------------------------------------------------------------------------</span>
0154 <span class="comment">% pause(.2)</span>
0155 <span class="comment">% if aa~=length(plotOPTS.Time_record)</span>
0156 <span class="comment">% delete(Plots(plotOPTS.figure).handles(:))</span>
0157 <span class="comment">% end</span>
0158 <span class="comment">%</span>
0159 <span class="comment">% end</span>
0160 <span class="comment">%</span>
0161 <span class="keyword">return</span>
0162</pre></div>
<hr><address>Generated on Thu 01-Feb-2018 09:49:00 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
0133 <span class="comment">% fprintf('Time step %i of %i\n',aa,length(plotOPTS.Time_record));</span>
0134 <span class="comment">% %plot map</span>
0135 <span class="comment">% hold on</span>
0136 <span class="comment">% Plots(plotOPTS.figure).handles=patch('Vertices',[lat,lon],'Faces',plotOPTS.mesh.tri,...</span>
0137 <span class="comment">% 'Cdata',squeeze(FVCOM.(plotOPTS.var_plot)(:,plotOPTS.nz_plot,aa)),'edgecolor','interp','facecolor','interp');</span>
0138 <span class="comment">%</span>
0139 <span class="comment">% caxis(plotOPTS.clims)</span>
0140 <span class="comment">% colorbar</span>
0141 <span class="comment">% % check if mesh elements are required</span>
0142 <span class="comment">% if plotOPTS.do_mesh</span>
0143 <span class="comment">% %plot vertices</span>
0144 <span class="comment">% Plots(plotOPTS.figure).mesh=patch('Vertices',[lon,lat],'Faces',plotOPTS.mesh.tri,...</span>
0145 <span class="comment">% 'EdgeColor',[0.6 0.6 0.6],'FaceColor','none');hold on</span>
0146 <span class="comment">%</span>
0147 <span class="comment">% end</span>
0148 <span class="comment">% setm(gca,'MapLatLimit',plotOPTS.range_lat,'MapLonLimit',[plotOPTS.range_lon])</span>
0149 <span class="comment">% %-----------------------------------------------------------------------</span>
0150 <span class="comment">% % Only in my case it needs adding 6 because proj automatically determines a</span>
0151 <span class="comment">% % reference long in strides of 6deg while m_map doesn't</span>
0152 <span class="comment">% %------------------------------------------------------------------------</span>
0153 <span class="comment">% pause(.2)</span>
0154 <span class="comment">% if aa~=length(plotOPTS.Time_record)</span>
0155 <span class="comment">% delete(Plots(plotOPTS.figure).handles(:))</span>
0156 <span class="comment">% end</span>
0157 <span class="comment">%</span>
0158 <span class="comment">% end</span>
0159 <span class="comment">%</span>
0160 <span class="keyword">return</span>
0161</pre></div>
<hr><address>Generated on Wed 20-Feb-2019 16:06:01 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>
\ No newline at end of file
......@@ -110,6 +110,6 @@ This function is called by:
0066
0067 <span class="keyword">return</span>
0068</pre></div>
<hr><address>Generated on Thu 01-Feb-2018 09:49:00 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
<hr><address>Generated on Wed 20-Feb-2019 16:06:01 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>
\ No newline at end of file
......@@ -118,6 +118,6 @@ This function is called by:
0074
0075 <span class="keyword">return</span>
0076</pre></div>
<hr><address>Generated on Thu 01-Feb-2018 09:49:00 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
<hr><address>Generated on Wed 20-Feb-2019 16:06:01 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>
\ No newline at end of file
......@@ -186,6 +186,6 @@ This function is called by:
0111 plot(X,Y,<span class="string">'r'</span>,<span class="string">'LineWidth'</span>,2.5)
0112 <span class="keyword">end</span>
0113</pre></div>
<hr><address>Generated on Thu 01-Feb-2018 09:49:00 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
<hr><address>Generated on Wed 20-Feb-2019 16:06:01 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>
\ No newline at end of file
......@@ -194,6 +194,6 @@ This function is called by:
0120 delete(Plots(plotOPTS.figure).handles(:))
0121 <span class="keyword">end</span>
0122 <span class="keyword">end</span></pre></div>
<hr><address>Generated on Thu 01-Feb-2018 09:49:00 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
<hr><address>Generated on Wed 20-Feb-2019 16:06:01 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>
\ No newline at end of file
......@@ -238,6 +238,6 @@ This function is called by:
0164 <span class="comment">% delete(Plots(plotOPTS.figure).handles(:))</span>
0165 <span class="comment">% end</span>
0166 <span class="comment">% end</span></pre></div>
<hr><address>Generated on Thu 01-Feb-2018 09:49:00 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
<hr><address>Generated on Wed 20-Feb-2019 16:06:01 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>
\ No newline at end of file
......@@ -224,6 +224,6 @@ This function is called by:
0150 <span class="comment">% delete(Plots(plotOPTS.figure).handles(:))</span>
0151 <span class="comment">% end</span>
0152 <span class="comment">% end</span></pre></div>
<hr><address>Generated on Thu 01-Feb-2018 09:49:00 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
<hr><address>Generated on Wed 20-Feb-2019 16:06:01 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>
\ No newline at end of file
......@@ -184,6 +184,6 @@ This function is called by:
0120 <span class="comment">% vector plots requires correct lat and lon for u and v positions: FVCOM.xc</span>
0121 <span class="comment">% and FVCOM.yc. Remember to extract them</span>
0122 PLotoutV=<a href="do_vector_plot.html" class="code" title="function [Plots]=do_vector_plot(plotOPTS,FVCOM)">do_vector_plot</a>(plotOPTS,FVCOM)</pre></div>
<hr><address>Generated on Thu 01-Feb-2018 09:49:00 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
<hr><address>Generated on Wed 20-Feb-2019 16:06:01 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>
\ No newline at end of file
......@@ -16,4 +16,5 @@ digraph m2html {
do_vector_plot_MatlabMap [URL="do_vector_plot_MatlabMap.html"];
do_vector_plot_MatlabMapC [URL="do_vector_plot_MatlabMapC.html"];
example_surface_plot [URL="example_surface_plot.html"];
show_max_CFL [URL="show_max_CFL.html"];
}
\ No newline at end of file
......@@ -34,6 +34,6 @@
</map>
</center>
<hr><address>Generated on Thu 01-Feb-2018 09:48:59 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
<hr><address>Generated on Wed 20-Feb-2019 16:06:00 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>
\ No newline at end of file
......@@ -19,7 +19,7 @@
<h2>Matlab files in this directory:</h2>
<table>
<tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="README.html">README</a></td><td>README for FVCOM_postproc </td></tr><tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="do_residual.html">do_residual</a></td><td>DO_RESIDUAL Takes the u and v vectors of a model output and calculates </td></tr><tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="do_residual_plot.html">do_residual_plot</a></td><td>Take the output of do_residual and plot as a vector figure. Summarises a </td></tr><tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="do_surface_plot.html">do_surface_plot</a></td><td> </td></tr><tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="do_surface_plotMatlabMap.html">do_surface_plotMatlabMap</a></td><td> </td></tr><tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="do_surface_plotVel.html">do_surface_plotVel</a></td><td>reads image and plots tracks or stations </td></tr><tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="do_surface_plotVelMatlabMap.html">do_surface_plotVelMatlabMap</a></td><td>reads image and plots tracks or stations </td></tr><tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="do_transect_plot.html">do_transect_plot</a></td><td> </td></tr><tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="do_vector_plot.html">do_vector_plot</a></td><td> </td></tr><tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="do_vector_plot_MatlabMap.html">do_vector_plot_MatlabMap</a></td><td> </td></tr><tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="do_vector_plot_MatlabMapC.html">do_vector_plot_MatlabMapC</a></td><td> </td></tr><tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="example_surface_plot.html">example_surface_plot</a></td><td>Sample script to extract and generate m_map contours of tracer variables </td></tr></table>
<tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="README.html">README</a></td><td>README for FVCOM_postproc </td></tr><tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="do_residual.html">do_residual</a></td><td>DO_RESIDUAL Takes the u and v vectors of a model output and calculates </td></tr><tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="do_residual_plot.html">do_residual_plot</a></td><td>Take the output of do_residual and plot as a vector figure. Summarises a </td></tr><tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="do_surface_plot.html">do_surface_plot</a></td><td> </td></tr><tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="do_surface_plotMatlabMap.html">do_surface_plotMatlabMap</a></td><td> </td></tr><tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="do_surface_plotVel.html">do_surface_plotVel</a></td><td>reads image and plots tracks or stations </td></tr><tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="do_surface_plotVelMatlabMap.html">do_surface_plotVelMatlabMap</a></td><td>reads image and plots tracks or stations </td></tr><tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="do_transect_plot.html">do_transect_plot</a></td><td> </td></tr><tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="do_vector_plot.html">do_vector_plot</a></td><td> </td></tr><tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="do_vector_plot_MatlabMap.html">do_vector_plot_MatlabMap</a></td><td> </td></tr><tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="do_vector_plot_MatlabMapC.html">do_vector_plot_MatlabMapC</a></td><td> </td></tr><tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="example_surface_plot.html">example_surface_plot</a></td><td>Sample script to extract and generate m_map contours of tracer variables </td></tr><tr><td><img src="../matlabicon.gif" alt="" border="">&nbsp;<a href="show_max_CFL.html">show_max_CFL</a></td><td>SHOW_MAX-CFL Function to find the max CFL encountered in each mesh element during an FVCOM model </td></tr></table>
<h2>Subsequent directories:</h2>
......@@ -29,6 +29,6 @@
<ul style="list-style-image:url(../simulinkicon.gif)">
<li>View the <a href="graph.html">Graph</a>.</li>
</ul>
<hr><address>Generated on Thu 01-Feb-2018 09:48:59 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
<hr><address>Generated on Wed 20-Feb-2019 16:06:00 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>
\ No newline at end of file
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
"http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
<title>Description of show_max_CFL</title>
<meta name="keywords" content="show_max_CFL">
<meta name="description" content="SHOW_MAX-CFL Function to find the max CFL encountered in each mesh element during an FVCOM model">
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
<meta name="generator" content="m2html v1.5 &copy; 2003-2005 Guillaume Flandin">
<meta name="robots" content="index, follow">
<link type="text/css" rel="stylesheet" href="../m2html.css">
</head>
<body>
<a name="_top"></a>
<div><a href="../index.html">Home</a> &gt; <a href="index.html">fvcom_postproc</a> &gt; show_max_CFL.m</div>
<!--<table width="100%"><tr><td align="left"><a href="../index.html"><img alt="<" border="0" src="../left.png">&nbsp;Master index</a></td>
<td align="right"><a href="index.html">Index for fvcom_postproc&nbsp;<img alt=">" border="0" src="../right.png"></a></td></tr></table>-->
<h1>show_max_CFL
</h1>
<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="box"><strong>SHOW_MAX-CFL Function to find the max CFL encountered in each mesh element during an FVCOM model</strong></div>
<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="box"><strong>function [ maxCFLs, fighandle ] = show_max_CFL(grdFile, depFile, ncFile, extTS, coordtype, fig_flag) </strong></div>
<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="fragment"><pre class="comment">SHOW_MAX-CFL Function to find the max CFL encountered in each mesh element during an FVCOM model
run. NB need enough RAM to load all u, v and zeta data from the given
ncFile.
Inputs: grdFile, depFile: casename_grd.dat and casename_dep.dat files
from an FVCOM model
ncFile : Output from FVCOM. Must include u, v, and zeta.
extTS : External timestep at which to calculate CFL, in
seconds.
coordtype: 'lonlat' or 'cartesian'. If cartesian, coordinates
should be in metres. If lonlat, they should be in
degrees.
fig_flag: Should a figure be plotted? true or false.</pre></div>
<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../matlabicon.gif)">
</ul>
This function is called by:
<ul style="list-style-image:url(../matlabicon.gif)">
</ul>
<!-- crossreference -->
<h2><a name="_subfunctions"></a>SUBFUNCTIONS <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<ul style="list-style-image:url(../matlabicon.gif)">
<li><a href="#_sub1" class="code">function [distm]=haversine(lat1,lon1,lat2,lon2)</a></li></ul>
<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function [ maxCFLs, fighandle ] = show_max_CFL(grdFile, depFile, ncFile, extTS, coordtype, fig_flag)</a>
0002 <span class="comment">%SHOW_MAX-CFL Function to find the max CFL encountered in each mesh element during an FVCOM model</span>
0003 <span class="comment">%run. NB need enough RAM to load all u, v and zeta data from the given</span>
0004 <span class="comment">%ncFile.</span>
0005 <span class="comment">% Inputs: grdFile, depFile: casename_grd.dat and casename_dep.dat files</span>
0006 <span class="comment">% from an FVCOM model</span>
0007 <span class="comment">% ncFile : Output from FVCOM. Must include u, v, and zeta.</span>
0008 <span class="comment">% extTS : External timestep at which to calculate CFL, in</span>
0009 <span class="comment">% seconds.</span>
0010 <span class="comment">% coordtype: 'lonlat' or 'cartesian'. If cartesian, coordinates</span>
0011 <span class="comment">% should be in metres. If lonlat, they should be in</span>
0012 <span class="comment">% degrees.</span>
0013 <span class="comment">% fig_flag: Should a figure be plotted? true or false.</span>
0014 <span class="comment">%</span>
0015
0016 <span class="comment">% The formulation of CFL used here is taken from the MIKE 3 manual</span>
0017
0018 <span class="comment">% Simon Waldman (Marine Scotland Science / Heriot-Watt University), March 2018</span>
0019
0020 <span class="comment">%% Check inputs</span>
0021
0022 assert( nargin == 6, <span class="string">'Wrong number of arguments.'</span>);
0023 assert( exist(grdFile, <span class="string">'file'</span>)==2, <span class="string">'Cannot find file %s'</span>, grdFile );
0024 assert( exist(depFile, <span class="string">'file'</span>)==2, <span class="string">'Cannot find file %s'</span>, depFile );
0025 assert( exist(ncFile, <span class="string">'file'</span>)==2, <span class="string">'Cannot find file %s'</span>, ncFile );
0026 assert( ischar(coordtype), <span class="string">'Coord type must be ''lonlat'' or ''cartesian''.'</span> );
0027 assert( islogical( fig_flag ), <span class="string">'fig_flag should be logical.'</span> );
0028
0029 <span class="keyword">switch</span> coordtype
0030 <span class="keyword">case</span> <span class="string">'lonlat'</span>
0031 lonlat=true;
0032 <span class="keyword">case</span> <span class="string">'latlon'</span>
0033 lonlat=true;
0034 <span class="keyword">case</span> <span class="string">'cartesian'</span>
0035 lonlat=false;
0036 <span class="keyword">otherwise</span>
0037 error(<span class="string">'Coord type must be ''lonlat'' or ''cartesian''.'</span>);
0038 <span class="keyword">end</span>
0039
0040 <span class="comment">% does the ncFile have velocities &amp; water levels?</span>
0041
0042 info = ncinfo(ncFile);
0043 assert( any(strcmp(<span class="string">'u'</span>, {info.Variables.Name})) &amp;&amp; any(strcmp(<span class="string">'v'</span>, {info.Variables.Name})) &amp;&amp; <span class="keyword">...</span>
0044 any(strcmp(<span class="string">'zeta'</span>, {info.Variables.Name})), <span class="string">'netCDF file must include the fields u, v and zeta.'</span> );
0045
0046 <span class="comment">%% Do the stuff.</span>
0047
0048 disp(<span class="string">'Reading mesh &amp; bathymetry...'</span>);
0049 M = read_fvcom_mesh( grdFile ); <span class="comment">%NB this function puts lon and lat in the x and y fields.</span>
0050 M.h = read_fvcom_bath( depFile );
0051 <span class="comment">% calculate element depths from mean of nodes</span>
0052 M.hc = mean(M.h(M.tri),2);
0053
0054 disp(<span class="string">'Reading U velocities...'</span>);
0055 U = ncread(ncFile, <span class="string">'u'</span>);
0056 disp(<span class="string">'Reading V velocities...'</span>);
0057 V = ncread(ncFile, <span class="string">'v'</span>);
0058 disp(<span class="string">'Reading surface elevations...'</span>);
0059 Z = ncread(ncFile, <span class="string">'zeta'</span>);
0060
0061 g = 9.81;
0062
0063 NumTS = size(U, 3);
0064 NumEl = size(U, 1);
0065
0066 <span class="comment">% Find water depths for each element at each TS. First calculate this for</span>
0067 <span class="comment">% nodes, then average them for elements.</span>
0068 disp(<span class="string">'Calculating depths...'</span>);
0069 NodeDepths = repmat(M.h, 1, NumTS) + Z; <span class="comment">%repmat because Z has a time dimension, M.h doesn't.</span>
0070 <span class="keyword">for</span> n = 1:3
0071 tmp(:,:,n) = NodeDepths(M.tri(:,n), :);
0072 <span class="keyword">end</span>
0073 ElDepths = mean(tmp, 3); <span class="comment">%this has dimensions of element x TS.</span>
0074 clear Z NodeDepths tmp; <span class="comment">%save some memory</span>
0075
0076 <span class="comment">% Find minimum characteristic length for each element. This is approximated by the</span>
0077 <span class="comment">% minimum edge length. It could be shorter with really long thin triangles,</span>
0078 <span class="comment">% but the 30 degree internal angle minimum for FVCOM means we shouldn't have those.</span>
0079
0080 disp(<span class="string">'Calculating triangle sizes...'</span>);
0081 CharLen = nan(NumEl, 1);
0082 <span class="keyword">for</span> e = 1:NumEl <span class="comment">%for each element</span>
0083 xv = M.x(M.tri(e,:)); <span class="comment">%vertices.</span>
0084 yv = M.y(M.tri(e,:));
0085 <span class="comment">%close the triangle by copying the first vertex to the end</span>
0086 xv(4) = xv(1);
0087 yv(4) = yv(1);
0088 <span class="comment">%find the edge lengths</span>
0089 <span class="keyword">if</span> lonlat
0090 <span class="keyword">for</span> a = 1:3
0091 edges(a) = <a href="#_sub1" class="code" title="subfunction [distm]=haversine(lat1,lon1,lat2,lon2)">haversine</a>(yv(a), xv(a), yv(a+1), xv(a+1));
0092 <span class="keyword">end</span>
0093 <span class="keyword">else</span>
0094 edges = sqrt( diff(xv).^2 + diff(yv).^2 );
0095 <span class="keyword">end</span>
0096 CharLen(e) = min(edges);
0097 <span class="keyword">end</span>
0098
0099 <span class="comment">% For each element and TS, find the layer with the highest U and V magnitudes.</span>
0100 <span class="comment">% Technically doing this with each component rather than the vector magnitude</span>
0101 <span class="comment">% is wrong, but it'll usually be close and where wrong, it'll overestimate the CFL, so it's safe.</span>
0102 maxU = squeeze(max(abs(U), [], 2));
0103 maxV = squeeze(max(abs(V), [], 2));
0104
0105 <span class="comment">% Find CFL for each element at each TS</span>
0106 CFL = ( 2 .* sqrt( g .* ElDepths ) + maxU + maxV ) .* repmat( (extTS ./ CharLen), 1, NumTS );
0107 <span class="comment">%This is based on equation 6.1 on pg 33 of the MIKE hydrodynamic module</span>
0108 <span class="comment">%manual (modified for using a single characteristic length rather than</span>
0109 <span class="comment">%deltaX/deltaY)</span>
0110
0111 <span class="comment">% find the max over time for each element</span>
0112 MaxCFL = max(CFL, [], 2);
0113
0114 [val, I] = max(MaxCFL);
0115 fprintf(<span class="string">'Max CFL reached with an external timestep of %.2f secs is approx. %.3f, in Element %i.\n'</span>, extTS, val, I);
0116
0117 <span class="comment">% find how long the timestep probably could go. We set the CFL to 0.8 and</span>
0118 <span class="comment">% apply the same formula backwards.</span>
0119 TargetCFL = 0.8;
0120 MaxTSs = repmat( (TargetCFL .* CharLen), 1, NumTS ) ./ ( 2 .* sqrt( g .* ElDepths ) + maxU + maxV );