Commit 5eb0e887 authored by Pierre Cazenave's avatar Pierre Cazenave

Merge branch 'dev' of github.com:pwcazenave/fvcom-toolbox into geoff.

parents 03df2384 9894493d
This diff is collapsed.
......@@ -24,10 +24,13 @@ function [Mobj] = add_obc_nodes_list(Mobj,Nlist,ObcName,ObcType,plotFig)
% Author(s):
% Geoff Cowles (University of Massachusetts Dartmouth)
% Pierre Cazenave (Plymouth Marine Laboratory)
%
% Karen Amoudry (National Oceanography Centre, Liverpool)
%
% Revision history:
% 2012-11-26 Add ability to turn off the figures.
% 2013-01-02 KJA bug fix: amended usage of 'unique' in line 53 to
% prevent it from sorting the values it returns. Amended by Pierre to
% support pre-2012 versions of MATLAB whilst giving the same result.
%
%==========================================================================
subname = 'add_obc_nodes';
......@@ -45,7 +48,10 @@ end
%--------------------------------------------------------------------------
% Get a unique list and make sure they are in the range of node numbers
%--------------------------------------------------------------------------
Nlist = unique(Nlist);
% Make this works in versions of MATLAB older than 2012a (newer versions
% can just use unique(A, 'stable'), but checking versions is a pain).
[~, Nidx] = unique(Nlist);
Nlist = Nlist(sort(Nidx));
if(max(Nlist) > Mobj.nVerts);
fprintf('your open boundary node number exceed the total number of nodes in the domain\n');
......
......@@ -20,12 +20,16 @@ function [Mobj] = add_river_nodes_list(Mobj,Nlist,RiverName)
%
% Author(s):
% Geoff Cowles (University of Massachusetts Dartmouth)
% Karen Amoudry (National Oceanography Centre, Liverpool)
%
% Note:
% Uses ginput2 which allows zooming before selecting points and displays
% clicked points realtime
%
% Revision history
% 2013-01-02 KJA bug fix: amended usage of 'unique' to prevent it from
% sorting the values it returns. Amended by Pierre to support pre-2012
% versions of MATLAB whilst giving the same result.
%
%==========================================================================
subname = 'add_river_nodes_list';
......@@ -38,7 +42,10 @@ end
%--------------------------------------------------------------------------
% Get a unique list and make sure they are in the range of node numbers
%--------------------------------------------------------------------------
Nlist = unique(Nlist);
% Make this works in versions of MATLAB older than 2012a (newer versions
% can just use unique(A, 'stable'), but checking versions is a pain).
[~, Nidx] = unique(Nlist);
Nlist = Nlist(sort(Nidx));
if max(Nlist) > Mobj.nVerts
fprintf('your river node number(s) exceed the total number of nodes in the domain\n');
......
......@@ -22,7 +22,13 @@ function [spongeRadius] = calc_sponge_radius(Mobj,Nlist)
% spongeRadius = calc_sponge_radius(Mobj,Nlist)
%
% Author(s)
% Karen Thurston (National Oceanography Centre, Liverpool)
% Karen Amoudry (National Oceanography Centre, Liverpool)
% Pierre Cazenave (Plymouth Marine Laboratory)
%
% Revision history:
% 2013-01-02 KJA bug fix: amended usage of 'unique' to prevent it from
% sorting the values it returns. Amended by Pierre to support pre-2012
% versions of MATLAB whilst giving the same result.
%
%==========================================================================
subname = 'calc_sponge_radius';
......@@ -35,7 +41,10 @@ end
%--------------------------------------------------------------------------
% Get a unique list and make sure they are in the range of node numbers
%--------------------------------------------------------------------------
Nlist = unique(Nlist);
% Make this works in versions of MATLAB older than 2012a (newer versions
% can just use unique(A, 'stable'), but checking versions is a pain).
[~, Nidx] = unique(Nlist);
Nlist = Nlist(sort(Nidx));
spongeRadius = 100000+zeros(size(Nlist));
......@@ -43,7 +52,9 @@ spongeRadius = 100000+zeros(size(Nlist));
for i =1:length(Nlist)
% Find the neighbouring nodes
[r,c]=find(Mobj.tri==Nlist(i));
neighbours = unique(Mobj.tri(r,:));
neighbours = Mobj.tri(r,:);
[~,neighidx] = unique(Mobj.tri(r,:));
neighbours = neighbours(sort(neighidx));
% Remove the node of interest from the neighbours list
n = find(neighbours~=Nlist(i));
......
function [metvar,X_send,Y_send] = extract_mesoscale(fname,ndays)
%load in the mesoscale data from mesoscale operational POLCOMS met input data and convert to cs3 operational surge model met data grid
% Script from JMB
A= load (fullfile('/work/kthurs/Met_data_processing_scripts/MESOSCALE/OUTPUT/',fname));
tint=24/3;%Time interval 3hrs
days=ndays+2-1/tint; % number of days in month +2-1 to include extra day of output up to 21:00, the output before midnight of the next day
% day 1 = 00:00 of the start day, day 1.5 = 12:00 of the start day
%POLCOMS meso met grid
dx=0.11; dy=0.11;
x1=-13;
y1=48.39;
nx=218;
ny=136;
x2=(nx-1)*dx+x1;
y2=(ny-1)*dy+y1;
[X,Y]=meshgrid(x1:dx:x2,y1:dy:y2);%grid
X_send = x1:dx:x2;
Y_send = y1:dy:y2;
C=zeros(ny*days*tint,nx);%ygrid*number of days * outputs per day, xgrid
jj=1;
for ii=1:(ny*nx):length(A) %total length every 3 hours output
% B = reshape(A(ii,ii+29647),136,218);
C(jj:jj+ny-1,:) = reshape(A(ii:ii+nx*ny-1),nx,ny)';
jj=jj+ny;
end
% Reshape C to have form ny x nx x days*tint
metvar = zeros(ny,nx,days*tint);
jj=1;
for kk = 1:ny:length(C)
metvar(:,:,jj) = C(kk:kk+ny-1,:);
jj = jj+1;
end
%%
% % Load coastline
% load '~/POLCOMS/POLCOMS_matlab/N_Atlantic_coast_m.dat';
% x=N_Atlantic_coast_m(:,1);
% y=N_Atlantic_coast_m(:,2);
%
% %
% for kk=1:ny:length(C)
% figure(3)
% pcolor(X,Y,C(kk:kk+ny-1,:))
% shading interp
% hold on
% colorbar
% % %Plot coastline
% % plot(x,y,'k')
% % drawnow
% end
function [Mobj] = get_AMM(Mobj,StartDate,EndDate,ModelFolder)
% Extract boundary forcing information from NOC Operational Tide Surge
% Model output output
% Model output.
%
% function get_AMM(Mobj,StartDate,EndDate,ModelFolder)
%
......
......@@ -44,6 +44,7 @@ function Mobj = get_FVCOM_rivers(Mobj, dist_thresh)
%
% Author(s):
% Pierre Cazenave (Plymouth Marine Laboratory)
% Karen Amoudry (National Oceanography Centre, Liverpool)
%
% Revision history:
% 2013-03-27 - First version.
......
......@@ -49,14 +49,14 @@ if isunix % Unix?
metfname = ['/bank/jane/met/',datestr(inputConf.startDate,'YYYY'),...
'/',lower(datestr(inputConf.startDate,'mmmYY')),'nae10R.dat'];
comprfname = '/login/jane/NAE2/metintco.cs3x.nae2.compress.2';
setupfname = '/work/jane/cs3x/prep/setupcs3xSGIL.uda';
setupfname = '/work/jane/cs3x/prep/setupcs3xSGIl.uda';
% elevfname = ['/bank/jane/cs3x/sarray.uda.',...
% datestr(inputConf.startDate,'YYYY')];
elseif ispc % Or Windows?
metfname = ['\\store\bank\jane\met\',datestr(inputConf.startDate,'YYYY'),...
'\',lower(datestr(inputConf.startDate,'mmmYY')),'nae10R.dat'];
comprfname = '\\store\kthurs\from_Jane\metintco.cs3x.nae2.compress.2';
setupfname = '\\store\work\jane\cs3x\prep\setupcs3xSGIL.uda';
setupfname = '\\store\work\jane\cs3x\prep\setupcs3xSGIl.uda';
% elevfname = ['\\store\bank\jane\cs3x\sarray.uda.',...
% datestr(inputConf.startDate,'YYYY')];
end
......
......@@ -66,11 +66,21 @@ function data = get_NCEP_forcing(Mobj, modelTime)
% the toolbox's need. Also, we're not actually using 'pevpr' for the
% calculation of evaporation since we're estimating that from the latent
% heat net flux ('lhtfl'), so it's superfluous anyway.
% 2013-06-28 Changed the way the Matlab version is determiend. Now using
% release date rather then version number. For example version 7.13 >
% verion 7.7 but 7.13 is not greater than 7.7.
%
%==========================================================================
subname = 'get_NCEP_forcing';
% Define date that matlab version 7.14 was released.
% OPeNDAP was included in version 7.14
% see http://en.wikipedia.org/wiki/MATLAB and
% https://publicwiki.deltares.nl/display/OET/OPeNDAP+access+with+Matlab
version_7_14_date = datenum(2012,3,1);
%version_7_13_date = datenum(2011,9,1);
global ftbverbose;
if ftbverbose
fprintf('\nbegin : %s\n', subname)
......@@ -111,6 +121,10 @@ ncep.lhtfl = ['http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanaly
ncep.shtfl = ['http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis/surface_gauss/shtfl.sfc.gauss.',num2str(year),'.nc'];
% ncep.pevpr = ['http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis/surface_gauss/pevpr.sfc.gauss.',num2str(year),'.nc'];
% Possible future data to use?
% Skin temperature
% ncep.skt = ['http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis/surface_gauss/skt.sfc.gauss.',num2str(year),'.nc'];
% The fields below can be used to create the net shortwave and longwave
% fluxes if the data you're using don't include net fluxes. Subtract the
% downward from upward fluxes to get net fluxes.
......@@ -132,7 +146,7 @@ for aa = 1:length(fields)
% libraries to load the OPeNDAP data, otherwise we need the relevant
% third-party toolbox.
out = ver('MATLAB');
if str2double(out.Version) > 7.13
if datenum(out.Date) > version_7_14_date % Look at the date rather than the version number
%ncid_info = ncinfo(ncep.(fields{aa}));
ncid = netcdf.open(ncep.(fields{aa}));
......@@ -208,7 +222,7 @@ for aa = 1:length(fields)
if iscell(index_lon)
data.(fields{aa}).lon = data_lon.lon(cat(1,index_lon{:}));
if str2double(out.Version) > 7.13
if datenum(out.Date) > version_7_14_date % Look at the date rather than the version number
% varidlon = netcdf.inqVarID(ncid,'lon');
% varidtime = netcdf.inqVarID(ncid,'time');
% varidlat = netcdf.inqVarID(ncid,'lat');
......@@ -282,7 +296,7 @@ for aa = 1:length(fields)
% We have a straightforward data extraction
data.(fields{aa}).lon = data_lon.lon(index_lon);
if str2double(out.Version) > 7.13
if datenum(out.Date) > version_7_14_date % Look at the date rather than the version number
varid = netcdf.inqVarID(ncid,(fields{aa}));
% [varname,xtype,dimids,natts] = netcdf.inqVar(ncid,varid);
% [~,length1] = netcdf.inqDim(ncid,dimids(1))
......@@ -327,6 +341,14 @@ end
% Now we have some data, we need to create some additional parameters
% required by FVCOM.
% FVCOM's sign convention is the opposite of the NCEP data for heat fluxes
% (FVCOM: positive = downward flux = ocean heating, negative = upward flux
% = ocean cooling. NCEP: positive = upward flux = ocean cooling, negative =
% downward flux = ocean heating). So, rather than do the corrections in
% create_files.m or wherever, do them here instead.
% data.nlwrs.data = -data.nlwrs.data;
% data.nswrs.data = -data.nswrs.data;
% Convert precipitation from kg/m^2/s to m/s (required by FVCOM) by
% dividing by freshwater density (kg/m^3).
data.prate.data = data.prate.data/1000;
......
......@@ -90,7 +90,7 @@ for ii = 1:todo
ncdata.(getVar).data = data;
else
if ndims(data) < 3
if strcmpi(varlist{var}, 'time')
if strcmpi(getVar, 'time')
% If the dimension is time, we need to be a bit more
% clever since we'll need a concatenated time series
% (in which values are continuous and from which we
......
......@@ -120,12 +120,6 @@ for t = 1:nt
% Get the current 3D array of PML POLCOMS-ERSEM results.
pctemp3 = pc.ETWD.data(:, :, :, t);
pcsalt3 = pc.x1XD.data(:, :, :, t);
% Flip the vertical layer dimension to make the POLCOMS data go from
% surface to seabed to match its depth data and to match how FVCOM
% works.
pctemp3 = flipdim(pctemp3, 3);
pcsalt3 = flipdim(pcsalt3, 3);
% Preallocate the intermediate results arrays.
itempz = nan(nf, nz);
......@@ -298,12 +292,12 @@ end
%%
% % Plot a vertical profile for a boundary node (for my Irish Sea case, this
% % is one of the ones along the Celtic Sea boundary). Also plot the
% % distribution of interpolated values over the POLCOMS data. Add the
% % location of the vertical profile (both FVCOM and POLCOMS) to the plot.
% Plot a vertical profile for a boundary node (for my Irish Sea case, this
% is one of the ones along the Celtic Sea boundary). Also plot the
% distribution of interpolated values over the POLCOMS data. Add the
% location of the vertical profile (both FVCOM and POLCOMS) to the plot.
% nn = 55; % open boundary index
% tt = 1; % time index
% tt = 1; % time index
%
% % Get the corresponding indices for the POLCOMS data
% [~, xidx] = min(abs(lon(1, :) - fvlon(nn)));
......@@ -322,7 +316,7 @@ end
% subplot(2,2,2)
% % Although POLCOMS stores its temperature values from seabed to surface,
% % the depths are stored surface to seabed. Nice.
% plot(squeeze(pc.ETWD.data(xidx, yidx, :, 1)), flipud(squeeze(pc.depth.data(xidx, yidx, :, 1))), 'rx-')
% plot(squeeze(pc.ETWD.data(xidx, yidx, :, 1)), squeeze(pc.depth.data(xidx, yidx, :, 1)), 'rx-')
% xlabel('Temperature (^{\circ}C)')
% ylabel('Depth (m)')
% title('POLCOMS')
......@@ -349,15 +343,15 @@ end
% dx = mean(diff(pc.lon.data));
% dy = mean(diff(pc.lat.data));
% pcolor(pc.lon.data - (dx / 2), pc.lat.data - (dy / 2), ...
% squeeze(pc.ETWD.data(:, :, end, tt))')
% squeeze(pc.ETWD.data(:, :, 1, tt))')
% shading flat
% axis('equal', 'tight')
% daspect([1.5, 1, 1])
% hold on
% % Add the interpolated surface data (first sigma layer)
% scatter(Mobj.lon(oNodes), Mobj.lat(oNodes), repmat(40, size(Mobj.lon(oNodes))), Mobj.temperature(:, 1, tt), 'filled', 'MarkerEdgeColor', 'k')
% scatter(Mobj.lon(oNodes), Mobj.lat(oNodes), 40, Mobj.temperature(:, 1, tt), 'filled', 'MarkerEdgeColor', 'k')
% axis([min(Mobj.lon(oNodes)), max(Mobj.lon(oNodes)), min(Mobj.lat(oNodes)), max(Mobj.lat(oNodes))])
% caxis([6, 12])
% caxis([6, 20])
% plot(lon(yidx, xidx), lat(yidx, xidx), 'rs') % polcoms is all backwards
% plot(Mobj.lon(oNodes(nn)), Mobj.lat(oNodes(nn)), 'wo')
% colorbar
This diff is collapsed.
......@@ -117,12 +117,10 @@ if ftbverbose
fprintf('%s : interpolate POLCOMS onto FVCOM''s vertical grid... ', subname)
end
% Permute the arrays to be x by y rather than y by x. Also flip the
% vertical layer dimension to make the POLCOMS data go from surface to
% seabed to match its depth data and to match how FVCOM works.
temperature = flipdim(permute(squeeze(pc.ETWD.data(:, :, :, tidx)), [2, 1, 3]), 3);
salinity = flipdim(permute(squeeze(pc.x1XD.data(:, :, :, tidx)), [2, 1, 3]), 3);
density = flipdim(permute(squeeze(pc.rholocalD.data(:, :, :, tidx)), [2, 1, 3]), 3);
% Permute the arrays to be x by y rather than y by x.
temperature = permute(squeeze(pc.ETWD.data(:, :, :, tidx)), [2, 1, 3]);
salinity = permute(squeeze(pc.x1XD.data(:, :, :, tidx)), [2, 1, 3]);
density = permute(squeeze(pc.rholocalD.data(:, :, :, tidx)), [2, 1, 3]);
depth = permute(squeeze(pc.depth.data(:, :, :, tidx)), [2, 1, 3]);
mask = depth(:, :, end) >= 0; % land is positive.
......@@ -221,33 +219,32 @@ if ftbverbose
end
%% Debugging figure
%
% close all
%
% tidx = 1; % time step to plot
% ri = 85; % column index
% ci = 95; % row index
%
% [~, idx] = min(sqrt((Mobj.lon - lon(ri, ci)).^2 + (Mobj.lat - lat(ri, ci)).^2));
%
% % Vertical profiles
% figure
% clf
%
% % The top row shows the temperature/salinity values as plotted against
% % index (i.e. position in the array). Since POLCOMS stores the seabed
% % first, its values appear at the bottom i.e. the profile is the right way
% % up. Just to make things interesting, the depths returned from the NetCDF
% % files are stored the opposite way (surface is the first value in the
% % array). So, if you plot temperature/salinity against depth, the profile
% % is upside down.
% %
% % Thus, the vertical distribution of temperature/salinity profiles should
% % match in the top and bottom rows. The temperature/salinity data are
% % flipped in those figures (either directly in the plot command, or via the
% % flipped arrays (temperature, salinity)).
% % index (i.e. position in the array). I had thought POLCOMS stored its data
% % from seabed to sea surface, but looking at the NetCDF files in ncview,
% % it's clear that the data are in fact stored surface to seabed (like
% % FVCOM). As such, the two plots in the top row should be upside down (i.e.
% % surface values at the bottom of the plot). The two lower rows should have
% % three lines which all match: the raw POLCOMS data, the POLCOMS data for
% % the current time step (i.e. that in 'temperature' and 'salinity') and the
% % interpolated FVCOM data against depth.
% %
% % Furthermore, the pc.*.data have the rows and columns flipped, so (ci, ri)
% % in pc.*.data and (ri, ci) in 'temperature', 'salinity' and 'depth'.
% % Needless to say, the two lines in the lower plots should overlap.
% % Also worth noting, the pc.*.data have the rows and columns flipped, so
% % (ci, ri) in pc.*.data and (ri, ci) in 'temperature', 'salinity' and
% % 'depth'. Needless to say, the two lines in the lower plots should
% % overlap.
%
% subplot(2,2,1)
% plot(squeeze(pc.ETWD.data(ci, ri, :, tidx)), 1:size(depth, 3), 'rx:')
......@@ -264,30 +261,36 @@ end
% subplot(2,2,3)
% % Although POLCOMS stores its temperature values from seabed to surface,
% % the depths are stored surface to seabed. Nice. Flip the
% % temperature/salinity data accordingly.
% plot(flipud(squeeze(pc.ETWD.data(ci, ri, :, tidx))), squeeze(pc.depth.data(ci, ri, :, tidx)), 'rx-')
% % temperature/salinity data accordingly. The lines here should match one
% % another.
% plot(squeeze(pc.ETWD.data(ci, ri, :, tidx)), squeeze(pc.depth.data(ci, ri, :, tidx)), 'rx-')
% hold on
% plot(squeeze(temperature(ri, ci, :)), squeeze(depth(ri, ci, :)), '.:')
% % Add the equivalent interpolated FVCOM data point
% plot(fvtemp(idx, :), Mobj.siglayz(idx, :), 'g.-')
% xlabel('Temperature (^{\circ}C)')
% ylabel('Depth (m)')
% title('Depth Temperature')
% legend('pc', 'temp', 'location', 'north')
% legend('pc', 'temp', 'fvcom', 'location', 'north')
% legend('boxoff')
%
% subplot(2,2,4)
% % Although POLCOMS stores its temperature values from seabed to surface,
% % the depths are stored surface to seabed. Nice. Flip the
% % temperature/salinity data accordingly.
% plot(flipud(squeeze(pc.x1XD.data(ci, ri, :, tidx))), squeeze(pc.depth.data(ci, ri, :, tidx)), 'rx-')
% % temperature/salinity data accordingly. The lines here should match one
% % another.
% plot(squeeze(pc.x1XD.data(ci, ri, :, tidx)), squeeze(pc.depth.data(ci, ri, :, tidx)), 'rx-')
% hold on
% plot(squeeze(salinity(ri, ci, :)), squeeze(depth(ri, ci, :)), '.:')
% % Add the equivalent interpolated FVCOM data point
% plot(fvsalt(idx, :), Mobj.siglayz(idx, :), 'g.-')
% xlabel('Salinity')
% ylabel('Depth (m)')
% title('Depth Salinity')
% legend('pc', 'salt', 'location', 'north')
% legend('pc', 'salt', 'fvcom', 'location', 'north')
% legend('boxoff')
%
% % Plot the sample location
% %% Plot the sample location
% figure
% dx = mean(diff(pc.lon.data));
% dy = mean(diff(pc.lat.data));
......
This diff is collapsed.
function [Mobj] = read_grid_mesh(varargin)
% Read .grid mesh files into Matlab mesh object
%
% [Mobj] = function read_grid_mesh(varargin)
%
% DESCRIPTION:
% Read NOCL .grid file
% Store in a matlab mesh object
%
% INPUT [keyword pairs]:
% 'grid' = NOCL .grid file (e.g. UK_grid.grid)
% [optional] 'coordinate' = coordinate system for output data [spherical; cartesian (default)]
% [optional] 'input_coord' = coordinate system for input data [spherical; cartesian (default)]
% [optional] 'project' = generate (x,y) coordinates if input is (lon,lat)
% generate (lon,lat) coordinates if input is (x,y)
% ['true' ; 'false' (default)], see myproject.m
% [optional] 'zone' = specify UTM zone for projection
% [optional] 'addCoriolis' = calculate Coriolis param (f), requires [lon,lat]
%
% OUTPUT:
% Mobj = matlab structure containing mesh data
%
% EXAMPLE USAGE
% Mobj = read_grid_mesh('grid','UK_grid.grid','coordinate','spherical','addCoriolis','true')
%
% Author(s):
% Geoff Cowles (University of Massachusetts Dartmouth)
% Pierre Cazenave (Plymouth Marine Laboratory)
% Karen Thurston (National Oceanography Centre Liverpool)
%
% Revision history (KJT)
% 2012-11-13 Adapted 'read_sms_mesh.m' to read NOCL .grid files
% 2012-11-19 Added input/output coordinate functionality
%
%==============================================================================
subname = 'read_grid_mesh';
global ftbverbose;
if(ftbverbose);
fprintf('\n')
fprintf(['begin : ' subname '\n'])
end;
have_grid = false;
have_bath = false;
have_lonlat = false;
have_xy = false;
userproject = false;
haveUTM = false;
addCoriolis = false;
%------------------------------------------------------------------------------
% Create a blank mesh object
%------------------------------------------------------------------------------
Mobj = make_blank_mesh();
%------------------------------------------------------------------------------
% Parse input arguments
%------------------------------------------------------------------------------
if(mod(length(varargin),2) ~= 0)
error(['incorrect usage of ',subname,', use keyword pairs'])
end;
for i=1:2:length(varargin)-1
keyword = lower(varargin{i});
if( ~ischar(keyword) )
error(['incorrect usage of ',subname,', check keywords'])
end;
switch(keyword(1:3))
case 'gri'
grid_fname = varargin{i+1};
have_grid = true;
case 'coo'
coord = varargin{i+1};
if(coord(1:3)=='sph')
coordinate = 'spherical';
else
coordinate = 'cartesian';
end;
case 'in_'
in_coord = varargin{i+1};
if(in_coord(1:3)=='sph')
in_coordinate = 'spherical';
have_lonlat = true;
else
in_coordinate = 'cartesian';
have_xy = true;
end;
case 'pro'
val = varargin{i+1};
if( val )
userproject = true;
else
userproject = false;
end;
case 'zon'
fullzone = varargin{i+1};
UTMzone = regexpi(fullzone,'\ ','split');
UTMzone=str2double(char(UTMzone{1}(1)));
haveUTM = true;
case 'add'
val = varargin{i+1};
if( val )
addCoriolis = true;
else
addCoriolis = false;
end;
otherwise
error(['Can''t understand property:' char(varargin{i+1})]);
end; %switch(keyword)
end;
%------------------------------------------------------------------------------
% Read the mesh from the .grid file
%------------------------------------------------------------------------------
fid = fopen(grid_fname,'r');
if(fid < 0)
error(['file: ' grid_fname ' does not exist']);
end;
% Count number of elements and vertices
if(ftbverbose);
fprintf(['reading from: ' grid_fname '\n'])
end;
lin = fgetl(fid); % ignore first line, it's the mesh name
c=textscan(fid,'%u %u',1); % get nElems and nVerts
nElems = c{1};
nVerts = c{2};
clear c
if(ftbverbose);
fprintf('nVerts: %d\n',nVerts);
fprintf('nElems: %d\n',nElems);
fprintf('reading in connectivity and grid points\n')
end;
% allocate memory to hold mesh and connectivity
tri = zeros(nElems,3);
x = zeros(nVerts,1);
y = zeros(nVerts,1);
h = zeros(nVerts,1);
lon = zeros(nVerts,1);
lat = zeros(nVerts,1);
ts = zeros(nVerts,1);
c = textscan(fid, '%u %f %f %f ', nVerts);
x = c{2};
y = c{3};
h = c{4};
clear c
c = textscan(fid, '%u %u %u %u %u', nElems);
tri(:,1) = c{3};
tri(:,2) = c{4};
tri(:,3) = c{5};
clear c
% Make sure we have bathymetry
if sum(h)==0
have_bath=false;
else
have_bath=true;
end
% Make sure we have positive depths
if sum(h>0) < sum(h<0)
h = -h;
end
% Build array of all the nodes in the open boundaries
c = textscan(fid, '%u %*[^\n]',1);
nOpenSeg = c{1}; % number of open boundary segments
clear c
lin=fgetl(fid); % skip the next line
c = textscan(fid, '%u %*[^\n]',1);
nOpenNodes = c{1}; % number of open boundary nodes
clear c
% Initialise SegCount variable
SegCount = [0,0];
for i=1:nOpenSeg % for each open boundary segment
c = textscan(fid, '%u %*[^\n]',1); % how many nodes in this segment?
SegCount(1) = 1+SegCount(2);
SegCount(2) = SegCount(1)+c{1}-1;
clear c
c = textscan(fid,'%u %*[^\n]',(SegCount(2)-SegCount(1)+1)); % get all the nodes in this segment
allNodes{i} = c{1};
clear c
end
if have_lonlat == true
lon = x;
lat = y;
x = x*0.0;
y = y*0.0;
% Just do a double check on the coordinates to make sure we don't
% actually have cartesian
if max(lon)>360
warning('You''ve specified spherical coordinates, but your upper longitude value exceeds 360 degrees. Are you sure you have spherical data?')
end
else
have_xy = true;
end;
%------------------------------------------------------------------------------
% Project if desired by user
%------------------------------------------------------------------------------
if(userproject)
if (in_coordinate(1:3)=='car')
fprintf('inverse projecting to get (lon,lat)\n')
utmZones=cellfun(@(x) repmat(x,length(x),1),fullzone,'uni',false);
[lon,lat] = utm2deg(x,y,utmZones{1});
Mobj.have_lonlat = true;
elseif (in_coordinate(1:3)=='sph')
fprintf('forward projecting to get (x,y)\n')
[x,y] = wgs2utm(lat,lon,UTMzone,'N');
have_xy = true;
end
end
%------------------------------------------------------------------------------
% Transfer to Mesh structure
%------------------------------------------------------------------------------
Mobj.nVerts = nVerts;
Mobj.nElems = nElems;
Mobj.nativeCoords = coordinate;
if have_lonlat==true
Mobj.have_lonlat = have_lonlat;
end;
if have_xy==true
Mobj.have_xy = have_xy;
end;
Mobj.have_bath = have_bath;
Mobj.read_obc_nodes = allNodes;
Mobj.x = x;
Mobj.y = y;
Mobj.ts = ts;
Mobj.lon = lon;
Mobj.lat = lat;
Mobj.h = h;
Mobj.tri = tri;
if addCoriolis==true
Mobj.have_cor = true;
Mobj = add_coriolis(Mobj,'uselatitude');
end
if(ftbverbose);
fprintf(['end : ' subname '\n'])
end;
fclose(fid);
......@@ -25,24 +25,24 @@ function write_FVCOM_forcing(Mobj, fileprefix, data, infos, fver)
%
% The fields in data may be called any of:
% - 'u10', 'v10', 'uwnd', 'vwnd' - wind components
% - 'slp' - sea level pressure
% - 'pevpr' - evaporation
% - 'prate' - precipitation
% - 'slp' - sea level pressure
% - 'nlwrs' - net longwave radiation*,**
% - 'nswrs' - net shortwave radiation*,**
% - 'shtfl' - sensible heat net flux*,**
% - 'lhtfl' - latent heat net flux*,**
% - 'lon' - longitude (vector)
% - 'lat' - latitude (vector)