function [Mbc] = bias_correct(Od,Md,fit)
%BIAS_CORRECT Bias corrects the model input according to the given
% observational data. Bias correction options are for normal distributions
% (useful for variables like temperature) and two-parameter gamma 
% distributions (useful for variables like precipitation). This function
% assumes that the two fields have the same spatial dimension but not time
% dimension.
%   [Mbc] = bias_correct(Od,Md,fit)
%   Where 
%       Od[space,time] = observational data
%       Md[space,time] = model data (to be corrected)
%       fit = type of distribution fit, 'n' = normal, 'g' = gamma
%
%       Mbc[space,time] = model bias-corrected data
%
% Nathan Steiger, Sept 2016, LDEO

if ndims(Od)>2 || ndims(Md)>2
    error('Only use one spatial dimension...')
end


Mbc=NaN(size(Md));

if fit == 'n'
    
    for i=1:size(Mbc,1)
        
        % Get distribution parameters
        [o_mu,o_sig]=normfit(Od(i,:));
        [m_mu,m_sig]=normfit(Md(i,:));
        
        % Convert model data to CDF
        Mcdf=normcdf(Md(i,:),m_mu,m_sig);
        
        % Invert CDF to new corrected distribution
        Mbc(i,:)=norminv(Mcdf,o_mu,o_sig);
        
    end
    
elseif fit == 'g'
    
    for i=1:size(Mbc,1)
        
        % Get distribution parameters
        o_prms=gamfit(Od(i,:));
        m_prms=gamfit(Md(i,:));
        
        % Convert model data to CDF
        Mcdf=gamcdf(Md(i,:),m_prms(1),m_prms(2));
        
        % Invert CDF to new corrected distribution
        Mbc(i,:)=gaminv(Mcdf,o_prms(1),o_prms(2));
        
    end
    
else
    error('Wrong fit distribution option...')
end


end

