% loc2mods2
%
% This script converts an LOC netCDF file created by HydroBase to
% a MODS formatted ASCII file.

%-------------------------------------------------------------------------------
% Set run parameters.
%
%    basename.....Base name for input/output files.
%    mincnt.......Minimum acceptible number of data points that went into
%                 creating a climatology point.
%    in_suff......Suffix for input files.
%    out_suff.....Suffix for input files.
%    ttl..........Descriptive title for MODS formatted output file.
%    time.........Time to write to Climatology file. [yr mon day hr min sec]
%    julian_off...Off set to subtract from Julian dates.
%-------------------------------------------------------------------------------

%basename = ['7305'; ...
%            '7306'; ...
%            '7307'; ...
%            '7405'; ...
%            '7406'; ...
%            '7407'];
%basename = ['7306'; ...
%            '7307'; ...
%            '7405'; ...
%            '7406'; ...
%            '7407'];
%basename = ['7307'; ...
%            '7405'; ...
%            '7406'; ...
%            '7407'];
basename = ['7306'; ...
            '7307'; ...
            '7406'; ...
            '7407'];
%basename = '7306';

%basename = ['aug_all_GOM.7306'; ...
%            'aug_all_GOM.7307'; ...
%            'aug_all_GOM.7406'; ...
%            'aug_all_GOM.7407'];

%basename = 'massbay_april_m7_stab_top';

%mincnt = 5;
%mincnt = 3;
mincnt = 1;

in_suff  = '.nc';
out_suff = '.mods';
if (mincnt>1),
   out_suff = ['.',num2str(mincnt),out_suff];
end;

out_suff = ['.p25',out_suff];

%ttl = 'LOC 2 degree climatology, March';
%ttl = 'LOC 1 degree climatology, August';
%ttl = 'LOC 1/4 degree climatology, Mar';
%ttl = 'LOC 1/4 degree climatology, Jan';
%ttl = 'LOC 1/4 degree climatology, May';
ttl = 'LOC 1/4 degree climatology, Jul-Sep';
if (mincnt>1),
   ttl = [ttl,', ',num2str(mincnt),' data vals. min'];
end;

time = [2006 8 15 12 0 0];

julian_off = 2440000;

%-------------------------------------------------------------------------------
% Set julian date for casts.
%-------------------------------------------------------------------------------

jdate = julian(time) - julian_off + 0.5;

clear time;

%-------------------------------------------------------------------------------
% Set netCDF parameters.
%-------------------------------------------------------------------------------

status = mexcdf ('setopts',0);

nc_nowrite = mexcdf ('parameter','nc_nowrite');

%-------------------------------------------------------------------------------
% Loop over all files.
%-------------------------------------------------------------------------------

[nfile nchar] = size(basename);

