% Daily Transport Process v91. (2/15/2005)
% 
% Filename: daily_transport_process_v9.m (version 9)
% by: Chris Meinen & Rigo Garcia
%
% Revised version: 9/14/2005
% Changes: 12/12/2005
% Changes: 1/9/2006 Include a version to use the NI system
% Changes: 1/26/2006 Include to work with NI system voltage data after
%          remove the UW system
% Changes: 11/16/2006 Include the correction to the bad voltage period of
%          April 12 to July 12, 2006.
% Changes: 1/2/2008 Include the correction to the period of 2007, where
%          the transport estimated by the dropsonde were allways higher
%          than the cable estiamted...
% Changes: 8/12/2008 Include the history of how many dropsonde was used in
%          the past to calibrate the cable/voltage/transport...
%
% Changes: 10/14/2008 Include correction to the period of 2008, from
%          6/12/2008, when a new earth grond wire was installed.
%
% This script is the automated program that does the daily Florida Cable
% processing and produces the daily mean transport estimates.  It should be
% run in a Matlab batch mode using the system crontab.  
%
% This script is designed to:
% 1) Read the voltage for a particular day,
% 2) Calibrates the voltage from the Florida Current submarine 
%    cable (West Palm Beach, FL to Eight Mile Rock, Grand Bahama
%    Island) into transport,
% 3) Despike the raw transport time series,
% 4) Remove the tide and magnetic signals from the transport time series,
% 5) Compute the mean daily transport
% 6) Update the data tables and email the results to the group
%
% Input:  none
% Output: - Daily file that reports the program progess
%         - File that is sent to the web every day with a new transport
%           estimation (if it can be estimated). This file is for the
%           complete year.
%         - Figure that is sent to the web every day, showing the transport
%           estimated for each day in the context of the time series.
%
% Subroutines called: lowpassfilt.m - for low pass filtering
%

clear all
%%
%% The voltage from the Florida Current submarine cable (West Palm Beach,
%% FL to Eight Mile Rock, Grand Bahama Island), is the basic time series
%% that we need in order to estimate the transport in the Florida Straits. 
%% This time series is saved in a yearly database (DB) file called 
%% voltage_data_'year'.mat.  A small file is also saved with the date of 
%% the last day that the program successfully read in the voltage; the 
%% date value is saved in the file voltage_last_day_read.mat. 
%% 
%load voltage_last_day_read     % variable -> last_day (UW system)
load voltage_ni_last_day_read  % variable -> last_day (NI_system)
voltage_last_day = last_day;
clear last_day

%%
%% An additional small file is used to monitor the last day when the 
%% program was able to calculate the transport.  Next we load in the date 
%% on which the processing was last successful. 
%%
load transport_last_day_processed  % variable -> last_day

%% The cable time series shows voltage changes that are due to changes in
%% transport through the Florida Straits, however there are also very
%% significant voltage changes that are due purely to changes in the
%% geomagnetic fields.  At high frequencies (periods shorter than a few
%% days) there are large ionospheric signals, resulting from solar winds
%% and flares, and those ionospheric signals propagate down and affect the
%% magnetic field that is ambient around the submarine cable.  Jimmy Larsen
%% attempted to remove the geomagnetic signals from the cable data in order
%% to extract the high frequency transport variations when he was running
%% the cable project.  We have determined that we do not believe that we
%% can remove the geomagnetic signals accurately enough to make the high
%% frequency transport data useful.  It is known, based on coincident
%% magnetic and electric field measurements, that the geomagnetic
%% variations only affect the electric field/voltage measurements at
%% periods of three days and less.  Therefore we have chosen to remove the
%% magnetic signals by applying a three-day lowpass-filter to the voltage
%% data.  This requires us to have three days after the day we are
%% processing so that the value for the day in question is not negatively
%% impacted by edge effects in the filter.  Because the voltage data we are
%% downloading is two days behind (due to electronics issues), this means
%% we are processing data with a five day lag.  
%%
%% The number of days that we need to process are:
%%
day_to_process = (voltage_last_day - 3) - last_day ; % We need three days of voltage after the day that 
                                                     % we need process
%%
%% If we have day(s) to process,
%%
if day_to_process > 0  

%%
%% We loop through the individual days, processing them one at a time:
%%
for last_day = last_day + 1 : (voltage_last_day - 3);
    day_vec = datevec(last_day) ; 
    year = day_vec(1) ;                 
    stryear = num2str(year);
    month = day_vec(2);
    strmonth = num2str(month);
    day = day_vec(3);
    strday = num2str(day);
    jday = datetojulian(year,month,day);
    strjd = num2str(jday);
%%
%% The month must have two digits
%%
  if length(strmonth)==1
      strmonth=['0' strmonth];
  end
%%
%% The day must have two digits
%%
  if length(strday)==1
      strday=['0' strday];
  end
%%
%% The julian day must have three digits
%%
  if length(strjd)==1
      strjd=['00' strjd];
  elseif length(strjd)==2
          strjd=['0' strjd];
  end

%%
day_num=datenum(year,month,day);
%%  
%% We open the DB that has the records of voltage: voltage_data.mat 
%% We ask for the records of voltage that belong to that julian day and the
%% 26 days before and three days after (we want to apply the lowpass-filter
%% to a 30 day chunk of data).
%% 
% load voltage_data % -> datenumbers, voltage
load voltage

ind_days=find(datenumbers>=day_num-26 & datenumbers<day_num+4);

voltage_raw=voltage(ind_days); % Voltage for the complete days periods
datenumber=datenumbers(ind_days);
    
clear voltage datenumbers

%% We ajust here the voltage for the UW "bad" periods that were corrected:
ind1=find(datenumber>=datenum(2002,7,3,0,0,0) & datenumber<datenum(2002,8,12,23,59,59));
ind2=find(datenumber>=datenum(2002,8,14,0,0,0) & datenumber<datenum(2002,9,2,23,59,59));
ind3=find(datenumber>=datenum(2002,10,8,0,0,0) & datenumber<datenum(2003,5,7,23,59,59));
ind4=find(datenumber>=datenum(2003,5,11,0,0,0) & datenumber<datenum(2003,8,7,23,59,59));
ind5=find(datenumber>=datenum(2003,8,9,0,0,0) & datenumber<datenum(2003,9,8,23,59,59));
ind6=find(datenumber>=datenum(2003,9,13,0,0,0) & datenumber<datenum(2003,10,6,23,59,59));
%% Here we do the correction to the NI bad period during April, 2006 to July
%% 12, 2006 (Created in November 14, 2006)
ind7=find(datenumber>=datenum(2006,4,12,0,0,0) & datenumber<datenum(2006,7,11));

