function [p_d,p_lat,p_lon] = make_vslite_pseudos_2(mdl_prior,fmeth)

%MAKE_VSLITE_PSEUDOS Make VS-lite based pseudoproxies for a specified
% length of time

% Get model information
%[~,gcm_syr,gcm_eyr]=model_info(mdl_prior);
[gcm_syr,gcm_eyr]=modelyrs(mdl_prior);

% Use all available data for prior
p_yrs=(gcm_syr:gcm_eyr)-gcm_syr+1;

if mdl_prior=='e'
    % Surface temperature
    %load('/d3/nsteiger/echam5wiso/e5_latlon.mat')
    %load('/d3/nsteiger/echam5wiso/tas_mon_e5.mat')
    disp('USING ECHAM5 COARSENED TO SPEEDY RESOLUTION')
    load('/d3/nsteiger/echam5wiso/echam5_coarsened2speedy_vars.mat')
    tas_mon=tas_crsn2speedy+273.15;
    lat=clat;lon=clon;
    if fmeth=='P'
        % Precipitation []
        %load('/d3/nsteiger/echam5wiso/pr_mon_e5.mat')
        %m_mon=pr_mon;
        %clear pr_mon
        m_mon=pr_crsn2speedy;
    else
        % Soil moisture
        %load('/home/disk/humid4/nathanjs/echamruns/mrso_r_mon.mat')
        %m_mon=mrso_r_mon;
        %clear mrso_r_mon
        error('Missing option...')
    end
    
elseif strcmp(mdl_prior,'cesm_lme009')
    
    % Surface temperature
    ens='009';%disp(['LME #' ens])
    varpth=['/d6/haibo/lme/cam5/BLMTRC5CN/ALLFORCEINGS/' ens];
    Xv=ncread([varpth '/TREFHT/b.e11.BLMTRC5CN.f19_g16.' ens '.cam.h0.TREFHT.085001-184912.nc'],'TREFHT');
    
    tas_mon=permute(Xv,[2 1 3]);
    lat=ncread([varpth '/TREFHT/b.e11.BLMTRC5CN.f19_g16.' ens '.cam.h0.TREFHT.085001-184912.nc'],'lat');
    lon=ncread([varpth '/TREFHT/b.e11.BLMTRC5CN.f19_g16.' ens '.cam.h0.TREFHT.085001-184912.nc'],'lon');
    
    if fmeth=='P'
        % Precipitation [mm/s]; kg/m^2/s - the same as mm/s
        Xv=ncread([varpth '/RAIN/b.e11.BLMTRC5CN.f19_g16.' ens '.clm2.h0.RAIN.085001-184912.nc'],'RAIN');
        m_mon=permute(Xv,[2 1 3]);
    else
        error('Soil moisture not available')
    end
    
elseif strcmp(mdl_prior,'cesm_lme010')
    
    % Surface temperature
    ens='010';%disp(['LME #' ens])
    varpth=['/d6/haibo/lme/cam5/BLMTRC5CN/ALLFORCEINGS/' ens];
    Xv=ncread([varpth '/TREFHT/b.e11.BLMTRC5CN.f19_g16.' ens '.cam.h0.TREFHT.085001-184912.nc'],'TREFHT');
    
    tas_mon=permute(Xv,[2 1 3]);
    lat=ncread([varpth '/TREFHT/b.e11.BLMTRC5CN.f19_g16.' ens '.cam.h0.TREFHT.085001-184912.nc'],'lat');
    lon=ncread([varpth '/TREFHT/b.e11.BLMTRC5CN.f19_g16.' ens '.cam.h0.TREFHT.085001-184912.nc'],'lon');
    
    if fmeth=='P'
        % Precipitation [mm/s]; kg/m^2/s - the same as mm/s
        Xv=ncread([varpth '/RAIN/b.e11.BLMTRC5CN.f19_g16.' ens '.clm2.h0.RAIN.085001-184912.nc'],'RAIN');
        %Xvr=ncread([varpth '/RAIN/b.e11.BLMTRC5CN.f19_g16.' ens '.clm2.h0.RAIN.085001-184912.nc'],'RAIN');
        %Xvs=ncread([varpth '/SNOW/b.e11.BLMTRC5CN.f19_g16.' ens '.clm2.h0.RAIN.085001-184912.nc'],'SNOW');
        %Xv=Xvr+Xvs; % add rain and snow
        m_mon=permute(Xv,[2 1 3]);
    else
        error('Soil moisture not available')
    end
    
