function [ B, context ] = SHAMDB0D( context, RLAT,FLON,T,RZ )
%SHAMDB0D COMPUTES THE HOURLY VALUES OF B0 FROM A SET OF SH COEFFICIENTS
%
%	SUBROUTINE  SHAMDB0D (RLAT,FLON,T,RZ,B)
%-------------------------------------------------------------------
%	COMPUTES THE HOURLY VALUES OF B0 FROM A SET OF SH COEFFICIENTS
%	IN A POINT OF A GIVEN GEOCENTRIC LATITUDE AND LONGITUDE
%	OF THE EARTH'S SURFACE FOR A GIVEN MONTH AND A GIVEN SUSPOT NUMER
%
% INPUT:	RLAT    The geogrphic latitude on the meridian given by 
%					the local time (FLON), where the modified dip
%                   latitude is the same as of the orginal site.
%			FLON	=LONGITUDE+15.*UT(hours)
%			T		Month as a REAL number (1.0 to 12.0)
%			RZ		12-month running mean
% OUTOUT	B		=B0
%
%  Blanch E., D. Arrazola, D. Altadill, D. Buresova, M. Mosert, 
%     Adv. Space Res. 39, 701-710, 2007.
%  Altadill, D., D. Arrazola, E. Blanch, D. Buresova, 
%     Adv. Space Res. 42, 610-616, 2008.
%  Altadill, D., J.M. Torta, and E. Blanch, 
%     Adv. Space Res. 43,1825-1834, 2009.
%-------------------------------------------------------------------
  persistent IBO JBO KDIM LDIM L RE TZERO IFIT IB KINT LINT ...
    KEXT LEXT IE GANM GBNM HANM HBNM KMAX KT;
  if isempty(IBO)
    IBO=0;
    JBO=1;
    KDIM=6;
    LDIM=4;
    L=-1;
%      DIMENSION   FN(0:KDIM,0:KDIM), CONST(0:KDIM,0:KDIM)
%      DIMENSION   GNM(0:KDIM,0:KDIM,1-IBO-JBO:LDIM),
%     *			HNM(0:KDIM,0:KDIM,1-IBO-JBO:LDIM)
%	  DIMENSION   GANM(0:KDIM,0:KDIM,1-IBO-JBO:LDIM),
%     *			HANM(0:KDIM,0:KDIM,1-IBO-JBO:LDIM),
%     *            GBNM(0:KDIM,0:KDIM,1-IBO-JBO:LDIM),
%     *			HBNM(0:KDIM,0:KDIM,1-IBO-JBO:LDIM)
%      DIMENSION   BINT(0:KDIM,0:KDIM,1-IBO-JBO:LDIM),
%     *            BEXT(0:KDIM,0:KDIM,1-IBO-JBO:LDIM)
%      CHARACTER*1 IE       
%      COMMON  BINT,BEXT,RE,TZERO,IFIT,IB,KINT,LINT,KEXT,
%     *              LEXT,KMAX,FN
%c     ,CONST

    %THETA = 180.;
    %RE = 6371.2;
    RE=WGS84.equatorialRadius/1000.0; % km
    TZERO = 1.0;
    IFIT = -1;
    %ICEN = 0;
    %IREF = 0;
    IB = 2;
    KINT = 6;
    LINT = 4;
    KEXT = 0;
    LEXT = -1;