if length(ind1)~=0
    voltage_raw(ind1)=voltage_raw(ind1)+0.4272;
end
if length(ind2)~=0
    voltage_raw(ind2)=voltage_raw(ind2)-0.3753;
end
if length(ind3)~=0
    voltage_raw(ind3)=voltage_raw(ind3)+0.3179;
end
if length(ind4)~=0
    voltage_raw(ind4)=voltage_raw(ind4)+0.6271;
end
if length(ind5)~=0
    voltage_raw(ind5)=voltage_raw(ind5)+0.3884;
end
if length(ind6)~=0
    voltage_raw(ind6)=voltage_raw(ind6); % No ajustment in this period
end

% New instructions by Rigo Garcia in 11/16/2006
if length(ind7)~=0 % In this period we detrened the voltage data and are saved in the 
    %                matlab file bad_period_2006.mat. How we create this
    %                file is explained in the matlab scripts:
    %                best_day_to_fix_bad_period_May_Jul2006.m, and
    %                correct_voltage_bad_May2006_Jul2006.m
    %
    load bad_period_2006 % -> bad_dates, voltage_dt
    ind_bad=find(bad_dates>=datenumber(ind7(1)) & bad_dates<=datenumber(ind7(end)));
    voltage_corr=voltage_dt(ind_bad);
    voltage_raw(ind7)=voltage_corr;
end

%%
%% Create a Voltage Report variable
%%
voltage=reshape(voltage_raw,24*60,length(voltage_raw)/(24*60));
voltage_mean=nanmean(voltage);
voltage_nan=nansum(isnan(voltage));

Day=27;
voltage_day_mean=voltage_mean(Day);
voltage_day_minute_record=24*60-voltage_nan(Day);
voltage_day_bad_minute_record=0;

voltage_report=[voltage_day_mean;voltage_day_minute_record;voltage_day_bad_minute_record];

clear voltage voltage_day_* voltage_mean

%%
%% We must now to process the data, obtain the transport, and despike
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% New version created in 2/14/2005, 3 day-low pass filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% The raw transport is obtained using
%%
transport_raw=(-24.42*voltage_raw)-1;

%%
%% Now we define the periods where the UW system have been used:
%% From April 19, 2002 to 9/3/2004 and from 1/19/2005 to 8/3/2005
%% 12/8/2005 the UW system was connected at the same time that the new NI
%% system. "Good data" was obtained from 12/10/2005
%%
ind_uw_1=find(datenumber>=datenum(2002,4,19) & datenumber<=datenum(2004,9,4));
ind_uw_2=find(datenumber>datenum(2005,1,19) & datenumber<datenum(2005,8,3));
ind_uw_3=find(datenumber>=datenum(2005,12,9) & datenumber<datenum(2006,1,1));

%%
%% Now we define the periods where the "new BBIS"system have been used:
%% From Octuber, 2004 to 1/18/2005 and from 8/4/2005 to 12/8/2005
%%
ind_bbis_1=find(datenumber>=datenum(2004,10,29) & datenumber<=datenum(2005,1,18));
ind_bbis_2=find(datenumber>=datenum(2005,8,4) & datenumber<=datenum(2005,12,8));

%%
%% From 1/1/2006 we decide to use the new NI system to estimate the Florida
%% Transport. We disconnect the UW system the day January 24, 2006, but we
%% leave the first three days using the same way to estimate the transport.
%%
ind_ni_1=find(datenumber>=datenum(2006,1,1) & datenumber<datenum(2006,1,26));  % NI system
ind_ni_2=find(datenumber>=datenum(2006,1,27) & datenumber<datenum(2007,1,1));
ind_ni_3=find(datenumber>=datenum(2007,1,1) & datenumber<datenum(2007,12,5));
ind_ni_35=find(datenumber>=datenum(2007,12,5) & datenumber<=datenum(2007,12,31));
ind_ni_4=find(datenumber>=datenum(2008,1,1) & datenumber<=datenum(2008,3,12));
ind_ni_5=find(datenumber>=datenum(2008,3,13) & datenumber<=datenum(2008,6,12));
ind_ni_6=find(datenumber>=datenum(2008,6,13) & datenumber<datenum(2014,1,23)); % Adjusted in 5/9/2014
ind_ni_7=find(datenumber>=datenum(2014,1,23) & datenumber<datenum(2014,3,13)); % Add in 5/9/2014 for NIC period in 2014
ind_ni_8=find(datenumber>=datenum(2014,3,13)); % Add in 5/9/2014
%%
%% For those UW periods we apply an additional correction to take into
%% account the fact that due to a low resistance within the UW voltage
%% recording system (built by Rex Johnson at the Univ. of Washington) that
%% appears to put a load on the Florida cable itself, lowering the voltage
%% by 0.1 volts.  As a result we apply a correction of 2.4 volts, which
%% corresponds to 0.1 volt multiplied by the calibration coefficient of
%% 24.42 Sv/volt that was derived by Jimmy Larsen from comparison between
%% cable and Pegasus data in the early 1980s.  
%%
transport_raw(ind_uw_1)=(-24.42*voltage_raw(ind_uw_1)) - 1 + 2.4 ;
transport_raw(ind_uw_2)=(-24.42*voltage_raw(ind_uw_2)) - 1 + 2.4 ;
transport_raw(ind_uw_3)=(-24.42*voltage_raw(ind_uw_3)) - 1 + 2.4 ;

%%
%% For those "new BBIS" periods we apply an additional correction to take into
%% account the fact that due to ..... Wait for Chris's comments....
%% (9/14/2005)
%%
transport_raw(ind_bbis_1)=(-24.42*voltage_raw(ind_bbis_1)) - 1 - 8 ;
transport_raw(ind_bbis_2)=(-24.42*voltage_raw(ind_bbis_2)) - 1 - 8 ;