elseif strcmp(mdl_prior,'gfdlesm2m_pi')
    
    % Surface temperature
    Xv=ncread('/d1/nsteiger/climate-data/gfdl-esm2m-pi/tas_Amon_GFDL-ESM2M_piControl_r1i1p1_000101-050012.nc','tas');
    
    tas_mon=permute(Xv,[2 1 3]);
    lat=ncread('/d1/nsteiger/climate-data/gfdl-esm2m-pi/tas_Amon_GFDL-ESM2M_piControl_r1i1p1_000101-050012.nc','lat');
    lon=ncread('/d1/nsteiger/climate-data/gfdl-esm2m-pi/tas_Amon_GFDL-ESM2M_piControl_r1i1p1_000101-050012.nc','lon');
    
    if fmeth=='P'
        % Precipitation [mm/s]; kg/m^2/s - the same as mm/s
        Xv=ncread('/d1/nsteiger/climate-data/gfdl-esm2m-pi/pr_Amon_GFDL-ESM2M_piControl_r1i1p1_000101-050012.nc','pr');
        m_mon=permute(Xv,[2 1 3]);
    else
        error('Soil moisture not available')
    end
elseif strcmp(mdl_prior,'gfdlesm2g_pi')
    
    % Surface temperature
    Xv=ncread('/d1/nsteiger/climate-data/gfdl-esm2g-pi/tas_Amon_GFDL-ESM2G_piControl_r1i1p1_000101-050012.nc','tas');
    
    tas_mon=permute(Xv,[2 1 3]);
    lat=ncread('/d1/nsteiger/climate-data/gfdl-esm2g-pi/tas_Amon_GFDL-ESM2G_piControl_r1i1p1_000101-050012.nc','lat');
    lon=ncread('/d1/nsteiger/climate-data/gfdl-esm2g-pi/tas_Amon_GFDL-ESM2G_piControl_r1i1p1_000101-050012.nc','lon');
    
    if fmeth=='P'
        % Precipitation [mm/s]; kg/m^2/s - the same as mm/s
        Xv=ncread('/d1/nsteiger/climate-data/gfdl-esm2g-pi/pr_Amon_GFDL-ESM2G_piControl_r1i1p1_000101-050012.nc','pr');
        m_mon=permute(Xv,[2 1 3]);
    else
        error('Soil moisture not available')
    end    
    
elseif mdl_prior=='a' % icam5 historical
    
    % Surface temperature
    Xv=ncread('/home/disk/humid4/nathanjs/icam/tas_Amon_iCAM_historical_185001-201412.nc','TREFHT');
    tas_mon=flip(permute(Xv,[2 1 3]),1);
    load('/home/disk/humid4/nathanjs/icam/latlon_icam.mat')
    
    if fmeth=='P'
        % Precipitation [m/s]; kg/m^2/s - the same as mm/s
        Xv=ncread('/home/disk/humid4/nathanjs/icam/pr_Amon_iCAM_historical_185001-201412.nc','pr');
        Xv=Xv*1000;% convert to mm/s or kg/(m^2*s)
        m_mon=flip(permute(Xv,[2 1 3]),1);
    else
        error('Soil moisture not available')
    end
    
elseif mdl_prior=='c'
    
    % Surface temperature
    load('/d1/nsteiger/climate-data/ccsm4_last_mil/tas_past1000.mat')
    tas_mon=tas;
    clear tas
    load('/d1/nsteiger/climate-data/ccsm4_last_mil/latlon.mat')
    if fmeth=='P'
        % Precipitation
        load('/d1/nsteiger/climate-data/ccsm4_last_mil/pr_mon_ccsm4_past1000.mat');
        m_mon=pr_mon;
        clear pr_mon
    else
        error('Soil moisture not available')
    end
elseif mdl_prior=='y'
    
    % Surface temperature
    %load('/d1/nsteiger/speedy/SPEEDY_Temp2m_R80_mon_1000_2005.mat');
    load('/d1/nsteiger/climate-data/speedy/speedy_corrected_T_to_echam5wiso.mat');
    tas_mon=permute(s_trans_T(1:end-1,:,:),[2 3 1]);
    p_yrs=p_yrs(1:end-1); disp('SPEEDY DATA MISSING LAST YEAR...')
    load('/d1/nsteiger/climate-data/speedy/latlon_speedy.mat')
    if fmeth=='P'
        % Precipitation
        %load('/d1/nsteiger/speedy/SPEEDY_Precip_R80_mon_1000_2005.mat');
        load('/d1/nsteiger/climate-data/speedy/speedy_corrected_PR_to_echam5wiso.mat');
        %m_mon=Precip;
        m_mon=permute(s_trans_P(1:end-1,:,:),[2 3 1]);
    else
        error('Soil moisture not available')
    end
