function SSP=ssp_from_gdem(gdm_root,out_root,BATHY,iih_sort)
%
% Call SSP=ssp_from_gdem(gdm_root,out_root,BATHY,iih_sort);
%
% Converts GDEM data in GDM files to SSP file.
%
% Example call: SSP=ssp_from_gdem('hce02',[GUI.adi_dir,'HCE02_GDEM'],BATHY,1)
%
%
% COPYRIGHT (c) 2008
% 
% Applied Research Laboratories, The University of Texas at Austin
% (ARL:UT)
% All rights reserved.
% 
% This material may be reproduced by or for the U.S. Government
% only. This notice must appear in all copies of this file and its
% derivatives. Further distribution outside the U.S. Government is
% prohibited without the explicit written authorization from the
% Director of the Environmental Sciences Laboratory, ARL:UT.
% 
% THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND,
% EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO THE
% IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
% PURPOSE, OR NONINFRINGEMENT.
% 
% IN NO EVENT WILL ARL:UT BE LIABLE FOR ANY DAMAGES WHATSOEVER ARISING
% OUT OF THE USE OR INABILITY TO USE THIS SOFTWARE, EVEN IF ARL:UT HAS
% BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.
if nargin<4, iih_sort=0; end
%
D=dir([gdm_root,'*.gdm']); nfile=length(D);
%
iidiag=0; nprof=0;
Jday=zeros(nfile,1); Jhour=Jday; lat=Jday; long=Jday; hw=Jday; file_names=cell(nfile,1);
disp(['Processing ',int2str(nfile),' GDEM Files ...']);
for jfile=1:nfile,
    disp(D(jfile).name);
    fid=fopen(D(jfile).name); hdone=0; nline_skip=0;
	eof=0; 
% Read in first line of new header for a new profile and detect EOF:
	[line]=fgetl(fid);
    if isempty(line), line=fgetl(fid); end
	if ~isstr(line) | isempty(line), eof=1; end
	while ~eof
		nprof=nprof + 1;
        file_names{nprof}=D(jfile).name;
% Extract latitude, longitude, and number of points in profile:
		lat(nprof)=sscanf(line(5:13),'%g');
		long(nprof)=sscanf(line(19:28),'%g');
		day1=sscanf(line(38:40),'%g'); day2=sscanf(line(42:44),'%g');
		Jday(nprof)=round((day1+day2)/2);
		Jhour(nprof)=0;
		[line]=fgetl(fid);
		nz=sscanf(line(19:25),'%g');
% Read in matrix of (z,c):
		[cz_arr,count]=fscanf(fid,'%g',[4,nz]); cz_arr=cz_arr';
        ssp{nprof}=cz_arr(:,[1 4]); if iidiag, disp(nz), disp(ssp{nprof}), end
        hw(nprof)=min(ssp{nprof}(end,1),interp2(BATHY.lat,BATHY.long,BATHY.hw.',lat(nprof),long(nprof)));
		[line]=fgetl(fid);
        if isempty(line), line=fgetl(fid); end
		if ~isstr(line) | isempty(line), disp('EOF: '), eof=1; [line]=fgetl(fid); end
	end
	fclose(fid);
end
%
hmin_gdem=1;
hw=floor(hw*2)/2; hwmax=max(hw);
z_out=unique([0:2:50,55:5:100,110:10:200,225:25:min(hwmax,500),600:100:hwmax,hw(:)']);
ii=find(z_out<=hwmax); z_out=z_out(ii); nz_out=length(z_out);
data=NaN*ones(nz_out,nprof); nz=zeros(1,nprof); svp_bad=zeros(1,nprof);
for j=1:nprof,
    ii=find(ssp{j}(:,1)>=hmin_gdem);
    data(:,j)=interp1(ssp{j}(ii,1),ssp{j}(ii,2),z_out(:),'linear','extrap');
    ii=find(z_out>hw(j)); data(ii,j)=NaN;
    if isempty(ii), nz(j)=nz_out; else, nz(j)=ii(1)-1; end
%    ii=find(data(:,j)<1460 | data(:,j)>1570); data(ii,j)=NaN;
%    iibad=find(isnan(data(1:nz(j),j))); nbad=length(iibad);
%    if nbad>0, svp_bad(j)=1; end
end
%
abs_time=Jday + Jhour/24;
[dummy,iiok]=sort(abs_time); nok=length(iiok);
%
SSP.ndim=2;
if ~iih_sort,
	SSP.labels={'Depth - m','SSP Number'};
	SSP.scales={z_out, [1:nok]};
else
% Construct second set of scales:
    [h_sort,ii_sort]=sort(hw(iiok));
    iiok=iiok(ii_sort);
	SSP.labels={'Depth - m','Water Depth'};
	SSP.scales={z_out, h_sort};
end
SSP.dims=[nz_out, nok];
SSP.data=data(:,iiok); SSP.dlab='Sound Speed - m/s';
SSP.hw=hw(iiok);
SSP.nz=nz(iiok); %SSP.ssp_no=ssp_no(iiok);
SSP.lat=lat(iiok); SSP.long=long(iiok);
SSP.Jday=Jday(iiok); SSP.Jhour=Jhour(iiok);
SSP.file_names=file_names(iiok);
SSP.title=out_root;
SSP.clim=[min(SSP.data(:)),max(SSP.data(:))];
%
str=['save ',out_root,'_SSP SSP']; disp(str); eval(str);
