From d8fd1dbd6c354f898efa4e0c7be418d72576f459 Mon Sep 17 00:00:00 2001
From: jpolton <jelt@noc.ac.uk>
Date: Fri, 26 Apr 2019 16:17:02 +0100
Subject: [PATCH] update IC, BC, rivers
---
SCRIPTS/COPERNICUS_INTERP.m | 134 +++++++++++++++++
SCRIPTS/COPERNICUS_INTERP_SSH.m | 100 +++++++++++++
SCRIPTS/Coast_finder.m | 72 +++++++++
SCRIPTS/bdy_obc.m | 100 +++++++++++++
SCRIPTS/nemo_IC_salinity.m | 143 ++++++++++++++----
SCRIPTS/nemo_IC_temperature.m | 143 ++++++++++++++----
SCRIPTS/obc_SSH.m | 228 +++++++++++++++++++++++++++++
SCRIPTS/obc_Sal.m | 233 ++++++++++++++++++++++++++++++
SCRIPTS/obc_Temp.m | 233 ++++++++++++++++++++++++++++++
SCRIPTS/obc_U.m | 249 ++++++++++++++++++++++++++++++++
SCRIPTS/obc_V.m | 249 ++++++++++++++++++++++++++++++++
SCRIPTS/river_maker.m | 125 ++++++++++++++++
12 files changed, 1953 insertions(+), 56 deletions(-)
create mode 100644 SCRIPTS/COPERNICUS_INTERP.m
create mode 100644 SCRIPTS/COPERNICUS_INTERP_SSH.m
create mode 100644 SCRIPTS/Coast_finder.m
create mode 100644 SCRIPTS/bdy_obc.m
create mode 100644 SCRIPTS/obc_SSH.m
create mode 100644 SCRIPTS/obc_Sal.m
create mode 100644 SCRIPTS/obc_Temp.m
create mode 100644 SCRIPTS/obc_U.m
create mode 100644 SCRIPTS/obc_V.m
create mode 100644 SCRIPTS/river_maker.m
diff --git a/SCRIPTS/COPERNICUS_INTERP.m b/SCRIPTS/COPERNICUS_INTERP.m
new file mode 100644
index 0000000..e1f8789
--- /dev/null
+++ b/SCRIPTS/COPERNICUS_INTERP.m
@@ -0,0 +1,134 @@
+function [Variable_obc,z_obc] = COPERNICUS_INTERP(fname_data,z_data,variable_name,nbidta,nbjdta)
+
+
+
+
+%variable_name='so';
+%fname_data = 'Sal_16_0000.nc';
+%z_data = 'domain_cfg.nc';
+
+
+
+
+
+% This is the NEMO grid
+
+xout=double(ncread(z_data,'nav_lon'));
+yout=double(ncread(z_data,'nav_lat'));
+dzout=ncread(z_data,'e3t_0');
+
+% make x and y 3D array - this is because griddata needs it in this form
+xout=repmat(xout,1,1,size(dzout,3));
+yout=repmat(yout,1,1,size(dzout,3));
+
+
+
+% Here we read in the copernicus data
+
+Variable_in=ncread(fname_data,variable_name);
+xin=double(ncread(fname_data,'longitude'));
+yin=double(ncread(fname_data,'latitude'));
+zin=double(ncread(fname_data,'depth'));
+Variable_in=permute(Variable_in,[2 1 3]);
+
+
+
+
+% Build z (e3t_0 just gives dz)
+z=zeros(size(dzout,1),size(dzout,2),size(dzout,3));
+
+for i=1:size(dzout,3)
+
+
+ if i==1
+ z(:,:,i)=dzout(:,:,i);
+ else
+ z(:,:,i)=z(:,:,i-1)+dzout(:,:,i);
+ end
+
+
+
+end
+
+zin(1)=0;
+%zin(end)=max(z(:))+1;
+
+
+
+% make COPERNICUS coordinates 3D arrays (for griddata function)
+[Xin,Yin,Zin]=meshgrid(xin,yin,zin);
+
+
+% We only want to feed non-nan numbers from the copernicus grid.
+Xlr=Xin(~isnan(Variable_in));
+Ylr=Yin(~isnan(Variable_in));
+Zlr=Zin(~isnan(Variable_in));
+Variable_lr=Variable_in(~isnan(Variable_in));
+
+
+% Interpolate Copernicus grid (without nans) on to same grid - basically
+% "flood filling", but by interpolation, not extrapolation.
+
+Variable_in=griddata(Xlr,Ylr,Zlr,Variable_lr,Xin,Yin,Zin);
+
+clear Variable_lr Xlr Ylr Zlr
+
+% nans will still exists around the edges - let's find the nearest
+% neighours to fill them in.
+
+ind_nnan=~isnan(Variable_in(:));
+ind_nan = isnan(Variable_in(:));
+
+Xlnnan=Xin(ind_nnan);
+Ylnnan=Yin(ind_nnan);
+Zlnnan=Zin(ind_nnan);
+coordnnan=[Xlnnan Ylnnan Zlnnan];
+
+Xlnan=Xin(ind_nan);
+Ylnan=Yin(ind_nan);
+Zlnan=Zin(ind_nan);
+coordnan=[Xlnan Ylnan Zlnan];
+
+
+clear Xin Yin Zin Xlnan Ylnan Zlnan Xlnnan Ylnnan Zlnnan
+
+% I found you!
+ind=knnsearch(coordnnan,coordnan);
+
+clear coordnan coordnnan
+
+% tricky to read, but these 3 lines replace nans with nearest neighbours.
+% Trust me.
+
+Variable_TEMP=Variable_in(ind_nnan);
+Variable_TEMP=Variable_TEMP(ind);
+Variable_in(isnan(Variable_in))=Variable_TEMP;
+
+clear Variable_TEMP ind ind_nnan ind_nan
+
+
+
+% let's put our treated copernicus variable in NEMO space.
+
+Variable_interp=interp3(xin,yin,zin,Variable_in,xout,yout,z);
+
+clear Variable_in xin yin zin
+
+for i = 1:length(nbidta)
+
+z_obc(i,1,:) = squeeze(z(nbidta(i),nbjdta(i),:));
+Variable_obc(i,1,:) = squeeze( Variable_interp(nbidta(i),nbjdta(i),:));
+
+end
+
+
+% ncid=netcdf.open(fname,'WRITE');
+% varid = netcdf.inqVarID(ncid,'vosaline');
+%
+% netcdf.putVar(ncid,varid,Variable_interp)
+
+
+
+end
+
+
diff --git a/SCRIPTS/COPERNICUS_INTERP_SSH.m b/SCRIPTS/COPERNICUS_INTERP_SSH.m
new file mode 100644
index 0000000..f01c377
--- /dev/null
+++ b/SCRIPTS/COPERNICUS_INTERP_SSH.m
@@ -0,0 +1,100 @@
+function [Variable_interp,Variable_obc] = COPERNICUS_INTERP_SSH(fname_data,z_data,variable_name,nbidta,nbjdta)
+
+
+
+
+%variable_name='so';
+%fname_data = 'Sal_16_0000.nc';
+%z_data = 'domain_cfg.nc';
+
+
+
+
+
+% This is the NEMO grid
+
+xout=double(ncread(z_data,'nav_lon'));
+yout=double(ncread(z_data,'nav_lat'));
+
+
+% Here we read in the copernicus data
+
+Variable_in=ncread(fname_data,variable_name);
+xin=double(ncread(fname_data,'longitude'));
+yin=double(ncread(fname_data,'latitude'));
+
+
+
+% make COPERNICUS coordinates 3D arrays (for griddata function)
+[Xin,Yin]=meshgrid(xin,yin);
+
+
+% We only want to feed non-nan numbers from the copernicus grid.
+Xlr=Xin(~isnan(Variable_in));
+Ylr=Yin(~isnan(Variable_in));
+Variable_lr=Variable_in(~isnan(Variable_in));
+
+
+% Interpolate Copernicus grid (without nans) on to same grid - basically
+% "flood filling", but by interpolation, not extrapolation.
+
+Variable_in=griddata(Xlr,Ylr,Variable_lr,Xin,Yin);
+
+clear Variable_lr Xlr Ylr Zlr
+
+% nans will still exists around the edges - let's find the nearest
+% neighours to fill them in.
+
+ind_nnan=~isnan(Variable_in(:));
+ind_nan = isnan(Variable_in(:));
+
+Xlnnan=Xin(ind_nnan);
+Ylnnan=Yin(ind_nnan);
+coordnnan=[Xlnnan Ylnnan];
+
+Xlnan=Xin(ind_nan);
+Ylnan=Yin(ind_nan);
+coordnan=[Xlnan Ylnan];
+
+
+clear Xin Yin Zin Xlnan Ylnan Zlnan Xlnnan Ylnnan Zlnnan
+
+% I found you!
+ind=knnsearch(coordnnan,coordnan);
+
+clear coordnan coordnnan
+
+% tricky to read, but these 3 lines replace nans with nearest neighbours.
+% Trust me.
+
+Variable_TEMP=Variable_in(ind_nnan);
+Variable_TEMP=Variable_TEMP(ind);
+Variable_in(isnan(Variable_in))=Variable_TEMP;
+
+clear Variable_TEMP ind ind_nnan ind_nan
+
+
+
+% let's put our treated copernicus variable in NEMO space.
+
+Variable_interp=interp2(xin,yin,Variable_in,xout,yout);
+
+clear Variable_in xin yin zin
+
+for i = 1:length(nbidta)
+
+Variable_obc(i,1) = squeeze( Variable_interp(nbidta(i),nbjdta(i)));
+
+end
+
+
+% ncid=netcdf.open(fname,'WRITE');
+% varid = netcdf.inqVarID(ncid,'vosaline');
+%
+% netcdf.putVar(ncid,varid,Variable_interp)
+
+
+
+end
+
+
diff --git a/SCRIPTS/Coast_finder.m b/SCRIPTS/Coast_finder.m
new file mode 100644
index 0000000..46b0369
--- /dev/null
+++ b/SCRIPTS/Coast_finder.m
@@ -0,0 +1,72 @@
+function [lon_coast,lat_coast,coastal] = Coast_finder(bathy_fname)
+
+%bathy_fname = '../DATA/bathy_meter.nc';
+
+B=ncread(bathy_fname,'Bathymetry');
+lon=ncread(bathy_fname,'nav_lon');
+lat=ncread(bathy_fname,'nav_lat');
+
+% build a land sea mask.
+
+B(isnan(B))=0;
+B(B>0)=1;
+B=-B+1;
+
+% we define a costal point as a wet point that has at least one dry point
+% around it
+
+coastal = zeros(size(B,1),size(B,2));
+
+for i = 1:size(B,1)
+ for j = 1:size(B,2)
+
+
+
+ if B(i,j)==0
+
+
+ adj_points(1:4)=0;
+
+ if i<size(B,1)
+ adj_points(1)=B(i+1,j);
+ end
+
+ if i>1
+ adj_points(2)=B(i-1,j);
+ end
+
+ if j<size(B,2)
+ adj_points(3)=B(i,j+1);
+ end
+
+ if j>1
+ adj_points(4)=B(i,j-1);
+ end
+
+
+
+ if sum(adj_points)>0
+ coastal(i,j)=1;
+ end
+
+
+
+
+ end
+
+
+
+
+
+
+ end
+end
+
+coastal=logical(coastal);
+
+lon_coast=lon(coastal);
+lat_coast=lat(coastal);
+
+
+end
+
\ No newline at end of file
diff --git a/SCRIPTS/bdy_obc.m b/SCRIPTS/bdy_obc.m
new file mode 100644
index 0000000..d2354ef
--- /dev/null
+++ b/SCRIPTS/bdy_obc.m
@@ -0,0 +1,100 @@
+clear all
+clc
+
+
+fname = 'initcd_vosaline.nc';
+fname_data = 'Sal_16_0000.nc';
+z_data = 'domain_cfg.nc';
+
+Sin=ncread(fname_data,'so');
+xin=double(ncread(fname_data,'longitude'));
+yin=double(ncread(fname_data,'latitude'));
+zin=double(ncread(fname_data,'depth'));
+Sin=permute(Sin,[2 1 3]);
+
+
+xout=double(ncread(z_data,'nav_lon'));
+yout=double(ncread(z_data,'nav_lat'));
+dzout=ncread(z_data,'e3t_0');
+
+xout=repmat(xout,1,1,size(dzout,3));
+yout=repmat(yout,1,1,size(dzout,3));
+
+
+
+
+z=zeros(size(dzout,1),size(dzout,2),size(dzout,3));
+
+for i=1:size(dzout,3)
+
+
+ if i==1
+ z(:,:,i)=dzout(:,:,i);
+ else
+ z(:,:,i)=z(:,:,i-1)+dzout(:,:,i);
+ end
+
+
+
+end
+
+zin(1)=0;
+%zin(end)=max(z(:))+1;
+
+
+
+
+[Xin,Yin,Zin]=meshgrid(xin,yin,zin);
+
+
+
+Xlr=Xin(~isnan(Sin));
+Ylr=Yin(~isnan(Sin));
+Zlr=Zin(~isnan(Sin));
+Slr=Sin(~isnan(Sin));
+
+
+
+Sin=griddata(Xlr,Ylr,Zlr,Slr,Xin,Yin,Zin);
+
+
+ind_nnan=~isnan(Sin(:));
+ind_nan = isnan(Sin(:));
+
+Xlnnan=Xin(ind_nnan);
+Ylnnan=Yin(ind_nnan);
+Zlnnan=Zin(ind_nnan);
+coordnnan=[Xlnnan Ylnnan Zlnnan];
+
+Xlnan=Xin(ind_nan);
+Ylnan=Yin(ind_nan);
+Zlnan=Zin(ind_nan);
+coordnan=[Xlnan Ylnan Zlnan];
+
+ind=knnsearch(coordnnan,coordnan);
+
+S_TEMP=Sin(ind_nnan);
+S_TEMP=S_TEMP(ind);
+
+
+Sin(isnan(Sin))=S_TEMP;
+
+
+
+
+
+
+
+
+Sinterp=interp3(xin,yin,zin,Sin,xout,yout,z);
+
+
+
+
+
+ncid=netcdf.open(fname,'WRITE');
+varid = netcdf.inqVarID(ncid,'vosaline');
+
+netcdf.putVar(ncid,varid,Sinterp)
+
+
diff --git a/SCRIPTS/nemo_IC_salinity.m b/SCRIPTS/nemo_IC_salinity.m
index de693b5..41756f4 100644
--- a/SCRIPTS/nemo_IC_salinity.m
+++ b/SCRIPTS/nemo_IC_salinity.m
@@ -1,66 +1,151 @@
clear all
clc
-addpath('/projectsa/ETI/Tracer_work/Matlab_progs')
-fname = 'initcd_vosaline.nc';
-fname_data = 'Sal_16_0000.nc';
-z_data = 'domain_cfg.nc';
+% This programme takes Temperature from an input grid, treats it so all
+% nans are replaced by interpolated or nearest neighbour values and then
+% interpolates this on to the NEMO grid
-Sin=ncread(fname_data,'so');
+% Note, interpolation around land is a dark art - this may not give a
+% stable solution, but it at least works for me.
+
+% 24/04/19 - Ashley Brereton
+
+
+fname = 'initcd_vosaline.nc'; % name of output file
+fname_data = 'Sal_16_0000.nc'; % name of input file
+z_data = 'domain_cfg.nc'; % grid file
+
+% The input data is a regular grid. The co-ordinates are 1D vectors and the
+% variable is F(x,y,z,t). Later in this code, we make our co-ordinates each
+% a function of F(x,y,z) with meshgrid. If your data is in NEMO format,
+% replicate how (xout,yout,zout) are made for (xin,yin,zin).
+
+% Read in input variables
+
+Sin=ncread(fname_data,'so');
xin=double(ncread(fname_data,'longitude'));
yin=double(ncread(fname_data,'latitude'));
zin=double(ncread(fname_data,'depth'));
Sin=permute(Sin,[2 1 3]);
-Sout=double(ncread(fname,'vosaline'));
-xout=double(ncread(fname,'x'));
-yout=double(ncread(fname,'y'));
+
+% Read output grid information
+
+xout=double(ncread(z_data,'nav_lon'));
+yout=double(ncread(z_data,'nav_lat'));
dzout=ncread(z_data,'e3t_0');
-xout=repmat(xout,1,1,size(Sout,3));
-yout=repmat(yout,1,1,size(Sout,3));
+% need all grid variables as F(x,y,z) as this is what is needed by griddata
+xout=repmat(xout,1,1,size(dzout,3));
+yout=repmat(yout,1,1,size(dzout,3));
+% pre-allocate matrix
+zout=zeros(size(dzout,1),size(dzout,2),size(dzout,3));
-z=zeros(size(Sout,1),size(Sout,2),size(Sout,3));
-for i=1:size(Sout,3)
-
+
+% Build z from dz given in grid file
+
+for i=1:size(dzout,3)
if i==1
- z(:,:,i)=dzout(:,:,i);
+ zout(:,:,i)=dzout(:,:,i);
else
- z(:,:,i)=z(:,:,i-1)+dzout(:,:,i);
+ zout(:,:,i)=zout(:,:,i-1)+dzout(:,:,i);
end
-
-
-
+
end
+
+% fix surface depth to 0
+
zin(1)=0;
-zin(end)=max(z(:))+1;
-%
-% parfor i=1:size(Sin,3)
-% i
-% Sin(:,:,i)=smooth2(Sin(:,:,i),100,100);
-%
-% end
+% make all input grid information F(x,y,z)
+[Xin,Yin,Zin]=meshgrid(xin,yin,zin);
+
+
+% Remove all nan points.
+Xlr=Xin(~isnan(Sin));
+Ylr=Yin(~isnan(Sin));
+Zlr=Zin(~isnan(Sin));
+Tlr=Sin(~isnan(Sin));
+
+
+% Interpolate input on to whole input grid - getting rid of nans
+Sin=griddata(Xlr,Ylr,Zlr,Tlr,Xin,Yin,Zin);
+
+
+% All the nans that couldn't be filled on the input grid, we replace with
+% nearest neighbours.
+
+ind_nnan=~isnan(Sin(:));
+ind_nan = isnan(Sin(:));
+
+Xlnnan=Xin(ind_nnan);
+Ylnnan=Yin(ind_nnan);
+Zlnnan=Zin(ind_nnan);
+coordnnan=[Xlnnan Ylnnan Zlnnan];
+
+Xlnan=Xin(ind_nan);
+Ylnan=Yin(ind_nan);
+Zlnan=Zin(ind_nan);
+coordnan=[Xlnan Ylnan Zlnan];
+
+ind=knnsearch(coordnnan,coordnan);
+
+T_TEMP=Sin(ind_nnan);
+T_TEMP=T_TEMP(ind);
+
+
+Sin(isnan(Sin))=T_TEMP;
+
+
+
+% Then we interpolate fully treated input grid on to our domain.
+
+Sinterp=interp3(xin,yin,zin,Sin,xout,yout,zout);
+
+
+% converts days since 1900 to seconds since 1950...
+time=ncread( fname_data , 'time');
+time_counter = (time + 18262*24)*3600; % 18262 days between 1900 and 1950
+
+
+% write all the variables to file
+delete(fname)
+nccreate(fname,'nav_lon',...
+ 'Dimensions',{'x',size(xout,1),'y',size(yout,2)},...
+ 'Format','netcdf4_classic');
+ncwrite(fname,'nav_lon',nav_lon);
-Sinterp=interp3(xin,yin,zin,Sin,xout,yout,z);
+nccreate(fname,'nav_lat',...
+ 'Dimensions',{'x',size(xout,1),'y',size(yout,2)},...
+ 'Format','netcdf4_classic');
+ncwrite(fname,'nav_lat',nav_lat);
+nccreate(fname,'time_counter',...
+ 'Dimensions',{'time_counter',Inf},...
+ 'Format','netcdf4_classic');
+ncwrite(fname,'time_counter',time_counter);
+nccreate(fname,'gdept',...
+ 'Dimensions',{'x',size(zout,1),'y',size(zout,2),'z',size(zout,3)},...
+ 'Format','netcdf4_classic');
+ncwrite(fname,'gdept',z);
-ncid=netcdf.open(fname,'WRITE');
-varid = netcdf.inqVarID(ncid,'vosaline');
-netcdf.putVar(ncid,varid,Sinterp)
+nccreate(fname,'vosaline',...
+ 'Dimensions',{'x',size(zout,1),'y',size(zout,2),'z',size(zout,3),'time_counter',1},...
+ 'Format','netcdf4_classic');
+ncwrite(fname,'votemper',Sinterp);
diff --git a/SCRIPTS/nemo_IC_temperature.m b/SCRIPTS/nemo_IC_temperature.m
index 413ad8e..73ebe7d 100644
--- a/SCRIPTS/nemo_IC_temperature.m
+++ b/SCRIPTS/nemo_IC_temperature.m
@@ -1,64 +1,153 @@
clear all
clc
-fname = 'initcd_votemper.nc';
-fname_data = 'Temp_16_0000.nc';
-z_data = 'domain_cfg.nc';
+% This programme takes Temperature from an input grid, treats it so all
+% nans are replaced by interpolated or nearest neighbour values and then
+% interpolates this on to the NEMO grid
-Tin=ncread(fname_data,'thetao');
+% Note, interpolation around land is a dark art - this may not give a
+% stable solution, but it at least works for me.
+
+% 24/04/19 - Ashley Brereton
+
+
+fname = 'initcd_votemper.nc'; % name of output file
+fname_data = 'Temp_16_0000.nc'; % name of input file
+z_data = 'domain_cfg.nc'; % grid file
+
+
+% The input data is a regular grid. The co-ordinates are 1D vectors and the
+% variable is F(x,y,z,t). Later in this code, we make our co-ordinates each
+% a function of F(x,y,z) with meshgrid. If your data is in NEMO format,
+% replicate how (xout,yout,zout) are made for (xin,yin,zin).
+
+
+% Read in input variables
+
+Tin=ncread(fname_data,'thetao');
xin=double(ncread(fname_data,'longitude'));
yin=double(ncread(fname_data,'latitude'));
zin=double(ncread(fname_data,'depth'));
Tin=permute(Tin,[2 1 3]);
-Tout=double(ncread(fname,'votemper'));
-xout=double(ncread(fname,'x'));
-yout=double(ncread(fname,'y'));
+
+% Read output grid information
+
+xout=double(ncread(z_data,'nav_lon'));
+yout=double(ncread(z_data,'nav_lat'));
dzout=ncread(z_data,'e3t_0');
-xout=repmat(xout,1,1,size(Tout,3));
-yout=repmat(yout,1,1,size(Tout,3));
+% need all grid variables as F(x,y,z) as this is what is needed by griddata
+xout=repmat(xout,1,1,size(dzout,3));
+yout=repmat(yout,1,1,size(dzout,3));
+% pre-allocate matrix
+zout=zeros(size(dzout,1),size(dzout,2),size(dzout,3));
-z=zeros(size(Tout,1),size(Tout,2),size(Tout,3));
-for i=1:size(Tout,3)
-
+
+% Build z from dz given in grid file
+
+for i=1:size(dzout,3)
if i==1
- z(:,:,i)=dzout(:,:,i);
+ zout(:,:,i)=dzout(:,:,i);
else
- z(:,:,i)=z(:,:,i-1)+dzout(:,:,i);
+ zout(:,:,i)=zout(:,:,i-1)+dzout(:,:,i);
end
-
-
-
+
end
+
+% fix surface depth to 0
+
zin(1)=0;
-zin(end)=max(z(:))+1;
+
+
+% make all input grid information F(x,y,z)
+[Xin,Yin,Zin]=meshgrid(xin,yin,zin);
+
+
+% Remove all nan points.
+Xlr=Xin(~isnan(Tin));
+Ylr=Yin(~isnan(Tin));
+Zlr=Zin(~isnan(Tin));
+Tlr=Tin(~isnan(Tin));
+
+
+% Interpolate input on to whole input grid - getting rid of nans
+Tin=griddata(Xlr,Ylr,Zlr,Tlr,Xin,Yin,Zin);
+
+
+% All the nans that couldn't be filled on the input grid, we replace with
+% nearest neighbours.
+
+ind_nnan=~isnan(Tin(:));
+ind_nan = isnan(Tin(:));
+
+Xlnnan=Xin(ind_nnan);
+Ylnnan=Yin(ind_nnan);
+Zlnnan=Zin(ind_nnan);
+coordnnan=[Xlnnan Ylnnan Zlnnan];
+
+Xlnan=Xin(ind_nan);
+Ylnan=Yin(ind_nan);
+Zlnan=Zin(ind_nan);
+coordnan=[Xlnan Ylnan Zlnan];
+
+ind=knnsearch(coordnnan,coordnan);
+
+T_TEMP=Tin(ind_nnan);
+T_TEMP=T_TEMP(ind);
+
+
+Tin(isnan(Tin))=T_TEMP;
+
+
+
+% Then we interpolate fully treated input grid on to our domain.
+
+Tinterp=interp3(xin,yin,zin,Tin,xout,yout,zout);
+
+
+% converts days since 1900 to seconds since 1950...
+time=ncread( fname_data , 'time');
+time_counter = (time + 18262*24)*3600; % 18262 days between 1900 and 1950
-% parfor i=1:size(Tin,3)
-% i
-% Tin(:,:,i)=smooth2(Tin(:,:,i),100,100);
-%
-% end
+% write all the variables to file
+delete(fname)
+nccreate(fname,'nav_lon',...
+ 'Dimensions',{'x',size(xout,1),'y',size(yout,2)},...
+ 'Format','netcdf4_classic');
+ncwrite(fname,'nav_lon',nav_lon);
-Tinterp=interp3(xin,yin,zin,Tin,xout,yout,z);
+nccreate(fname,'nav_lat',...
+ 'Dimensions',{'x',size(xout,1),'y',size(yout,2)},...
+ 'Format','netcdf4_classic');
+ncwrite(fname,'nav_lat',nav_lat);
+nccreate(fname,'time_counter',...
+ 'Dimensions',{'time_counter',Inf},...
+ 'Format','netcdf4_classic');
+ncwrite(fname,'time_counter',time_counter);
+nccreate(fname,'gdept',...
+ 'Dimensions',{'x',size(zout,1),'y',size(zout,2),'z',size(zout,3)},...
+ 'Format','netcdf4_classic');
+ncwrite(fname,'gdept',z);
-ncid=netcdf.open(fname,'WRITE');
-varid = netcdf.inqVarID(ncid,'votemper');
-netcdf.putVar(ncid,varid,Tinterp)
+nccreate(fname,'votemper',...
+ 'Dimensions',{'x',size(zout,1),'y',size(zout,2),'z',size(zout,3),'time_counter',1},...
+ 'Format','netcdf4_classic');
+ncwrite(fname,'votemper',Tinterp);
diff --git a/SCRIPTS/obc_SSH.m b/SCRIPTS/obc_SSH.m
new file mode 100644
index 0000000..21ef517
--- /dev/null
+++ b/SCRIPTS/obc_SSH.m
@@ -0,0 +1,228 @@
+clear all
+clc
+
+% This programme reads in Salinity data, interpolates on to the NEMO grid
+% that you provide and outputs the points on the open boundaries, which you
+% will also provide.
+
+% Ingedients needed:
+
+% 1) domain_cfg.nc % Your NEMO grid
+% 2) bathy_meter.nc % Your bathymetry file
+% 3) coordinates.bdy.nc % Your coordinates file, generated when tides
+% you make tides
+
+% 4) Daily parent data. This has struction ./SAL16_dir/JAN/Sal_16_0000.nc
+
+% This code uses the parallel toolbox. This might have memory restrictions
+% based on the computer you are using. If you want this to run quickly, use
+% the following syntax: parpool(N) with N = 16 or N=31 if you have the
+% resources.
+
+% Note this is adapted only for Copernicus data but can easily adapted for
+% other datasets. Email a.brereton@liverpool.ac.uk for help.
+
+
+
+% we want to make the current directory visible so we can use subroutines
+addpath(pwd)
+
+% Control variables - make sure
+
+prefix='INDIAN_bdy_'; % prefix to output filename
+output_name='sossheig'; % output variable name
+variable_name='zos'; % variable from input file
+filename_prefix='SSH'; % e.g. 'INDIAN_bdy_TEMP...'
+Input_directory='SSH17_dir'; % Input file directory
+year = '2017'; % e.g. 'INDIAN_bdy_TEMP_y2017...'
+grid_data = 'domain_cfg.nc'; % NEMO grid information
+
+
+% e.g. 'INDIAN_bdy_SAL_y2017m01.nc'
+
+Input_month_suf{1}='01.nc';
+Input_month_suf{2}='02.nc';
+Input_month_suf{3}='03.nc';
+Input_month_suf{4}='04.nc';
+Input_month_suf{5}='05.nc';
+Input_month_suf{6}='06.nc';
+Input_month_suf{7}='07.nc';
+Input_month_suf{8}='08.nc';
+Input_month_suf{9}='09.nc';
+Input_month_suf{10}='10.nc';
+Input_month_suf{11}='11.nc';
+Input_month_suf{12}='12.nc';
+
+
+% Input directory month names e.g. ...../SAL16_dir/JAN/
+
+Input_month{1}='JAN';
+Input_month{2}='FEB';
+Input_month{3}='MAR';
+Input_month{4}='APR';
+Input_month{5}='MAY';
+Input_month{6}='JUN';
+Input_month{7}='JUL';
+Input_month{8}='AUG';
+Input_month{9}='SEP';
+Input_month{10}='OCT';
+Input_month{11}='NOV';
+Input_month{12}='DEC';
+
+
+
+% x co-ordinate indice of boundary
+
+nbidta = ncread('coordinates.bdy.nc','nbit'); % (yb, xb)
+
+
+
+
+
+% y co-ordinate indice of boundary
+
+nbjdta = ncread('coordinates.bdy.nc','nbjt'); % (yb,xb)
+
+
+
+
+
+% rim level, this code is only set up for a rimwidth of 1.
+
+nbrdta = nbidta*0+1; % (yb,xb)
+
+
+
+% nav_lon - we have this in domain_cfg (y,x)
+% nav_lat - we have this in domain_cfg (y,x)
+
+nav_lon = ncread(grid_data,'nav_lon');
+nav_lat = ncread(grid_data,'nav_lat');
+
+
+
+
+% create a land-sea mask from bathymetry
+
+B=ncread('bathy_meter.nc','Bathymetry');
+
+B(B>0)=1; % 1 is Water
+B(isnan(B))=0; % 0 is Land
+
+bdy_msk=B;
+
+
+
+
+
+
+for month=1:12
+
+
+
+% e.g. 'INDIAN_bdy_SAL_y2017m01.nc'
+
+output_filename = [prefix , filename_prefix , '_y',year,'m',Input_month_suf{month}];
+
+
+
+% cd into the data directory
+cd(Input_directory)
+cd(Input_month{month})
+filenames=dir('*.nc');
+
+
+
+% This loops through daily files, interpolating onto NEMO grid and
+% getting boundary points for the data.
+
+
+
+
+parfor j=1:length(filenames)
+
+
+ [~,Variable_obc] = COPERNICUS_INTERP_SSH(filenames(j).name,grid_data,variable_name,nbidta,nbjdta);
+
+
+
+
+ NEMO_variable_obc(:,:,j)=Variable_obc;
+
+
+
+
+ % deptht - this is the depths at the boundary points (z, yb, xb)
+ % votemper (etc) - (time_counter, z, yb, xb)
+
+ time(j)=ncread( filenames(j).name , 'time');
+ time_counter(j) = (time(j) + 18262*24)*3600; % 18262 days between 1900 and 1950
+
+end
+
+
+
+
+
+
+cd ../..
+
+
+
+
+% Write out all the relevant variables to monthly file
+
+
+delete(output_filename)
+
+
+nccreate(output_filename,'nbrdta',...
+ 'Dimensions',{'xb',length(nbidta),'yb',1},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nbrdta',nbrdta);
+
+nccreate(output_filename,'nbidta',...
+ 'Dimensions',{'xb',length(nbidta),'yb',1},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nbidta',nbidta);
+
+nccreate(output_filename,'nbjdta',...
+ 'Dimensions',{'xb',length(nbidta),'yb',1},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nbjdta',nbjdta);
+
+nccreate(output_filename,'nav_lon',...
+ 'Dimensions',{'x',size(nav_lon,1),'y',size(nav_lon,2)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nav_lon',nav_lon);
+
+nccreate(output_filename,'nav_lat',...
+ 'Dimensions',{'x',size(nav_lat,1),'y',size(nav_lat,2)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nav_lat',nav_lat);
+
+nccreate(output_filename,'bdy_msk',...
+ 'Dimensions',{'x',size(bdy_msk,1),'y',size(bdy_msk,2)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'bdy_msk',bdy_msk);
+
+nccreate(output_filename,'time_counter',...
+ 'Dimensions',{'time_counter',Inf},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'time_counter',time_counter);
+
+nccreate(output_filename,output_name,...
+ 'Dimensions',{'xb',size(NEMO_variable_obc,1),'yb',1,'time_counter',size(NEMO_variable_obc,3)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,output_name,NEMO_variable_obc);
+
+
+end
\ No newline at end of file
diff --git a/SCRIPTS/obc_Sal.m b/SCRIPTS/obc_Sal.m
new file mode 100644
index 0000000..1308853
--- /dev/null
+++ b/SCRIPTS/obc_Sal.m
@@ -0,0 +1,233 @@
+clear all
+clc
+
+% This programme reads in Salinity data, interpolates on to the NEMO grid
+% that you provide and outputs the points on the open boundaries, which you
+% will also provide.
+
+% Ingedients needed:
+
+% 1) domain_cfg.nc % Your NEMO grid
+% 2) bathy_meter.nc % Your bathymetry file
+% 3) coordinates.bdy.nc % Your coordinates file, generated when tides
+% you make tides
+
+% 4) Daily parent data. This has struction ./SAL16_dir/JAN/Sal_16_0000.nc
+
+% This code uses the parallel toolbox. This might have memory restrictions
+% based on the computer you are using. If you want this to run quickly, use
+% the following syntax: parpool(N) with N = 16 or N=31 if you have the
+% resources.
+
+% Note this is adapted only for Copernicus data but can easily adapted for
+% other datasets. Email a.brereton@liverpool.ac.uk for help.
+
+
+
+% we want to make the current directory visible so we can use subroutines
+addpath(pwd)
+
+% Control variables - make sure
+
+prefix='INDIAN_bdy_'; % prefix to output filename
+output_name='vosaline'; % output variable name
+variable_name='so'; % variable from input file
+filename_prefix='SAL'; % e.g. 'INDIAN_bdy_SAL...'
+Input_directory='SAL17_dir'; % Input file directory
+year = '2017'; % e.g. 'INDIAN_bdy_SAL_y2017...'
+grid_data = 'domain_cfg.nc'; % NEMO grid information
+
+
+% e.g. 'INDIAN_bdy_SAL_y2017m01.nc'
+
+Input_month_suf{1}='01.nc';
+Input_month_suf{2}='02.nc';
+Input_month_suf{3}='03.nc';
+Input_month_suf{4}='04.nc';
+Input_month_suf{5}='05.nc';
+Input_month_suf{6}='06.nc';
+Input_month_suf{7}='07.nc';
+Input_month_suf{8}='08.nc';
+Input_month_suf{9}='09.nc';
+Input_month_suf{10}='10.nc';
+Input_month_suf{11}='11.nc';
+Input_month_suf{12}='12.nc';
+
+
+% Input directory month names e.g. ...../SAL16_dir/JAN/
+
+Input_month{1}='JAN';
+Input_month{2}='FEB';
+Input_month{3}='MAR';
+Input_month{4}='APR';
+Input_month{5}='MAY';
+Input_month{6}='JUN';
+Input_month{7}='JUL';
+Input_month{8}='AUG';
+Input_month{9}='SEP';
+Input_month{10}='OCT';
+Input_month{11}='NOV';
+Input_month{12}='DEC';
+
+
+
+% x co-ordinate indice of boundary
+
+nbidta = ncread('coordinates.bdy.nc','nbit'); % (yb, xb)
+
+
+
+
+
+% y co-ordinate indice of boundary
+
+nbjdta = ncread('coordinates.bdy.nc','nbjt'); % (yb,xb)
+
+
+
+
+
+% rim level, this code is only set up for a rimwidth of 1.
+
+nbrdta = nbidta*0+1; % (yb,xb)
+
+
+
+% nav_lon - we have this in domain_cfg (y,x)
+% nav_lat - we have this in domain_cfg (y,x)
+
+nav_lon = ncread(grid_data,'nav_lon');
+nav_lat = ncread(grid_data,'nav_lat');
+
+
+
+
+% create a land-sea mask from bathymetry
+
+B=ncread('bathy_meter.nc','Bathymetry');
+
+B(B>0)=1; % 1 is Water
+B(isnan(B))=0; % 0 is Land
+
+bdy_msk=B;
+
+
+
+
+
+
+for month=1:12
+
+
+
+% e.g. 'INDIAN_bdy_SAL_y2017m01.nc'
+
+output_filename = [prefix , filename_prefix , '_y',year,'m',Input_month_suf{month}];
+
+
+
+% cd into the data directory
+cd(Input_directory)
+cd(Input_month{month})
+filenames=dir('*.nc');
+
+
+
+% This loops through daily files, interpolating onto NEMO grid and
+% getting boundary points for the data.
+
+parfor j=1:length(filenames)
+
+ if j==1
+
+ [Variable_obc,z_obc] = COPERNICUS_INTERP(filenames(j).name,grid_data,variable_name,nbidta,nbjdta);
+
+ deptht(:,:,:,j) = z_obc;
+
+
+ else
+
+ [Variable_obc,~] = COPERNICUS_INTERP(filenames(j).name,grid_data,variable_name,nbidta,nbjdta);
+
+ end
+
+
+ NEMO_variable_obc(:,:,:,j)=Variable_obc;
+
+
+ % Get time information. NEMO wants seconds since 1950. Copernicus data
+ % comes in hours since 1900.
+
+ time(j)=ncread( filenames(j).name , 'time');
+ time_counter(j) = (time(j) + 18262*24)*3600; % 18262 days between 1900 and 1950
+
+end
+
+
+
+
+cd ../..
+
+
+
+% Write out all the relevant variables to monthly file
+
+
+delete(output_filename)
+
+nccreate(output_filename,'nbrdta',...
+ 'Dimensions',{'xb',length(nbidta),'yb',1},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nbrdta',nbrdta);
+
+nccreate(output_filename,'nbidta',...
+ 'Dimensions',{'xb',length(nbidta),'yb',1},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nbidta',nbidta);
+
+nccreate(output_filename,'nbjdta',...
+ 'Dimensions',{'xb',length(nbidta),'yb',1},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nbjdta',nbjdta);
+
+nccreate(output_filename,'nav_lon',...
+ 'Dimensions',{'x',size(nav_lon,1),'y',size(nav_lon,2)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nav_lon',nav_lon);
+
+nccreate(output_filename,'nav_lat',...
+ 'Dimensions',{'x',size(nav_lat,1),'y',size(nav_lat,2)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nav_lat',nav_lat);
+
+nccreate(output_filename,'bdy_msk',...
+ 'Dimensions',{'x',size(bdy_msk,1),'y',size(bdy_msk,2)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'bdy_msk',bdy_msk);
+
+nccreate(output_filename,'time_counter',...
+ 'Dimensions',{'time_counter',Inf},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'time_counter',time_counter);
+
+nccreate(output_filename,'deptht',...
+ 'Dimensions',{'xb',size(deptht,1),'yb',1,'z',size(deptht,3)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'deptht',deptht);
+
+nccreate(output_filename,output_name,...
+ 'Dimensions',{'xb',size(NEMO_variable_obc,1),'yb',1,'z',size(NEMO_variable_obc,3),'time_counter',size(NEMO_variable_obc,4)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,output_name,NEMO_variable_obc);
+
+
+end
diff --git a/SCRIPTS/obc_Temp.m b/SCRIPTS/obc_Temp.m
new file mode 100644
index 0000000..78eee55
--- /dev/null
+++ b/SCRIPTS/obc_Temp.m
@@ -0,0 +1,233 @@
+clear all
+clc
+
+% This programme reads in Salinity data, interpolates on to the NEMO grid
+% that you provide and outputs the points on the open boundaries, which you
+% will also provide.
+
+% Ingedients needed:
+
+% 1) domain_cfg.nc % Your NEMO grid
+% 2) bathy_meter.nc % Your bathymetry file
+% 3) coordinates.bdy.nc % Your coordinates file, generated when tides
+% you make tides
+
+% 4) Daily parent data. This has struction ./SAL16_dir/JAN/Sal_16_0000.nc
+
+% This code uses the parallel toolbox. This might have memory restrictions
+% based on the computer you are using. If you want this to run quickly, use
+% the following syntax: parpool(N) with N = 16 or N=31 if you have the
+% resources.
+
+% Note this is adapted only for Copernicus data but can easily adapted for
+% other datasets. Email a.brereton@liverpool.ac.uk for help.
+
+
+
+% we want to make the current directory visible so we can use subroutines
+addpath(pwd)
+
+% Control variables - make sure
+
+prefix='INDIAN_bdy_'; % prefix to output filename
+output_name='votemper'; % output variable name
+variable_name='thetao'; % variable from input file
+filename_prefix='TEMP'; % e.g. 'INDIAN_bdy_TEMP...'
+Input_directory='TEMP17_dir'; % Input file directory
+year = '2017'; % e.g. 'INDIAN_bdy_TEMP_y2017...'
+grid_data = 'domain_cfg.nc'; % NEMO grid information
+
+
+% e.g. 'INDIAN_bdy_SAL_y2017m01.nc'
+
+Input_month_suf{1}='01.nc';
+Input_month_suf{2}='02.nc';
+Input_month_suf{3}='03.nc';
+Input_month_suf{4}='04.nc';
+Input_month_suf{5}='05.nc';
+Input_month_suf{6}='06.nc';
+Input_month_suf{7}='07.nc';
+Input_month_suf{8}='08.nc';
+Input_month_suf{9}='09.nc';
+Input_month_suf{10}='10.nc';
+Input_month_suf{11}='11.nc';
+Input_month_suf{12}='12.nc';
+
+
+% Input directory month names e.g. ...../SAL16_dir/JAN/
+
+Input_month{1}='JAN';
+Input_month{2}='FEB';
+Input_month{3}='MAR';
+Input_month{4}='APR';
+Input_month{5}='MAY';
+Input_month{6}='JUN';
+Input_month{7}='JUL';
+Input_month{8}='AUG';
+Input_month{9}='SEP';
+Input_month{10}='OCT';
+Input_month{11}='NOV';
+Input_month{12}='DEC';
+
+
+
+% x co-ordinate indice of boundary
+
+nbidta = ncread('coordinates.bdy.nc','nbit'); % (yb, xb)
+
+
+
+
+
+% y co-ordinate indice of boundary
+
+nbjdta = ncread('coordinates.bdy.nc','nbjt'); % (yb,xb)
+
+
+
+
+
+% rim level, this code is only set up for a rimwidth of 1.
+
+nbrdta = nbidta*0+1; % (yb,xb)
+
+
+
+% nav_lon - we have this in domain_cfg (y,x)
+% nav_lat - we have this in domain_cfg (y,x)
+
+nav_lon = ncread(grid_data,'nav_lon');
+nav_lat = ncread(grid_data,'nav_lat');
+
+
+
+
+% create a land-sea mask from bathymetry
+
+B=ncread('bathy_meter.nc','Bathymetry');
+
+B(B>0)=1; % 1 is Water
+B(isnan(B))=0; % 0 is Land
+
+bdy_msk=B;
+
+
+
+
+
+
+for month=1:12
+
+
+
+% e.g. 'INDIAN_bdy_SAL_y2017m01.nc'
+
+output_filename = [prefix , filename_prefix , '_y',year,'m',Input_month_suf{month}];
+
+
+
+% cd into the data directory
+cd(Input_directory)
+cd(Input_month{month})
+filenames=dir('*.nc');
+
+
+
+% This loops through daily files, interpolating onto NEMO grid and
+% getting boundary points for the data.
+
+parfor j=1:length(filenames)
+
+ if j==1
+
+ [Variable_obc,z_obc] = COPERNICUS_INTERP(filenames(j).name,grid_data,variable_name,nbidta,nbjdta);
+
+ deptht(:,:,:,j) = z_obc;
+
+
+ else
+
+ [Variable_obc,~] = COPERNICUS_INTERP(filenames(j).name,grid_data,variable_name,nbidta,nbjdta);
+
+ end
+
+
+ NEMO_variable_obc(:,:,:,j)=Variable_obc;
+
+
+ % Get time information. NEMO wants seconds since 1950. Copernicus data
+ % comes in hours since 1900.
+
+ time(j)=ncread( filenames(j).name , 'time');
+ time_counter(j) = (time(j) + 18262*24)*3600; % 18262 days between 1900 and 1950
+
+end
+
+
+
+
+cd ../..
+
+
+
+% Write out all the relevant variables to monthly file
+
+
+delete(output_filename)
+
+nccreate(output_filename,'nbrdta',...
+ 'Dimensions',{'xb',length(nbidta),'yb',1},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nbrdta',nbrdta);
+
+nccreate(output_filename,'nbidta',...
+ 'Dimensions',{'xb',length(nbidta),'yb',1},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nbidta',nbidta);
+
+nccreate(output_filename,'nbjdta',...
+ 'Dimensions',{'xb',length(nbidta),'yb',1},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nbjdta',nbjdta);
+
+nccreate(output_filename,'nav_lon',...
+ 'Dimensions',{'x',size(nav_lon,1),'y',size(nav_lon,2)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nav_lon',nav_lon);
+
+nccreate(output_filename,'nav_lat',...
+ 'Dimensions',{'x',size(nav_lat,1),'y',size(nav_lat,2)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nav_lat',nav_lat);
+
+nccreate(output_filename,'bdy_msk',...
+ 'Dimensions',{'x',size(bdy_msk,1),'y',size(bdy_msk,2)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'bdy_msk',bdy_msk);
+
+nccreate(output_filename,'time_counter',...
+ 'Dimensions',{'time_counter',Inf},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'time_counter',time_counter);
+
+nccreate(output_filename,'deptht',...
+ 'Dimensions',{'xb',size(deptht,1),'yb',1,'z',size(deptht,3)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'deptht',deptht);
+
+nccreate(output_filename,output_name,...
+ 'Dimensions',{'xb',size(NEMO_variable_obc,1),'yb',1,'z',size(NEMO_variable_obc,3),'time_counter',size(NEMO_variable_obc,4)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,output_name,NEMO_variable_obc);
+
+
+end
diff --git a/SCRIPTS/obc_U.m b/SCRIPTS/obc_U.m
new file mode 100644
index 0000000..9df650c
--- /dev/null
+++ b/SCRIPTS/obc_U.m
@@ -0,0 +1,249 @@
+clear all
+clc
+
+% This programme reads in Salinity data, interpolates on to the NEMO grid
+% that you provide and outputs the points on the open boundaries, which you
+% will also provide.
+
+% Ingedients needed:
+
+% 1) domain_cfg.nc % Your NEMO grid
+% 2) bathy_meter.nc % Your bathymetry file
+% 3) coordinates.bdy.nc % Your coordinates file, generated when tides
+% you make tides
+
+% 4) Daily parent data. This has struction ./U016_dir/JAN/Sal_16_0000.nc
+
+% This code uses the parallel toolbox. This might have memory restrictions
+% based on the computer you are using. If you want this to run quickly, use
+% the following syntax: parpool(N) with N = 16 or N=31 if you have the
+% resources.
+
+% Note this is adapted only for Copernicus data but can easily adapted for
+% other datasets. Email a.brereton@liverpool.ac.uk for help.
+
+
+
+% we want to make the current directory visible so we can use subroutines
+addpath(pwd)
+
+% Control variables - make sure
+
+prefix='INDIAN_bdy_'; % prefix to output filename
+output_name='vozocrtx'; % output variable name (3D)
+mean_output_name='vobtcrtx'; % output variable name (2D)
+input_name='uo'; % variable from input file
+filename_prefix='U'; % e.g. 'INDIAN_bdy_U...'
+Input_directory='U017_dir'; % Input file directory
+year = '2017'; % e.g. 'INDIAN_bdy_U_y2017...'
+grid_data = 'domain_cfg.nc'; % NEMO grid information
+
+
+% e.g. 'INDIAN_bdy_U_y2017m01.nc'
+
+Input_month_suf{1}='01.nc';
+Input_month_suf{2}='02.nc';
+Input_month_suf{3}='03.nc';
+Input_month_suf{4}='04.nc';
+Input_month_suf{5}='05.nc';
+Input_month_suf{6}='06.nc';
+Input_month_suf{7}='07.nc';
+Input_month_suf{8}='08.nc';
+Input_month_suf{9}='09.nc';
+Input_month_suf{10}='10.nc';
+Input_month_suf{11}='11.nc';
+Input_month_suf{12}='12.nc';
+
+
+% Input directory month names e.g. ...../SAL16_dir/JAN/
+
+Input_month{1}='JAN';
+Input_month{2}='FEB';
+Input_month{3}='MAR';
+Input_month{4}='APR';
+Input_month{5}='MAY';
+Input_month{6}='JUN';
+Input_month{7}='JUL';
+Input_month{8}='AUG';
+Input_month{9}='SEP';
+Input_month{10}='OCT';
+Input_month{11}='NOV';
+Input_month{12}='DEC';
+
+
+
+% x co-ordinate indice of boundary
+
+nbidta = ncread('coordinates.bdy.nc','nbit'); % (yb, xb)
+
+
+
+
+
+% y co-ordinate indice of boundary
+
+nbjdta = ncread('coordinates.bdy.nc','nbjt'); % (yb,xb)
+
+
+
+
+
+% rim level, this code is only set up for a rimwidth of 1.
+
+nbrdta = nbidta*0+1; % (yb,xb)
+
+
+
+% nav_lon - we have this in domain_cfg (y,x)
+% nav_lat - we have this in domain_cfg (y,x)
+
+nav_lon = ncread(grid_data,'nav_lon');
+nav_lat = ncread(grid_data,'nav_lat');
+
+
+
+
+% create a land-sea mask from bathymetry
+
+B=ncread('bathy_meter.nc','Bathymetry');
+
+B(B>0)=1; % 1 is Water
+B(isnan(B))=0; % 0 is Land
+
+bdy_msk=B;
+
+
+
+
+
+
+for month=1:12
+
+
+
+% e.g. 'INDIAN_bdy_SAL_y2017m01.nc'
+
+output_filename = [prefix , filename_prefix , '_y',year,'m',Input_month_suf{month}];
+
+
+
+% cd into the data directory
+cd(Input_directory)
+cd(Input_month{month})
+filenames=dir('*.nc');
+
+
+
+% This loops through daily files, interpolating onto NEMO grid and
+% getting boundary points for the data.
+
+parfor j=1:length(filenames)
+
+ if j==1
+
+ [Variable_obc,z_obc] = COPERNICUS_INTERP(filenames(j).name,grid_data,variable_name,nbidta,nbjdta);
+
+ deptht(:,:,:,j) = z_obc;
+
+
+ else
+
+ [Variable_obc,~] = COPERNICUS_INTERP(filenames(j).name,grid_data,variable_name,nbidta,nbjdta);
+
+ end
+
+
+ NEMO_variable_obc(:,:,:,j)=Variable_obc;
+
+
+ % Get time information. NEMO wants seconds since 1950. Copernicus data
+ % comes in hours since 1900.
+
+ time(j)=ncread( filenames(j).name , 'time');
+ time_counter(j) = (time(j) + 18262*24)*3600; % 18262 days between 1900 and 1950
+
+end
+
+
+
+
+cd ../..
+
+
+
+
+for j=1:size(NEMO_variable_obc,4)
+ for i=1:size(NEMO_variable_obc,1)
+ zz=squeeze(deptht(i,1,:));
+ UU=squeeze( NEMO_variable_obc(i,1,:,j));
+
+ mean_NEMO_variable_obc(i,1,j) = trapz(zz,UU)/max(zz);
+ end
+end
+
+
+
+nccreate(output_filename,'nbrdta',...
+ 'Dimensions',{'xb',length(nbidta),'yb',1},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nbrdta',nbrdta);
+
+nccreate(output_filename,'nbidta',...
+ 'Dimensions',{'xb',length(nbidta),'yb',1},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nbidta',nbidta);
+
+nccreate(output_filename,'nbjdta',...
+ 'Dimensions',{'xb',length(nbidta),'yb',1},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nbjdta',nbjdta);
+
+nccreate(output_filename,'nav_lon',...
+ 'Dimensions',{'x',size(nav_lon,1),'y',size(nav_lon,2)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nav_lon',nav_lon);
+
+nccreate(output_filename,'nav_lat',...
+ 'Dimensions',{'x',size(nav_lat,1),'y',size(nav_lat,2)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nav_lat',nav_lat);
+
+nccreate(output_filename,'bdy_msk',...
+ 'Dimensions',{'x',size(bdy_msk,1),'y',size(bdy_msk,2)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'bdy_msk',bdy_msk);
+
+nccreate(output_filename,'time_counter',...
+ 'Dimensions',{'time_counter',Inf},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'time_counter',time_counter);
+
+nccreate(output_filename,'deptht',...
+ 'Dimensions',{'xb',size(deptht,1),'yb',1,'z',size(deptht,3)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'deptht',deptht);
+
+nccreate(output_filename,output_name,...
+ 'Dimensions',{'xb',size(NEMO_variable_obc,1),'yb',1,'z',size(NEMO_variable_obc,3),'time_counter',size(NEMO_variable_obc,4)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,output_name,NEMO_variable_obc);
+
+
+nccreate(output_filename,mean_output_name,...
+ 'Dimensions',{'xb',size(NEMO_variable_obc,1),'yb',1,'time_counter',size(NEMO_variable_obc,4)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,mean_output_name,mean_NEMO_variable_obc);
+
+
+
+end
\ No newline at end of file
diff --git a/SCRIPTS/obc_V.m b/SCRIPTS/obc_V.m
new file mode 100644
index 0000000..3254268
--- /dev/null
+++ b/SCRIPTS/obc_V.m
@@ -0,0 +1,249 @@
+clear all
+clc
+
+% This programme reads in Salinity data, interpolates on to the NEMO grid
+% that you provide and outputs the points on the open boundaries, which you
+% will also provide.
+
+% Ingedients needed:
+
+% 1) domain_cfg.nc % Your NEMO grid
+% 2) bathy_meter.nc % Your bathymetry file
+% 3) coordinates.bdy.nc % Your coordinates file, generated when tides
+% you make tides
+
+% 4) Daily parent data. This has struction ./U016_dir/JAN/Sal_16_0000.nc
+
+% This code uses the parallel toolbox. This might have memory restrictions
+% based on the computer you are using. If you want this to run quickly, use
+% the following syntax: parpool(N) with N = 16 or N=31 if you have the
+% resources.
+
+% Note this is adapted only for Copernicus data but can easily adapted for
+% other datasets. Email a.brereton@liverpool.ac.uk for help.
+
+
+
+% we want to make the current directory visible so we can use subroutines
+addpath(pwd)
+
+% Control variables - make sure
+
+prefix='INDIAN_bdy_'; % prefix to output filename
+output_name='vomecrty'; % output variable name (3D)
+mean_output_name='vobtcrty'; % output variable name (2D)
+input_name='vo'; % variable from input file
+filename_prefix='V'; % e.g. 'INDIAN_bdy_U...'
+Input_directory='V017_dir'; % Input file directory
+year = '2017'; % e.g. 'INDIAN_bdy_U_y2017...'
+grid_data = 'domain_cfg.nc'; % NEMO grid information
+
+
+% e.g. 'INDIAN_bdy_U_y2017m01.nc'
+
+Input_month_suf{1}='01.nc';
+Input_month_suf{2}='02.nc';
+Input_month_suf{3}='03.nc';
+Input_month_suf{4}='04.nc';
+Input_month_suf{5}='05.nc';
+Input_month_suf{6}='06.nc';
+Input_month_suf{7}='07.nc';
+Input_month_suf{8}='08.nc';
+Input_month_suf{9}='09.nc';
+Input_month_suf{10}='10.nc';
+Input_month_suf{11}='11.nc';
+Input_month_suf{12}='12.nc';
+
+
+% Input directory month names e.g. ...../SAL16_dir/JAN/
+
+Input_month{1}='JAN';
+Input_month{2}='FEB';
+Input_month{3}='MAR';
+Input_month{4}='APR';
+Input_month{5}='MAY';
+Input_month{6}='JUN';
+Input_month{7}='JUL';
+Input_month{8}='AUG';
+Input_month{9}='SEP';
+Input_month{10}='OCT';
+Input_month{11}='NOV';
+Input_month{12}='DEC';
+
+
+
+% x co-ordinate indice of boundary
+
+nbidta = ncread('coordinates.bdy.nc','nbit'); % (yb, xb)
+
+
+
+
+
+% y co-ordinate indice of boundary
+
+nbjdta = ncread('coordinates.bdy.nc','nbjt'); % (yb,xb)
+
+
+
+
+
+% rim level, this code is only set up for a rimwidth of 1.
+
+nbrdta = nbidta*0+1; % (yb,xb)
+
+
+
+% nav_lon - we have this in domain_cfg (y,x)
+% nav_lat - we have this in domain_cfg (y,x)
+
+nav_lon = ncread(grid_data,'nav_lon');
+nav_lat = ncread(grid_data,'nav_lat');
+
+
+
+
+% create a land-sea mask from bathymetry
+
+B=ncread('bathy_meter.nc','Bathymetry');
+
+B(B>0)=1; % 1 is Water
+B(isnan(B))=0; % 0 is Land
+
+bdy_msk=B;
+
+
+
+
+
+
+for month=1:12
+
+
+
+% e.g. 'INDIAN_bdy_SAL_y2017m01.nc'
+
+output_filename = [prefix , filename_prefix , '_y',year,'m',Input_month_suf{month}];
+
+
+
+% cd into the data directory
+cd(Input_directory)
+cd(Input_month{month})
+filenames=dir('*.nc');
+
+
+
+% This loops through daily files, interpolating onto NEMO grid and
+% getting boundary points for the data.
+
+parfor j=1:length(filenames)
+
+ if j==1
+
+ [Variable_obc,z_obc] = COPERNICUS_INTERP(filenames(j).name,grid_data,variable_name,nbidta,nbjdta);
+
+ deptht(:,:,:,j) = z_obc;
+
+
+ else
+
+ [Variable_obc,~] = COPERNICUS_INTERP(filenames(j).name,grid_data,variable_name,nbidta,nbjdta);
+
+ end
+
+
+ NEMO_variable_obc(:,:,:,j)=Variable_obc;
+
+
+ % Get time information. NEMO wants seconds since 1950. Copernicus data
+ % comes in hours since 1900.
+
+ time(j)=ncread( filenames(j).name , 'time');
+ time_counter(j) = (time(j) + 18262*24)*3600; % 18262 days between 1900 and 1950
+
+end
+
+
+
+
+cd ../..
+
+
+
+
+for j=1:size(NEMO_variable_obc,4)
+ for i=1:size(NEMO_variable_obc,1)
+ zz=squeeze(deptht(i,1,:));
+ UU=squeeze( NEMO_variable_obc(i,1,:,j));
+
+ mean_NEMO_variable_obc(i,1,j) = trapz(zz,UU)/max(zz);
+ end
+end
+
+
+
+nccreate(output_filename,'nbrdta',...
+ 'Dimensions',{'xb',length(nbidta),'yb',1},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nbrdta',nbrdta);
+
+nccreate(output_filename,'nbidta',...
+ 'Dimensions',{'xb',length(nbidta),'yb',1},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nbidta',nbidta);
+
+nccreate(output_filename,'nbjdta',...
+ 'Dimensions',{'xb',length(nbidta),'yb',1},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nbjdta',nbjdta);
+
+nccreate(output_filename,'nav_lon',...
+ 'Dimensions',{'x',size(nav_lon,1),'y',size(nav_lon,2)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nav_lon',nav_lon);
+
+nccreate(output_filename,'nav_lat',...
+ 'Dimensions',{'x',size(nav_lat,1),'y',size(nav_lat,2)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'nav_lat',nav_lat);
+
+nccreate(output_filename,'bdy_msk',...
+ 'Dimensions',{'x',size(bdy_msk,1),'y',size(bdy_msk,2)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'bdy_msk',bdy_msk);
+
+nccreate(output_filename,'time_counter',...
+ 'Dimensions',{'time_counter',Inf},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'time_counter',time_counter);
+
+nccreate(output_filename,'deptht',...
+ 'Dimensions',{'xb',size(deptht,1),'yb',1,'z',size(deptht,3)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,'deptht',deptht);
+
+nccreate(output_filename,output_name,...
+ 'Dimensions',{'xb',size(NEMO_variable_obc,1),'yb',1,'z',size(NEMO_variable_obc,3),'time_counter',size(NEMO_variable_obc,4)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,output_name,NEMO_variable_obc);
+
+
+nccreate(output_filename,mean_output_name,...
+ 'Dimensions',{'xb',size(NEMO_variable_obc,1),'yb',1,'time_counter',size(NEMO_variable_obc,4)},...
+ 'Format','netcdf4_classic');
+
+ncwrite(output_filename,mean_output_name,mean_NEMO_variable_obc);
+
+
+
+end
\ No newline at end of file
diff --git a/SCRIPTS/river_maker.m b/SCRIPTS/river_maker.m
new file mode 100644
index 0000000..726d511
--- /dev/null
+++ b/SCRIPTS/river_maker.m
@@ -0,0 +1,125 @@
+clear all
+clc
+
+% Programme to find nearest coastal wet point to river source.
+
+% River data is
+
+river_fname = '../DATA/coastal-stns-Vol-monthly.updated-oct2007.nc';
+
+% NEMO bathymetry file is...
+
+bathy_fname='../DATA/bathy_meter.nc';
+coord_fname='../DATA/coordinates.nc';
+
+sc= (1000^3)/(365*86400);
+
+flow=ncread(river_fname,'vol_stn');
+m2s=ncread(river_fname,'ratio_m2s');
+
+flw=flow.*m2s.*sc;
+lonr=ncread(river_fname,'lon_mou');
+latr=ncread(river_fname,'lat_mou');
+
+
+
+lon_nemo=ncread(bathy_fname,'nav_lon');
+lat_nemo=ncread(bathy_fname,'nav_lat');
+
+e1t=ncread(coord_fname,'e1t');
+e2t=ncread(coord_fname,'e2t');
+
+
+
+
+border=0.5; % don't accept rivers around the border of a certain width
+
+
+
+ind=lonr < (min(lon_nemo(:))+border);
+
+lonr=lonr(~ind);
+latr=latr(~ind);
+flw=flw(~ind);
+
+
+
+ind= lonr > (max(lon_nemo(:))-border);
+
+lonr=lonr(~ind);
+latr=latr(~ind);
+flw=flw(~ind);
+
+
+
+ind=latr < (min(lat_nemo(:))+border);
+
+lonr=lonr(~ind);
+latr=latr(~ind);
+flw=flw(~ind);
+
+
+
+ind= latr > (max(lat_nemo(:))-border);
+
+lonr=lonr(~ind);
+latr=latr(~ind);
+flw=flw(~ind);
+
+
+
+[lon_coast,lat_coast,coastal]=Coast_finder(bathy_fname);
+coastal=find(coastal);
+
+
+[id,D]=knnsearch([lon_coast lat_coast],[lonr latr]);
+
+
+% Let's make an upper limit of a distance of 5 degrees
+
+upper_limit=5;
+
+id=id(D<upper_limit);
+coastal=coastal(D<upper_limit);
+
+
+
+rho0=1000;
+nx=size(lon_nemo,1);
+ny=size(lon_nemo,2);
+rorunoff=zeros(nx,ny);
+
+
+rorunoff(coastal(1))=flw(1)*rho0/(e1t(coastal(1))*e2t(coastal(1)));
+
+for i=2:length(id)
+
+ runoff=flw(i)*rho0/(e1t(coastal(i))*e2t(coastal(i)));
+
+ if coastal(i) == coastal(i-1)
+
+ rorunoff(coastal(i))=runoff+rorunoff(coastal(i-1));
+
+ else
+
+ rorunoff(coastal(i))=runoff;
+
+ end
+
+
+end
+
+rorunoff=repmat(rorunoff,1,1,12);
+
+
+delete('river_test.nc');
+output='river_test.nc';
+nccreate(output,'time_counter','Dimensions',{'time_counter',Inf},'Format','netcdf4_classic');
+nccreate(output,'lat','Dimensions',{'lat',ny},'Format','netcdf4_classic');
+nccreate(output,'lon','Dimensions',{'lon',nx},'Format','netcdf4_classic');
+nccreate(output,'rorunoff','Dimensions',{'lon',nx,'lat',ny,'time_counter',12},'Format','netcdf4_classic');
+ncwrite(output,'lon',(1:nx));
+ncwrite(output,'lat',(1:ny));
+ncwrite(output,'rorunoff',rorunoff);
+ncwrite(output,'time_counter',(1:12));
+
--
GitLab