else
    error('Wrong model type...')
end


% FIND PRIOR MONTHS DATA
p_monf=12*p_yrs;
p_mono=12*(p_yrs-1)+1;

m1=1;m2=12;
p_mons=zeros(12*length(p_yrs),1);
for i=1:length(p_yrs)
    p_mons(m1:m2)=p_mono(i):p_monf(i);
    m1=m1+12;m2=m2+12;
end


% get full size
[rows,cols,numMonths]=size(tas_mon);
yrs=numMonths/12;

%------- Get PAGES2k-based locations of tree rings ----------
%[~,tlat,tlon,LObs]=trw(lat,lon,[1890 2000]);
%------- Get expanded, denser NCDC network --------
% load('/d1/nsteiger/proxy-data/trees_all_NCDC_and_breitenmoser.mat')
% convert lon values so they are like GCM values and keep all proxy
% locations even if they overlap grid points b/c they can have different
% vs-lite growth parameters...
% [~,~,~,tlat,tlon,LObs] = proxy_latlon(lat,lon,plat,plon);

% get matching locations (comment above)
%load('/d1/nsteiger/proxy-data/psm-pseudo-proxies/speedy_trees_pages2k_meanparams_109.mat')
%load('/d1/nsteiger/psm-pseudo-proxies/SPEEDY_trees_pages2k_mparams_101_matchECHAM.mat')
%load('/d1/nsteiger/psm-pseudo-proxies/SPEEDY_trees_pages2k_meanparams_103.mat')
%load('/d1/nsteiger/proxy-data/psm-pseudo-proxies/cesm_lme010_trees_pages2k_m.mat')
%load('/d1/nsteiger/proxy-data/psm-pseudo-proxies/gfdlesm2g_pi_trees_pages2k_m.mat','p_lat','p_lon')


%load('/home/nsteiger/MATLAB/vs-lite/vs_lite_historical_params_T1upper75_calib_1901_2000.mat','p_lat','p_lon','T1','T2','M1','M2');
load('pseudo_vs_lite_historical_params_T1upper75_calib_1901_2000.mat','T1','T2','M1','M2','p_lat','p_lon')

% %disp('Removing tree ring sites for gfdl models...') % comment out below if not gfdl
%rmvtr=[42,157,221,222,261,268]; % for gfdl models
%rmvtr=[38,41,156,220,221,260,267]; % for cesm mean parameters
%T1(rmvtr)=[];T2(rmvtr)=[];M1(rmvtr)=[];M2(rmvtr)=[];
%p_lat(rmvtr)=[];p_lon(rmvtr)=[];

[~,~,~,tlat,tlon,LObs] = proxy_latlon(lat,lon,p_lat,p_lon);
%disp('Using GFDL locations...')

%----- Just take unique locations -------
[lObs,id_u,id_c]=unique(LObs,'stable');
tlat=tlat(id_u);tlon=tlon(id_u);
% Average the tree ring parameters
t1=zeros(length(id_u),1);t2=zeros(length(id_u),1);
m1=zeros(length(id_u),1);m2=zeros(length(id_u),1);
for i=1:length(id_u)
   t1(i)=mean(T1(find(id_c==i)));
   t2(i)=mean(T2(find(id_c==i))); 
   m1(i)=mean(M1(find(id_c==i))); 
   m2(i)=mean(M2(find(id_c==i))); 
end
T1=t1;T2=t2;M1=m1;M2=m2;

% rmvtr=[29]; % for cesm mean parameters
% T1(rmvtr)=[];T2(rmvtr)=[];M1(rmvtr)=[];M2(rmvtr)=[];
% tlat(rmvtr)=[];tlon(rmvtr)=[];lObs(rmvtr)=[];

%disp([num2str(length(tlat)) ' original trees...'])
%----- Keep all proxy locations ------
%lObs=LObs;

% Reshape temp and precip to feed into the tree ring model
tv_mon=zeros(rows*cols,numMonths);
for n=1:numMonths
    tv_mon(:,n)=reshape(tas_mon(:,:,n),rows*cols,1);