%%
%% When we connect the NI system together with the UW system 
%% account the fact that due to a low resistance within the UW voltage
%% recording system (built by Rex Johnson at the Univ. of Washington) that
%% appears to put a load on the Florida cable itself, lowering the voltage
%% by 0.2 volts in the UW.  As a result we apply a correction of 4.8 volts, which
%% corresponds to 0.2 volt multiplied by the calibration coefficient of
%% 24.42 Sv/volt that was derived by Jimmy Larsen from comparison between
%% cable and Pegasus data in the early 1980s.  
%%
transport_raw(ind_ni_1)=(-24.42*voltage_raw(ind_ni_1)) - 1 + 2.4 ;
transport_raw(ind_ni_2)=(-24.42*voltage_raw(ind_ni_2)) - 1 + 2.4 - 0.2*24.42 + 4;
transport_raw(ind_ni_3)=(-24.42*voltage_raw(ind_ni_3)) - 1 + 2.4 - 0.2*24.42 + 4 + 1.96;
transport_raw(ind_ni_35)=(-24.42*voltage_raw(ind_ni_35)) - 1 + 2.4 - 0.2*24.42 + 4;
transport_raw(ind_ni_4)=(-24.42*voltage_raw(ind_ni_4)) - 1 + 2.4 - 0.2*24.42 + 4 + 1.25; 
transport_raw(ind_ni_5)=(-24.42*voltage_raw(ind_ni_5)) - 1 + 2.4 - 0.2*24.42 + 4 - 0.5;
% transport_raw(ind_ni_6)=(-24.42*voltage_raw(ind_ni_6)) - 1 + 2.4 - 0.2*24.42 + 4 - 0.5 + 2.2;
transport_raw(ind_ni_6)=(-24.42*voltage_raw(ind_ni_6)) - 1 + 2.4 - 0.2*24.42 + 4 - 0.5 + 2.2 - 1.1;
transport_raw(ind_ni_7)=(-24.42*voltage_raw(ind_ni_7)) - 1 + 2.4 - 0.2*24.42 + 4 - 0.5 + 2.2 - 1.1 + 2.2; % Add 5/9/2014
transport_raw(ind_ni_8)=(-24.42*voltage_raw(ind_ni_8)) - 1 + 2.4 - 0.2*24.42 + 4 - 0.5 + 2.2 - 1.1; % Add 5/9/2014 (Same as ind_ni_6)

