function [Pt,data] = LT3d(x,y,zs,U,V,varargin)
%
%  This is a 3-D precipitation solution for the linear orographic precipitation model (Smith and Barstad , 2004).
%
%  Usage:
%
%  >> P = LT3d(x,y,z,U,V)
% 
%  >> [P,data] = LT3d(x,y,z,param1,value1,param2,value2, ... )
%
% Where:
%
%  x = 2d array of x coordinates in meters
%  x = 2d array of y coordinates in meters
%  z = 2d array of terrain heights 
%  U = Wind velocity in X direction
%  V = Wind velocity in Y direction
%  param1, value1 = model parameter/value pairs (see below)
%
%  P     = the simulated precipitation (mm/hr)
%  data  = a data structure containing several different model outputs
%
% Known model parameters, and their default settings, are as follows:
%
% 'z0'    = 0;            % at the height level where you want the precip
% 'Ndsq'  = 3.6e-5;       % dry BV freq sq.
% 'hw'    = 2500;         % scale of water vapor (see formula in SB04,appendix)
% 'Vt'    = 1.0;          % Integrated fall speed in m/s (=1 original code)
% 'cw'    = 1.1;          % sensitivity paper (set to 1.1 in here, see paper SB04)
% 'cwqv'  = cw*0.01;      % sensitivity times q_vs (set to: 0.01 kg/kg, see paper SB04)
% 'tauf'  = 1000;         % BS'11: =zg/Vt ( =hw/Vt for moisture level: around 500s (->750s) is normally good)
% 'tauc'  = 1000;         % probably around 500s, normally shorter for warmer condition
% 'Pbg'   = 0;            % Background precip (if not sure, use =0) add with units [mm/hr]
%
% Some LM papers :SB'04  :: Smith and Barstad (2004); JAS
%                 SE'07  :: Smith and Evans (2007); JHM
%                 BGS'07 :: Barstad et al. 2007; J Hydrol.
%                 BS'11  :: Barstad and Schuller (2011) Nov, JAS
%
%  Example 1: Gaussian mountain (theoretical)
%
%     hm = 500;         % mnt height
%     xm = 300.0e3;     % mountain location in domain 
%     ym = 300.0e3;
%     am = 15000;       % mountain half-width (m)
%     dx = 500;         % grid spacing
%
%     [x,y] = meshgrid( 0:dx:(xm*2), 0:dx:(xm*2) );
%     zs = hm*exp(-( ((x-xm)./am).^2 + ((y-ym)./am).^2 ));  % Gaussian
%     P  = LT3d(x,y,zs,10,0);
%
% Example 2: Realworld example with terrain data for the Nauelbuta mountains
%
%     zz = srtmread([-74 -72],[-38.5 -36.5]);
%     [xx,yy] = meshgrid( (0:(size(zz,2)-1))*100 , (0:(size(zz,2)-1))*100 );
%
%     P = LT3d(xx,yy,zz,+10,-5,'hw',3000);
%
%     clf;imagesc(P);hold on;contour(zz,[1,500,1000,1500],'w');caxis([0 3]);colorbar;set(gca,'ydir','normal');
%
%

%---------------------------------------------------------------------------------
% S-term paramters (see how to calculate them in SB'04):
%
% The defaut settings try to reproduce the fresult shoen in figure 4c of SB'04
%
z0    = 0;            % at the height level where you want the precip
Ndsq  = 3.6e-5;       % dry BV freq sq.
hw    = 2500;         % scale of water vapor (see formula in SB04,appendix)
Vt    = 1.0;          % Fall speed
cw    = 1.1;          % sensitivity paper (set to 1.1 in here, see paper SB04)
cwqv  = cw*0.01;      % sensitivity times q_vs (set to: 0.01 kg/kg, see paper SB04)
tauf  = [];           % BS'11: =zg/Vt ( =hw/Vt for moisture level: around 500s (->750s) is normally good)
tauc  = 1000;         % probably around 500s, normally shorter for warmer condition
Pbg   = 0;            % Background precip (if not sure, use =0) add with units [mm/hr]
dx    = [];
dy    = [];
%---------------------------------------------------------------------------------

% Variable settings
for iv = 1:2:length(varargin)
	        eval([varargin{iv},' = varargin{iv+1};']);
end

if isempty(dx)
	dx    = x(1,2)-x(1,1);
	dy    = dx;           % Must have same spacing in x and y
end
Nx    = size(zs,2);   % no of grid pts in X
Ny    = size(zs,1);   % no of grid pts in Y
Lx    = max(x);
Ly    = max(y);

% Compute solution in Fourier space
% -first the input to Precip-equation

Fzs  = fftshift(fft2(zs))/Nx;  % For now the input must be square!!
[k,l] = meshgrid( linspace(-pi/dx,pi/dx,Nx), linspace(-pi/dy,pi/dy,Ny) );
sig   = U*k + V*l;

% nonhydrostatic version, needed for smaller-scale mountains(?)
%
% discrim=Ndsq/U^2-k.^2;
%  if (discrim(ik) > 0)
%    m = sign(sig).*sqrt(discrim);
%  else
%    m = sqrt(discrim);                      % Check before use
%  end

% Get S and P terms
if isempty(tauf)
    taufN  = hw/(Vt);
else
    taufN = tauf;
end
m      = sign(sig) .* sqrt( (Ndsq./sig.^2).*(k.^2 + l.^2) );
%*********************************************************************************
m(isnan(m)) = 0;  % --> critical case when sig is 0.  In this case we set m to 
                  %     zero but there might be a better way to handle this...
%********************************************************************************
FSterm = i.*cwqv.*Fzs.*exp(i.*m.*z0).*sig.*exp(-z0.*(1-i.*m.*hw)/hw) ./ (1-i.*m.*hw);   % simple, hydrostatic solution
FSraw  = i.*cwqv.*Fzs.*sig.*exp(-z0./hw);                                               % the raw upslope (~Smith'79)
FPterm = FSterm./((1+i.*sig.*tauc).*(1+i.*sig.*taufN));                                 % this also include time delays and dynamics
w      = i.*sig.*Fzs.*exp(i.*m.*1000);

% bring it back from Fourier space:
S      = Nx*ifft2(ifftshift(FSterm));
Sraw   = Nx*ifft2(ifftshift(FSraw));
P      = Nx*ifft2(ifftshift(FPterm));
W      = Nx*ifft2(ifftshift(w));

% add background precip :
Pt = P+Pbg/3600;                                % to get it in [kg/(m^2 s)]

% trucate negative values as unphysical, and return precipitation in mm/hr
Pt = max( 0, real(Pt)*3600);

% Full output dataset
data.x  = x;
data.y  = y;
data.zs = zs;
data.Pt = Pt;
data.P  = max( 0, real(P));
data.W  = W;
data.S  = S;
data.Sraw = Sraw;

