write_FVCOM_tsobc

 example file for dumping a file to force temperature and salinity at the open b.

 function write_FVCOM_tsobc(basename,time,nSiglay,in_temp,in_salt)

 DESCRIPTION:
    Setup an FVCOM hydrographic open boundary forcing file. Supply either
    uniform values for temperature and salinity or 3D arrays (node,
    sigma_layers, time).

 INPUT
    Model case name (to find the bathymetry and open boundary .dat files).
    Time
    Number of sigma layers
    Boundary temperature (Celcius)
    Boundary salinity (psu)
    Mobj (optional)

 OUTPUT:
    FVCOM hydrographic open boundary file

 Author(s):
    Geoff Cowles (University of Massachusetts Dartmouth)
    Pierre Cazenave (Plymouth Marine Laboratory)

 Revision history
    2012-06-15 Added support for native MATLAB NetCDF routines. Requires
    MATLAB 2010a or higher.
    2012-07-16 Removed hard-coded nSiglay and nSiglev and instead moved to
    arguments list.
    2012-10-08 Updated help to reflect the fact nSiglev is calculated as
    nSiglay+1.
    2012-11-09 Add new arguments to use user defined temperature and
    salinity.
    2013-01-09 Add support for 3D input temperature and salinity (such as
    might be generated with get_POLCOMS_tsobc.m.
    KJT: Add better check for the size of the input arrays (works with 
    scalars).

==============================================================================

0001 function write_FVCOM_tsobc(basename,time,nSiglay,in_temp,in_salt,Mobj)
0002 % example file for dumping a file to force temperature and salinity at the open b.
0003 %
0004 % function write_FVCOM_tsobc(basename,time,nSiglay,in_temp,in_salt)
0005 %
0006 % DESCRIPTION:
0007 %    Setup an FVCOM hydrographic open boundary forcing file. Supply either
0008 %    uniform values for temperature and salinity or 3D arrays (node,
0009 %    sigma_layers, time).
0010 %
0011 % INPUT
0012 %    Model case name (to find the bathymetry and open boundary .dat files).
0013 %    Time
0014 %    Number of sigma layers
0015 %    Boundary temperature (Celcius)
0016 %    Boundary salinity (psu)
0017 %    Mobj (optional)
0018 %
0019 % OUTPUT:
0020 %    FVCOM hydrographic open boundary file
0021 %
0022 % Author(s):
0023 %    Geoff Cowles (University of Massachusetts Dartmouth)
0024 %    Pierre Cazenave (Plymouth Marine Laboratory)
0025 %
0026 % Revision history
0027 %    2012-06-15 Added support for native MATLAB NetCDF routines. Requires
0028 %    MATLAB 2010a or higher.
0029 %    2012-07-16 Removed hard-coded nSiglay and nSiglev and instead moved to
0030 %    arguments list.
0031 %    2012-10-08 Updated help to reflect the fact nSiglev is calculated as
0032 %    nSiglay+1.
0033 %    2012-11-09 Add new arguments to use user defined temperature and
0034 %    salinity.
0035 %    2013-01-09 Add support for 3D input temperature and salinity (such as
0036 %    might be generated with get_POLCOMS_tsobc.m.
0037 %    KJT: Add better check for the size of the input arrays (works with
0038 %    scalars).
0039 %
0040 %==============================================================================
0041 
0042 if nargin == 5
0043     warning(['Assuming uniform terrain-following sigma coordinates. ',...
0044         'To specify different sigma coordintes, supply a MATLAB mesh ',...
0045         'structure with fields siglay and siglev.'])
0046 end
0047 
0048 subname = 'write_FVCOM_tsobc';
0049 global ftbverbose;
0050 if(ftbverbose);
0051   fprintf('\n')
0052   fprintf(['begin : ' subname '\n'])
0053 end;
0054 
0055 fvcom_bathy = [basename, '_dep.dat'];
0056 fvcom_obc   = [basename, '_obc.dat'];
0057 tsOBCFile = [basename, '_tsobc.nc'];
0058 
0059 %------------------------------------------------------------------------------
0060 % read in the FVCOM open boundary node data (need node numbers and dimension)
0061 %------------------------------------------------------------------------------
0062 fid = fopen(fvcom_obc,'r');
0063 if(fid  < 0)
0064   error(['file: ' fvcom_obc ' does not exist']);
0065 end;
0066 C = textscan(fid, '%s %s %s %s %d', 1);
0067 nObc = C{5};
0068 obc_nodes = zeros(nObc,1);
0069 if(ftbverbose); fprintf('reading obc file\n'); end;
0070 if(ftbverbose); fprintf('# nodes %d\n',nObc); end;
0071 for i=1:nObc
0072   C = textscan(fid, '%d %d %d', 1);
0073   obc_nodes(i) = C{2};
0074 end;
0075 
0076 if(ftbverbose); fprintf('obc reading complete\n');end;
0077 
0078 %------------------------------------------------------------------------------
0079 % read in the FVCOM bathymetry data (need bathymetry on open boundary nodes)
0080 %------------------------------------------------------------------------------
0081 fid = fopen(fvcom_bathy,'r');
0082 if(fid  < 0)
0083   error(['file: ' fvcom_bathy ' does not exist']);
0084 end;
0085 C = textscan(fid, '%s %s %s %d', 1);
0086 Nverts = C{4};
0087 h = zeros(Nverts,1);
0088 if(ftbverbose); fprintf('reading bathymetry file\n');end;
0089 if(ftbverbose); fprintf('# nodes %d\n',Nverts);end;
0090 for i=1:Nverts
0091   C = textscan(fid, '%f %f %f', 1);
0092   h(i) = C{3};
0093 end;
0094 if(ftbverbose); fprintf('min depth %f max depth %f\n',min(h),max(h));end;
0095 if(ftbverbose); fprintf('bathymetry reading complete\n');end;
0096 fclose(fid);
0097 
0098 %--------------------------------------------------------------
0099 % Generate the requisite data
0100 %--------------------------------------------------------------
0101 
0102 % extract bathymetry at open boundary nodes
0103 obc_h = h(obc_nodes);
0104 
0105 % time
0106 % time = 0:1:31.;
0107 nTimes = numel(time);
0108 
0109 nSiglev = nSiglay + 1;
0110 
0111 % Create or process the temperature and salinity arrays.
0112 if max(size(in_temp)) == 1
0113     inc = 1/real(nSiglay);
0114     siglev = 0:-inc:-1;
0115     siglay = nan(1, nSiglay);
0116     for i=1:nSiglay
0117         siglay(i) = mean(siglev(i:i+1));
0118     end
0119     % initialize temperature/salinity arrays
0120     temp = zeros(nObc,nSiglay,nTimes);
0121     salt = zeros(nObc,nSiglay,nTimes);
0122 
0123     % set a constant temperature and salinity
0124     obc_temp = repmat(in_temp, 1, nTimes);
0125     obc_salt = repmat(in_salt, 1, nTimes);
0126 
0127     % set variable temperature and salinity
0128     % for i=1:nTimes
0129     %     obc_temp(i) = 18. + 2.*real(i-1)/nTimes;
0130     %     obc_salt(i) = 30. - 5.*real(i-1)/nTimes;
0131     % end
0132 
0133     % Create 3D array from three 1D arrays
0134     % temp = repmat(obc_temp, [nObc, nSiglay, 1]);
0135     % salt = repmat(obc_salt, [nObc, nSiglay, 1]);
0136     for i=1:nObc
0137         for j=1:nSiglay
0138             temp(i,j,:) = obc_temp;
0139             salt(i,j,:) = obc_salt;
0140         end
0141     end
0142 else
0143     % We have a 3D array already so we just need a couple of stats.
0144     temp = in_temp;
0145     salt = in_salt;
0146 
0147     if nargin == 6 && isfield(Mobj, 'siglay') && isfield(Mobj, 'siglev')
0148         siglev = Mobj.siglev;
0149         siglay = Mobj.siglay;
0150     else
0151         warning('Assuming uniform terrain-following sigma coordinates')
0152         inc = 1/real(nSiglay);
0153         siglev = 0:-inc:-1;
0154         siglay = nan(1, nSiglay);
0155     end
0156 
0157     if nSiglev ~= size(in_temp, 2) + 1 || length(siglev) ~= size(in_temp, 2) + 1 || length(siglev) ~= size(in_salt, 2) + 1
0158         error('Specified number sigma levels does not match supplied data')
0159     end
0160     if nSiglay ~= size(in_temp, 2) || length(siglay) ~= size(in_temp, 2) || length(siglay) ~= size(in_salt, 2)
0161         error('Specified number of sigma layers does not match supplied data')
0162     end
0163 end
0164 
0165 %--------------------------------------------------------------
0166 % set NetCDF variables and dump to file
0167 %--------------------------------------------------------------
0168 
0169 % open boundary forcing
0170 nc = netcdf.create(tsOBCFile, 'clobber');
0171 
0172 % define global attributes
0173 netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'type','FVCOM RIVER FORCING FILE')
0174 netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'title','simple open boundary hydrography test')
0175 netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'type','FVCOM TIME SERIES OBC TS FILE')
0176 netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'history','generated using the fvcom-toolbox')
0177 
0178 
0179 % define dimensions
0180 nobc_dimid=netcdf.defDim(nc,'nobc',nObc);
0181 datestrlen_dimid=netcdf.defDim(nc,'Datestrln',26);
0182 time_dimid=netcdf.defDim(nc,'time',netcdf.getConstant('NC_UNLIMITED'));
0183 siglay_dimid=netcdf.defDim(nc,'siglay',nSiglay);
0184 siglev_dimid=netcdf.defDim(nc,'siglev',nSiglev);
0185 
0186 % variables
0187 % nc{'river_names'} = ncchar('rivers', 'namelen');
0188 
0189 time_varid=netcdf.defVar(nc,'time','NC_FLOAT',time_dimid);
0190 netcdf.putAtt(nc,time_varid,'long_name','time');
0191 netcdf.putAtt(nc,time_varid,'units','days since 1858-11-17 00:00:00');
0192 netcdf.putAtt(nc,time_varid,'time_zone','UTC');
0193 
0194 itime_varid=netcdf.defVar(nc,'Itime','NC_INT',time_dimid);
0195 netcdf.putAtt(nc,itime_varid,'units','days since 1858-11-17 00:00:00');
0196 netcdf.putAtt(nc,itime_varid,'format','modified julian day (MJD)');
0197 netcdf.putAtt(nc,itime_varid,'time_zone','UTC');
0198 
0199 itime2_varid=netcdf.defVar(nc,'Itime2','NC_INT',time_dimid);
0200 netcdf.putAtt(nc,itime2_varid,'units','msec since 00:00:00');
0201 netcdf.putAtt(nc,itime2_varid,'time_zone','UTC');
0202 
0203 nobc_varid=netcdf.defVar(nc,'obc_nodes','NC_INT',nobc_dimid);
0204 netcdf.putAtt(nc,nobc_varid,'long_name','Open Boundary Node Number');
0205 netcdf.putAtt(nc,nobc_varid,'grid','obc_grid');
0206 netcdf.putAtt(nc,nobc_varid,'type','data');
0207 
0208 obc_h_varid=netcdf.defVar(nc,'obc_h','NC_FLOAT',nobc_dimid);
0209 netcdf.putAtt(nc,obc_h_varid,'long_name','Open Boundary Depth');
0210 netcdf.putAtt(nc,obc_h_varid,'units','m');
0211 netcdf.putAtt(nc,obc_h_varid,'grid','obc_grid');
0212 netcdf.putAtt(nc,obc_h_varid,'type','data');
0213 
0214 obc_siglev_varid=netcdf.defVar(nc,'siglev','NC_FLOAT',siglev_dimid);
0215 netcdf.putAtt(nc,obc_siglev_varid,'long_name','ocean_sigma/general_coordinate');
0216 netcdf.putAtt(nc,obc_siglev_varid,'grid','obc_grid');
0217 
0218 obc_siglay_varid=netcdf.defVar(nc,'siglay','NC_FLOAT',siglay_dimid);
0219 netcdf.putAtt(nc,obc_siglay_varid,'long_name','ocean_sigma/general_coordinate');
0220 netcdf.putAtt(nc,obc_siglay_varid,'grid','obc_grid');
0221 
0222 obc_temp_varid=netcdf.defVar(nc,'obc_temp','NC_FLOAT',[nobc_dimid,siglay_dimid,time_dimid]);
0223 netcdf.putAtt(nc,obc_temp_varid,'long_name','sea_water_temperature');
0224 netcdf.putAtt(nc,obc_temp_varid,'units','Celcius');
0225 netcdf.putAtt(nc,obc_temp_varid,'grid','obc_grid');
0226 
0227 obc_salinity_varid=netcdf.defVar(nc,'obc_salinity','NC_FLOAT',[nobc_dimid,siglay_dimid,time_dimid]);
0228 netcdf.putAtt(nc,obc_salinity_varid,'long_name','sea_water_salinity');
0229 netcdf.putAtt(nc,obc_salinity_varid,'units','PSU');
0230 netcdf.putAtt(nc,obc_salinity_varid,'grid','obc_grid');
0231 
0232 % end definitions
0233 netcdf.endDef(nc);
0234 
0235 % write data
0236 netcdf.putVar(nc,nobc_varid,obc_nodes);
0237 netcdf.putVar(nc,obc_h_varid,obc_h);
0238 netcdf.putVar(nc,obc_siglev_varid,siglev);
0239 netcdf.putVar(nc,obc_siglay_varid,siglay);
0240 netcdf.putVar(nc,time_varid,0,numel(time),time);
0241 netcdf.putVar(nc,itime_varid,floor(time));
0242 netcdf.putVar(nc,itime2_varid,0,numel(time),mod(time,1)*24*3600*1000);
0243 
0244 netcdf.putVar(nc,obc_temp_varid,temp);
0245 netcdf.putVar(nc,obc_salinity_varid,salt);
0246 
0247 % close file
0248 netcdf.close(nc);
0249 
0250 if(ftbverbose); fprintf(['end   : ' subname '\n']);end;