%
% From 1/23/2006 the UW system was disconnected and only the NI system is
% taking the voltage data. In this case we need to drecrece that value by:
% 0.2 * 24.42 volts = 4.8840 Sv
% Note: In March 31, 2006 we recalculate the transport doing a correction
% of 4 Sv (the last term in transport_raw equation, due to we realize three
% dropsonde cruise and compare with the cable transport estimate at those
% days. The mean of the difference if the 4 Sv correction.
% The days of the cruise were:
% 2/2/2006  -> Dropsonde: 27.6 Sv , Cable: 26.0 Sv
% 3/14/2006 -> Dropsonde: 33.5 Sv , Cable: 27.3 Sv
% 3/27/2006 -> Dropsonde: 28.3 Sv , Cable: 23.8 Sv
% The NI system was connected alone in January 24, 2006 but I don't want to
% reprossed the data the first three days, 24,25, and 26 because is the
% time that the NI system have to adjust itself, and start capturing good
% data...
%
% (1/2/2008) Note:
% The last year (2007) the dropsonde transpors estimated report values
% higher than the cable transport in the first seven cruise of that year.
% In December 5th, 2007 Uli change the cable connection and in the last
% cruise (12/7/2007) we obtain the transport estimated by the cable higher
% than the dropsonde for that day.
% To fix this problem, we calculate the mean of the difference between the
% dropsonde and the cable for that days (1.9571 Sv) and add this quantity
% to the calculation of cable transport in that period of time.
% The days of the cruise were:
% 3/29/2007  -> Dropsonde: 31.6 Sv; Cable: 29.8 Sv
% 6/19/2007  -> Dropsonde: 28.7 Sv; Cable: 27.2 Sv
% 7/10/2007  -> Dropsonde: 34.1 Sv; Cable: 32.7 Sv
% 9/27/2007  -> Dropsonde: 30.3 Sv; Cable: 28.3 Sv
% 10/4/2007  -> Dropsonde: 30.2 Sv; Cable: 27.4 Sv
% 11/6/2007  -> Dropsonde: 26.8 Sv; Cable: 24.9 Sv
% 11/28/2007 -> Dropsonde: 36.9 Sv; Cable: 34.6 Sv
%
% In January 3th, 2008 a new system recording was instaled in Bahamas: NIB
%
% 3/17/2008 Note:
% In March 6, 2008 the NIB recording system in Bahamas fail. On March 13th,
% 2008 a new recording system was instaled: NIC.
% In the period that NIB was used, two dropsonde cruise were done to
% compare with the transpot estimated by the NIB system. In both cruise the
% transport estimated using the dropsonde was higher than the cable.
% To adjust the cable data, we calculate the mean of the difference between the
% dropsonde and the cable for that days (1.25 Sv) and add this quantity
% to the calculation of cable transport in the period from 1/1/2008 to
% 3/12/2008.
% The days of the cruises were:
% 1/23/2008 -> Dropsonde: 29.5 Sv; Cable: 28.9 Sv
% 1/29/2008 -> Dropsonde: 28.4 Sv; Cable: 26.5 Sv

% 5/12/2008 Note:
% From March 13, 2008 we need to adjust the calibration according to the
% results of two more dropsonde (until we change the wire in the near
% future). Subtract the mean of the difference (0.5 Sv) between the cable and dropsonde
% for that days, to the calibration equation used in this period...
% The days of the cruises were:
% 04/22/2008 -> Dropsonde: 26.6 Sv; Cable: 26.3 Sv
% 05/07/2008 -> Dropsonde: 24.5 Sv; Cable: 25.8 Sv

% 8/12/2008 Note:
% We determinate how many dropsonde have been used in the past to calibrate
% the cable/voltage/transport in the FS. Some of then has been documented
% above. Before 2004, when Chris M. and I (RG) start working on this, we
% don't know for sure how man dropsonde were used to calibrated the cable
% transport estimated. We assume that at some point ALL of them were used
% to do so. After 2004, several ajustment (calibration) we make to the
% cable, due mainly to the fact that we change a lot the recording systems
% from 2004 to final of 2005/begining of 2006. In January of 2006 we
% instaled a new generation of recording system (NI system) but again, we
% have change a lot in these three year periods (until today Aug. 2008).
% Anyway, the dropsonde that have been used to calibrate the cable from
% 2004 until the present are:
% 01/09/2004 -> Dropsonde: 33.3 Sv
% 05/07/2004 -> Dropsonde: 29.3 Sv
% 05/24/2004 -> Dropsonde: 32.1 Sv
% 06/07/2004 -> Dropsonde: 28.9 Sv
% 06/11/2004 -> Dropsonde: 30.8 Sv
% 08/24/2004 -> Dropsonde: 35.1 Sv  We don't know the value that the cable have before ajustment
% 11/29/2004 -> Dropsonde: 32.8 Sv
% 02/17/2005 -> Dropsonde: 29.0 Sv
% 05/18/2005 -> Dropsonde: 31.5 Sv
% 08/31/2005 -> Dropsonde: 30.3 Sv
% 02/02/2006 -> Dropsonde: 27.6 Sv; Cable: 26.0 Sv
% 03/14/2006 -> Dropsonde: 33.5 Sv; Cable: 27.3 Sv
% 03/27/2006 -> Dropsonde: 28.3 Sv; Cable: 23.8 Sv
% 03/29/2007 -> Dropsonde: 31.6 Sv; Cable: 29.8 Sv
% 06/19/2007 -> Dropsonde: 28.7 Sv; Cable: 27.2 Sv
% 07/10/2007 -> Dropsonde: 34.1 Sv; Cable: 32.7 Sv
% 09/27/2007 -> Dropsonde: 30.3 Sv; Cable: 28.3 Sv
% 10/04/2007 -> Dropsonde: 30.2 Sv; Cable: 27.4 Sv
% 11/06/2007 -> Dropsonde: 26.8 Sv; Cable: 24.9 Sv
% 11/28/2007 -> Dropsonde: 36.9 Sv; Cable: 34.6 Sv
% 01/23/2008 -> Dropsonde: 29.5 Sv; Cable: 28.9 Sv
% 01/29/2008 -> Dropsonde: 28.4 Sv; Cable: 26.5 Sv
% 04/22/2008 -> Dropsonde: 26.6 Sv; Cable: 26.3 Sv
% 05/07/2008 -> Dropsonde: 24.5 Sv; Cable: 25.8 Sv

% 10/14/2008 Note of change in the script...
% On 6/12/2008 we have installed a new wire to the earth ground.
% From that date untill today, three dropsondes cruise have been taken
% place in the Florida Straits. The mean of the difference between
% Dropsonde and the Cable transport is 2.23 Sv. 
% Chris say that we use 2.2 Sv
% To adjust the cable results, we add this quantity (2.2 Sv) to the
% equation that calculate the cable transport from 6/13/2008 to the
% present.
% The dropsonde cruise's dates are:
% 07/14/2008 -> Dropsonde: 31.4 Sv , Cable: 30.1 Sv
% 10/02/2008 -> Dropsonde: 29.5 Sv , Cable: 27.8 Sv
% 10/09/2008 -> Dropsonde: 29.1 Sv , Cable: 25.4 Sv
% The cable transport estimates were re-processed on 10/14/2008...

% 2/18/2009 Note of change in the script to calculate the transport ...
% After the adjustment to the cable (see notes above) we have made two
% more dropsonde cruise. Using the adjusted cable result (after we apply
% the correction above), now we have:
% Dropsonde dates:
% 07/14/2008 -> Dropsonde: 31.4 Sv , Cable: 32.3 Sv
% 10/02/2008 -> Dropsonde: 29.5 Sv , Cable: 30.0 Sv
% 10/09/2008 -> Dropsonde: 29.1 Sv , Cable: 27.6 Sv
% 12/05/2008 -> Dropsonde: 35.4 Sv , Cable: 37.2 Sv
% 02/13/2009 -> Dropsonde: 35.4 Sv , Cable: 38.0 Sv
% 03/17/2009 -> Dropsonde: 33.6 Sv , Cable: 35.8 Sv
%
% The mean(dropsonde - cable) = -1.08 (-1.1) Sv.
% To adjust the cable results, we add this quantity (-1.1 Sv) to the
% equation that calculate the cable transport from 6/13/2008 to the
% present.
% IMPORTANT NOTE: We will wait until we have another dropsonde cruise to
% make the adjustment to the cable ..... (after a conversation with Chris
% Meinen).... Todas estas notas estan pendientes de aplicacion hasta que
% hagamos un nuevo crucero de dropsonde (2/18/2009)....
% IMPORTANT TASK: On March 23, 2009 we add the last dropsonde cruise taken 
% on 3/17/2009. Now we make the ajustment to the cable calibration
% equation...(see notes above)
% The cable transport estimates were re-processed on 03/23/2009...10:23 am
%
% 5/9/2014 Note: We want to ajust the voltage-transport equation because at
% the begining of 2014 several recording system were used. In particular,
% during the period of 1/23/2014 to 3/12/2014 the system NIC was used and
% two dropsonde cruises take place during that period.
% Dropsonde dates:
% 01/28/2014 -> Dropsonde: 33.8 Sv , Cable: 30.4 Sv
% 02/06/2014 -> Dropsonde: 30.9 Sv , Cable: 29.9 Sv
%
% The mean (dropsonde-cable) = 2.2 Sv
% For that period (1/23/2014 to 3/12/2014) we will add this quantity to the
% transport equation...
% From 3/13/2014 (when the NIB system was installed) to the future, we will
% continue using the same voltage-transport equation as before Jan. 23,
% 2014...
% New values of cable transport obtained (re-processed) in 5/9/2014 at
% 09:56 am

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% We now despike the transport time series
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  This is a very simple despike routine for the Florida Cable transport.  
%  Looking at several chunks of the data it was determined that the standard 
%  deviation of the transport time series was about 3-4 Sv, while there were 
%  occasional spikes (noise) that were much larger.  This routine simply cuts 
%  any values of transport larger than 3 standard deviations from the median 
%  value of transport for that day.  This routine requires the data to be
%  in units of even days (i.e. 3 days or 300 days but not 2.37 days).  Also
%  it assumes the data is minute data, and will crash if applied to
%  something that doesn't have a whole number multiple of 1440 samples.  
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  First toss out the obvious outliers (defined as a few Sverdrups larger 
%  or smaller than the max-min values observed during 1982-1998).  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dummy=transport_raw;

bad=find(dummy>50 | dummy<15);
dummy(bad)=NaN;
n_bad = length(bad);
clear bad

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Now do a more detailed cleanup using a 3 std criterion
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

numdays=length(dummy)/1440;
dummyRS=reshape(dummy,1440,numdays); % Reshape the transport array 

dummymedian=nanmedian(dummyRS);   % Calculate the median
dummymedian=(dummymedian'*ones(1,1440))';  % Reshape the data 

dummydiff=abs(dummyRS-dummymedian);  % Calculate the difference from daily median

[bad1,bad2]=find(dummydiff>(3*4));  %  Look for outliers beyond 3 std times 4 Sv

output=dummyRS;  % Replace outliers with NaNs
for i=1:length(bad1)
  output(bad1(i),bad2(i))=NaN;
end
n_bad1 = length(bad1);
clear i bad1 bad2

%%
%% Create some variables that will be needed for the Transport Report
%%

transport_nan_partial=nansum(isnan(output));
transport_day_nan=transport_nan_partial(Day);

transport_report=transport_day_nan-voltage_nan(Day);

%%

output=reshape(output,1440*numdays,1);  % Reshape the output array

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

transport_good=output;

clear dummy output  voltage_raw transport_raw


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Now we need to remove the known bad periods of voltage data.  These are
%% periods when we believe that there were serious electronics problems
%% with the voltage recording system, and hence we belive the voltage
%% measurements are not useful.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% These are the bad periods:
% 9/15/2001 (one day)
% 11/22/2001 - 1/26/2002
% 7/2/2002 - 10/8/2002
% 5/5/2003 - 10/13/2003
% 7/25/2005 - 8/3/2005 Note: The UW system fail and the "new" BBIS system
%                            was instaled
% 5/1/2006 - 7/12/2006: The wire connection faile and Uli Rivero change that
%                       connection (included in 7/20/2006. See index 11)
%
% Indices with problem:
ind_1=find(datenumber>=datenum(2001,9,15,0,0,0) & datenumber<datenum(2001,9,15,23,59,59)); ind_1=ind_1(:);
ind_2=find(datenumber>=datenum(2001,11,22,0,0,0) & datenumber<datenum(2002,1,26,23,59,59)); ind_2=ind_2(:);
%ind_3=find(datenumber>=datenum(2002,7,2,0,0,0) & datenumber<datenum(2002,10,8,23,59,59)); ind_3=ind_3(:);
ind_3=find(datenumber>=datenum(2002,7,2,0,0,0) & datenumber<datenum(2002,7,2,23,59,59)); ind_3=ind_3(:);
ind_4=find(datenumber>=datenum(2002,8,13,0,0,0) & datenumber<datenum(2002,8,13,23,59,59)); ind_4=ind_4(:);
ind_5=find(datenumber>=datenum(2002,9,3,0,0,0) & datenumber<datenum(2002,10,7,23,59,59)); ind_5=ind_5(:);
ind_6=find(datenumber>=datenum(2003,5,8,0,0,0) & datenumber<datenum(2003,5,10,23,59,59)); ind_6=ind_6(:);
ind_7=find(datenumber>=datenum(2003,8,8,0,0,0) & datenumber<datenum(2003,8,8,23,59,59)); ind_7=ind_7(:);
ind_8=find(datenumber>=datenum(2003,9,9,0,0,0) & datenumber<datenum(2003,9,12,23,59,59)); ind_8=ind_8(:);
ind_9=find(datenumber>=datenum(2003,10,7,0,0,0) & datenumber<datenum(2003,10,13,23,59,59)); ind_9=ind_9(:);
%ind_4=find(datenumber>=datenum(2003,5,5,0,0,0) & datenumber<datenum(2003,10,13,23,59,59)); ind_4=ind_4(:);
%ind_5=find(datenumber>=datenum(2005,7,25,0,0,0) & datenumber<datenum(2005,8,3,23,59,59)); ind_5=ind_5(:);
ind_10=find(datenumber>=datenum(2005,7,26,0,0,0) & datenumber<datenum(2005,8,3,23,59,59)); ind_10=ind_10(:);
%ind_11=find(datenumber>=datenum(2006,5,1,0,0,0) & datenumber<datenum(2006,7,12,23,59,59)); ind_11=ind_11(:);
ind_11=find(datenumber>=datenum(2006,7,10,0,0,0) & datenumber<datenum(2006,7,13,23,59,59)); ind_11=ind_11(:);
ind_12=find(datenumber>=datenum(2019,6,29,0,0,0) & datenumber<datenum(2019,7,31,23,59,59)); ind_12=ind_12(:);

% Group the indices and order them
ind=sort([ind_1;ind_2;ind_3;ind_4;ind_5;ind_6;...
          ind_7;ind_8;ind_9;ind_10;ind_11;ind_12]);

clear ind_*

% Replace the voltage values in those bad time periods with NaNs. 
transport_good(ind)=NaN;

%%
%% Reshape the DB of transport and the dates, and obtain the numbers of
%% days that we are processing
%%
n=length(transport_good);
minute_day=24*60;
number_day=n/minute_day;

dates=reshape(datenumber,minute_day,number_day);
transport=reshape(transport_good,minute_day,number_day);

transport_nan_sum=sum(isnan(transport)); % Row vector of daily number of NaNs
transport_nan=transport_nan_sum(:);
ind=find(transport_nan(:)<=minute_day/2); % Days with at least 720 good minutes

ind_nan=find(transport_nan(:)>minute_day/2); % Days with less than 720 good minutes
transport_nan(:)=0;
transport_nan(ind_nan)=1;

clear datenumber transport_good

%%
%% Calculate the daily mean of the raw transport
%%
dates_mean=mean(dates);
dates_mean=dates_mean';

transport_raw_daily_mean=NaN*ones(size(transport_nan));
flag=NaN*ones(size(transport_nan));

transport_raw_daily_mean(ind)=nanmean(transport(:,ind));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Obtain the transport using daily raw mean or three-day lowpass filter
%%
%% We have basically three possibilities in our processing scheme; 
%%
%%  1) There is less than 720 good minutes of data on a given day.  In this
%%     case we will report a NaN for the transport on that day because we
%%     do not think we can provide a reasonable estimate of the daily mean.
%%  2) There are more than 720 good minutes of data in a given day, but
%%     there is insufficient data in at least one of the days within three
%%     days before or after the processing day.  In this case we must use
%%     simply the daily mean of the transport because we cannot apply the
%%     three-day lowpass filter.  
%%  3) There are more than 720 good minutes of data on a given day, and
%%     there is sufficient data during both the three days before and after
%%     the processing day.  In this case we apply the three-day lowpass
%%     filter and then daily average the minute values of the filtered data
%%     to yield the final daily mean.  
%%
%%  There is a "flag" variable, called flag, that indicates which of the
%%  three above points is true for any given day.  

Day=27;  % Day that we are working on

transport_day=transport(:,Day);
transportcorr=NaN*ones(size(transport_day)); % If we can't do the filter we have NaNs

if transport_nan(Day)>0
            transport_daily_mean=transport_raw_daily_mean(Day); % NaN because of reason (1) above
            flag(Day)=2;
    elseif transport_nan(Day)<=0 & sum(transport_nan(Day-5:Day+3))>0 % sum(transport_nan)>0 % 
            transport_daily_mean=transport_raw_daily_mean(Day);
            flag(Day)=1; % We have data but we don't do filter (reason 2 above)
            
    elseif transport_nan(Day)<=0 & sum(transport_nan(Day-5:Day+3))<=0 % sum(transport_nan)<=0 %
            
            transport_days=transport(:);
            ind_nan_days=find(isnan(transport_days));
            
            transport_LP=lowpassfilt(interpNaN(transport_days),minute_day*3,1);
            transport_LP(ind_nan_days)=NaN;
            
            transport_LPr=reshape(transport_LP,minute_day,length(transport_nan));            
            
            transport_daily_mean=nanmean(transport_LPr(:,end-3));
            transportcorr=transport_LPr(Day);
            flag(Day)=0;
 end

flag=flag(Day);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% We obtain the Transport Daily Mean value only if we have more than a
%% half day of minute data (1400 minutes a day / 2)
%%

transport_nan = transport_nan_sum(Day);
%%
%% We also calculate the number of good point and the parameter that give
%% the "accuracy" of the daily mean transport in that day
%%
ind_good_minutes=(1440-transport_nan)/1440 ;

%%
%% We form the transport report that is sent by email every day to the
%% Cable Group
%%
transport_report = [transport_report;transport_nan; transport_daily_mean];
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Now we save the daily transport values in the corresponding Data Base
%transport_tide=transport_day;    % Provisional solution NEED CHANGE !!!!
transport_tide=NaN*ones(size(transport_day));    % We are no longer calculating the tide directly

%%
eval(['load transport_data_',stryear,';'])
eval(['ind = find (transport_data_',stryear,'(:,1)==jday);'])
eval(['transport_data_',stryear,'(ind,8)=transport_day ;'])
eval(['transport_data_',stryear,'(ind,9)=transport_tide ;'])
eval(['transport_data_',stryear,'(ind,10)=transportcorr ;'])
eval(['save transport_data_',stryear,' transport_data_',stryear,';'])
eval(['clear transport_data_',stryear,';'])
%%
eval(['load transport_processed_dates_',stryear,';'])
eval(['transport_processed_dates_',stryear,'(jday,3)=1 ;'])

eval(['save transport_processed_dates_',stryear,' transport_processed_dates_',stryear,';'])
eval(['clear transport_processed_dates_',stryear,';'])
%%
eval(['load transport_daily_mean_',stryear,';'])
eval(['transport_daily_mean_',stryear,'(jday,5)=transport_daily_mean ;'])
eval(['transport_daily_mean_',stryear,'(jday,6)=flag ;'])

eval(['save transport_daily_mean_',stryear,' transport_daily_mean_',stryear,';'])
%%
     eval(['temp = transport_daily_mean_',stryear,';'])
     dt = datenum(temp(:,2),temp(:,3),temp(:,4),12,0,0);
     T = temp(:,5); % Transport daily mean
     F = temp(:,6); % Flag columns
     Output = [ temp(:,2) temp(:,3) temp(:,4) T F];
     
%%
%% We save the last value of date and Transport, and Flag
T_d=transport_daily_mean ;
dt_d=dt(jday);
F_d=flag;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Open a file FC_transport_daily_mean.mat where all the Daily Mean Transport
%% from 2000 until the present are.
%%
load FC_transport_daily_mean  % Variables: datenumbers (datenum format); 
%                               transport_all (the complete list of the
%                               daily mean values of transport from 2000 to
%                               the present)
%                               flag: 0 - Good transport daily mean data,
%                                     1 - Transport daily mean estimated
%                                         without filter
%                                     2 - No enougth transport minute data
%

%%
%% We add now the new values of date and Transport daily mean
%% 
datenumbers=[datenumbers;dt_d];
transport_all=[transport_all;T_d];
flag=[flag;F_d];

[Y,I]=sort(datenumbers);
transport_all=transport_all(I);
flag=flag(I);

J=find(diff(Y)==0);
Y(J)=[];
transport_all(J)=[];
flag(J)=[];

% Redefining the variables
datenumbers=Y;

% Save the result
save FC_transport_daily_mean datenumbers transport_all flag

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Here we open and save a file that need to be send to the web
%%
 date_completed=['Table updated: ' date];
 