end
clear tas_mon
mv_mon=zeros(rows*cols,numMonths);
for n=1:numMonths
    mv_mon(:,n)=reshape(m_mon(:,:,n),rows*cols,1);
end
clear m_mon


% Convert temps from K to deg C and take only proxy locations
tp_mon=tv_mon(lObs,:)-273.15;

% Convert if necessary and take only proxy locations
if fmeth=='P'
    % Convert precipitation to monthly accumulation
    if strcmp(mdl_prior,'y');
        opt=2; % SPEEDY = mm/day
    else
        opt=1; % kg*m-2*s-1
    end
    mp_mon = pcp_conv(mv_mon(lObs,:),1:yrs,opt);
else
    mp_mon=mv_mon(lObs,:);
end

% free up space
clear tv_mon mv_mon

% disp('Using median estimated parameters for vs-lite (from 2013 paper)')
% T1=6*ones(size(T1)); % no growth below this
% T2=18*ones(size(T2)); % not sensitive to temp. above this
% M1=0.035*ones(size(M1)); % no growth below this
% M2=0.8*ones(size(M2)); % not sensitive to moisture above this

% PRECIPITATION/MOISTURE BASED FORWARD MODEL
%[p_d,p_lat,p_lon] = trw_model(...
%    tp_mon,mp_mon,fmeth,p_mons,tlat,tlon,lObs,param);
[p_d,p_lat,p_lon] = trw_model(...
    tp_mon,mp_mon,fmeth,p_mons,tlat,tlon,lObs,T1,T2,M1,M2);

%THROW OUT EXTREME OUTLIER TREE RINGS
maxp=max(abs(p_d'));
outl=find(maxp > 10); % throw out if greater than 5 std units

if ~isempty(outl)
    disp([num2str(length(outl)) ' outlier tree ring(s) removed '])
    save('tr_outliers.mat','outl')
end

%T1(outl)=[];
%T2(outl)=[];
%M1(outl)=[];
%M2(outl)=[];
p_d(outl,:)=[];
p_lat(outl)=[];
p_lon(outl)=[];



%save('output_test9_2.mat','T1','T2','M1','M2','p_d','p_lat','p_lon')

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   FORWARD MODEL FOR TREE RING WIDTH
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [trw_hxb,tlat,tlon,lObs] = trw_model(...
    tp_mon,mp_mon,fmeth,p_mons,tlat,tlon,lObs,T1,T2,M1,M2)
%TRW_MODEL Compute tree ring forward model

obsNum=length(lObs);
pyrs=length(p_mons)/12;

% Take only required years
t_mon_p=tp_mon(:,p_mons);m_mon_p=mp_mon(:,p_mons);

% Reshape calibration data for calibrating TRW ensemble
k=1;l=12;
t_mon_p_r=zeros(obsNum,12,pyrs);
m_mon_p_r=zeros(obsNum,12,pyrs);
for n=1:pyrs
    t_mon_p_r(:,:,n)=t_mon_p(:,k:l);
    m_mon_p_r(:,:,n)=m_mon_p(:,k:l);
    k=k+12; l=l+12;
end

% Compute tree ring widths for H(X_b)
trw_hxb=zeros(obsNum,pyrs);
gT=zeros(obsNum,12,pyrs); % growth from T and M
gM=zeros(obsNum,12,pyrs);%,gT(q,:,:),gM(q,:,:)
for r=1:obsNum
        
    [trw_hxb(r,:),gT(r,:,:),gM(r,:,:)]=VSLite_v2_3(1,pyrs,tlat(r),T1(r),T2(r),M1(r),M2(r),...
        squeeze(t_mon_p_r(r,:,:)),squeeze(m_mon_p_r(r,:,:)),...
        'hydroclim',fmeth);
end

%save('vsl_output.mat','-v7.3')

% Check for NaNs in growth model output, leave those out
tr_nanids=find(isnan(trw_hxb(:,2))); % 2nd year of all

if ~isempty(tr_nanids)
    disp([num2str(length(tr_nanids)) ' tree rings removed '...
        'because temperatures too cold'])
    trw_hxb(tr_nanids,:)=[];
    %disp('Saving removed tree locations...')
    save('tr_nanids.mat','tr_nanids')
end

tlat(tr_nanids)=[];
tlon(tr_nanids)=[];
lObs(tr_nanids)=[];


end






