Commit d57ddf68 authored by Ricardo Torres's avatar Ricardo Torres 💬

Added scripts to read polcoms gcoms files (physics and ersem) to generate boundary conditions.

parent d8dd703c
function write_FVCOM_river_ERSEM(RiverFile,RiverName,time,flux,temp,salt,n1p,n3n,n4n,n5s,RiverInfo1,RiverInfo2)
% Write FVCOM 3.x NetCDF river file
%
% function write_FVCOM_river(RiverFile,RiverName,time,flux,temp,salt,RiverInfo1,RiverInfo2)
%
% DESCRIPTION:
% Write river flux, temperature, and salinity to an FVCOM river file.
% Flux, temperature and salinity must be calculated prior to being given
% here as the raw values in the arrays are simply written out as is to
% the NetCDF file.
%
% INPUT
% RiverFile : FVCOM 3.x NetCDF river forcing file
% RiverName : Name of the actual River
% time : Timestamp array in modified Julian day
% flux : Total river flux in m^3/s (dimensions [time, nRivernodes])
% temp : Temperature in C (dimensions [time, nRivernodes])
% salt : Salinity in PSU (dimensions [time, nRivernodes])
% RiverInfo1 : Global attribute title of file
% RiverInfo2 : Global attribute info of file
%
% OUTPUT:
% FVCOM NetCDF river file with flux, temperature and salinity.
%
% EXAMPLE USAGE
% write_FVCOM_river('tst_riv.nc', {'Penobscot'}, time, flux, temp, ...
% salt, 'Penobscot Flux', 'source: USGS')
%
% Author(s):
% Geoff Cowles (University of Massachusetts Dartmouth)
% Pierre Cazenave (Plymouth Marine Laboratory)
%
% Revision history
% 2013-03-21 Modified to take a list of river nodes rather than a single
% river spread over multiple nodes. This means you have to scale your
% inputs prior to using this function. This also means I have broken
% backwards compatibility with the old way of doing it (i.e. this
% function previously wrote only a single river's data but spread over a
% number of nodes). I removed the sediment stuff as the manual makes no
% mention of this in the river input file. Also added support for writing
% to NetCDF using MATLAB's native tools.
% 2013-03-21 Transpose the river data arrays to the correct shape for the
% NetCDF file.
%
%==========================================================================
subname = 'write_FVCOM_river';
global ftbverbose
if ftbverbose
fprintf('\nbegin : %s\n', subname)
end
[nTimes, nRivnodes] = size(flux);
if ftbverbose
fprintf('# of river nodes: %d\n', nRivnodes);
fprintf('# of time frames: %d\n', nTimes);
end
[year1, month1, day1, ~, ~, ~] = mjulian2greg(time(1));
[year2, month2, day2, ~, ~, ~] = mjulian2greg(time(end));
if ftbverbose
fprintf('time series begins at:\t%04d %02d %02d\n', year1, month1, day1);
fprintf('time series ends at:\t%04d %02d %02d\n', year2, month2, day2);
end
clear year? month? day?
%--------------------------------------------------------------------------
% dump to netcdf file
%--------------------------------------------------------------------------
% river node forcing
nc = netcdf.create(RiverFile, 'clobber');
% global variables
netcdf.putAtt(nc, netcdf.getConstant('NC_GLOBAL'), 'type', 'FVCOM RIVER FORCING FILE')
netcdf.putAtt(nc, netcdf.getConstant('NC_GLOBAL'), 'title', RiverInfo1)
netcdf.putAtt(nc, netcdf.getConstant('NC_GLOBAL'), 'info', RiverInfo2)
netcdf.putAtt(nc, netcdf.getConstant('NC_GLOBAL'), 'history', 'File created using write_FVCOM_river.m from the MATLAB fvcom-toolbox')
% dimensions
namelen_dimid = netcdf.defDim(nc, 'namelen', 80);
rivers_dimid = netcdf.defDim(nc, 'rivers', nRivnodes);
time_dimid = netcdf.defDim(nc, 'time', netcdf.getConstant('NC_UNLIMITED'));
date_str_len_dimid = netcdf.defDim(nc, 'DateStrLen', 26);
% variables
river_names_varid = netcdf.defVar(nc, 'river_names', 'NC_CHAR', [namelen_dimid, rivers_dimid]);
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, 'format', 'modified julian day (MJD)');
netcdf.putAtt(nc, time_varid, 'time_zone', 'UTC');
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');
times_varid = netcdf.defVar(nc,'Times','NC_CHAR',[date_str_len_dimid, time_dimid]);
netcdf.putAtt(nc, times_varid, 'time_zone','UTC');
river_flux_varid = netcdf.defVar(nc, 'river_flux', 'NC_FLOAT', [rivers_dimid, time_dimid]);
netcdf.putAtt(nc, river_flux_varid, 'long_name', 'river runoff volume flux');
netcdf.putAtt(nc, river_flux_varid, 'units', 'm^3s^-1');
river_temp_varid = netcdf.defVar(nc, 'river_temp', 'NC_FLOAT', [rivers_dimid, time_dimid]);
netcdf.putAtt(nc, river_temp_varid, 'long_name', 'river runoff temperature');
netcdf.putAtt(nc, river_temp_varid, 'units', 'Celsius');
river_salt_varid = netcdf.defVar(nc, 'river_salt', 'NC_FLOAT', [rivers_dimid, time_dimid]);
netcdf.putAtt(nc, river_salt_varid, 'long_name', 'river runoff salinity');
netcdf.putAtt(nc, river_salt_varid, 'units', 'PSU');
river_n1p_varid = netcdf.defVar(nc, 'N1p', 'NC_FLOAT', [rivers_dimid, time_dimid]);
netcdf.putAtt(nc, river_n1p_varid, 'long_name', 'river phosphate concentrations');
netcdf.putAtt(nc, river_n1p_varid, 'units', 'mmolm^-3');
river_n3n_varid = netcdf.defVar(nc, 'N3n', 'NC_FLOAT', [rivers_dimid, time_dimid]);
netcdf.putAtt(nc, river_n3n_varid, 'long_name', 'river Nitrate concentrations');
netcdf.putAtt(nc, river_n3n_varid, 'units', 'mmolm^-3');
river_n4n_varid = netcdf.defVar(nc, 'N4n', 'NC_FLOAT', [rivers_dimid, time_dimid]);
netcdf.putAtt(nc, river_n4n_varid, 'long_name', 'river ammonium concentrations');
netcdf.putAtt(nc, river_n4n_varid, 'units', 'mmolm^-3');
river_n5s_varid = netcdf.defVar(nc, 'N5s', 'NC_FLOAT', [rivers_dimid, time_dimid]);
netcdf.putAtt(nc, river_n5s_varid, 'long_name', 'river silicate concentrations');
netcdf.putAtt(nc, river_n5s_varid, 'units', 'mmolm^-3');
% end definitions
netcdf.endDef(nc);
% river names (must be 80 character strings)
rString = char();
for i = 1:nRivnodes
% Left-aligned 80 character string.
rString = [rString, sprintf('%-80s', RiverName{i})];
end
netcdf.putVar(nc, river_names_varid, rString);
% dump dynamic data
netcdf.putVar(nc, time_varid, 0, nTimes, time);
netcdf.putVar(nc, itime_varid, 0, nTimes, floor(time));
netcdf.putVar(nc, itime2_varid, 0, nTimes, mod(time, 1)*24*3600*1000);
netcdf.putVar(nc, river_flux_varid, flux');
netcdf.putVar(nc, river_temp_varid, temp');
netcdf.putVar(nc, river_salt_varid, salt');
netcdf.putVar(nc, river_n1p_varid, n1p');
netcdf.putVar(nc, river_n3n_varid, n3n');
netcdf.putVar(nc, river_n4n_varid, n4n');
netcdf.putVar(nc, river_n5s_varid, n5s');
%n1p,n3n,n4n,n5s
% build the time string and output to NetCDF.
nStringOut = char();
[nYr, nMon, nDay, nHour, nMin, nSec] = mjulian2greg(time);
for tt = 1:nTimes
nDate = [nYr(tt), nMon(tt), nDay(tt), nHour(tt), nMin(tt), nSec(tt)];
nStringOut = [nStringOut, sprintf('%04i/%02i/%02i %02i:%02i:%02i ', nDate)];
end
netcdf.putVar(nc, times_varid, nStringOut);
netcdf.close(nc);
if ftbverbose
fprintf(['end : ' subname '\n'])
end
function write_FVCOM_tsobcERSEM(basename,time,nSiglay,Mobj,Nutrients)
% Export temperature and salinity forcing at the open boundary.
%
% function write_FVCOM_tsobc(basename,time,nSiglay,in_temp,in_salt)
%
% DESCRIPTION:
% Setup an FVCOM hydrographic open boundary forcing file. Supply either
% uniform values for temperature and salinity or 3D arrays (node,
% sigma_layers, time).
%
% INPUT
% Model case name (to find the bathymetry and open boundary .dat files).
% Time
% Number of sigma layers
% Mobj (optional)
% Nutrients
%
% OUTPUT:
% FVCOM ERSEM open boundary file. No temp or salt added to file!!
%
% Author(s):
% Geoff Cowles (University of Massachusetts Dartmouth)
% Pierre Cazenave (Plymouth Marine Laboratory)
% Karen Amoudry (National Oceanography Centre, Liverpool)
% Ricardo Torres (Plymouth Marine Laboratory)
% PWC Revision history
% 2012-06-15 Added support for native MATLAB NetCDF routines. Requires
% MATLAB 2010a or higher.
% 2012-07-16 Removed hard-coded nSiglay and nSiglev and instead moved to
% arguments list.
% 2012-10-08 Updated help to reflect the fact nSiglev is calculated as
% nSiglay+1.
% 2012-11-09 Add new arguments to use user defined temperature and
% salinity.
% 2013-01-09 Add support for 3D input temperature and salinity (such as
% might be generated with get_POLCOMS_tsobc.m.
%
% KJA Revision history
% Undated - Add better check for the size of the input arrays (works with
% scalars).
% 2013-08-16 - Updated output of Itime2 to avoid rounding errors
% when converting from double to single format.
% RJT Revision history
% 2013-12-05 Added functionality to output ERSEM nutrients
%==============================================================================
if nargin == 4
warning(['Assuming uniform terrain-following sigma coordinates. ',...
'To specify different sigma coordintes, supply a MATLAB mesh ',...
'structure with fields siglay and siglev.'])
disp(['No Nutrient fields provided! no nutrients to be output!!! Aborting'])
return
end
if nargin == 5
NNuts=length(Nutrients)./3
if NNuts==0
disp(['Nutrient fields are empty! no nutrients to be output!!!'])
else
disp(['We are adding nutrients to the BC file!!!',Nutrients{1:NNuts}])
tsOBCFile = [basename, '_Nutsobc.nc'];
end
end
subname = 'write_FVCOM_Nutsobc';
global ftbverbose;
if(ftbverbose);
fprintf('\n')
fprintf(['begin : ' subname '\n'])
end;
fvcom_bathy = [basename, '_dep.dat'];
fvcom_obc = [basename, '_obc.dat'];
%------------------------------------------------------------------------------
% read in the FVCOM open boundary node data (need node numbers and dimension)
%------------------------------------------------------------------------------
fid = fopen(fvcom_obc,'r');
if(fid < 0)
error(['file: ' fvcom_obc ' does not exist']);
end;
C = textscan(fid, '%s %s %s %s %d', 1);
nObc = C{5};
obc_nodes = zeros(nObc,1);
if(ftbverbose); fprintf('reading obc file\n'); end;
if(ftbverbose); fprintf('# nodes %d\n',nObc); end;
for i=1:nObc
C = textscan(fid, '%d %d %d', 1);
obc_nodes(i) = C{2};
end;
if(ftbverbose); fprintf('obc reading complete\n');end;
%------------------------------------------------------------------------------
% read in the FVCOM bathymetry data (need bathymetry on open boundary nodes)
%------------------------------------------------------------------------------
fid = fopen(fvcom_bathy,'r');
if(fid < 0)
error(['file: ' fvcom_bathy ' does not exist']);
end;
C = textscan(fid, '%s %s %s %d', 1);
Nverts = C{4};
h = zeros(Nverts,1);
if(ftbverbose); fprintf('reading bathymetry file\n');end;
if(ftbverbose); fprintf('# nodes %d\n',Nverts);end;
for i=1:Nverts
C = textscan(fid, '%f %f %f', 1);
h(i) = C{3};
end;
if(ftbverbose); fprintf('min depth %f max depth %f\n',min(h),max(h));end;
if(ftbverbose); fprintf('bathymetry reading complete\n');end;
fclose(fid);
%--------------------------------------------------------------
% Generate the requisite data
%--------------------------------------------------------------
% extract bathymetry at open boundary nodes
obc_h = h(obc_nodes);
% time
% time = 0:1:31.;
nTimes = numel(time);
nSiglev = nSiglay + 1;
for nuts=1:NNuts
in_temp=Mobj.(Nutrients{nuts});
% Create or process the temperature and salinity arrays.
if max(size(in_temp)) == 1
inc = 1/real(nSiglay);
siglev = 0:-inc:-1;
siglay = nan(1, nSiglay);
for i=1:nSiglay
siglay(i) = mean(siglev(i:i+1));
end
eval([Nutrients{nuts},'= zeros(nObc,nSiglay,nTimes);'])
eval(['obc_',Nutrients{nuts},' = repmat(Mobj.(',Nutrients{nuts},'), 1, nTimes);'])
for i=1:nObc
for j=1:nSiglay
eval([Nutrients{nuts},'(i,j,:)=obc_',Nutrients{nuts},';'])
end
end
else
% We have a 3D array already so we just need a couple of stats.
eval([Nutrients{nuts},'= Mobj.',Nutrients{nuts},';'])
if nargin >= 4 && isfield(Mobj, 'siglay') && isfield(Mobj, 'siglev')
siglev = Mobj.siglev;
siglay = Mobj.siglay;
else
warning('Assuming uniform terrain-following sigma coordinates')
inc = 1/real(nSiglay);
siglev = 0:-inc:-1;
siglay = nan(1, nSiglay);
end
in_test= Mobj.(Nutrients{nuts});
if nSiglay ~= size(in_test, 2) || length(siglay) ~= size(in_test, 2)
error('Specified number of sigma layers does not match supplied data')
end
end
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'),'type','FVCOM RIVER FORCING FILE')
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'title','simple open boundary hydrography test')
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'type','FVCOM TIME SERIES OBC BIO FILE') %% THIS IS IMPORTANT
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'history','File generated using write_FVCOM_tsobcERSEM.m from the MATLAB fvcom-toolbox')
% define dimensions
nobc_dimid=netcdf.defDim(nc,'nobc',nObc);
datestrlen_dimid=netcdf.defDim(nc,'Datestrln',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
% nc{'river_names'} = ncchar('rivers', 'namelen');
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');
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');
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',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',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
varidN{nuts}=['obc_',Nutrients{nuts},'_varid'];
eval([varidN{nuts},'=netcdf.defVar(nc,''',Nutrients{nuts},''',''NC_FLOAT'',[nobc_dimid,siglay_dimid,time_dimid]);'])
eval(['netcdf.putAtt(nc,',varidN{nuts},',''long_name'',''',Nutrients{NNuts+nuts},''');'])
eval(['netcdf.putAtt(nc,',varidN{nuts},',''units'',''',Nutrients{NNuts*2+nuts},''');'])
eval(['netcdf.putAtt(nc,',varidN{nuts},',''grid'',''obc_grid'');'])
% 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');
% netcdf.putAtt(nc,obc_salinity_varid,'grid','obc_grid');
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);
netcdf.putVar(nc,time_varid,0,numel(time),time);
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));
for nuts=1:NNuts
eval(['netcdf.putVar(nc,',varidN{nuts},',',Nutrients{nuts},');'])
end
% close file
netcdf.close(nc);
if(ftbverbose); fprintf(['end : ' subname '\n']);end;
function [polcoms]=calc_scoord(polcoms)
% %n=18;iesub=63;jesub=36;
%dir_path='/users/modellers/rito/research/Projects/MERSEA/dataWCH/WCH2_rito/data2003/';
%filenameb= 'WCH.bathy_WCH2'
%filenameiu= 'WCH.ipexu_WCH2'
%filenameib= 'WCH.ipexb_WCH2'
% program to calculate the scoordinates from model run
%dir_data = '/data/milkyway/rito/lolo/mpi_H_POLCOMS_VIGO_3D_OK_2/setups/data2002_03/'
% Model parameters
%n=18;iesub=50;jesub=70;
% S coordinate parameters
hs = polcoms.bathy;
% Read ipexb
polcoms.ipexbM= zeros(polcoms.iesub,polcoms.jesub);
for aa=1:length(polcoms.isea)
polcoms.ipexbM(polcoms.isea(aa),polcoms.jsea(aa))=1;
end
polcoms.ipexuM= zeros(polcoms.iesub,polcoms.jesub);
for aa=1:length(polcoms.npusea)
polcoms.ipexuM(polcoms.iusea(aa),polcoms.jusea(aa))=1;
end
% Evenly spaced S values
%
dsc = 1.0d0/(polcoms.params.n-2);
sval(1) = -1.0d0;
for k=2:polcoms.params.n-2
sval(k) = sval(k-1)+dsc;
end
sval(polcoms.params.n-1) = 0.0d0;
%JPO START
sval(polcoms.params.n) = 0.0d0;
%JPO END
for j=1:polcoms.jesub
for i=1:polcoms.iesub
sigo(1,i,j) = -1.0d0;
for k=2:polcoms.params.n-2
% c for identical to sigma on shelf
% c
if (hs(i,j)>polcoms.scoord.hc)
ffh=(hs(i,j)-polcoms.scoord.hc)/hs(i,j);
cs = (1.0d0-polcoms.scoord.bb)*(sinh(polcoms.scoord.theta*sval(k)))/sinh(polcoms.scoord.theta)+...
polcoms.scoord.bb*(tanh(polcoms.scoord.theta*(sval(k)+0.5d0))...
-tanh(0.5d0*polcoms.scoord.theta))/...
(2.*tanh(0.5d0*polcoms.scoord.theta));
sigo(k,i,j) = sval(k)+ffh*(cs-sval(k));
else
sigo(k,i,j) = sval(k);
end
end
sigo(polcoms.params.n-1,i,j) = 0.0d0;
% CJPO START
sigo(polcoms.params.n,i,j) = 0.0d0;
% CJPO END
%
end
%
end
%
% Define coordinates on B points
%
ds = 0.0;
for j=1:polcoms.jesub
for i=1:polcoms.iesub
for k=1:polcoms.params.n-2
ds(k,i,j) = sigo(k+1,i,j)-sigo(k,i,j);
end
%c
%c Set surface level to ensure correct sum for ds
%c
sds = 0.0d0;
for k=1:polcoms.params.n-3
sds = sds+ds(k,i,j);
end
ds(polcoms.params.n-2,i,j) = 1.0d0-sds;
%c
dsu(1,i,j) = ds(1,i,j);
for k=2:polcoms.params.n-2
dsu(k,i,j) = 0.5d0*(ds(k,i,j)+ds(k-1,i,j));
end
dsu(polcoms.params.n-1,i,j) = ds(polcoms.params.n-2,i,j);
%c
%sig(0,i,j) = -1.0d0;
sig(1,i,j) = -1.0d0+0.5d0*ds(1,i,j);
for k=2:polcoms.params.n-2
sig(k,i,j) = sig(k-1,i,j)+dsu(k,i,j);
end
sig(polcoms.params.n-1,i,j) = 0.0d0;
%CJPO START
sig(polcoms.params.n,i,j) = 0.0d0;
%CJPO END
end
end
sigov = nan*ones(size(sig));
dsv = nan*ones(size(dsu));
dsuv = nan*ones(size(dsu));
%c
%c Average coordinates onto U points
%c
for j=2:polcoms.jesub
for i=2:polcoms.iesub
if ( polcoms.ipexbM(i,j)~=0 | polcoms.ipexuM(i,j)~=0 )
sigov(1,i,j) = -1.0d0;
for k=2:polcoms.params.n-2
sigov(k,i,j) = 0.25d0*...
(sigo(k,i ,j )+sigo(k,i-1,j )...
+sigo(k,i-1,j-1)+sigo(k,i ,j-1));
end
sigov(polcoms.params.n-1,i,j) = 0.0d0;
%c
end
%c
for k=1:polcoms.params.n-2
dsv(k,i,j) = sigov(k+1,i,j)-sigov(k,i,j);
end
%c
%cjth Set surface level to ensure correct sum for DS
%c
sds=0.;
for k=1:polcoms.params.n-3
sds=sds+dsv(k,i,j);
end
dsv(polcoms.params.n-2,i,j) = 1.0d0-sds;
%c
dsuv(1,i,j) = dsv(1,i,j);
dsuv(polcoms.params.n-1,i,j) = dsv(polcoms.params.n-2,i,j);
for k=2:polcoms.params.n-2
dsuv(k,i,j) = 0.5d0*(dsv(k,i,j)+dsv(k-1,i,j));
end
%c
% sigv(0,i,j) = -1.0d0;
sigv(1,i,j) = -1.0d0+0.5d0*dsv(1,i,j);
for k=2:polcoms.params.n-2
sigv(k,i,j) = sigv(k-1,i,j)+dsuv(k,i,j);
end
sigv(polcoms.params.n-1,i,j) = 0.0d0;
%CJPO START
sigv(polcoms.params.n,i,j) = 0.0d0;
%CJPO END
end
end
sigv(:,1,:)=sigo(:,1,:);sigv(:,:,1)=sigo(:,:,1);
polcoms.sig=sig;
polcoms.sigv=sigv;
polcoms.sigo=sigo;
polcoms.sigov=sigov;
polcoms.ds=ds;
polcoms.dsu=dsu;
polcoms.dsv=dsv;
polcoms.dsuv=dsuv;
return
function Mobj=distance_to_coast(Mobj)
CD=pwd;
% setup HPJ solver to calculate the distance function for the SMS mesh
[~,~,~,bnd] = connectivity([Mobj.x,Mobj.y],Mobj.tri);
pv = [Mobj.x(find(bnd)) Mobj.y(find(bnd));Mobj.x(find(bnd(1))) Mobj.y(find(bnd(1)))];
% remove nodestring from coast.
BCnodes=Mobj.read_obc_nodes{1};
coast_ind=find(bnd);
coast_ind(1:max(BCnodes))=[];
% % calculate distance function
myParam.Max_Tri_per_Star = 20;
myParam.NumGaussSeidel = 40;
myParam.INF_VAL = 1000000;
myParam.TOL = 1e-12;
% in this case, we will assume the standard Euclidean metric
myMetric = [];
myTM.Vtx=[Mobj.x(:),Mobj.y(:)];
myTM.DoFmap=[Mobj.tri];
myTM.NegMask=false(size(bnd));
myBdy.Nodes=coast_ind(:);
myBdy.Data=zeros(size(myBdy.Nodes));
%
cd /users/modellers/rito/matlab/HJB_Solver_Package/HJB_Solver_Package
%
%
SEmex = SolveEikonalmex(myTM,myBdy,myParam,myMetric);
tic
Mobj.dist = SEmex.Compute_Soln;
cd(CD)
end
function Mobj = find_relaxation_boundary(Mobj)
% Find the elements which fall along the boundary.
%
% Mobj = find_boundary_elements(Mobj)
%
% DESCRIPTION:
% Find the elements which are bounded by the open boundaries described by
% the nodes in Mobj.read_obc_nodes.
%
% INPUT:
% Mobj - required fields:
% - read_obc_nodes
% - obc_nodes
% - tri
%
% OUTPUT:
% Mobj - new field of a cell array read_obc_elements which contains the
% IDs of the elements which fall on the model open boundaries and
% nObcElements which is the total number of boundary elements
% along each boundary.