update IC, BC, rivers

SCRIPTS/COPERNICUS_INTERP.m

+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
+
+

SCRIPTS/COPERNICUS_INTERP_SSH.m

+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
+
+

SCRIPTS/Coast_finder.m

+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

SCRIPTS/bdy_obc.m

+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)
+
+

SCRIPTS/nemo_IC_salinity.m

 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);
 
 

SCRIPTS/nemo_IC_temperature.m

 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);
 
 

SCRIPTS/obc_SSH.m

+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

SCRIPTS/obc_Sal.m

+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

SCRIPTS/obc_Temp.m

+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

SCRIPTS/obc_U.m

+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

SCRIPTS/obc_V.m

+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

SCRIPTS/river_maker.m

+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));
+