Commit 53a714f0 authored by Pierre Cazenave's avatar Pierre Cazenave

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

parents 73647a8c 1468701e
......@@ -103,7 +103,7 @@ end
'Cdata',squeeze(FVCOM.(plotOPTS.var_plot)(:,plotOPTS.nz_plot,aa)),'edgecolor','interp','facecolor','interp');
end
caxis(plotOPTS.clims)
colorbar h
colorbar('peer',gca,'EastOutside')
% check if mesh elements are required
if plotOPTS.do_mesh
%plot vertices
......
......@@ -141,7 +141,7 @@ netcdf.endDef(nc2read);
netcdf.putVar(nc2read, cweights_varid, nest.weight_cell);
netcdf.putVar(nc2read, nweights_varid, nest.weight_node');
netcdf.putVar(nc2read, nweights_varid, nest.weight_node);
catch e
fprintf(e.message)
error('Adding variable %s failed - does the variable already exist?', 'weight_cell')
......
......@@ -243,9 +243,11 @@ for obc_idx = 1:Mobj.nObs
% setdiff
if ~isempty(truecandidate)
for dd=1:length(truecandidate)
Nested.read_obc_elems{end} = cat(2,Nested.read_obc_elems{end},truecandidate(dd));
Nested.weight_cell{end} = cat(2,Nested.weight_cell{end},Nested.weight_cell{end}(end));
end
Nested.read_obc_elems{end} = cat(2,Nested.read_obc_elems{end},truecandidate(dd));
if conf.Nested_type(obc_idx) ~= 1
Nested.weight_cell{end} = cat(2,Nested.weight_cell{end},Nested.weight_cell{end}(end));
end
end
end
if ftbverbose
fprintf('\n')
......
......@@ -155,7 +155,7 @@ Mobj.siglevz = repmat(Mobj.h, 1, nlev) .* Mobj.siglev;
Mobj.siglayz = repmat(Mobj.h, 1, nlev-1) .* Mobj.siglay;
Mobj.siglevzc = repmat(Mobj.hc, 1, nlev) .* Mobj.siglevc;
Mobj.siglayzc = repmat(Mobj.hc, 1, nlev-1) .* Mobj.siglayc;
Mobj.Hmin = Hmin;
if ftbverbose
fprintf('done.\n')
fprintf('end : %s\n', subname)
......
function [ outputDTs ] = Times2Matlab( Times )
%TIMES2DATETIME Converts the 'Times' variable from a FVCOM output file to
%a Matlab Datetime array. This isn't gonna work if the FVCOM model isn't using
% Julian dates.
% Input: 'Times' variable, already read from a FVCOM netCDF output file.
% Output: Array of Matlab Datetime objects, one for each timestep.
% NB time zones are not considered.
% Simon Waldman / Marine Scotland Science 2018.
assert(ischar(Times), 'Input should be a character array.');
outputDTs = datetime( Times', 'InputFormat', 'yyyy-MM-dd''T''HH:mm:ss.SSSSSS');
end
function [Riv]=read_FVCOM_river_file(river_root,Mobj)
% function [Riv]=read_river_file(river_root,Mobj)
% Reads FVCOM river files
%
% function [Riv]=read_river_file(river_root,Mobj)
%
% DESCRIPTION:
% Reads both *.nc and *nml files associated with FVCOM river input files
%
% INPUT
% river_root = full address of river filename without extension
% Mobj = needs bathymetry, nodes, elements and triangulation table.
%
% OUTPUT:
% Riv = structure variable with all variables in the file
%
% EXAMPLE USAGE
% [Riv]=read_river_file('~/tapas_v0_riv_new',Mobj)
%
% Author(s):
% Ricardo Torres (Plymouth Marine Laboratory)
%
% Revision history
%
% 2018-07-22 First version.
%
%==============================================================================
[~, subname] = fileparts(mfilename('fullpath'));
global ftbverbose
if ftbverbose
fprintf('\nbegin : %s \n', subname)
end
% Check if they are in the same order... (a couple of manual test suggest
ncRivers = [river_root '.nc'];
fidnml = fopen ([river_root '.nml']);
% get variables from netcdf file
info = ncinfo(ncRivers);
% extract all variables
varnames={};
for nn =1:length(info.Variables)
varnames{nn}=info.Variables(nn).Name;
end
Riv=[];nn=1;
for vv=varnames
Riv.(vv{1}).data = ncread(ncRivers,vv{1});
if isfloat(Riv.(vv{1}).data)
Riv.(vv{1}).data =double(Riv.(vv{1}).data);
% extract units
for rr = 1:length(info.Variables(nn).Attributes)
if startsWith(info.Variables(nn).Attributes(rr).Name ,'units')
Riv.(vv{1}).units = info.Variables(nn).Attributes(rr).Value
end
end
end
nn=nn+1;
end
pos =0;
while ~feof(fidnml)
pos = pos+1;
lines=cell(3,1);
for dd=1:3
lines{dd} = fgetl(fidnml);
% fprintf(fout,'%s\n',line)
end
% keep name
[river_name]= extractAfter(lines{2},'= ');
Riv.river_name(pos,1)={river_name(1:end-1)};
% read grid location
line = fgetl(fidnml);
% extract grid location
t =strfind(line,'=');
Riv.IDX(pos,1) = str2double(line(t+1:end));
% check if IDX is in the list to remove
line1 = fgetl(fidnml);
line2 = fgetl(fidnml);
end
fclose (fidnml);
% extract river positions
for rr=1:length(Riv.IDX)
Riv.lon.data(rr) = Mobj.lon(Riv.IDX(rr,1));
Riv.lat.data(rr) = Mobj.lat(Riv.IDX(rr,1));
Riv.x.data(rr) = Mobj.x(Riv.IDX(rr,1));
Riv.y.data(rr) = Mobj.y(Riv.IDX(rr,1));
end
% %% for each river estimate the annual cycle mean and river
% % flow/concentration relationship
% % rnut = Riv2.N1_p;
% % read tamar nutrient observations.
% ncfile = fullfile('/data/euryale4/backup/mbe/Experiments/Rosa_rivers/Analysis/tamar_nutrient_nc','gunnislake_dayavg_%s.nc')
% vars = {'si','amm','phs','nitrate','o2'}
% varsR = {'N5_s','N4_n','N1_p','N3_n','O2_o'}
% tamar.time = ncread(sprintf(ncfile,vars{1}),'time');
% for nn=1:length(vars)
% tamar.([varsR{nn},'T']) = ncread(sprintf(ncfile,vars{nn}),'time');
% tamar.(varsR{nn}) = ncread(sprintf(ncfile,vars{nn}),vars{nn});
% end
%%
Riv.mtime.data = (Riv.time.data+ 678942);
% rr=10;
\ No newline at end of file
......@@ -196,10 +196,12 @@ try
Itime.ID=netcdf.inqVarID(nc,'Times');
Itime.sData=ncread(file_netcdf,'Times')
Itime.Data(1) = datenum(Itime.sData(:,1)','yyyy-mm-ddTHH:MM:SS');
Itime.Data(2) = datenum(Itime.sData(:,end-1)','yyyy-mm-ddTHH:MM:SS');
% Itime.Data(2) = datenum(Itime.sData(:,end-1)','yyyy-mm-ddTHH:MM:SS');
Itime.Data(2) = datenum(Itime.sData(:,end)','yyyy-mm-ddTHH:MM:SS');
start_date= Itime.Data(1);
end_date = Itime.Data(2);
var_time = datenum(Itime.sData(:,1:end-1)','yyyy-mm-ddTHH:MM:SS');
% var_time = datenum(Itime.sData(:,1:end-1)','yyyy-mm-ddTHH:MM:SS');
var_time = datenum(Itime.sData(:,1:end)','yyyy-mm-ddTHH:MM:SS');
% Itime.idx=find(strcmpi(vars,'Itime'));
% Itime.ID=netcdf.inqVarID(nc,'Itime');
% Itime.Data(1) = netcdf.getVar(nc,Itime.ID,0,1,'int32');
......@@ -225,7 +227,8 @@ catch me
Itime.idx=find(strcmpi(vars,'time'));
Itime.ID=netcdf.inqVarID(nc,'time');
Itime.Data(1) = netcdf.getVar(nc,Itime.ID,0,1,'double');
Itime.Data(2) = netcdf.getVar(nc,Itime.ID,last_entry-1,1,'double');
% Itime.Data(2) = netcdf.getVar(nc,Itime.ID,last_entry-1,1,'double');
Itime.Data(2) = netcdf.getVar(nc,Itime.ID,last_entry,1,'double');
[start_date,end_date] = deal(Itime.Data(1)+time_offset,Itime.Data(end)+time_offset);
DeltaT=(end_date-start_date)./last_entry;
var_time = start_date:DeltaT:(end_date-DeltaT);
......@@ -505,17 +508,20 @@ for aa=1:length(varnames)
read_stride(do_time)=stride.(cc_names{do_time});
eval([varnames{aa},'=netcdf.getVar(nc,varID,read_start,read_count,read_stride,''double'');'])
else % we are looking at stations or depth layers
for cc=1:length(start.(cc_names{(do_restrict)}))
for cc=1:length(start.(cc_names{find(do_restrict)}))
read_start(find(do_restrict))=start.(cc_names{find(do_restrict)})(cc);
read_count(find(do_restrict))=count.(cc_names{find(do_restrict)})(cc);
read_stride(find(do_restrict))=stride.(cc_names{find(do_restrict)});
var_dump=netcdf.getVar(nc,varID,read_start,read_count,read_stride,'double');
switch dimName(find(do_restrict))
switch dimName{find(do_restrict)}
case 'node' | 'nele'
eval([varnames{aa},'(cc,:,:)=var_dump;'])
case 'siglay' | 'siglev'
eval([varnames{aa},'(:,cc,:)=var_dump;'])
case 'time' % this only happens if we only have one timestamp
eval([varnames{aa},'=var_dump;'])
end
clear var_dump
end
......
function FVCOM = restrict_spatial_indices(FVCOM,mask_nodes,mask_elems);
% Eliminates the FVCOM nodes and elements in the lists mask_nodes and
% mask_elems.
%
% function FVCOM = restrict_spatial_indices(FVCOM,mask_nodes,mask_elems);
%
% DESCRIPTION:
% Loops through all variables in FVCOM and restricts spatial dimensions
%
%
% INPUT:
% FVCOM = Structure variable with FVCOM output data
% mask_nodes = list of node indices to remove
% mask_elems = list of elements indices to remove
%
% OUTPUT:
% FVCOM with fewer nodes and elements!
%
% EXAMPLE USAGE:
% FVCOM.temp = temperature field (node,levels,times)
% FVCOM.u = velocity field (elements,levels,times)
% mask_nodes = (1:300) can be the boundary/nesting nodes
% mask_elems = (1:400) can be the boundary/nesting elements
% FVCOM = restrict_spatial_indices(FVCOM,mask_nodes,mask_elems);
%
% Author(s):
% Ricardo Torres (Plymouth Marine Laboratory)
%
% Revision history:
% 2018-09-17 First version
%
%==========================================================================
[~, subname] = fileparts(mfilename('fullpath'));
global ftbverbose
if ftbverbose
fprintf('\nbegin : %s\n', subname)
end
vnames = fields (FVCOM);
if isfield(FVCOM,'x');
nodes = length(FVCOM.x);
elseif isfield(FVCOM,'lon')
nodes = length(FVCOM.lon);
else
warning('No easily identifiable variable with node dimensions positions... e.g. x/y or lon/lat and I will continue')
end
if isfield(FVCOM,'xc');
elems = length(FVCOM.xc);
elseif isfield(FVCOM,'lonc')
elems = length(FVCOM.lonc);
else
warning('No easily identifiable variable with element dimensions positions... e.g. xc/yc or lonc/latc and I cannot continue')
end
if exist('nodes','var')
else
nodes=0;
end
if exist('elems','var')
else
elems=0;
end
for vv=vnames'
switch size(FVCOM.(vv{1}),1) % In FVCOM variable structure, the first dimension is always the spatial dimension if it is present
case nodes
disp(['Clipping variable FVCOM.',vv{1}])
switch ndims(FVCOM.(vv{1}))
case 1
FVCOM.(vv{1})(mask_nodes)=[];
case 2
FVCOM.(vv{1})(mask_nodes,:)=[];
case 3
FVCOM.(vv{1})(mask_nodes,:,:)=[];
end
case elems
disp(['Clipping variable FVCOM.',vv{1}])
switch ndims(FVCOM.(vv{1}))
case 1
FVCOM.(vv{1})(mask_elems)=[];
case 2
FVCOM.(vv{1})(mask_elems,:)=[];
case 3
FVCOM.(vv{1})(mask_elems,:,:)=[];
end
end
end
function dist = sigma_tanh(nlev,dl,du)
function [Mobj] = sigma_tanh(nlev,dl,du,Mobj,sigma_file)
% Generate a tanh sigma coordinate distribution.
%
% Mobj = sigma_tanh(nlev, dl, du)
% Mobj = sigma_tanh(nlev, dl, du,Mobj)
%
% DESCRIPTION:
% Generate a tanh vertical sigma coordinate distribution.
......@@ -12,12 +12,15 @@ function dist = sigma_tanh(nlev,dl,du)
% coordinates are parallel with uniform thickness.
% du: The upper depth boundary from the surface, up to which the
% coordinates are parallel with uniform thickness.
% Mobj: [optional] mesh object file
%
% OUTPUT:
% dist: Tanh vertical sigma coordinate distribution.
% Mobj.sigma: Tanh vertical sigma coordinate distribution.
% Mobj.sig... All the sigma layers variables such as siglev, siglevz,
% siglayz, etc...
%
% EXAMPLE USAGE:
% Mobj = read_sigma(nlev, dl, du)
% Mobj = read_sigma(nlev, dl, du,Mobj)
%
% Author(s):
% Geoff Cowles (University of Massachusetts Dartmouth)
......@@ -37,3 +40,41 @@ for k = 1:nlev-1
x3 = x2+tanh(du);
dist(k+1) = (x1+x2)/x3-1.0;
end
if nargin >= 4
conf.nlev = nlev;
Mobj.sigma = dist;
Mobj.siglev = zeros(Mobj.nVerts,conf.nlev);
Mobj.siglevc = zeros(Mobj.nElems,conf.nlev);
Mobj.siglayc = zeros(Mobj.nElems,conf.nlev-1);
Mobj.siglev = repmat(Mobj.sigma,Mobj.nVerts,1);
Mobj.siglay = Mobj.siglev(:,1:end-1) + (diff(Mobj.siglev,1,2)/2);
for zz = 1:conf.nlev-1
Mobj.siglevc(:, zz) = nodes2elems(Mobj.siglev(:, zz), Mobj);
Mobj.siglayc(:, zz) = nodes2elems(Mobj.siglay(:, zz), Mobj);
end
% An extra level compared with layers.
Mobj.siglevc(:, zz + 1) = nodes2elems(Mobj.siglev(:, zz + 1), Mobj);
% Finally, make some [depth, sigma] arrays.
Mobj.siglevz = repmat(Mobj.h, 1, conf.nlev) .* Mobj.siglev;
Mobj.siglayz = repmat(Mobj.h, 1, conf.nlev-1) .* Mobj.siglay;
if isfield(Mobj, 'hc')
Mobj.siglevzc = repmat(Mobj.hc, 1, conf.nlev) .* Mobj.siglevc;
Mobj.siglayzc= repmat(Mobj.hc, 1, conf.nlev-1) .* Mobj.siglayc;
end
else
Mobj = dist;
end
% generate sigma file
% Save to the given file name.
if nargin==5
fout = fopen(sigma_file, 'wt');
assert(fout >= 0, 'Error opening sigma file: %s', sigma_file)
fprintf(fout, 'NUMBER OF SIGMA LEVELS = %d\n', nlev);
fprintf(fout, 'SIGMA COORDINATE TYPE = TANH\n');
fprintf(fout, 'DU = %4.1f\n', du);
fprintf(fout, 'DL = %4.1f\n', dl);
fprintf(fout,'\n');
fclose(fout);
end
return
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment