function TMINf = fill_missing_temperature_days(TMIN, PRCP, PRCPf, tcal)
% Fill missing temperature values via conditional sampling, designed to be
% run after "stochastic_filling" of rainfall. Full calendar years are
% assumed.
%
% INPUTS:
% TMIN(nrec,nsta):  data to fill, one column per station, 
%                   with NaN missing data flag
% PRCP(nrec,nsta):  precip, with missing data flagged with NaN
% PRCPf(nrec,nsta): precip, with missing data already FILLED using SWG
% tcal(nrec,6):     matlab date vector 
%
% OUTPUT:
% TMINf(nrec,nsta): output temperature matrix with missing data filled
%
% 10-mar-2009; AWR

% transposition so function arguments match matlab convention!
TMIN=TMIN'; PRCP=PRCP'; PRCPf=PRCPf'; 
% first insert missing data records for Feb 29s
% "2" fields are 366-day year
[nsta nrec]=size(TMIN);
nyr=floor(nrec/365);

TMIN2=NaN*ones(nsta,366*nyr); 
PRCP2=NaN*ones(nsta,366*nyr); PRCPf2=NaN*ones(nsta,366*nyr);
i=0; idum=zeros(1,366*nyr);
for t=1:nrec
  i=i+1;
  idum(i)=0;
  TMIN2(:,i)=TMIN(:,t); 
  PRCP2(:,i)=PRCP(:,t); PRCPf2(:,i)=PRCPf(:,t);
  if tcal(t,2)==2 && tcal(t,3)==28 && tcal(t+1,3)==1
    i=i+1;
    TMIN2(:,i)=NaN; 
    PRCP2(:,i)=NaN; PRCPf2(:,i)=NaN;
    idum(i)=1; % tag
  end
end
TMIN2=reshape(TMIN2,nsta,366,nyr); 
PRCP2=reshape(PRCP2,nsta,366,nyr); PRCPf2=reshape(PRCPf2,nsta,366,nyr);
% calendar
cmon=[31 29 31 30 31 30 31 31 30 31 30 31];
cmonE=cumsum(cmon); % month endings
cmonB=cumsum(cmon)-cmon+1; % month starts

%% construct calendar-month subsets of wet & dry days & fill
TMIN2f=NaN*ones(nsta,366,nyr);
for s=1:nsta % station 1: Tmin is all missing for non-missing precip
  for m=1:12
    pre=PRCP2(s,cmonB(m):cmonE(m),:); % select cal month
    pref=PRCPf2(s,cmonB(m):cmonE(m),:); % filled precip
    minT=TMIN2(s,cmonB(m):cmonE(m),:);
    TMINdry=minT(~isnan(minT) & pre==0); % subset of non-missing dry days
    TMINwet=minT(~isnan(minT) & pre>0);  % subset of non-missing wet days

    if ~isempty(TMINdry)
      % missing days that are dry in filled-precip
      Nsmpl=length(find(isnan(minT) & pref==0)); % sample size
      minT(isnan(minT) & pref==0)=randsample(TMINdry,Nsmpl,'true');
    else
      disp(['TMINdry is empty for s=' num2str(s) ' and m=' num2str(m)]);
    end
    if ~isempty(TMINwet) 
      % missing days that are wet in filled-precip
      Nsmpl=length(find(isnan(minT) & pref>0)); % sample size
      minT(isnan(minT) & pref>0)=randsample(TMINwet,Nsmpl,'true'); 
    else
      disp(['TMINwet is empty for s=' num2str(s) ' and m=' num2str(m)]);
    end     
    % NB: missing TMIN days that have missing precip
    % don't exist because precip data has been filled with SWG   
    TMIN2f(s,cmonB(m):cmonE(m),:)=minT;
  end
end
% remove Feb 29s in non-leap years
TMIN2f=reshape(TMIN2f,nsta,366*nyr);
TMINf=TMIN2f(:,idum==0)';
