function Tm = gwmt_d(dlat,dlon,hgt,dmjd, hod, mode)

%
% input data
% ----------
% dlat: latitude in degree scale
% dlon: longitude in degree scale
% hgt: orthometric height in m   scale
% dmjd: modified julian date
% hod: hour of day
% mode: computation mode

% if mode is 'day', dlat, dlon and dmjd are
% scaler, but the hgt and hod are a 1*N matrix;
%
% if mode is 'bl', dlat, dlon and hgt are 3D matrix with the same size,
% dmjd and hour are scaler;
%
% output data
% -----------
% Tm: atmospheric-weighted-mean temperature in GPS atmospheric sounding (K)
% 'hgt'mode: Tm is the size of length(hod)*length(hgt)
%

% check the inputs
if  nargin < 5
    error('MATLAB:gwmt_d:WrongComputation','insufficient inputs');
end

if (nargin<6)
    mode = 'day';
end

%read parameters
load('grid_parameters.mat');
if ~exist('parameters','var') %!=
    error('MATLAB:gwmt_d:WrongComputation','Do not find parameters matrix');
end
hr = [0, 6, 12, 18, 24]';
standard_h= [0 2000 5000 9000]';
% if hgt<=2000
%     k = 1;
% elseif hgt<=5000
%     k = 2;
% elseif hgt<=9000
%     k = 3;
% else
%     error('MATLAB:gwmt_d:WrongComputation','hgt is not between 0 and 9 km');
% end
switch mode
    case 'day'
        %% hgt mode
        [rowd cold] = size(dmjd);
        [rowh colh] = size(hod);
        [rowH colH] = size(hgt);
        if rowd~=1 || cold ~= 1 || rowh*colh~=max(rowh,colh) || colh<1 || rowH~=rowh || colH~=colh%!=
            error('MATLAB:gwmt_d:InputError', 'fail to compute Tm.');
        end
        if colh>1 %&& colh == 1
            hod = hod';
            hgt = hgt';
        end
        len = max(rowh, colh);
        % change the reference epoch to January 1 2010
        dmjd = dmjd-55197;
        Tm_grid = nan(2,2,4,5);
        Tm_hr = nan(len,5);
        Tm = nan(size(hod));
        %annual and semi-annual variations
        cosfy = cos(2*pi*dmjd/365.25);
        sinfy =  sin(2*pi*dmjd/365.25);
        coshy = cos(4*pi*dmjd/365.25);
        sinhy =  sin(4*pi*dmjd/365.25);
        
        i = floor((dlat - 90)/(-5)) + 1;
        if dlat == -90
            i = 36;
        end
        j = floor((dlon + 180)/5) + 1;
        if dlon == 180
            j = 72;
        end
        [L, B] = meshgrid([5*(j-1)-180, 5*j-180], [-5*(i-1)+90, -5*i+90]);
        %UTC 00 06 12 18 24
        for t = 1:5
            %t = 5 for UTC 24 (the next day)
            if t ==5
                cosfy = cos(2*pi*(dmjd+1)/365.25);
                sinfy =  sin(2*pi*(dmjd+1)/365.25);
                coshy = cos(4*pi*(dmjd+1)/365.25);
                sinhy =  sin(4*pi*(dmjd+1)/365.25);
            end
            for k = 1:4%0 2 5 9 km
                Tm_grid(:,:,k, t) = [ parameters{i,j,k,mod(t-1,4)+1}*[ones(size(dmjd')), cosfy',sinfy',coshy',sinhy']',  parameters{i,j+1,k,mod(t-1,4)+1}*[ones(size(dmjd')), cosfy',sinfy',coshy',sinhy']';...
                    parameters{i+1,j,k,mod(t-1,4)+1}*[ones(size(dmjd')), cosfy',sinfy',coshy',sinhy']',  parameters{i+1,j+1,k,mod(t-1,4)+1}*[ones(size(dmjd')), cosfy',sinfy',coshy',sinhy']';];
            end
            %interp
            Tm_interp = [interp2(L,B, Tm_grid(:,:,1,t), dlon, dlat, 'linear'), interp2(L,B, Tm_grid(:,:,2,t), dlon, dlat, 'linear'),...
                interp2(L,B, Tm_grid(:,:,3,t), dlon, dlat, 'linear'), interp2(L,B, Tm_grid(:,:,4,t), dlon, dlat, 'linear')]';
            Tm_hr(:,t) =  interp1(standard_h, Tm_interp, hgt,'linear','extrap');
            
        end
        % UTC 00 06 12 18 24
        for n = 1:len
            Tm(n) = interp1(hr, Tm_hr(n,:), hod(n), 'spline');
        end
        return;
        %% time
    case 'globe'
        if ~isequal(size(dlat), size(dlon), size(hgt)) || ~isscalar(hod) || ~isscalar(dmjd)%hod mjd
            error('MATLAB:gwmt_d_vec:InputError', 'fail to compute Tm.');
        end
        [row, col] = size(dlat);
        % change the reference epoch to January 1 2000
        dmjd = dmjd-mjuliandate(2000,1,1);
        Tm_grid = nan(37,73,4,5);
        Tm_interp = nan(row, col, 4,5);
        Tm_hr = nan(1,5);
        Tm = nan([row, col]);
        
        for i = 1:37
            for j = 1:73
                for k = 1:4
                    Tm_grid(i,j,k,:) = [parameters{i,j,k,1}*[1, cos(2*pi*dmjd/365.25),sin(2*pi*dmjd/365.25),cos(4*pi*dmjd/365.25),sin(4*pi*dmjd/365.25)]',...
                        parameters{i,j,k,2}*[1, cos(2*pi*dmjd/365.25),sin(2*pi*dmjd/365.25),cos(4*pi*dmjd/365.25),sin(4*pi*dmjd/365.25)]',...
                        parameters{i,j,k,3}*[1, cos(2*pi*dmjd/365.25),sin(2*pi*dmjd/365.25),cos(4*pi*dmjd/365.25),sin(4*pi*dmjd/365.25)]',...
                        parameters{i,j,k,4}*[1, cos(2*pi*dmjd/365.25),sin(2*pi*dmjd/365.25),cos(4*pi*dmjd/365.25),sin(4*pi*dmjd/365.25)]',...
                        parameters{i,j,k,1}*[1, cos(2*pi*(dmjd+1)/365.25),sin(2*pi*(dmjd+1)/365.25),cos(4*pi*(dmjd+1)/365.25),sin(4*pi*(dmjd+1)/365.25)]']; %the last one is for UTC 24 (the next day)
                end
            end
        end
        %        toc
        % 2D grid
        resolution = 5;
        [L, B] = meshgrid(-180:resolution:180, 90:-resolution:-90);
        %UTC 00 06 12 18 24
        % horizontal interp
        %tic;
        for layer = 1:4
            for t = 1:5
                %
                Tm_interp(:,:,layer,t) = interp2(L,B, Tm_grid(:,:,layer,t), dlon, dlat, 'linear');
            end
        end
        %toc;
        % vertical intp
        %tic
        for i = 1:row
            for j = 1:col
                Tm_hr = [interp1(standard_h, reshape(Tm_interp(i,j,:,1), [],1), hgt(i,j),'linear','extrap'), ...
                    interp1(standard_h, reshape(Tm_interp(i,j,:,2), [],1), hgt(i,j),'linear','extrap'), ...
                    interp1(standard_h, reshape(Tm_interp(i,j,:,3), [],1), hgt(i,j),'linear','extrap'), ...
                    interp1(standard_h, reshape(Tm_interp(i,j,:,4), [],1), hgt(i,j),'linear','extrap'),...
                    interp1(standard_h, reshape(Tm_interp(i,j,:,5), [],1), hgt(i,j),'linear','extrap')];
                % UTC 00 06 12 18 24
                Tm(i,j) = interp1(hr, Tm_hr, hod, 'spline');
            end
        end
        %toc
        return;
    otherwise
        error('MATLAB:gwmt_d:WrongComputation','unknown mode');
end