% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   SAMPLE MATLAB CODE TO READ OPENDAP FILES FROM topaz.nersc.no/thredds/
%
% Prerequisites: loaddap matlab routine, available from http://www.opendap.org/download/ml-structs.html
%
% See the matlab help of loaddap for more info.
%
% Written: 30th March 2006, Knut Liseter
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Check if loaddap.m is installed
S=which('loaddap.m');
if (isempty(S))
   disp('ERROR: loaddap.m is not in your path or is not installed!!!')
   disp('adjust matlab path, or install the ml-structs tools, available from:')
   disp('http://www.opendap.org/download/ml-structs.html')
else

   disp('***********************************************************************')
   disp('Retrieving latest TOPAZ forecast for the Arctic')


   %Get data set attributes - this retrieves attributes for the Arctic "Best Estimate" data set.
   %Other data sets are described (via a web browser) on http://topaz.nersc.no/thredds/
   OpenDAP_URL='http://topaz.nersc.no/thredds/dodsC/topaz/mersea-ipv2/arctic/tmipv2a-class1-b-be';
   att=loaddap('-A',OpenDAP_URL);

   % Note that "att" now contains info on all fields in the data set, explore it if you want to plot
   % other fields than those below

   % The topaz datasets uses offset and scaling factors to save disk space, the following retrieves
   % offset and scaling factors for uice, vice (east/north ice velocities) and fice (ice concentration)
   fice_sc=att.fice.scale_factor ; % fice scaling factor
   fice_os=att.fice.add_offset   ; % fice offset
   %
   uice_sc=att.uice.scale_factor ; % uice scaling factoror
   uice_os=att.uice.add_offset   ; % uice offset
   %
   vice_sc=att.vice.scale_factor ; % vice scaling factoror
   vice_os=att.vice.add_offset   ; % vice offset

   % The attributes contains a lot of useful information on the variables.
   % Here we retrieve the size of the time variable, since we are interested in the latest forecast
   % below.
   numtimes=att.time.DODS_ML_Size;


   %Encode last index into a text variable, used when calling loaddap
   last_index=numtimes-1;             % count starts from zero, so this is the last index
   clast=num2str(last_index,'%4.4i'); % Text version of last index


   % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Loaddap can be called something like this:
   % loaddap('http://topaz.nersc.no:80/thredds/dodsC/topaz/mersea-ip/arctic/tmipa-class1-be?fice')
   % which retrieves ice concentration (?fice) from the opendap url. This gets all ice concentration
   % fields for all times. To collect for a single time:
   % loaddap('http://topaz.nersc.no:80/thredds/dodsC/topaz/mersea-ip/arctic/tmipa-class1-be?fice[0:0]')
   % Similar logic applies when selecting for sub-regions and for specific depth levels. See help file for loaddap
   % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


   % Get latest forecast of ice concentration - 1 - From thredds
   % Dont be confused with all the string concatenation, the loaddap call is in principle the same as above
   loaddap([ OpenDAP_URL '?fice[' clast ':' clast ']'])
   loaddap([ OpenDAP_URL '?uice[' clast ':' clast ']'])
   loaddap([ OpenDAP_URL '?vice[' clast ':' clast ']'])

   % This is in principle all we need, since this call returns structures "fice"  "uice" and "vice" which
   % contain coordinate variables as well. This includes time and positions on grid. However, for the arctic topaz
   % fields, the grid is a map projection, so its nice to get the longitude and latitude corresponding to the grid 
   % positions as well:
   loaddap([ OpenDAP_URL '?longitude'])  % 2D field, so no time index require
   loaddap([ OpenDAP_URL '?latitude'])   % 2D field, so no time index require



   % Get fields from structures
   flduice=uice.uice;
   fldvice=vice.vice;
   fldfice=fice.fice;

   %% Get Time in serial days since 1 1 1970 (seconds since 1970 1  1 Default on thredds server)
   fldtime=fice.time/86400;

   % Time in serial days since 1 1 0
   fldtime=fldtime+datenum(1970,1,1);

   %Date string(s) - nice to have
   ctime=datestr(fldtime);
   ctime2=datestr(ctime,1);

   % Get lon, lat from structures
   fldlon=longitude.longitude;
   fldlat=latitude.latitude;
    
    
   % Scale stuff, using scale factors retreived from the attributes of the data set
   flduice=flduice*uice_sc+uice_os;
   fldvice=fldvice*vice_sc+vice_os;
   fldfice=fldfice*fice_sc+fice_os;

   % This is not strictly neccesary, but it sets missing field values to zero,
   % which produces nicer default plots
   I=find(fldfice<=0); 
   fldfice(I)=nan;
   flduice(I)=nan;
   fldvice(I)=nan;

   % Thats it, plot the sucker
   P=pcolor(fldfice); set(P,'LineStyle','none')

   disp('')
   disp('TOPAZ data retrieval successful')
   disp('***********************************************************************')
end