%     *     /180.,6371.2,1.0,  -1,  0,   0,   2, 6,   4,   0,   -1/
%   	cnst = ...
%       	 [1.,1.,1.,1.,1.73205,0.866025,1.,2.44949,1.93649,0.7905691,1., ...
%       	  3.16228,3.35410,2.09165,0.739510,1.,3.87298,5.12348, ...
%       	  4.18330,2.21853,0.701561,1.,4.58258,7.24569,7.24569,4.96078, ...
%             2.32681,0.671693];
%     CONST = zeros(KDIM+1,KDIM+1);
%     K = 1;
%     for N=0:KDIM
%       for M=0:N
%         CONST(N+1,M+1) = cnst(K);
%         K = K + 1;
%       end
%     end

    IE = 'I';

    dat = ...
       [176.746, 0.233,   0.000, 0.000, ...
       -109.413, 0.072,   0.000, 0.000, -66.000,-0.004,   0.000, 0.000, ...
         36.874,-0.018,   0.000, 0.000, -19.515, 0.040,   0.000, 0.000, ...
         94.998, 0.724,   0.000, 0.000, ...
       -116.900,-0.971,   0.000, 0.000, -93.849,-0.590,   0.000, 0.000, ...
         80.579, 0.425,   0.000, 0.000, -19.205,-0.220,   0.000, 0.000, ...
        -94.824, 0.115,   6.055, 0.265, ...
         84.720,-0.161,  -7.101,-0.374,  35.200,-0.138,   1.043,-0.350, ...
        -23.960, 0.109,  -2.478, 0.133,  25.550,-0.049,  -3.143, 0.003, ...
        -29.418,-0.823,   0.000, 0.000, ...
         40.070, 0.800,   0.000, 0.000,  24.019, 0.478,   0.000, 0.000, ...
        -13.175,-0.265,   0.000, 0.000,   8.799, 0.090,   0.000, 0.000, ...
        -31.099,-1.117,  -1.906, 0.498, ...
         43.807, 1.406,  -3.216,-0.520,  37.957, 0.902,  -0.717,-0.570, ...
        -40.232,-0.583,  12.171, 0.244,  11.595, 0.241,  -4.890, 0.054, ...
        -87.665, 1.635, 117.581,-1.765, ...
        123.444,-2.119,-146.917, 2.131,  81.137,-1.173, -99.063, 1.548, ...
        -42.646, 0.681,  61.263,-0.811,  17.550,-0.408, -24.374, 0.260, ...
         54.538,-0.170,   0.000, 0.000, ...
        -71.552, 0.361,   0.000, 0.000, -50.565,-0.077,   0.000, 0.000, ...
         36.653,-0.071,   0.000, 0.000, -10.816, 0.236,   0.000, 0.000, ...
        -31.138, 1.156,  37.307,-1.407, ...
         40.390,-1.390, -34.573, 1.730,  41.597,-0.835, -41.318, 1.550, ...
        -19.779, 0.404,  15.954,-0.696,  -1.706,-0.220,   5.084, 0.040, ...
        -57.671, 0.045,  42.683,-0.800, ...
         71.491, 0.048, -49.441, 0.980,  47.893,-0.037, -36.191, 0.562, ...
        -26.638,-0.029,  20.346,-0.384,   9.998, 0.067,  -6.787, 0.213, ...
         90.187,-1.198, -66.032,-0.056, ...
       -119.148, 1.428,  81.202, 0.022, -63.375, 0.754,  53.070, 0.149, ...
         39.249,-0.436, -30.898,-0.052, -27.293, 0.301,  12.838,-0.067, ...
       -110.261, 1.509,   0.000, 0.000, ...
        164.956,-1.761,   0.000, 0.000, 103.699,-1.005,   0.000, 0.000, ...
        -55.127, 0.569,   0.000, 0.000,  25.376,-0.315,   0.000, 0.000, ...
       -104.655, 1.341, 109.057,-1.367, ...
        139.129,-1.730,-127.325, 1.532,  88.526,-1.068,-106.461, 1.397, ...
        -38.306, 0.508,  56.240,-0.798,  17.239,-0.267,  -7.766, 0.058, ...
         -6.494,-1.253,   5.714, 0.132, ...
          3.547, 1.545,  -5.372,-0.106,  -4.343, 1.103,  -3.393,-0.017, ...
          2.454,-0.626,  -3.297,-0.025,   5.871, 0.160,   2.040,-0.036, ...
         50.814,-0.230, -25.094, 0.817, ...
        -65.502, 0.304,  32.267,-1.075, -44.176, 0.019,  14.606,-0.605, ...
         27.869,-0.009,  -5.147, 0.387, -11.041, 0.131,   5.922,-0.225, ...
         77.825,-0.728, 128.501,-0.810, ...
        -87.685, 0.838,-164.016, 1.103, -74.431, 0.807, -95.539, 0.498, ...
         40.631,-0.454,  49.950,-0.292,  -4.229, 0.000, -29.666, 0.272, ...
        152.380,-1.232,   0.000, 0.000, ...
       -192.098, 1.514,   0.000, 0.000,-132.417, 1.370,   0.000, 0.000, ...
         82.894,-0.709,   0.000, 0.000, -28.162, 0.050,   0.000, 0.000, ...
        -12.633, 1.192,  47.246,-1.193, ...
         -5.488,-1.387, -67.206, 1.486,  -9.917,-0.914, -34.438, 0.552, ...
         13.185, 0.477,  21.225,-0.387,   0.586,-0.208, -15.426, 0.419, ...
         -4.478,-0.118,  17.908, 0.175, ...
         -0.417, 0.067, -27.047,-0.241,   7.636, 0.028, -10.075,-0.109, ...
        -10.582, 0.005,  14.496, 0.086,   0.421, 0.001, -12.200,-0.041, ...
         16.086, 0.321,  47.044,-0.126, ...
        -24.823,-0.280, -62.615, 0.210, -12.030,-0.136, -44.003,-0.023, ...
          4.929, 0.137,  28.340,-0.009,  -4.688,-0.057,  -9.315, 0.103, ...
         28.023,-0.031, -21.535, 0.115, ...
        -31.946, 0.011,  24.143,-0.180, -21.019,-0.057,  24.108,-0.116, ...
         13.969, 0.004, -13.823, 0.042,  -6.860, 0.031,   0.546,-0.035, ...
       0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000, ...
       0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000, ...
        -31.994, 0.409,   0.000, 0.000, ...
         18.217,-0.458,   0.000, 0.000,  39.280,-0.754,   0.000, 0.000, ...
        -20.453, 0.324,   0.000, 0.000,  -8.111, 0.139,   0.000, 0.000, ...
         28.765,-0.477, -28.368, 0.516, ...
        -50.604, 0.751,  25.725,-0.471, -23.444, 0.283,  29.966,-0.558, ...
          2.759,-0.146, -10.824, 0.341,  -7.419, 0.206,  -3.711, 0.056, ...
         42.429,-0.415,   1.993, 0.117, ...
        -53.162, 0.555,   7.229,-0.246, -19.307, 0.039,  -8.028, 0.028, ...
          9.849,-0.035,   6.834, 0.033, -17.010, 0.272,   4.668,-0.129, ...
          4.546,-0.359, -57.796, 0.359, ...
          0.738, 0.343,  73.027,-0.423,  -7.421, 0.420,  56.067,-0.327, ...
          5.093,-0.279, -37.581, 0.226,   3.636,-0.041,  10.910,-0.059, ...
         88.440,-0.393, -69.598, 0.643, ...
       -109.481, 0.532,  82.266,-0.765, -59.229, 0.182,  55.279,-0.580, ...
         28.514,-0.057, -30.282, 0.326, -22.924, 0.164,  11.602,-0.073, ...
       0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000, ...
       0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000, ...
       0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000, ...
       0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000];
    %  DATA (((GANM(N,M,J),GBNM(N,M,J),HANM(N,M,J),HBNM(N,M,J),
    % *    	J=0,LDIM), M=0,N), N=0,KDIM)
    GANM = zeros(KDIM+1,KDIM+1,LDIM-(1-IBO-JBO)+1);
    GBNM = zeros(KDIM+1,KDIM+1,LDIM-(1-IBO-JBO)+1);
    HANM = zeros(KDIM+1,KDIM+1,LDIM-(1-IBO-JBO)+1);
    HBNM = zeros(KDIM+1,KDIM+1,LDIM-(1-IBO-JBO)+1);
    K = 1;
    for N=0:KDIM
      for M=0:N
        for J=1-IBO-JBO:LDIM
          GANM(N+1,M+1,J-(1-IBO-JBO)+1) = dat(K);
          %GANM(M+1,N+1,J-(1-IBO-JBO)+1) = dat(K);
          GBNM(N+1,M+1,J-(1-IBO-JBO)+1) = dat(K+1);
          %GBNM(M+1,N+1,J-(1-IBO-JBO)+1) = dat(K+1);
          HANM(N+1,M+1,J-(1-IBO-JBO)+1) = dat(K+2);
          %HANM(M+1,N+1,J-(1-IBO-JBO)+1) = dat(K+2);
          HBNM(N+1,M+1,J-(1-IBO-JBO)+1) = dat(K+3);
          %HBNM(M+1,N+1,J-(1-IBO-JBO)+1) = dat(K+3);
          K = K + 4;
        end
      end
    end

    KMAX = max(KINT,KEXT);
    if KMAX > KDIM
      % stop
      fprintf(1,'ERROR - bad KINT(%d) or KEXT(%d) in SHAMDB0D (%d > %d)\n',KINT,KEXT,KMAX,KDIM);
      return;
    end
    KT = max(LINT,LEXT);
    if KT > LDIM
      % stop
      fprintf(1,'ERROR - bad LINT(%d) or LEXT(%d) in SHAMDB0D (%d > %d)\n',LINT,LEXT,KT,LDIM);
      return;
    end
  end
  context.RE = RE;
  context.TZERO = TZERO;
  context.IFIT = IFIT;
  context.IB = IB;
  context.KINT = KINT;
  context.LINT = LINT;
  context.KEXT = KEXT;
  context.LEXT = LEXT;
  context.KMAX = KMAX;
  context.FN = zeros(KDIM+1,KDIM+1);
  context.BINT = zeros(KDIM+1,KDIM+1,LDIM-(1-IBO-JBO)+1);
  context.BEXT = zeros(KDIM+1,KDIM+1,LDIM-(1-IBO-JBO)+1);
  GNM = zeros(KDIM+1,KDIM+1,LDIM-(1-IBO-JBO)+1);
  HNM = zeros(KDIM+1,KDIM+1,LDIM-(1-IBO-JBO)+1);

	for N=0:KMAX
    for M=0:N

      for J=1-IBO-JBO:KT
        GNM(N+1,M+1,J-(1-IBO-JBO)+1)=GANM(N+1,M+1,J-(1-IBO-JBO)+1) ...
                                    +GBNM(N+1,M+1,J-(1-IBO-JBO)+1)*RZ;
        HNM(N+1,M+1,J-(1-IBO-JBO)+1)=HANM(N+1,M+1,J-(1-IBO-JBO)+1) ...
                                    +HBNM(N+1,M+1,J-(1-IBO-JBO)+1)*RZ;
      end

      if IE == 'I'
        if N > KINT
          continue;
        end
        LJ = LINT;
      else
        if N > KEXT
          continue;
        end
        LJ = LEXT;
      end

      context.FN(N+1,M+1) = N;

      if M == 0
        for J=1-IBO-JBO:KT
          if IE == 'I'
            context.BINT(N+1,M+1,J-(1-IBO-JBO)+1)   = GNM(N+1,M+1,J-(1-IBO-JBO)+1);
          else
            context.BEXT(N+1,M+1,J-(1-IBO-JBO)+1)   = GNM(N+1,M+1,J-(1-IBO-JBO)+1);
          end
        end
      else
        for J=1-IBO-JBO:LJ
          if IE == 'I'
            context.BINT(N+1,M+1,J-(1-IBO-JBO)+1)   = GNM(N+1,M+1,J-(1-IBO-JBO)+1);
            context.BINT(M-1+1,N+1,J-(1-IBO-JBO)+1) = HNM(N+1,M+1,J-(1-IBO-JBO)+1);
          else
            context.BEXT(N+1,M+1,J-(1-IBO-JBO)+1)   = GNM(N+1,M+1,J-(1-IBO-JBO)+1);
            context.BEXT(M-1+1,N+1,J-(1-IBO-JBO)+1) = HNM(N+1,M+1,J-(1-IBO-JBO)+1);
          end
        end
      end
    end
  end
  %     **********************************************************
  %     SYNTHESIZES THE VALUE OF B0 FROM THE MODEL	
  %     **********************************************************
  [~,~,B,context] = SCHNEVPD(context,RZ,RLAT,FLON,0,T,L);

end

