function [X_obs,xlat,xlon,x_yrs] = load_obs_S(obs_dtp,mon_avg_o,mon_avg_f)
%LOAD_S  Load observational data for statistical PSMs
%   [X,xlat,xlon,id_X] = load_obs_S(model,state_tp,mon_avg_o,mon_avg_f,s_yrs)
%     Given a model code 'obs_dtp', state vector variable type 'state_tp',
%     start and end months for averaging, and the years from the model
%     'm_yrs', this function will load the data saved to the computer
%     'humid'. The outputs are the state vector 'X', the latitudes and
%     longitudes of the gridded data, and the locations of the state
%     variables 'id_X' which are in the same order as given in 'state_tp'.
%     The state variables are annually averaged according to the start
%     month and end moth inputs, unless these are empty in which case the
%     original monthly data is returned.
%     The possible data source codes are:
%     'cru_ts3_tmp'
%     'hadcrut4'
%     'bearth'
%     'mlost'
%     'cru_ts3_pre'
%     'soda_sst'
%     'soda_sss'
%     'my_spei'
%     
%
%   Nathan Steiger, Sept 2016


switch obs_dtp
    
        
    case 'bearth'
        
        
        tas=ncread('./Land_and_Ocean_LatLong1_1850_2016.nc','temperature');
        tas=tas(:,:,1:1992); % just take from 1850 to 2015
        %disp('Using un-infilled B. Earth temperature data...')
        %tas=ncread('./input-recon/b_earth_landocean_mon_1880_2014_infilled_RegEMiridge.nc','tas');
        %disp('Using infilled B. Earth temperature data...')
        tas=permute(tas,[2 1 3]);
        %lon = double(ncread('/home/disk/humid2/nathanjs/berkeley-earth/Land_and_Ocean_LatLong1.nc','longitude'));
        %lat = double(ncread('/home/disk/humid2/nathanjs/berkeley-earth/Land_and_Ocean_LatLong1.nc','latitude'));
        lon = double(ncread('./Land_and_Ocean_LatLong1_1850_2016.nc','longitude'));
        lat = double(ncread('./Land_and_Ocean_LatLong1_1850_2016.nc','latitude'));
        [tas,lon]=regrid2gcm(tas,lon);
        xlat=lat(:);xlon=lon(:);
        
        tas_inf=mon2ann(tas,mon_avg_o,mon_avg_f);
        tas_inf(tas_inf==0)=NaN; % make fill values nans
        
        %x_yrs=1880:(1880+size(tas_inf,3)-1); % for infilled data
        x_yrs=1850:(1850+size(tas_inf,3)-1);
        
        X_obs=reshape(tas_inf,size(tas_inf,1)*size(tas_inf,2),size(tas_inf,3));
        % Remove the time mean to be consistent with prior
        %X_obs=X_obs-repmat(nanmean(X_obs,2),[1,size(X_obs,2)]);
        
        
        
        
end

