Commit b8dc34cf authored by Pierre Cazenave's avatar Pierre Cazenave

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

parents 4387bdb1 64795aaf
% Create the necessary files for FVCOM from a given SMS input unstructured
% grid.
%
% This version uses the following forcings:
% - E-HYPE river climatology discharge based on CEH rivers
% - EA river temperature climatology
% - River salinity based on my river morphology classification
% - CFS surface forcing (wind, heat, preciptation/evaporation)
% - HYCOM open boundary temperature/salinity and initial conditions
% - TPXO predicted tidal elevations at the boundary
%
% This requires:
% - my fork of the fvcom-toolbox
% > https://gitlab.em.pml.ac.uk/pica/fvcom-toolbox
% - the Tide Model Driver toolbox and data
% > http://polaris.esr.org/ptm_index.html
%
% Author(s):
% Pierre Cazenave (Plymouth Marine Laboratory)
%
% Revision history:
%
% 2015-05-15 Poached the Irish Sea domain version of this script. Removed
% the version history as it was a bit long and pointless.
% 2015-08-14 Use ECMWF-ERA20C forcing instead of NCEP.
% matlabrc
close all
clc
cd('/users/modellers/rito/Models/git/fvcom/locate-implementation/loc_v0/grids')
global ftbverbose
ftbverbose = 1; % be noisy
addpath('/users/modellers/rito/Models/git/fvcom/fvcom-toolbox/utilities')
addpath('/users/modellers/rito/Models/git/fvcom/fvcom-toolbox/fvcom_prepro/')
% addpath('/users/modellers/pica/Code/MATLAB/func/')
addpath('/users/modellers/rito/matlab/air-sea')
conf.obc_tides.dirTMD = '/users/modellers/pica/Code/MATLAB/toolboxes/TMD/';
addpath(conf.obc_tides.dirTMD)
addpath(fullfile(conf.obc_tides.dirTMD, 'FUNCTIONS'))
addpath('scriptlets')
%%%------------------------------------------------------------------------
%%% INPUT CONFIGURATION
%%%
%%% Define the model input parameters here e.g. number of tidal components,
%%% type of boundary forcing, estimated current velocity, sponge layer
%%% coefficient and radius, forcing source etc.
%%%
%%%------------------------------------------------------------------------
conf.base = '/users/modellers/rito/Models/git/fvcom/locate-implementation/loc_v0/grids';
% Which version of FVCOM are we using (for the forcing file formats)?
conf.FVCOM_version = '3.2.2';
%%%------------------------------------------------------------------------
%%% Spatial stuff
%%%------------------------------------------------------------------------
% Case name for the model inputs and outputs
conf.casename = 'lyme_bay_v04';
conf.coordType = 'cartesian'; % 'cartesian' or 'spherical'
% Input grid UTM Zone (if applicable)
conf.utmZone = {'30 U'}; % syntax for utm2deg
% Give some names to the boundaries. This must match the number of node
% strings defined in SMS. Ideally, the order of the names should match the
% order in which the boundaries were made in SMS.
conf.boundaryNames = {'South'}; % Number is (relatively) important!
%%%------------------------------------------------------------------------
%%% Model constants and stuff
%%%------------------------------------------------------------------------
% Type of open boundary treatment (see Table 6.1 (?) in the manual).
% 1 - Active (ASL): sea level is specified at the open boundary.
% 2 - Clamped: zeta = 0 at the boundary (Beardsley and Haidvogel, 1981).
% 3 - Implicit Gravity Wave Radiation.
% 4 - Partial Clamped Gravity Wave Radiation (Blumberg and Kantha, 1985).
% 5 - Explicit Orlanski Radiation (Orlanski, 1976; Chapman, 1985)
conf.obc_type = 1;
% Sponge layer parameters
conf.sponge.radius = 15000; % in metres, or -1 for variable width.
conf.sponge.coeff = 0.001;
% z0 value in metres (for uniform) or 'variable' (requires a separate
% roughness (z0) distribution grid).
conf.bedRoughness = 0.035; % or 0.015, 0.025 or 0.03 - Davies and Jones (1990) shelf model
% Estimated velocity (m/s) and tidal range (m) for time step estimate
conf.estVel = 4;
conf.estRange = 14;
%%%------------------------------------------------------------------------
%%% END OF INPUT CONFIGURATION
%%%------------------------------------------------------------------------
%% Process the files needed for all months.
% Read the input mesh and bathymetry. Also creates the data necessary for
% the Coriolis correction in FVCOM.
Mobj = read_sms_mesh(...
'2dm', fullfile(conf.base,conf.casename, [conf.casename, '.2dm']),...
'coordinate', conf.coordType, 'addCoriolis', false);
% 'bath', fullfile(conf.base,conf.casename, [conf.casename, '.dat']),...
% Clean out some unneeded fields.
Mobj = rmfield(Mobj, 'riv_nodes');
% Convert the coordinates to whatever is missing.
if Mobj.have_lonlat == 0
utmZones = cellfun(@(x) repmat(x, length(Mobj.x), 1), conf.utmZone, 'uni', false);
[Mobj.lat, Mobj.lon] = utm2deg(Mobj.x, Mobj.y, utmZones{1});
Mobj.have_lonlat = true;
clear utmZones
end
if Mobj.have_xy == 0
tmpZone = regexpi(conf.utmZone, '\ ', 'split');
[Mobj.x, Mobj.y] = wgs2utm(Mobj.lon, Mobj.lon, str2double(char(tmpZone{1}(1))), 'N');
Mobj.have_xy = true;
clearvars tmpZone
end
% Add grid metrics.
Mobj = setup_metrics(Mobj);
% write ADMESH 14 format file
% Create a Coriolis file from the bathy which varies with latitude. Given
% the size of the domain, this is probably necessary.
if ~Mobj.have_cor
Mobj = add_coriolis(Mobj, 'uselatitude');
end
% Parse the open boundary nodes and add accordingly.
if Mobj.have_strings
for i = 1:size(Mobj.read_obc_nodes, 2)
nodeList = double(cell2mat(Mobj.read_obc_nodes(i)));
Mobj = add_obc_nodes_list(Mobj, nodeList, conf.boundaryNames{i}, conf.obc_type);
end
clear nodeList
end
% Create a sponge layer
if Mobj.have_strings
for i = 1:size(Mobj.read_obc_nodes, 2)
% Get the list of nodes in this boundary
nodeList = double(cell2mat(Mobj.read_obc_nodes(i)));
if conf.sponge.radius < 0 % if we want a variable sponge radius
if i==1
% Create an array to store the radii
Mobj.sponge_rad = zeros(size(Mobj.sponge_nodes));
end
% calculate the sponge radius
conf.spongeRadius = calc_sponge_radius(Mobj, nodeList);
% Add the sponge nodes to the list
Mobj = add_sponge_nodes_list(Mobj,nodeList,...
[conf.boundaryNames{i}, ' sponge'], conf.sponge.radius,...
conf.sponge.coeff);
else
Mobj = add_sponge_nodes_list(Mobj,nodeList,...
[conf.boundaryNames{i}, ' sponge'], conf.sponge.radius,...
conf.sponge.coeff);
end
clear nodeList
end
end
[Mobj] = write_admesh_mesh(Mobj,'output_directory',fullfile(conf.base,conf.casename),'native_coord','spherical');
% [Mobj] = write_admesh_mesh(Mobj,varargin)
save locate_sms_v0 Mobj
\ No newline at end of file
% convert admesh edge matlab file from cartesian to spherical
base_dir = '~/Models/git/fvcom/locate-implementation/loc_v0/grids/lymebay/admesh/'
addpath '/users/modellers/rito/Code/fvcom-toolbox/utilities/'
admesh_file = 'Locate_500m'
Mobj=read_admesh_mesh('msh',fullfile(base_dir,[admesh_file,'.14']),'coordinate','spherical');
write_SMS_2dm(fullfile(base_dir,[admesh_file,'.2dm']),Mobj.tri,Mobj.lon,Mobj.lat,Mobj.h,[])
z1 = utmzone(50,-4);
[ellipsoid,estr] = utmgeoid(z1);
utmstruct = defaultm('utm');
utmstruct.zone = z1;
utmstruct.geoid = ellipsoid;
utmstruct = defaultm(utmstruct);
%coast = load([base_dir 'coast_rosa.mat']);
[Mobj.x,Mobj.y]=mfwdtran(utmstruct,Mobj.lat,Mobj.lon);
Mobj.have_xy = 1;
Mobj.have_lonlat = 0;
[Mobj] = write_admesh_mesh(Mobj,'output_directory',base_dir,'native_coord','cartesian');
write_SMS_2dm(fullfile(base_dir,[admesh_file,'admesh.2dm']),Mobj.tri,Mobj.x,Mobj.y,Mobj.h,[])
% read cartesian 2dm mesh from sms to add y-shift appearing after admesh
% processing
Mobj = read_sms_mesh(...
'2dm', fullfile('~/Models/git/fvcom/locate-implementation/loc_v0/grids/', ['lyme_bay_v05_post_admesh', '.2dm']),...
'coordinate', 'cartesian', 'addCoriolis', false);
% 'bath', fullfile(conf.base,conf.casename, [conf.casename, '.dat']),...
% Read nodestring of buffer area from ROSA big mesh
base_dir = '/users/modellers/rito';
% base_dir = '~/Models/FVCOM/fvcom-projects/tapas/grids/lymebay/admesh/'
project_dir = fullfile(base_dir, '/Models/FVCOM/fvcom-projects/conway/');
cd(fullfile(project_dir, 'matlab'))
addpath(fullfile(base_dir, 'Code/fvcom-toolbox/utilities'))
addpath(fullfile(base_dir, 'Code/fvcom-toolbox/fvcom_prepro/'))
%%%------------------------------------------------------------------------
%%% INPUT CONFIGURATION
%%%
%%% Define the model input parameters here e.g. number of tidal components,
%%% type of boundary forcing, estimated current velocity, sponge layer
%%% coefficient and radius, forcing source etc.
%%%
%%%------------------------------------------------------------------------
conf.base = fullfile(project_dir, 'run');
%%%------------------------------------------------------------------------
%%% Time stuff
%%%------------------------------------------------------------------------
% Case name for the model inputs and outputs
conf.casename = 'conway_v0_aqua_v16_band_nest_large_single_ND';
conf.coordinates = 'cartesian'; % 'cartesian' or 'spherical'
% conf.run_coordinates = conf.coordinates;
% Input grid UTM Zone (if applicable)
Mobj = read_sms_mesh(...
'2dm', fullfile(conf.base, '..', 'grids', [conf.casename, '.2dm']),...
'coordinate', conf.coordinates, 'addCoriolis', false);
Boundary.nodes = [Mobj.read_obc_nodes{:}];
Boundary.x = Mobj.x(Boundary.nodes);
Boundary.y = Mobj.y(Boundary.nodes);
write_SMS_cst('../grids/buffer_nodes_conway_aqua_v16.cst',Boundary.x,Boundary.y)
%% Now convert the high res subdomain nodestring to build the buffer area between coarse domain and highres domain
cd(fullfile(project_dir, 'matlab'))
conf.base = fullfile(project_dir, 'run');
% Case name for the model inputs and outputs
conf.casename = 'conway_v1';
conf.coordinates = 'cartesian'; % 'cartesian' or 'spherical'
Mobj = read_sms_mesh(...
'2dm', fullfile(conf.base, '..', 'grids', [conf.casename, '.2dm']),...
'coordinate', conf.coordinates, 'addCoriolis', false);
Boundary.nodes = [Mobj.read_obc_nodes{:}];
Boundary.x = Mobj.x(Boundary.nodes);
Boundary.y = Mobj.y(Boundary.nodes);
write_SMS_cst('../grids/boundary_nodes_conway_v1.cst',Boundary.x,Boundary.y)
This diff is collapsed.
function [field] = smoothfield(fieldin,Mobj,nLoops,SmoothPts) function [field] = smoothfield2(fieldin,Mobj,nLoops,SmoothPts)
% Smooth a vertex-based field using minimum value of surrounding nodes % Smooth a vertex-based field using minimum value of surrounding nodes
% %
...@@ -25,8 +25,7 @@ function [field] = smoothfield(fieldin,Mobj,nLoops,SmoothPts) ...@@ -25,8 +25,7 @@ function [field] = smoothfield(fieldin,Mobj,nLoops,SmoothPts)
% Revision history % Revision history
% %
%============================================================================== %==============================================================================
fprintf('\nbegin : %s\n', subname) global ftbverbose
endglobal ftbverbose
[~, subname] = fileparts(mfilename('fullpath')); [~, subname] = fileparts(mfilename('fullpath'));
if ftbverbose if ftbverbose
fprintf('\nbegin : %s\n', subname) fprintf('\nbegin : %s\n', subname)
...@@ -43,8 +42,8 @@ end; ...@@ -43,8 +42,8 @@ end;
if(exist('SmoothPts')) if(exist('SmoothPts'))
nPts = length(SmoothPts); nPts = length(SmoothPts);
else else
nPts = Mobj.nVerts; nPts = Mobj.nElems;
SmoothPts = 1:Mobj.nVerts; SmoothPts = 1:Mobj.nElems;
end; end;
if(~Mobj.have_mets) if(~Mobj.have_mets)
......
...@@ -10,6 +10,7 @@ function write_FVCOM_z0(z0,filename,mytitle,cbcmin) ...@@ -10,6 +10,7 @@ function write_FVCOM_z0(z0,filename,mytitle,cbcmin)
% INPUT % INPUT
% z0 = user defined roughness field (m) % z0 = user defined roughness field (m)
% roughness is defined on the elements % roughness is defined on the elements
% expect values between 3 10^-3 (gravel) and .2 10^-3 (i.e. 0.0002) for mud
% filename = filename to dump to % filename = filename to dump to
% mytitle = title of the case (set as global attribute) % mytitle = title of the case (set as global attribute)
% cbcmin = minimum value of CBC (optional, defaults to 0.0018 if % cbcmin = minimum value of CBC (optional, defaults to 0.0018 if
......
...@@ -48,7 +48,12 @@ function [Plots]=do_surface_plotMatlabMap(plotOPTS,FVCOM) ...@@ -48,7 +48,12 @@ function [Plots]=do_surface_plotMatlabMap(plotOPTS,FVCOM)
% adds m_map to matlab paths. file is in utilities directory. % adds m_map to matlab paths. file is in utilities directory.
% amend according to your m_map installation paths % amend according to your m_map installation paths
figure(plotOPTS.figure);clf figure(plotOPTS.figure);
try delete(plotOPTS.PlotoutS(plotOPTS.figure).handles)
catch
clf
end
% generate figure with correct projection lat and lon range ellipsoid and % generate figure with correct projection lat and lon range ellipsoid and
% zone. % zone.
if isfield(plotOPTS,'Lontick') if isfield(plotOPTS,'Lontick')
...@@ -68,8 +73,12 @@ axesm('mercator','MapLatLimit',plotOPTS.range_lat,'MapLonLimit',[plotOPTS.range_ ...@@ -68,8 +73,12 @@ axesm('mercator','MapLatLimit',plotOPTS.range_lat,'MapLonLimit',[plotOPTS.range_
% add coastline if present % add coastline if present
if (isfield(plotOPTS,'coastline_file') && ~isempty(plotOPTS.coastline_file) ) if (isfield(plotOPTS,'coastline_file') && ~isempty(plotOPTS.coastline_file) )
coast=load(plotOPTS.coastline_file); if isfield (plotOPTS,'PlotoutS') && ~isempty(plotOPTS.PlotoutS(plotOPTS.figure).handles)
geoshow([coast.ncst(:,2)],[coast.ncst(:,1)],'Color','black') else
% coast=load(plotOPTS.coastline_file);
geoshow(plotOPTS.coastline_file,'Color','black')
end
end end
%------------------------------------------------------------------------------ %------------------------------------------------------------------------------
......
...@@ -47,6 +47,10 @@ function [Plots]=do_vector_plot_MatlabMap(plotOPTS,FVCOM) ...@@ -47,6 +47,10 @@ function [Plots]=do_vector_plot_MatlabMap(plotOPTS,FVCOM)
% adds m_map to matlab paths. file is in utilities directory. % adds m_map to matlab paths. file is in utilities directory.
% ammend according to your m_map installation paths % ammend according to your m_map installation paths
figure(plotOPTS.figure); figure(plotOPTS.figure);
try delete(plotOPTS.PlotoutV(plotOPTS.figure).handles{1})
catch
clf
end
% generate figure with correct projection lat and lon range ellipsoid and % generate figure with correct projection lat and lon range ellipsoid and
% zone. % zone.
if isfield(plotOPTS,'Lontick') if isfield(plotOPTS,'Lontick')
...@@ -61,15 +65,18 @@ else ...@@ -61,15 +65,18 @@ else
end end
axesm('mercator','MapLatLimit',plotOPTS.range_lat,'MapLonLimit',[plotOPTS.range_lon],'MeridianLabel','on',...
'ParallelLabel','on','MLineLocation',MerTick,'PLineLocation',ParTick,'LabelUnits','dm')
H=axesm('mercator','MapLatLimit',plotOPTS.range_lat,'MapLonLimit',[plotOPTS.range_lon],'MeridianLabel','on',...
'ParallelLabel','on','MLineLocation',MerTick,'PLineLocation',ParTick,'LabelUnits','dm')
% add coastline if present % add coastline if present
if (isfield(plotOPTS,'coastline_file') && ~isempty(plotOPTS.coastline_file) ) if (isfield(plotOPTS,'coastline_file') && ~isempty(plotOPTS.coastline_file) )
coast=load(plotOPTS.coastline_file); if isfield (plotOPTS,'PlotoutV') && ~isempty(plotOPTS.PlotoutV(plotOPTS.figure).handles)
geoshow([coast.ncst(:,2)],[coast.ncst(:,1)],'Color','black') else
% coast=load(plotOPTS.coastline_file);
geoshow(plotOPTS.coastline_file,'Color','black','LineWidth',2)
setm (H,'MapLatLimit',plotOPTS.range_lat,'MapLonLimit',[plotOPTS.range_lon])
end
end end
%----------------------------------------------------------------------- %-----------------------------------------------------------------------
% Convert element positions from FVCOM to lat and lon using m_ll2ll.m from % Convert element positions from FVCOM to lat and lon using m_ll2ll.m from
......
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