Commit 00c69b0e authored by Pierre Cazenave's avatar Pierre Cazenave
Browse files

Merge branch 'dev' of gitlab.ecosystem-modelling.pml.ac.uk:fvcom/fvcom-toolbox into dev

parents e8e27f01 a307264a
......@@ -175,6 +175,7 @@ nemo.o = (nemo.o / 16) *1000 ./ tmp; % Nemo oxygen concentrations are for O rath
nemo.p = (nemo.p / 35.5)*1000 ./ tmp;%g/s to mmol/m3
nemo.sio3 = (nemo.sio3 / 28) *1000./ tmp;%g/2 to mmol/m3
nemo.bioalk= nemo.bioalk./ tmp / 1000; % bioalk is in umol/s need umol/kg
nemo.dic= nemo.dic./12./ tmp *1000; % dic is in gC/s need mmol/m3
% total alkalinity is already in umol/Kg as expected by ERSEM.
clear tmp
......
......@@ -193,45 +193,45 @@ Times_varid = netcdf.defVar(nc, 'Times' ,'NC_CHAR', ...
[datestrlen_dimid, time_dimid]);
netcdf.putAtt(nc, Times_varid, 'time_zone', 'UTC');
x_varid = netcdf.defVar(nc, 'x', 'double', ...
x_varid = netcdf.defVar(nc, 'x', 'NC_FLOAT', ...
node_dimid);
netcdf.putAtt(nc, x_varid, 'units', 'meters');
netcdf.putAtt(nc, x_varid, 'long_name', 'nodal x-coordinate');
y_varid = netcdf.defVar(nc, 'y', 'double', ...
y_varid = netcdf.defVar(nc, 'y', 'NC_FLOAT', ...
node_dimid);
netcdf.putAtt(nc, y_varid, 'units', 'meters');
netcdf.putAtt(nc, y_varid, 'long_name', 'nodal y-coordinate');
xc_varid = netcdf.defVar(nc, 'xc', 'double', ...
xc_varid = netcdf.defVar(nc, 'xc', 'NC_FLOAT', ...
elem_dimid);
netcdf.putAtt(nc, xc_varid, 'units', 'meters');
netcdf.putAtt(nc, xc_varid, 'long_name', 'zonal x-coordinate');
yc_varid = netcdf.defVar(nc, 'yc', 'double', ...
yc_varid = netcdf.defVar(nc, 'yc', 'NC_FLOAT', ...
elem_dimid);
netcdf.putAtt(nc, yc_varid, 'units', 'meters');
netcdf.putAtt(nc, yc_varid, 'long_name', 'zonal y-coordinate');
lon_varid = netcdf.defVar(nc, 'lon', 'double', ...
lon_varid = netcdf.defVar(nc, 'lon', 'NC_FLOAT', ...
node_dimid);
netcdf.putAtt(nc, lon_varid, 'units', 'degrees_east');
netcdf.putAtt(nc, lon_varid, 'standard_name', 'longitude');
netcdf.putAtt(nc, lon_varid, 'long_name', 'nodal longitude');
lat_varid = netcdf.defVar(nc, 'lat', 'double', ...
lat_varid = netcdf.defVar(nc, 'lat', 'NC_FLOAT', ...
node_dimid);
netcdf.putAtt(nc, lat_varid, 'units', 'degrees_north');
netcdf.putAtt(nc, lat_varid, 'standard_name', 'latitude');
netcdf.putAtt(nc, lat_varid, 'long_name', 'nodal latitude');
lonc_varid = netcdf.defVar(nc, 'lonc', 'double', ...
lonc_varid = netcdf.defVar(nc, 'lonc', 'NC_FLOAT', ...
elem_dimid);
netcdf.putAtt(nc, lonc_varid, 'units', 'degrees_east');
netcdf.putAtt(nc, lonc_varid, 'standard_name', 'longitude');
netcdf.putAtt(nc, lonc_varid, 'long_name', 'zonal longitude');
latc_varid = netcdf.defVar(nc, 'latc', 'double', ...
latc_varid = netcdf.defVar(nc, 'latc', 'NC_FLOAT', ...
elem_dimid);
netcdf.putAtt(nc, latc_varid, 'units', 'degrees_north');
netcdf.putAtt(nc, latc_varid, 'standard_name', 'latitude');
......
function write_FVCOM_tsobcERSEM(basename,time,nSiglay,Mobj,ERSEMdata,varargin)
function write_FVCOM_tsobcERSEM(tsOBCFile,FileExists,time,nSiglay,Mobj,ERSEMdata,varargin)
% Export temperature and salinity forcing at the open boundary.
%
% function write_FVCOM_tsobc(basename,time,nSiglay,in_temp,in_salt)
......@@ -9,17 +9,18 @@ function write_FVCOM_tsobcERSEM(basename,time,nSiglay,Mobj,ERSEMdata,varargin)
% sigma_layers, time).
%
% INPUT
% basename - Model case name (to find the bathymetry and open boundary
% .dat files).
% tsOBCFile - Filename to use as output
% FileExists - True if tsOBCFile already exists and we only want to add
% ERSEM variables to it.
% time - Time (Modified Julian Days)
% nSiglay - Number of sigma layers
% in_temp - Boundary temperature (Celsius)
% in_salt - Boundary salinity (psu)
% Mobj - Mesh Object with the following fields:
% - nObs - number of open boundaries
% - read_obc_nodes - open boundary node cell array (length = nObs)
% - siglay - sigma layer definitions
% - siglev - sigma level definitions
% ERSEMdata - struct with information to output to Netcdf File
% Optional keyword-argument pairs:
%
% 'strtime' = set to true to output the 'Times' variable
......@@ -68,23 +69,15 @@ function write_FVCOM_tsobcERSEM(basename,time,nSiglay,Mobj,ERSEMdata,varargin)
% when converting from double to single format.
% RJT Revision history
% 2013-12-05 Added functionality to output ERSEM nutrients
% 2017-04-10 Option to write to an existing file... i.e. a Nesting file
%==============================================================================
[~, subname] = fileparts(mfilename('fullpath'));
global ftbverbose
if ftbverbose
fprintf('\nbegin : %s\n', subname)
end
if nargin == 7
NNuts=length(ERSEMdata)
if NNuts==0
disp(['ERSEM fields are empty! no nutrients to be output!!!'])
else
disp(['We are adding nutrients to the BC file!!!',ERSEMdata(1:NNuts).name])
tsOBCFile = [basename, '_ERSEMobc.nc'];
end
fprintf('\nbegin : %s\n', subname)
end
NNuts=length(ERSEMdata)
% Default to string times as FVCOM looks for these first.
strtime = true;
inttime = false;
......@@ -137,67 +130,95 @@ end
%--------------------------------------------------------------------------
% Set netCDF variables and dump to file
%--------------------------------------------------------------------------
% open boundary forcing
nc = netcdf.create(tsOBCFile, 'clobber');
% define global attributes
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'title','Open boundary ERSEM nudging')
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'type','FVCOM TIME SERIES OBC FABM FILE')
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'history', sprintf('File created with %s from the MATLAB fvcom-toolbox', subname))
% define dimensions
nobc_dimid=netcdf.defDim(nc,'nobc',nObc);
datestrlen_dimid=netcdf.defDim(nc,'DateStrLen',26);
time_dimid=netcdf.defDim(nc,'time',netcdf.getConstant('NC_UNLIMITED'));
siglay_dimid=netcdf.defDim(nc,'siglay',nSiglay);
siglev_dimid=netcdf.defDim(nc,'siglev',nSiglev);
% variables
if strtime
Times_varid=netcdf.defVar(nc,'Times','NC_CHAR',[datestrlen_dimid, time_dimid]);
netcdf.putAtt(nc,Times_varid,'time_zone','UTC');
end
if floattime
time_varid=netcdf.defVar(nc,'time','NC_FLOAT',time_dimid);
netcdf.putAtt(nc,time_varid,'long_name','time');
netcdf.putAtt(nc,time_varid,'units','days since 1858-11-17 00:00:00');
netcdf.putAtt(nc,time_varid,'time_zone','UTC');
end
if inttime
itime_varid=netcdf.defVar(nc,'Itime','NC_INT',time_dimid);
netcdf.putAtt(nc,itime_varid,'units','days since 1858-11-17 00:00:00');
netcdf.putAtt(nc,itime_varid,'format','modified julian day (MJD)');
netcdf.putAtt(nc,itime_varid,'time_zone','UTC');
itime2_varid=netcdf.defVar(nc,'Itime2','NC_INT',time_dimid);
netcdf.putAtt(nc,itime2_varid,'units','msec since 00:00:00');
netcdf.putAtt(nc,itime2_varid,'time_zone','UTC');
if FileExists
% open boundary forcing
nc = netcdf.open(tsOBCFile, 'WRITE');
% read dimensions from the
% define dimensions
dimids = netcdf.inqDimIDs(nc);
for dd =1:length(dimids)
dimidname=netcdf.inqDim(nc,dimids(dd));
switch dimidname
case{'node'}
nobc_dimid=netcdf.inqDimID(nc,dimidname);
case{'DateStrLen'}
datestrlen_dimid=netcdf.inqDimID(nc,dimidname);
case{'time'}
time_dimid=netcdf.inqDimID(nc,dimidname);
% read time in case is different from data in file
file_times = netcdf.getVar(nc,netcdf.inqVarID(nc,'Times'))';
file_timem = datenum(file_times,'yyyy-mm-dd HH:MM:SS.FFF');
case{'siglay'}
siglay_dimid=netcdf.inqDimID(nc,dimidname);
case{'siglev'}
siglev_dimid=netcdf.inqDimID(nc,dimidname);
end
end
netcdf.reDef(nc)
else
% open boundary forcing
nc = netcdf.create(tsOBCFile, 'clobber');
% define global attributes
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'title','Open boundary ERSEM nudging')
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'type','FVCOM TIME SERIES OBC FABM FILE')
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'history', sprintf('File created with %s from the MATLAB fvcom-toolbox', subname))
% define dimensions
nobc_dimid=netcdf.defDim(nc,'nobc',nObc);
datestrlen_dimid=netcdf.defDim(nc,'DateStrLen',26);
time_dimid=netcdf.defDim(nc,'time',netcdf.getConstant('NC_UNLIMITED'));
siglay_dimid=netcdf.defDim(nc,'siglay',nSiglay);
siglev_dimid=netcdf.defDim(nc,'siglev',nSiglev);
% variables
if strtime
Times_varid=netcdf.defVar(nc,'Times','NC_CHAR',[datestrlen_dimid, time_dimid]);
netcdf.putAtt(nc,Times_varid,'time_zone','UTC');
end
if floattime
time_varid=netcdf.defVar(nc,'time','NC_FLOAT',time_dimid);
netcdf.putAtt(nc,time_varid,'long_name','time');
netcdf.putAtt(nc,time_varid,'units','days since 1858-11-17 00:00:00');
netcdf.putAtt(nc,time_varid,'time_zone','UTC');
end
if inttime
itime_varid=netcdf.defVar(nc,'Itime','NC_INT',time_dimid);
netcdf.putAtt(nc,itime_varid,'units','days since 1858-11-17 00:00:00');
netcdf.putAtt(nc,itime_varid,'format','modified julian day (MJD)');
netcdf.putAtt(nc,itime_varid,'time_zone','UTC');
itime2_varid=netcdf.defVar(nc,'Itime2','NC_INT',time_dimid);
netcdf.putAtt(nc,itime2_varid,'units','msec since 00:00:00');
netcdf.putAtt(nc,itime2_varid,'time_zone','UTC');
end
nobc_varid=netcdf.defVar(nc,'obc_nodes','NC_INT',nobc_dimid);
netcdf.putAtt(nc,nobc_varid,'long_name','Open Boundary Node Number');
netcdf.putAtt(nc,nobc_varid,'grid','obc_grid');
netcdf.putAtt(nc,nobc_varid,'type','data');
obc_h_varid=netcdf.defVar(nc,'obc_h','NC_FLOAT',nobc_dimid);
netcdf.putAtt(nc,obc_h_varid,'long_name','Open Boundary Depth');
netcdf.putAtt(nc,obc_h_varid,'units','m');
netcdf.putAtt(nc,obc_h_varid,'grid','obc_grid');
netcdf.putAtt(nc,obc_h_varid,'type','data');
obc_siglev_varid=netcdf.defVar(nc,'siglev','NC_FLOAT',[nobc_dimid,siglev_dimid]);
netcdf.putAtt(nc,obc_siglev_varid,'long_name','ocean_sigma/general_coordinate');
netcdf.putAtt(nc,obc_siglev_varid,'grid','obc_grid');
obc_siglay_varid=netcdf.defVar(nc,'siglay','NC_FLOAT',[nobc_dimid,siglay_dimid]);
netcdf.putAtt(nc,obc_siglay_varid,'long_name','ocean_sigma/general_coordinate');
netcdf.putAtt(nc,obc_siglay_varid,'grid','obc_grid');
end
nobc_varid=netcdf.defVar(nc,'obc_nodes','NC_INT',nobc_dimid);
netcdf.putAtt(nc,nobc_varid,'long_name','Open Boundary Node Number');
netcdf.putAtt(nc,nobc_varid,'grid','obc_grid');
netcdf.putAtt(nc,nobc_varid,'type','data');
obc_h_varid=netcdf.defVar(nc,'obc_h','NC_FLOAT',nobc_dimid);
netcdf.putAtt(nc,obc_h_varid,'long_name','Open Boundary Depth');
netcdf.putAtt(nc,obc_h_varid,'units','m');
netcdf.putAtt(nc,obc_h_varid,'grid','obc_grid');
netcdf.putAtt(nc,obc_h_varid,'type','data');
obc_siglev_varid=netcdf.defVar(nc,'siglev','NC_FLOAT',[nobc_dimid,siglev_dimid]);
netcdf.putAtt(nc,obc_siglev_varid,'long_name','ocean_sigma/general_coordinate');
netcdf.putAtt(nc,obc_siglev_varid,'grid','obc_grid');
obc_siglay_varid=netcdf.defVar(nc,'siglay','NC_FLOAT',[nobc_dimid,siglay_dimid]);
netcdf.putAtt(nc,obc_siglay_varid,'long_name','ocean_sigma/general_coordinate');
netcdf.putAtt(nc,obc_siglay_varid,'grid','obc_grid');
% nutrients here
for nuts=1:NNuts
......@@ -206,6 +227,9 @@ for nuts=1:NNuts
eval(['netcdf.putAtt(nc,',varidN{nuts},',''long_name'',''',ERSEMdata(nuts).long_name,''');'])
eval(['netcdf.putAtt(nc,',varidN{nuts},',''units'',''',ERSEMdata(nuts).units,''');'])
eval(['netcdf.putAtt(nc,',varidN{nuts},',''grid'',''obc_grid'');'])
% enable compression on the big variables.
eval(['netcdf.defVarDeflate(nc, ',varidN{nuts},', true, true, 7);'])
% obc_salinity_varid=netcdf.defVar(nc,'obc_salinity','NC_FLOAT',[nobc_dimid,siglay_dimid,time_dimid]);
% netcdf.putAtt(nc,obc_salinity_varid,'long_name','sea_water_salinity');
% netcdf.putAtt(nc,obc_salinity_varid,'units','PSU');
......@@ -219,34 +243,63 @@ end
% end definitions
netcdf.endDef(nc);
% write data
netcdf.putVar(nc,nobc_varid,obc_nodes);
netcdf.putVar(nc,obc_h_varid,obc_h);
netcdf.putVar(nc,obc_siglev_varid,siglev);
netcdf.putVar(nc,obc_siglay_varid,siglay);
if strtime
nStringOut = char();
[nYr, nMon, nDay, nHour, nMin, nSec] = mjulian2greg(time);
for i=1:nTimes
nDate = [nYr(i), nMon(i), nDay(i), nHour(i), nMin(i), nSec(i)];
nStringOut = [nStringOut, sprintf('%04i/%02i/%02i %02i:%02i:%09.6f', nDate)];
if ~FileExists
% grid and time information already exist in file...
% write data
netcdf.putVar(nc,nobc_varid,obc_nodes);
netcdf.putVar(nc,obc_h_varid,obc_h);
netcdf.putVar(nc,obc_siglev_varid,siglev);
netcdf.putVar(nc,obc_siglay_varid,siglay);
if strtime
nStringOut = char();
[nYr, nMon, nDay, nHour, nMin, nSec] = mjulian2greg(time);
for i=1:nTimes
nDate = [nYr(i), nMon(i), nDay(i), nHour(i), nMin(i), nSec(i)];
nStringOut = [nStringOut, sprintf('%04i/%02i/%02i %02i:%02i:%09.6f', nDate)];
end
netcdf.putVar(nc,Times_varid,[0, 0],[26, nTimes],nStringOut);
end
if floattime
netcdf.putVar(nc,time_varid,0,numel(time),time);
end
if inttime
netcdf.putVar(nc,itime_varid,floor(time));
% netcdf.putVar(nc,itime2_varid,0,numel(time),mod(time,1)*24*3600*1000); % PWC original
% KJA edit: avoids rounding errors when converting from double to single
% Rounds to nearest multiple of the number of msecs in an hour
netcdf.putVar(nc,itime2_varid,0,numel(time),round((mod(time,1)*24*3600*1000)/(3600*1000))*(3600*1000));
end
for nuts=1:NNuts
eval(['netcdf.putVar(nc,',varidN{nuts},',Mobj.(ERSEMdata(nuts).name));'])
disp(['Finished with variable, ',ERSEMdata(nuts).name])
end
else % file exist and time could be different... check and interpolate if necessary
if length(time) < length(file_timem)
for nuts=1:NNuts
data= Mobj.(ERSEMdata(nuts).name);
dataint = nan(size(data, 1),size(data, 2),length(file_timem));
for nn=1:size(data, 1)
[X1,Y1]=meshgrid(file_timem- 678942,siglay(nn,:));
[X,Y]=meshgrid(time,siglay(nn,:));
% interpolate ERSEMdata...
dataint(nn,:,:) = interp2(X,Y,squeeze(data(nn,:,:)),X1,Y1);
end
Nut_id= netcdf.inqVarID(nc,ERSEMdata(nuts).name);
netcdf.putVar(nc,Nut_id,dataint);
disp(['Finished with variable, ',ERSEMdata(nuts).name])
end
else
% everything is in the same time frequency
for nuts=1:NNuts
eval(['netcdf.putVar(nc,',varidN{nuts},',Mobj.(ERSEMdata(nuts).name));'])
disp(['We have reversed the Z dimension of variable , ',ERSEMdata(nuts).name])
disp(['Finished with variable, ',ERSEMdata(nuts).name])
end
end
netcdf.putVar(nc,Times_varid,[0, 0],[26, nTimes],nStringOut);
end
if floattime
netcdf.putVar(nc,time_varid,0,numel(time),time);
end
if inttime
netcdf.putVar(nc,itime_varid,floor(time));
% netcdf.putVar(nc,itime2_varid,0,numel(time),mod(time,1)*24*3600*1000); % PWC original
% KJA edit: avoids rounding errors when converting from double to single
% Rounds to nearest multiple of the number of msecs in an hour
netcdf.putVar(nc,itime2_varid,0,numel(time),round((mod(time,1)*24*3600*1000)/(3600*1000))*(3600*1000));
end
for nuts=1:NNuts
eval(['netcdf.putVar(nc,',varidN{nuts},',Mobj.(ERSEMdata(nuts).name));'])
end
% close file
netcdf.close(nc);
......
function [Plots]=do_vector_plot_MatlabMapC(plotOPTS,FVCOM,tt)
%
% Function to display vector maps of FVCOM currents (i.e. U,V)
%
% [Plots]=do_vector_plot(plotOPTS,FVCOM)
%
% DESCRIPTION:
% Generates vector maps of currents using m_map toolbox
%
% INPUT:
% plotOPTS = structure array with predefined options for generating the
% maps
% FVCOM = data extracted from FVCOM NC file. See read_netCDF_FVCOM for
% details
%
% plotOPTS.range_lat: [50.1000 50.4000]
% plotOPTS.range_lon: [-4.5000 -3.8500]
% plotOPTS.fig_name: 'co2_S5_slowleak'
% plotOPTS.mesh: [1x1 struct]
% plotOPTS.coastline_file: '../mat/tamar3_0coast.mat'
% plotOPTS.zone: 30
% plotOPTS.ell: 'grs80'
% plotOPTS.var_plot: 'PH'
% plotOPTS.clims: [6 8]
% plotOPTS.do_mesh: 0
% plotOPTS.nz_plot: 1
% plotOPTS.figure: 1
% plotOPTS.Time_record: 7.3271e+05
% plotOPTS.nz_plot_vec: 1
% plotOPTS.data_dec: 5
% plotOPTS.vel_sca: 5
% plotOPTS.pause: 0.5000
%
% OUTPUT:
% Plots = structure array with figure handles
%
% EXAMPLE USAGE
% [Plots]=do_vector_plot(plotOPTS,FVCOM)
%
% Author(s):
% Ricardo Torres and Pierre Cazenave (Plymouth Marine Laboratory)
%
% Revision history
%
%==============================================================================
% adds m_map to matlab paths. file is in utilities directory.
% ammend according to your m_map installation paths
figure(plotOPTS.figure);
% generate figure with correct projection lat and lon range ellipsoid and
% zone.
axesm('mercator','MapLatLimit',plotOPTS.range_lat,'MapLonLimit',[plotOPTS.range_lon],'MeridianLabel','on',...
'ParallelLabel','on','MLineLocation',[0.1],'PLineLocation',[0.1],'LabelUnits','dm')
% add coastline if present
if (isfield(plotOPTS,'coastline_file') && ~isempty(plotOPTS.coastline_file) )
coast=load(plotOPTS.coastline_file);
geoshow([coast.ncst(:,2)],[coast.ncst(:,1)],'Color','black')
end
%-----------------------------------------------------------------------
% Convert element positions from FVCOM to lat and lon using m_ll2ll.m from
% utilities directory. This accesses m_map functions.
% In my case it needs adding 6 because of discrepancies between proj and m_map.
% Proj automatically determines a
% reference long in strides of 6deg while m_map doesn't
%------------------------------------------------------------------------
% only plot vectors inside lat and lon range and ...
if ~isfield(plotOPTS.mesh,'latc')
try
plotOPTS.mesh.latc = nodes2elems(plotOPTS.mesh.lat,plotOPTS.mesh);
plotOPTS.mesh.lonc = nodes2elems(plotOPTS.mesh.lon,plotOPTS.mesh);
catch
disp(['We have no access to nodes2elems in the fvcom_prepro directory'])
end
end
igood = find ( plotOPTS.mesh.lonc < plotOPTS.range_lon(2) & plotOPTS.mesh.lonc > plotOPTS.range_lon(1) &...
plotOPTS.mesh.latc < plotOPTS.range_lat(2) & plotOPTS.mesh.latc > plotOPTS.range_lat(1));
% decimate positions. Plot every plotOPTS.data_dec position.
igood=igood(1:plotOPTS.data_dec:end);
%------------------------------------------------------------------------
% Select how many layers to plot
%------------------------------------------------------------------------
ND=ndims(FVCOM.(plotOPTS.var_plotu))
switch ND
case 2
nLayers=1;
colourSpec=[0 0 0];
case 3
if isfield(plotOPTS,'nz_plot_vec')
nLayers=size(plotOPTS.nz_plot_vec,2);
nLayersRange=plotOPTS.nz_plot_vec;
else
nLayers=size(plotOPTS.nz_plot,2);
nLayersRange=plotOPTS.nz_plot;
end
% choose colors for vectors
if nLayers==1
colourSpec=[0 0 0];
else
colourSpec=colormap(hsv(nLayers));
end
end
% Preallocate outputs
u=nan(size(igood,1),nLayers,1);
v=nan(size(igood,1),nLayers,1);
% Check if we're running
% for aa=1:length(plotOPTS.Time_record)
aa=tt;
fprintf('Time step %i of %i\n',aa,length(FVCOM.Time_record));
% Mesh goes underneath vectors.
% if plotOPTS.do_mesh
% % plot vertices
% Plots(plotOPTS.figure).mesh=patch('Vertices',[double(plotOPTS.mesh.lat),double(plotOPTS.mesh.lon)],'Faces',plotOPTS.mesh.tri,...
% 'EdgeColor',[0.6 0.6 0.6],'FaceColor','none');hold on
% end
switch ND
case 2
u(:,1)=squeeze(FVCOM.(plotOPTS.var_plotu)(igood,(aa)));
v(:,1)=squeeze(FVCOM.(plotOPTS.var_plotv)(igood,(aa)));
case 3
for ii=1:nLayers
u(:,ii)=squeeze(FVCOM.(plotOPTS.var_plotu)(igood,nLayersRange(ii),(aa)));
v(:,ii)=squeeze(FVCOM.(plotOPTS.var_plotv)(igood,nLayersRange(ii),(aa)));
end
end
for jj=1:nLayers
% [Plots(plotOPTS.figure).handles{jj}]=quiver(double(plotOPTS.mesh.lonc(igood)),double(plotOPTS.mesh.latc(igood)),u,v,plotOPTS.vel_sca,'k','filled');
[Plots(plotOPTS.figure).handles{jj}]=quiverwcolorbar(double(plotOPTS.mesh.lonc(igood)),double(plotOPTS.mesh.latc(igood)),u,v,...
plotOPTS.vel_sca,'bounds',plotOPTS.clims)
% [Plots(plotOPTS.figure).handles{jj}]=quiverm(double(FVCOM.latc(igood)),double(FVCOM.lonc(igood)),u./100,v./100,'k',0,'filled'),
% [Plots(plotOPTS.figure).handles(jj),~]=m_vec(plotOPTS.vel_sca,x(igood),y(igood),...
% u(:,jj),v(:,jj),colourSpec(jj,:),...
% 'shaftwidth',1,'headwidth',2);
% hold on
end
pause(plotOPTS.pause)
% if aa~=length(plotOPTS.Time_record)
% delete(Plots(plotOPTS.figure).handles(:))
% end
% end
function hh = quiverwcolorbar(varargin)
% Melissa Day 5/24/2013
% Upgrade of Andrew Diamond's quiverc2wcmap to generate a quiver plot
% with arrows colored according to vector magnitude.
% Functional differences from quiverc2wcmap:
% 1) Allows user to specify colormap boundaries using 'bound': changes
% colorbar axes AND corresponding vector coloring
% (much more useful for intercomparison of datasets)
% 2) Improved fidelity to magnitude colors in large datasets
% (at a cost of increased computation time...)
% Uses some improvements from DS's vfield_color
% 3) Small clarifications added
% Example:
% x = rand(1,50).*100;
% y = rand(1,50).*100;
% u = rand(1,50) .* 10;
% v = rand(1,50) .* 10;
% scale = 0;
% figure; quiverwcolorbar(x',y',u',v',scale); %compare to:
% figure; quiverwcolorbar(x',y',u',v',scale,'bounds',[0 10]);
%----------------
% function hh = quiverc2wcmap(varargin)
% Andrew Diamond 3/17/2005
% This is based off of Tobias Höffken which was based off of Bertrand Dano
% keeping with their main intent to show a color w/r to magnitude quiver plot
% while maintaining a large amount of quiver API backwards compatability.
% Functional differences from quiverc2:
% 1) This works under 6.5.1
% 2) It handles NaNs
% 3) It draws a colormap that is w/r to the quiver magnitudes (hard coded to
% 20 segments of the colormap as per quiverc2 - seems fine for a quiver).
% 4) Various bug fixes (I think)
% In order to do this I needed some small hacks on 6.5.1's quiver to take a
% color triplet. I have included as part of this file in a subfunction below.
%----------------
% Comments from quiverc2
% changed Tobias Höffken 3-14-05
% totally downstripped version of the former
% split input field into n segments and do a quiver qhich each of them
% Modified version of Quiver to plots velocity vectors as arrows
% with components (u,v) at the points (x,y) using the current colormap
% Bertrand Dano 3-3-03
% Copyright 1984-2002 The MathWorks, Inc.
% changed T. Höffken 14.03.05, for high data resolution
% using fixed color "spacing" of 20
%QUIVERC Quiver color plot.
% QUIVERC(X,Y,U,V) plots velocity vectors as arrows with components (u,v)
% at the points (x,y). The matrices X,Y,U,V must all be the same size
% and contain corresponding position and velocity components (X and Y
% can also be vectors to specify a uniform grid). QUIVER automatically
% scales the arrows to fit within the grid.
%
% QUIVERC(U,V) plots velocity vectors at equally spaced points in
% the x-y plane.
%
% QUIVERC(U,V,S) or QUIVER(X,Y,U,V,S) automatically scales the
% arrows to fit within the grid and then stretches them by S. Use
% S=0 to plot the arrows without the automatic scaling.
%
% QUIVERC(...,LINESPEC) uses the plot linestyle specified for
% the velocity vectors. Any marker in LINESPEC is drawn at the base
% instead of an arrow on the tip. Use a marker of '.' to specify
% no marker at all. See PLOT for other possibilities.