get_HYCOM_forcing

 Get mean flow, temperature and salinity data from HYCOM model outputs
 through their OPeNDAP server. 
 
 data = get_HYCOM_forcing(Mobj, modelTime)
 
 DESCRIPTION:
   Using OPeNDAP, extract the necessary parameters to create an FVCOM
   forcing file. Requires the OPeNDAP toolbox (see below for where to get
   it).
 
 INPUT: 
   Mobj - MATLAB mesh object
   modelTime - Modified Julian Date start and end times
 
 OUTPUT:
   data - struct of the data necessary to force FVCOM. These can be
   interpolated onto an unstructured grid in Mobj using
   grid2fvcom.m.
 
 The parameters which are obtained from the HYCOM data are:
     - temperature
     - salinity
     - u flow component
     - v flow component
     - sea surface height (ssh)
 
 REQUIRES:
   The OPeNDAP toolbox:
       http://www.opendap.org/pub/contributed/source/ml-toolbox/
 
 Author(s)
   Pierre Cazenave (Plymouth Marine Laboratory)
 
 Revision history:
   2013-01-31 First version.
 
==========================================================================

0001 function data = get_HYCOM_forcing(Mobj, modelTime)
0002 % Get mean flow, temperature and salinity data from HYCOM model outputs
0003 % through their OPeNDAP server.
0004 %
0005 % data = get_HYCOM_forcing(Mobj, modelTime)
0006 %
0007 % DESCRIPTION:
0008 %   Using OPeNDAP, extract the necessary parameters to create an FVCOM
0009 %   forcing file. Requires the OPeNDAP toolbox (see below for where to get
0010 %   it).
0011 %
0012 % INPUT:
0013 %   Mobj - MATLAB mesh object
0014 %   modelTime - Modified Julian Date start and end times
0015 %
0016 % OUTPUT:
0017 %   data - struct of the data necessary to force FVCOM. These can be
0018 %   interpolated onto an unstructured grid in Mobj using
0019 %   grid2fvcom.m.
0020 %
0021 % The parameters which are obtained from the HYCOM data are:
0022 %     - temperature
0023 %     - salinity
0024 %     - u flow component
0025 %     - v flow component
0026 %     - sea surface height (ssh)
0027 %
0028 % REQUIRES:
0029 %   The OPeNDAP toolbox:
0030 %       http://www.opendap.org/pub/contributed/source/ml-toolbox/
0031 %
0032 % Author(s)
0033 %   Pierre Cazenave (Plymouth Marine Laboratory)
0034 %
0035 % Revision history:
0036 %   2013-01-31 First version.
0037 %
0038 %==========================================================================
0039 
0040 subname = 'get_HYCOM_forcing';
0041 
0042 global ftbverbose;
0043 if(ftbverbose);
0044     fprintf('\n')
0045     fprintf(['begin : ' subname '\n'])
0046 end
0047 
0048 % Get the extent of the model domain (in spherical)
0049 if ~Mobj.have_lonlat
0050     error('Need spherical coordinates to extract the forcing data')
0051 else
0052     % Add a buffer of one grid cell in latitude and two in longitude to
0053     % make sure the model domain is fully covered by the extracted data.
0054     [dx, dy] = deal(2.5, 2.5); % NCEP resolution in degrees
0055     extents = [min(Mobj.lon(:))-(2*dx), max(Mobj.lon(:))+(2*dx), min(Mobj.lat(:))-dy, max(Mobj.lat(:))+dy];
0056 end
0057 
0058 if modelTime(end) - modelTime(1) > 365
0059     error('Can''t (yet) process more than a year at a time.')
0060 end
0061 
0062 [yearStart, monStart, dayStart, hStart, mStart, sStart] = mjulian2greg(modelTime(1));
0063 [yearEnd, monEnd, dayEnd, hEnd, mEnd, sEnd] = mjulian2greg(modelTime(end));
0064 
0065 if yearEnd ~= yearStart
0066     error('Can''t (yet) process across a year boundary.')
0067 else
0068     year = yearEnd;
0069 end
0070 
0071 if year < 2008 && monStart < 9
0072     error('Not using the legacy HYCOM model output. Select a start date from September 2008 onwards.')
0073 elseif (year >= 2008 && monStart >= 9) && (year < 2009 && monEnd <= 5)
0074     url = ['http://tds.hycom.org/thredds/dodsC/GLBa0.08/expt_90.6/', num2str(year)];
0075 elseif (year >= 2009 && monStart >= 5) && (year < 2011 && monEnd <= 1)
0076     url = ['http://tds.hycom.org/thredds/dodsC/GLBa0.08/expt_90.8/', num2str(year)];
0077 elseif (year >= 2011 && monStart >= 1)
0078     url = 'http://tds.hycom.org/thredds/dodsC/GLBa0.08/expt_90.9';
0079 else
0080     error('Date is bef?')
0081 end
0082 
0083 % Get the days of year for the start and end dates
0084 indStart = floor(datenum([yearStart, monStart, dayStart, hStart, mStart, sStart])) - datenum([yearStart, 1, 1]);
0085 indEnd = ceil(datenum([yearEnd, monEnd, dayEnd, hEnd, mEnd, sEnd])) - datenum([yearEnd, 1, 1]);
0086 % Create the necessary string for the time indices.
0087 tInd = sprintf('[%i:1:%i]', indStart, indEnd);
0088 
0089 % Set up a struct of the HYCOM data sets in which we're interested.
0090 hycom.temp  = [url, '?temperature'];
0091 hycom.salt  = [url, '?salinity'];
0092 hycom.u     = [url, '?u'];
0093 hycom.v     = [url, '?v'];
0094 hycom.ssh   = [url, '?ssh'];
0095 hycom.time  = [url, '?MT'];
0096 hycom.X     = [url, '?X'];
0097 hycom.Y     = [url, '?Y'];
0098 hycom.lon   = [url, '?Longitude'];
0099 hycom.lat   = [url, '?Latitude'];
0100 
0101 % % The HYCOM OPeNDAP server is spectacularly slow (at the moment?). As such,
0102 % % extracting even the entire lat/long arrays is painfully slow. Instead of
0103 % % directly interrogating the server, we'll assume a uniformly spaced 1/12th
0104 % % degree grid and use that to find the indices. We'll give ourself a bit of
0105 % % wiggle-room to accommodate the inevitable inaccuracy in this method.
0106 % [grid.X, grid.Y] = meshgrid(0:1/12:360, -90:1/12:90);
0107 % grid.ids = inbox(extents, grid.X, grid.Y);
0108 %
0109 % if extents(1) < 0 && extents(2) < 0
0110 %     % This is OK, we can just shunt the values by 360.
0111 %     extents(1) = extents(1) + 360;
0112 %     extents(2) = extents(2) + 360;
0113 %     index_lon = find(data_lon.lon > extents(1) & data_lon.lon < extents(2));
0114 % elseif extents(1) < 0 && extents(2) > 0
0115 %     % This is the tricky one. We'll do two passes to extract the
0116 %     % western chunk first (extents(1)+360 to 360), then the eastern
0117 %     % chunk (0-extents(2))
0118 %     index_lon{1} = find(data_lon.lon >= extents(1) + 360);
0119 %     index_lon{2} = find(data_lon.lon <= extents(2));
0120 % else
0121 %     % Dead easy, we're in the eastern hemisphere, so nothing too
0122 %     % strenuous here
0123 %     index_lon = find(data_lon.lon > extents(1) & data_lon.lon < extents(2));
0124 % end
0125 %
0126 % % Latitude is much more straightforward
0127 % data_lat = loaddap([ncep.(fields{aa}),'?lat']);
0128 % index_lat = find(data_lat.lat > extents(3) & data_lat.lat < extents(4));
0129 
0130 
0131 
0132 
0133 % Now get the latitude and longitude to calculate the indices representing
0134 % the model domain.
0135 data.X.idx = loaddap(hycom.X);
0136 data.Y.idx = loaddap(hycom.Y);
0137 xIdx = length(data.X.idx.X) - 1;
0138 yIdx = length(data.Y.idx.Y) - 1;
0139 data.lon.all = loaddap([hycom.lon, sprintf('[%i:1:%i]', 0, 0), sprintf('[%i:1:%i]', 0, xIdx)]);
0140 data.lat.all = loaddap([hycom.lat, sprintf('[%i:1:%i]', 0, yIdx), sprintf('[%i:1:%i]', 0, 0)]);
0141 
0142 fields = fieldnames(hycom);
0143 
0144 for aa = 1:length(fields)
0145     
0146     % Store the downloaded data in a struct with associated spatial and
0147     % temporal data.
0148     data.(fields{aa}).data = [];
0149     data.(fields{aa}).time = [];
0150     data.(fields{aa}).lat = [];
0151     data.(fields{aa}).lon = [];
0152     data_attributes.(fields{aa}) = [];
0153 
0154     % Get attributes from which to calculate length of various fields.
0155     data_attributes.(fields{aa}) = loaddap('-A', [hycom.(fields{aa})]);
0156     % Get the data time and convert to Modified Julian Day.
0157     data_time = loaddap(hycom.time);
0158 
0159     timevec = datevec((data_time.time)/24+365);
0160     data.time = greg2mjulian(timevec(:,1), timevec(:,2), timevec(:,3), ...
0161         timevec(:,4), timevec(:,5), timevec(:,6));
0162     % Clip the time to the given range
0163     data_time_mask = data.time >= modelTime(1) & data.time <= modelTime(end);
0164     data_time_idx = 1:size(data.time,1);
0165     data_time_idx = data_time_idx(data_time_mask);
0166     data.time = data.time(data_time_mask);
0167