function [ y ] = VDRIFT( xt,xl,param )
%VDRIFT CALCULATES EQUATORIAL VERTICAL DRIFT
%****************************************************************************
%
%        subroutine vdrift(xt,xl,param,y)
%-------------------------------------------------------------------
%       SUBROUTINE CALCULATES EQUATORIAL VERTICAL DRIFT AS DESCRIBED 
%       IN SCHERLIESS AND FEJER, JGR, 104, 6829-6842, 1999
%
%       INPUT:   XT: SOLAR LOCAL TIME  [h]
%                XL: GEOGRAPHIC LONGITUDE (+ EAST) [degrees]
%               
%           PARAM: 2-DIM ARRAY (DOY,F10.7CM)
%                  DOY     :Day of Year has to run from 1 to 365(366)
%                  F10.7cm : F10.7cm solar flux (daily value)
%             
%       OUTPUT:   Y: EQUATORIAL VERTICAL DRIFT [m/s]
%
%-------------------------------------------------------------------

%c        IMPLICIT REAL*8 (A-H,O-Z)
%        IMPLICIT REAL (A-H,O-Z)
%       
%c        real*8 param(2),coeff(624),coeff1(594),coeff2(30),funct(6)
%c        real*8 xt,xl,y
%c        real*8 bspl4,bspl4_time,bspl4_long
%
%        real param(2),coeff(624),coeff1(594),coeff2(30),funct(6)
%        real xt,xl,y
%        real bspl4,bspl4_time,bspl4_long
%
%        integer i,j,ind,il,kk
%        integer index_t,dim_t,index_l,dim_l,index,dim,nfunc
  persistent index_t index_l dim nfunc coeff;
  if isempty(index_t)
    index_t = 13;
    %dim_t = 78;
    index_l = 8;
    %dim_l = 48;
    %index = 104;
    dim = 624;
    nfunc = 6;

    coeff1 = [ ...
        -10.80592, -9.63722,-11.52666, -0.05716, -0.06288,  0.03564, ...
         -5.80962, -7.86988, -8.50888, -0.05194, -0.05798, -0.00138, ...
          2.09876,-19.99896, -5.11393, -0.05370, -0.06585,  0.03171, ...
        -10.22653, -3.62499,-14.85924, -0.04023, -0.01190, -0.09656, ...
         -4.85180,-26.26264, -6.20501, -0.05342, -0.05174,  0.02419, ...
        -13.98936,-18.10416, -9.30503, -0.01969, -0.03132, -0.01984, ...
        -18.36633,-24.44898,-16.69001,  0.02033, -0.03414, -0.02062, ...
        -20.27621,-16.95623,-36.58234,  0.01445, -0.02044, -0.08297, ...
          1.44450,  5.53004,  4.55166, -0.02356, -0.04267,  0.05023, ...
          5.50589,  7.05381,  1.94387, -0.03147, -0.03548,  0.01166, ...
          3.24165, 10.05002,  4.26218, -0.03419, -0.02651,  0.07456, ...
          7.02218,  0.06708,-11.31012, -0.03252, -0.01021, -0.09008, ...
         -3.47588, -2.82534, -4.17668, -0.03719, -0.01519,  0.06507, ...
         -4.02607,-11.19563,-10.52923, -0.00592, -0.01286, -0.00477, ...
        -11.47478, -9.57758,-10.36887,  0.04555, -0.02249,  0.00528, ...
        -14.19283,  7.86422, -8.76821,  0.05758, -0.02398, -0.04075, ...
         14.58890, 36.63322, 27.57497,  0.01358, -0.02316,  0.04723, ...
         12.53122, 29.38367, 21.40356, -0.00071, -0.00553,  0.01484, ...
         18.64421, 26.27327, 18.32704,  0.00578,  0.03349,  0.11249, ...
          4.53014,  6.15099,  7.41935, -0.02860, -0.00395, -0.08394, ...
         14.29422,  9.77569,  2.85689, -0.00107,  0.04263,  0.10739, ...
          7.17246,  4.40242, -1.00794,  0.00089,  0.01436,  0.00626, ...
          7.75487,  5.01928,  4.36908,  0.03952, -0.00614,  0.03039, ...
         10.25556,  8.82631, 24.21745,  0.05492, -0.02968,  0.00177, ...
         21.86648, 24.03218, 39.82008,  0.00490, -0.01281, -0.01715, ...
         19.18547, 23.97403, 34.44242,  0.01978,  0.01564, -0.02434, ...
         26.30614, 14.22662, 31.16844,  0.06495,  0.19590,  0.05631, ...
         21.09354, 25.56253, 29.91629, -0.04397, -0.08079, -0.07903, ...
         28.30202, 16.80567, 38.63945,  0.05864,  0.16407,  0.07622, ...
         22.68528, 25.91119, 40.45979, -0.03185, -0.01039, -0.01206, ...
         31.98703, 24.46271, 38.13028, -0.08738, -0.00280,  0.01322, ...
         46.67387, 16.80171, 22.77190, -0.13643, -0.05277, -0.01982, ...
         13.87476, 20.52521,  5.22899,  0.00485, -0.04357,  0.09970, ...
         21.46928, 13.55871, 10.23772, -0.04457,  0.01307,  0.06589, ...
         16.18181, 16.02960,  9.28661, -0.01225,  0.14623, -0.01570, ...
         18.16289, -1.58230, 14.54986, -0.00375, -0.00087,  0.04991, ...
         10.00292, 11.82653,  0.44417, -0.00768,  0.15940, -0.01775, ...
         12.15362,  5.65843, -1.94855, -0.00689,  0.03851,  0.04851, ...
         -1.25167,  9.05439,  0.74164,  0.01065,  0.03153,  0.02433, ...
        -15.46799, 18.23132, 27.45320,  0.00899, -0.00017,  0.03385, ...
          2.70396, -0.87077,  6.11476, -0.00081,  0.05167, -0.08932, ...
          3.21321, -1.06622,  5.43623,  0.01942,  0.05449, -0.03084, ...
         17.79267, -3.44694,  7.10702,  0.04734, -0.00945,  0.11516, ...
          0.46435,  6.78467,  4.27231, -0.02122,  0.10922, -0.03331, ...
         15.31708,  1.70927,  7.99584,  0.07462,  0.07515,  0.08934, ...
          4.19893,  6.01231,  8.04861,  0.04023,  0.14767, -0.04308, ...
          9.97541,  5.99412,  5.93588,  0.06611,  0.12144, -0.02124, ...
         13.02837, 10.29950, -4.86200,  0.04521,  0.10715, -0.05465, ...
          5.26779,  7.09019,  1.76617,  0.09339,  0.22256,  0.09222, ...
          9.17810,  5.27558,  5.45022,  0.14749,  0.11616,  0.10418, ...
          9.26391,  4.19982, 12.66250,  0.11334,  0.02532,  0.18919, ...
         13.18695,  6.06564, 11.87835,  0.26347,  0.02858,  0.14801, ...
         10.08476,  6.14899, 17.62618,  0.09331,  0.08832,  0.28208, ...
         10.75302,  7.09244, 13.90643,  0.09556,  0.16652,  0.22751, ...
          6.70338, 11.97698, 18.51413,  0.15873,  0.18936,  0.15705, ...
          5.68102, 23.81606, 20.65174,  0.19930,  0.15645,  0.08151, ...
         29.61644,  5.49433, 48.90934,  0.70710,  0.40791,  0.26325, ...
         17.11994, 19.65380, 44.88810,  0.45510,  0.41689,  0.22398, ...
          8.45700, 34.54442, 27.25364,  0.40867,  0.37223,  0.22374, ...
         -2.30305, 32.00660, 47.75799,  0.02178,  0.43626,  0.30187, ...
          8.98134, 33.01820, 33.09674,  0.33703,  0.33242,  0.41156, ...
         14.27619, 20.70858, 50.10005,  0.30115,  0.32570,  0.45061, ...
         14.44685, 16.14272, 45.40065,  0.37552,  0.31419,  0.30129, ...
          6.19718, 18.89559, 28.24927,  0.08864,  0.41627,  0.19993, ...
          7.70847, -2.36281,-21.41381,  0.13766,  0.05113, -0.11631, ...
         -9.07236,  3.76797,-20.49962,  0.03343,  0.08630,  0.00188, ...
         -8.58113,  5.06009, -6.23262,  0.04967,  0.03334,  0.24214, ...
        -27.85742,  8.34615,-27.72532, -0.08935,  0.15905, -0.03655, ...
          2.77234,  0.14626, -4.01786,  0.22338, -0.04478,  0.18650, ...
          5.61364, -3.82235,-16.72282,  0.26456, -0.03119, -0.08376, ...
         13.35847, -6.11518,-16.50327,  0.28957, -0.01345, -0.19223, ...
         -5.37290, -0.09562,-27.27889,  0.00266,  0.22823, -0.35585, ...
        -15.29676,-18.36622,-24.62948, -0.31299, -0.23832, -0.08463, ...
        -23.37099,-13.69954,-26.71177, -0.19654, -0.18522, -0.20679, ...
        -26.33762,-15.96657,-42.51953, -0.13575, -0.00329, -0.28355, ...
        -25.42140,-14.14291,-21.91748, -0.20960, -0.19176, -0.32593, ...
        -23.36042,-23.89895,-46.05270, -0.10336,  0.03030, -0.21839, ...
        -19.46259,-21.27918,-32.38143, -0.17673, -0.15484, -0.11226, ...
        -19.06169,-21.13240,-34.01677, -0.25497, -0.16878, -0.11004, ...
        -18.39463,-16.11516,-19.55804, -0.19834, -0.23271, -0.25699, ...
        -19.93482,-17.56433,-18.58818,  0.06508, -0.18075,  0.02796, ...
        -23.64078,-18.77269,-22.77715, -0.02456, -0.12238,  0.02959, ...
        -12.44508,-21.06941,-19.36011,  0.02746, -0.16329,  0.19792, ...
        -26.34187,-19.78854,-24.06651, -0.07299, -0.03082, -0.03535, ...
        -10.71667,-26.04401,-16.59048,  0.02850, -0.09680,  0.15143, ...
        -18.40481,-23.37770,-16.31450, -0.03989, -0.00729, -0.01688, ...
         -9.68886,-20.59304,-18.46657,  0.01092, -0.07901,  0.03422, ...
         -0.06685,-19.24590,-29.35494,  0.12265, -0.24792,  0.05978, ...
        -15.32341, -9.07320,-13.76101, -0.17018, -0.15122, -0.06144, ...
        -14.68939,-14.82251,-13.65846, -0.11173, -0.14410, -0.07133, ...
        -18.38628,-18.94631,-19.00893, -0.08062, -0.14481, -0.12949, ...
        -16.15328,-17.40999,-14.08705, -0.08485, -0.06896, -0.11583, ...
        -14.50295,-16.91671,-25.25793, -0.06814, -0.13727, -0.12213, ...
        -10.92188,-14.10852,-24.43877, -0.09375, -0.11638, -0.09053, ...
        -11.64716,-14.92020,-19.99063, -0.14792, -0.08681, -0.12085, ...
        -24.09766,-16.14519, -8.05683, -0.24065, -0.05877, -0.23726, ...
        -25.18396,-15.02034,-15.50531, -0.12236, -0.09610, -0.00529, ...
        -15.27905,-19.36708,-12.94046, -0.08571, -0.09560, -0.03544, ...
         -7.48927,-16.00753,-13.02842, -0.07862, -0.10110, -0.05807];
    coeff2 = [ ...
        -13.06383,-27.98698,-18.80004, -0.05875, -0.03737, -0.11214, ...
        -13.67370,-16.44925,-16.12632, -0.07228, -0.09322, -0.05652, ...
        -22.61245,-21.24717,-18.09933, -0.05197, -0.07477, -0.05235, ...
        -27.09189,-21.85181,-20.34676, -0.05123, -0.05683, -0.07214, ...
        -27.09561,-22.76383,-25.41151, -0.10272, -0.02058, -0.16720];
    coeff = zeros(dim,1);
    for i=1:dim-length(coeff2)
      coeff(i)=coeff1(i);
    end
    for i=1:length(coeff2)
      coeff(i+dim-length(coeff2))=coeff2(i);
    end
  end

  funct = IRI2012.G(param,xl);

  y=0.;

  for i=1:index_t
    for il=1:index_l
      kk=index_l*(i-1)+il;
      for j=1:nfunc
         ind=nfunc*(kk-1)+j;
         bspl4=IRI2012.BSPL4_TIME(i,xt)*IRI2012.BSPL4_LONG(il,xl);
         y=y+bspl4*funct(j)*coeff(ind);
      end
    end
  end

end