for m = 1:nfile,

   %=========================
   %=== Open netCDF file. ===
   %=========================

   ncfile =  [basename(m,:),in_suff];
   disp (' ');
   disp (['   Working on: ',ncfile]);
   ncid   = mexcdf ('ncopen',ncfile,nc_nowrite);

   if (ncid<0),
      disp (' ');
      disp ('***Error:  LOC2MODS - unable to open file:');
      disp (['           ',setstr(34),ncfile,setstr(34)]);
      disp (' ');
      return;
   end;

   %============================
   %=== Get file dimensions. ===
   %============================

   %----------------------------------
   %--- Number of latitude points. ---
   %----------------------------------

   name = 'lat';
   [nlat,status] = get_dim (ncid,name);

   if (status<0),
      disp (' ');
      disp (['***Error:  LOC2MODS - unable to read dimension ',setstr(34), ...
                                            name,setstr(34),' in file:']);
      disp (['           ',setstr(34),ncfile,setstr(34)]);
      disp (' ');
      return;
   end;

   %-----------------------------------
   %--- Number of longitude points. ---
   %-----------------------------------

   name = 'lon';
   [nlon,status] = get_dim (ncid,name);

   if (status<0),
      disp (' ');
      disp (['***Error:  LOC2MODS - unable to read dimension ',setstr(34), ...
                                            name,setstr(34),' in file:']);
      disp (['           ',setstr(34),ncfile,setstr(34)]);
      disp (' ');
      return;
   end;

   %-------------------------------
   %--- Number of depth levels. ---
   %-------------------------------

   name = 'de';
   [nlev,status] = get_dim (ncid,name);

   if (status<0),
      disp (' ');
      disp (['***Error:  LOC2MODS - unable to read dimension ',setstr(34), ...
                                            name,setstr(34),' in file:']);
      disp (['           ',setstr(34),ncfile,setstr(34)]);
      disp (' ');
      return;
   end;

   %------------------------------
   %--- Number of time levels. ---
   %------------------------------

   name = 'time';
   [ntim,status] = get_dim (ncid,name);

   if (status<0),
      disp (' ');
      disp (['***Error:  LOC2MODS - unable to read dimension ',setstr(34), ...
                                            name,setstr(34),' in file:']);
      disp (['           ',setstr(34),ncfile,setstr(34)]);
      disp (' ');
      return;
   end;

   %========================
   %=== Get fill values. ===
   %========================

   %-------------------------------
   %--- Temperature fill value. ---
   %-------------------------------

   name = 'te';
   [temp_fill,status] = get_att (ncid,name,'_FillValue');

   if (status<0),
      disp (' ');
      disp (['***Error:  LOC2MODS - unable to read attribute ',setstr(34), ...
                               name,':_FillValue',setstr(34),' in file:']);
      disp (['           ',setstr(34),ncfile,setstr(34)]);
      disp (' ');
      return;
   end;

   %----------------------------
   %--- Salinity fill value. ---
   %----------------------------

   name = 'sa';
   [salt_fill,status] = get_att (ncid,name,'_FillValue');

   if (status<0),
      disp (' ');
      disp (['***Error:  LOC2MODS - unable to read attribute ',setstr(34), ...
                               name,':_FillValue',setstr(34),' in file:']);
      disp (['           ',setstr(34),ncfile,setstr(34)]);
      disp (' ');
      return;
   end;

   %==============================================
   %=== Extract useful casts from netCDF file. ===
   %==============================================

   %----------------------------
   %--- Extract cast depths. ---
   %----------------------------

   name = 'de';
   [zlev,status] = get_var (ncid,name,0,nlev);

   if (status<0),
      disp (' ');
      disp (['***Error:  LOC2MODS - unable to read variable ',setstr(34), ...
                                           name,setstr(34),' in file:']);
      disp (['           ',setstr(34),ncfile,setstr(34)]);
      disp (' ');
      return;
   end;

   %--------------------------------------------
   %---       Extract ending latitude.       ---
   %--- yes, those bozos decrement latitude. ---
   %--------------------------------------------

   name = 'latmax';
   [latmax,status] = get_att (ncid,[],name);

   if (status<0),
      disp (' ');
      disp (['***Error:  LOC2MODS - unable to read attribute ',setstr(34), ...
                                        ':',name,setstr(34),' in file:']);
      disp (['           ',setstr(34),ncfile,setstr(34)]);
      disp (' ');
      return;
   end;

   %-----------------------------------
   %--- Extract latitude step size. ---
   %-----------------------------------

   name = 'latincr';
   [dlat,status] = get_att (ncid,[],name);

   if (status<0),
      disp (' ');
      disp (['***Error:  LOC2MODS - unable to read attribute ',setstr(34), ...
                                        ':',name,setstr(34),' in file:']);
      disp (['           ',setstr(34),ncfile,setstr(34)]);
      disp (' ');
      return;
   end;

   %-----------------------------------
   %--- Extract starting longitude. ---
   %-----------------------------------

   name = 'lonmin';
   [lonmin,status] = get_att (ncid,[],name);

   if (status<0),
      disp (' ');
      disp (['***Error:  LOC2MODS - unable to read attribute ',setstr(34), ...
                                        ':',name,setstr(34),' in file:']);
      disp (['           ',setstr(34),ncfile,setstr(34)]);
      disp (' ');
      return;
   end;

   %------------------------------------
   %--- Extract longitude step size. ---
   %------------------------------------

   name = 'lonincr';
   [dlon,status] = get_att (ncid,[],name);

   if (status<0),
      disp (' ');
      disp (['***Error:  LOC2MODS - unable to read attribute ',setstr(34), ...
                                        ':',name,setstr(34),' in file:']);
      disp (['           ',setstr(34),ncfile,setstr(34)]);
      disp (' ');
      return;
   end;

   %----------------------------------------
   %--- Loop over all positions & times. ---
   %----------------------------------------

   castid = 0;

   for n = 0:(ntim-1),
      for j = 0:(nlat-1),
         for i = 0:(nlon-1),

            zwk = zlev;

            %----------------------------
            %--- Extract temperature. ---
            %----------------------------

            name = 'te';
            [twk,status] = get_var (ncid,name,[n j i 0],[1 1 1 nlev]);

            if (status<0),
               disp (' ');
               disp (['***Error:  LOC2MODS - unable to read variable ', ...
                           setstr(34),name,setstr(34),' in file:']);
               disp (['           ',setstr(34),ncfile,setstr(34)]);
               disp (['           (i,j,n) = (',num2str(i),',', ...
                              num2str(j),',',num2str(n),')']);
               disp (' ');
               return;
            end;

            if (mincnt>1),
               name = [name,'_cnt'];
               [tcnt,status] = get_var (ncid,name,[n j i 0],[1 1 1 nlev]);

               if (status<0),
                  disp (' ');
                  disp (['***Error:  LOC2MODS - unable to read variable ', ...
                              setstr(34),name,setstr(34),' in file:']);
                  disp (['           ',setstr(34),ncfile,setstr(34)]);
                  disp (['           (i,j,n) = (',num2str(i),',', ...
                                 num2str(j),',',num2str(n),')']);
                  disp (' ');
                  return;
               end;
            end;

            %-------------------------
            %--- Extract salinity. ---
            %-------------------------

            name = 'sa';
            [swk,status] = get_var (ncid,name,[n j i 0],[1 1 1 nlev]);

            if (status<0),
               disp (' ');
               disp (['***Error:  LOC2MODS - unable to read variable ', ...
                           setstr(34),name,setstr(34),' in file:']);
               disp (['           ',setstr(34),ncfile,setstr(34)]);
               disp (['           (i,j,n) = (',num2str(i),',', ...
                               num2str(j),',',num2str(n),')']);
               disp (' ');
               return;
            end;

            if (mincnt>1),
               name = [name,'_cnt'];
               [scnt,status] = get_var (ncid,name,[n j i 0],[1 1 1 nlev]);

               if (status<0),
                  disp (' ');
                  disp (['***Error:  LOC2MODS - unable to read variable ', ...
                              setstr(34),name,setstr(34),' in file:']);
                  disp (['           ',setstr(34),ncfile,setstr(34)]);
                  disp (['           (i,j,n) = (',num2str(i),',', ...
                                 num2str(j),',',num2str(n),')']);
                  disp (' ');
                  return;
               end;
            end;

            %------------------------
            %--- Remove bad data. ---
            %------------------------

            if (mincnt>1),
               ind = find ( (zwk<0) | (twk==temp_fill) | (swk==salt_fill) | ...
                            (tcnt<mincnt) | (scnt<mincnt) );
             else
               ind = find ( (zwk<0) | (twk==temp_fill) | (swk==salt_fill) );
            end;

            if (~isempty(ind)),
               no_pt  = length(zwk);
               if (length(ind)==no_pt),
                  zwk = [];
                  twk = [];
                  swk = [];
                elseif (min(ind)~=1),
                  ind2 = min(ind):no_pt;
                  zwk(ind2) = [];
                  twk(ind2) = [];
                  swk(ind2) = [];
                  clear ind2;
                else,
                  dind = diff(ind);
                  nst  = min(find(dind>1));
                  ind2 = ind(nst+1):no_pt;
                  zwk(ind2) = [];
                  twk(ind2) = [];
                  swk(ind2) = [];
                  ind2 = ind(1:nst);
                  zwk(ind2) = [];
                  twk(ind2) = [];
                  swk(ind2) = [];
                  clear dind ind2 nst;
               end;
               clear no_pt;
            end;

            %------------------
            %--- Save cast. ---
            %------------------

            if (~isempty(zwk)),

               castid = castid + 1;

               hinfo(castid,1) = 3;
               hinfo(castid,2) = max(size(zwk));
               hinfo(castid,3) = castid;
               hinfo(castid,4) = lonmin + (i+0.5).*dlon;
               hinfo(castid,5) = latmax - (j+0.5).*dlat;

               name = 'bottom';
               [zbot,status] = get_var (ncid,name,[n j i],[1 1 1]);

               if (status<0),
                  disp (' ');
                  disp (['***Error:  LOC2MODS - unable to read variable ', ...
                                 setstr(34),name,setstr(34),' in file:']);
                  disp (['           ',setstr(34),ncfile,setstr(34)]);
                  disp (['           (i,j,n) = (',num2str(i),',',num2str(j), ...
                                                ',',num2str(n),')']);
                  disp (' ');
                  return;
               end;

               hinfo(castid,6) = zbot;
               hinfo(castid,7) = jdate;
               hinfo(castid,8:10) = [0.1 1e-3 1e-3];
               hinfo(castid,11) = 0;

               dummy = [setstr(39),'CTD: z t s',setstr(39)];
               ind = 1:max(size(dummy));
               htype(castid,ind) = dummy;

               ind = 1:max(size(zwk));
               z(castid,ind) = zwk';
               t(castid,ind) = twk';
               s(castid,ind) = swk';

               zwk = [];
               twk = [];
               swk = [];

            end;
         end;
      end;
   end;

   clear castid dlat dlon i j latmin lonmin n name nlat nlev nlon ntim ...
         salt_fill swk temp_fill twk zbot zwk zlev;

   %=====================================
   %=== Close input netCDF data file. ===
   %=====================================

   status = mexcdf ('ncclose',ncid);

   if (status<0),
      disp (' ');
      disp ('***Error:  LOC2MODS - unable to close file:');
      disp (['           ',setstr(34),ncfile,setstr(34)]);
      disp (' ');
      return;
   end;

   clear ncfile status ncid;

   %========================
   %=== Write cast file. ===
   %========================

   if (exist('z')==1),

      %---------------------------
      %--- Create file header. ---
      %---------------------------

      dummy = [' title = ',ttl,':  ',basename(m,:)];
      ind = 1:max(size(dummy));
      header (1,ind) = dummy;

      [nsta, mxpts] = size(z);
      dummy = [' stations = ',num2str(nsta)];
      ind = 1:max(size(dummy));
      header (2,ind) = dummy;

      dummy = [' str_time = ',num2str(jdate,16),', ',pjh_greg(jdate-0.5)];
      ind = 1:max(size(dummy));
      header (3,ind) = dummy;

      dummy = [' end_time = ',num2str(jdate,16),', ',pjh_greg(jdate-0.5)];
      ind = 1:max(size(dummy));
      header (4,ind) = dummy;

      dummy = [' Jday_offset = ',num2str(julian_off,16)];
      ind = 1:max(size(dummy));
      header (5,ind) = dummy;

      dummy = [' lng_min = ',num2str(min(hinfo(:,4)),16)];
      ind = 1:max(size(dummy));
      header (6,ind) = dummy;

      dummy = [' lng_max = ',num2str(max(hinfo(:,4)),16)];
      ind = 1:max(size(dummy));
      header (7,ind) = dummy;

      dummy = [' lat_min = ',num2str(min(hinfo(:,5)),16)];
      ind = 1:max(size(dummy));
      header (8,ind) = dummy;

      dummy = [' lat_max = ',num2str(max(hinfo(:,5)),16)];
      ind = 1:max(size(dummy));
      header (9,ind) = dummy;

      dummy = ' format = ascii, record interleaving';
      ind = 1:max(size(dummy));
      header (10,ind) = dummy;

      dummy = ' type = CTD';
      ind = 1:max(size(dummy));
      header (11,ind) = dummy;

      dummy = ' fields = depth, temperature, salinity';
      ind = 1:max(size(dummy));
      header (12,ind) = dummy;

      dummy = ' units  = meter, Celsius, PSU';
      ind = 1:max(size(dummy));
      header (13,ind) = dummy;

      dummy = [' creation_date = ',date_stamp];
      ind = 1:max(size(dummy));
      header (14,ind) = dummy;

      dummy = 'END';
      ind = 1:max(size(dummy));
      header (15,ind) = dummy;

      clear dummy ind mxpts nsta;

   %------------------------
   %--- Write cast file. ---
   %------------------------

      hydnam = [basename(m,:),out_suff];
      disp (['   ... Creating: ',hydnam]);
      status = whydro (hydnam,header,hinfo,htype,z,t,s);

      if (status<0),
         disp (' ');
         disp ('***Error:  LOC2MODS - unable to write to file:');
         disp (['           ',setstr(34),hydnam,setstr(34)]);
         disp (' ');
         return;
      end;

      clear hydnam header hinfo htype z t s status;

    else

      disp (['   ... Eliminated all data from: ',[basename(m,:),in_suff]]);

   end;

end;