%  eval(['fid=fopen(''FC_cable_transport_',stryear,'.dat'',''w'');'])
%  
%  fprintf(fid,'%%  Transport estimate of the Florida Current  \n');
%  fprintf(fid,'%%  using voltage data of the submarine cable  \n');
%  fprintf(fid,'%%  from West Palm Beach, FL to Eight Mile Rock, \n');
%  fprintf(fid,'%%  Grand Bahama Island \n');
%  fprintf(fid,'%% \n');
%  fprintf(fid,'%% \n');
%  if Output(1,1)==2006
%      fprintf(fid,'%%  PLEASE NOTE: During the period from April 12, 2006     \n');
%      fprintf(fid,'%%               through July 9, 2006 we experienced some  \n');
%      fprintf(fid,'%%               technical difficulties with the Florida   \n');
%      fprintf(fid,'%%               Current cable equipment. Transport data   \n'); 
%      fprintf(fid,'%%               from this period has been corrected using \n');
%      fprintf(fid,'%%               dropsonde section data and statistical    \n');
%      fprintf(fid,'%%               analysis. We believe this data is good    \n');
%      fprintf(fid,'%%               quality, however some caution should be   \n');
%      fprintf(fid,'%%               used with data from this period.          \n');
%      fprintf(fid,'%%                                                         \n');
%      fprintf(fid,'%%                                                         \n');
%  elseif Output(1,1)==2019 % Add on 10/2/2019
%      fprintf(fid,'%%  PLEASE NOTE: After August 1, 2019, a new cable voltage \n');
%      fprintf(fid,'%%               recording system is in use. There can be  \n');
%      fprintf(fid,'%%               a small constant voltage offset that can  \n');
%      fprintf(fid,'%%               translate to transport estimates. We will \n'); 
%      fprintf(fid,'%%               be able to calibrate and remove any       \n');
%      fprintf(fid,'%%               offset that may be present with the new   \n');
%      fprintf(fid,'%%               recording system, but this will take us a \n');
%      fprintf(fid,'%%               few months as we collect enough dropsonde \n');
%      fprintf(fid,'%%               cruises to evaluate and update our        \n');
%      fprintf(fid,'%%               calibration. Until then, please use this  \n');
%      fprintf(fid,'%%               data with caution.                        \n');
%      fprintf(fid,'%%                                                         \n');
%      fprintf(fid,'%%                                                         \n');
%  elseif Output(1,1)==2020 % Add on 1/3/2020 (same as 2019, until we calibrate the voltage data)
%      fprintf(fid,'%%  PLEASE NOTE: After August 1, 2019, a new cable voltage \n');
%      fprintf(fid,'%%               recording system is in use. There can be  \n');
%      fprintf(fid,'%%               a small constant voltage offset that can  \n');
%      fprintf(fid,'%%               translate to transport estimates. We will \n'); 
%      fprintf(fid,'%%               be able to calibrate and remove any       \n');
%      fprintf(fid,'%%               offset that may be present with the new   \n');
%      fprintf(fid,'%%               recording system, but this will take us a \n');
%      fprintf(fid,'%%               few months as we collect enough dropsonde \n');
%      fprintf(fid,'%%               cruises to evaluate and update our        \n');
%      fprintf(fid,'%%               calibration. Until then, please use this  \n');
%      fprintf(fid,'%%               data with caution.                        \n');
%      fprintf(fid,'%%                                                         \n');
%      fprintf(fid,'%%                                                         \n');
% elseif Output(1,1)==2023 | Output(1,1)==2024 % Add on 1/26/2024
%      fprintf(fid,'%%  PLEASE NOTE: Beginning in November 2023, a mechanical  \n');
%      fprintf(fid,'%%               issue with the Florida Straits cable has  \n');
%      fprintf(fid,'%%               resulted in an interruption of the daily  \n');
%      fprintf(fid,'%%               Florida Current transport time series.    \n');
%      fprintf(fid,'%%               We are working to resolve the problem as  \n');
%      fprintf(fid,'%%               soon as possible.                         \n');
%      fprintf(fid,'%%                                                         \n');
%      fprintf(fid,'%%                                                         \n');
%  end 
% % 
%  fprintf(fid,'%%  For further information please contact\n');
%  fprintf(fid,'%%  Dr. Denis Volkov\n');
%  fprintf(fid,'%%  Denis.Volkov@noaa.gov\n');
%  fprintf(fid,'%%  (305) 361-4344 \n');
%  fprintf(fid,'%% \n%% \n');
%  fprintf(fid,'%%  Flag: 0 - Good data\n');
%  fprintf(fid,'%%        1 - Estimated without filter\n');
%  fprintf(fid,'%%        2 - No data\n');
%  fprintf(fid,'%% \n%% \n');
%  fprintf(fid,['%%  ' date_completed '\n%% \n%%\n']);
%  fprintf(fid,'%%  Year  Month  Day   Transport   Flag\n');
%  fprintf(fid,'%%                       [Sv] \n');
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 eval(['fid=fopen(''FC_cable_transport_',stryear,'_v2.dat'',''w'');'])
 fprintf(fid,['%%  DAILY FLORIDA CURRENT TRANSPORT ESTIMATES FOR ' stryear '  \n']); 

 fprintf(fid,'%% \n');
 fprintf(fid,'%%  These estimates are correlated to voltage data continuously recorded from  \n');
 fprintf(fid,'%%  a submarine telecommunications cable spanning the Florida Straits between  \n');
 fprintf(fid,'%%  West Palm Beach, Florida and Eight Mile Rock, Grand Bahama Island. \n');
 fprintf(fid,'%% \n');
 fprintf(fid,'%%  This data product is produced and publicly distributed by NOAA/AOML''s \n');
 fprintf(fid,'%%  Western Boundary Time Series (WBTS) project.  Funding is provided by \n');
 fprintf(fid,'%%  NOAA/AOML and NOAA''s Global Ocean Monitoring and Observing (GOMO) Program. \n');
 fprintf(fid,'%% \n');
 fprintf(fid,['%%  FILENAME: FC_cable_transport_' stryear '_v2.dat \n']); % Tratar de cambian el 2024...
 fprintf(fid,'%% \n');
 fprintf(fid,'%%  NOTE: The voltage data used to generate the daily Florida Current \n');
 fprintf(fid,'%%  transport estimates are now being corrected for the secular change in the \n');
 fprintf(fid,'%%  earth''s geomagnetic field over time. See Volkov et al., 2024:\n');
 fprintf(fid,'%%  https://doi.org/10.1038/s41467-024-51879-5 \n');
 fprintf(fid,'%%  \n');
 if Output(1,1)>=2023 % Temporary statement (to be removed after the cable system be restored)
     fprintf(fid,'%%  NOTE: Beginning in November 2023, a mechanical issue with the Florida\n');
     fprintf(fid,'%%  Straits submarine telecommunications cable has resulted in an interruption\n');
     fprintf(fid,'%%  of this daily Florida Current transport time series.  We are working to \n');
     fprintf(fid,'%%  resolve the problem as soon as possible. \n');
     fprintf(fid,'%%  \n');
 end
 fprintf(fid,'%%  For further information, please visit the project web site at:\n');
 fprintf(fid,'%%  https://www.aoml.noaa.gov/western-boundary-time-series/ \n');
 fprintf(fid,'%%  or contact: \n');
 fprintf(fid,'%%  \n');
 fprintf(fid,'%%  Dr. Denis Volkov\n');
 fprintf(fid,'%%  Denis.Volkov@noaa.gov\n');
 fprintf(fid,'%%  (305) 361-4344 \n');
 fprintf(fid,'%%  \n%% \n');
 fprintf(fid,'%%  COLUMN 1 - Year \n');
 fprintf(fid,'%%  COLUMN 2 - Month \n');
 fprintf(fid,'%%  COLUMN 3 - Day \n');
 fprintf(fid,'%%  COLUMN 4 - Transport in Sverdrups (Sv), 1 Sv = 1,000,000 meters^3 / second \n');
 fprintf(fid,'%%  COLUMN 5 - Quality Flag\n');
 fprintf(fid,'%%  \n');
 fprintf(fid,'%%  Quality Flag:  0 - Cable data are good \n');
 fprintf(fid,'%%                 1 - Cable data are estimated without a filter \n');
 fprintf(fid,'%%                 2 - No cable data are available \n');
 fprintf(fid,'%% \n%% \n');
 fprintf(fid,['%%  ' date_completed '\n%% \n%%\n']);
 fprintf(fid,'%%  Year  Month  Day   Transport   Flag\n');
 fprintf(fid,'%%                       [Sv] \n');
