Commit fff744bc authored by Ricardo Torres's avatar Ricardo Torres

new function to read FVCOM river files

parent ef72117c
......@@ -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
......
......@@ -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 [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 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