% New on 9/6/2024
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ii = 1 : jday
    for jj = 1 : 5
        if jj==1 
           fprintf(fid,'%7.0f ',Output(ii,jj));
          elseif jj==2 
           fprintf(fid,'%5.0f ',Output(ii,jj));
          elseif jj==3 
           fprintf(fid,'%4.0f ',Output(ii,jj));
          elseif jj==4 
           fprintf(fid,'%9.1f ',Output(ii,jj));
          elseif jj==5 
           fprintf(fid,'%7.0f ',Output(ii,jj));
        end
     end
   fprintf(fid,'\n');
end
fclose(fid);

%%%%%%%%%%%%
eval(['clear transport_daily_mean_',stryear,';'])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Here we save the last day that the transport has been processed
 save transport_last_day_processed last_day
 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Now we send the Transport's table to the Florida Current Webpage:
%% www.aoml.noaa.gov/phod/floridacurrent
%%
table=send_transport_table_to_web(last_day);  % Puesto a prueba el 11/17/2008
clear table
% !./send_transport_table_to_web.sh    % Si falla, poner esta instruccion
                                       % en la function de arriba
                                       % (send_transport_tab....)
                                       % I include this sentence in the
                                       % matlab script:
                                       % send_transport_table_to_web.m
                                       %

% This statement was commented in 7/10/2006 due to we was experiencing some
% technical difficulties with the Florida Current cable equipment.
% See "new table" called FC_cable_transport_2006_partial.dat
% On 7/25/2006 we uncomment the statement and restarted the FC webpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Here we create the Daily Process Report
%%
%% Daily Report:
%%
strdate = [strmonth '_' strday '_' stryear];

%%
%% The email address of people that we send the daily transport report
%%
%mail_address = 'rigoberto.garcia@noaa.gov';
mail_address = 'cable.transport@noaa.gov';

%%

 eval(['fid=fopen(''Cable_report_',strdate,'.txt'',''w'');'])
 fprintf(fid,'Daily Cable Processing Report - File generated by daily_transport_process_v91.m \n\n');
 fprintf(fid,'This report describes the basic processing steps needed to obtain the \n');
 fprintf(fid,'transport estimates of the Florida Current, using voltage data measured \n');
 fprintf(fid,'on the submarine telephone cable from West Palm Beach, FL, to Eight Mile \n');
 fprintf(fid,'Rock, Grand Bahama Island. \n\n');
 
 fprintf(fid,'Step   I: Read the minute voltage values from the submarine telephone cable.\n');
 fprintf(fid,'          Should include 1440 records. Typical range: -1.8 to -1.1 volts \n\n');

 fprintf(fid,'Step  II: Estimate the Florida Current Transport, which is obtained by:\n');
 fprintf(fid,'          1) Scale voltage into transport: Transport = (-24.42*volt)-1\n');
 fprintf(fid,'          2) Removing large spikes \n');
 fprintf(fid,'          3) Apply a three-day low pass filter\n');
 fprintf(fid,'          4) Calculate the daily transport average \n\n');

 %% Day that the report is generated
 fprintf(fid,'Date: ');
 fprintf(fid,'%2.0f',month);
 fprintf(fid,'/%2.0f',day);
 fprintf(fid,'/%4.0f',year);
 
 %% Voltage report
 fprintf(fid,'\n\nStep I report: \n\n');
 
 if length(voltage_report) == 1
    fprintf(fid,'We can''t estimate the FC transport today: \n');
    fprintf(fid,'The voltage data are not available \n\n\n\n\n');
    fprintf(fid,'For more information, please contact Denis Volkov (x 4344) or Rigo Garcia (x 4348).\n');
    subject = ['WARNING_Cable_Process_Fail_' strdate] ;
    eval(['!mail -s ',subject,' ',mail_address,' < Cable_report_',strdate,'.txt ;'])
else
    fprintf(fid,'Daily Mean Voltage: ');
    fprintf(fid,'%25.2f ',voltage_report(1));
    fprintf(fid,'volts\n');
    fprintf(fid,'Total number of minute records: ');
    fprintf(fid,'%13.0f ',voltage_report(2));
    fprintf(fid,'\nNumber of records with bad time values*: ');
    fprintf(fid,'%4.0f ',voltage_report(3));
    
    %% Transport process report

    fprintf(fid,'\n\nStep II report: \n\n');
    fprintf(fid,'Number of minute transport values removed during despike: ');
    fprintf(fid,'%4.0f ',transport_report(1));
    transport_records = 1440 - transport_report(2);
    fprintf(fid,'\nFinal number of good minute transport estimates: ');
    fprintf(fid,'%13.0f ',transport_records);
    fprintf(fid,'\nDaily Mean Transport: ');
    fprintf(fid,'%40.1f ',transport_report(3)');
    fprintf(fid,'Sv');
    fprintf(fid,'\n\n\n');
    fprintf(fid,'* These observations erroneously list a time of 24:00:00. \n\n\n');
    
    fprintf(fid,'For more information, please contact Denis Volkov (x 4344) or Rigo Garcia (x 4348).\n');
    
    %% Send email
       subject = ['Cable_report_' strdate];
       eval(['!mail -s ',subject,' ',mail_address,' < Cable_report_',strdate,'.txt ;'])
%    end
 end
fclose(fid);


%% Put the daily process report in the correct directory
Report_directory = ['/home/cable/Florida_Current_Transport/Cable_Process/Cable_Report/Cable_Report_',stryear,'/'];
eval(['!cp Cable_report_',strdate,'.txt ',Report_directory,';'])
eval(['!rm Cable_report_',strdate,'.txt ;'])


%%
%% Remove all the report variables. It is important because we do not need to save the
%% every day variables (.mat)

% eval(['!rm ',voltage_report_directory,'voltage_report_',strjd,'_',stryear,'.mat;'])
% eval(['!rm ',geomagnetic_report_directory,'geomagnetic_report_bsl_',strjd,'_',stryear,'.mat;'])
% % eval(['!rm geomagnetic_report_sjg_',strjd,'_',stryear,'.mat;'])
% % eval(['!rm geomagnetic_report_frd_',strjd,'_',stryear,'.mat;'])
    
end % day_to_process loop

else   % if day_to_process > 0
    'The transport process is updated. We used the last day of voltage data.'
end
