/* mrs97.f -- translated by f2c (version 19980913).
   You must link the resulting object file with the libraries:
	-lf2c -lm   (in that order)
*/

#include "f2c.h"

/* Table of constant values */

static integer c__1 = 1;
static doublereal c_b29 = 10.;

/* Subroutine */ int mrs97x_(x, q, mode, upv, dnv, usea, dsea, str, chm, bot, 
	glu)
doublereal *x, *q;
integer *mode;
doublereal *upv, *dnv, *usea, *dsea, *str, *chm, *bot, *glu;
{
    /* Initialized data */

    static doublereal xmin = 1e-5;
    static doublereal xmax = 1.;
    static doublereal qsqmin = 1.25;
    static doublereal qsqmax = 1e7;

    /* Format strings */
    static char fmt_99[] = "(\002  WARNING:  Q^2 VALUE IS OUT OF RANGE   \
\002)";
    static char fmt_98[] = "(\002  WARNING:   X  VALUE IS OUT OF RANGE   \
\002)";

    /* Builtin functions */
    integer s_wsfe(), e_wsfe();

    /* Local variables */
    extern /* Subroutine */ int mrsr1_(), mrsr2_(), mrsr3_(), mrsr4_();
    static doublereal q2;

    /* Fortran I/O blocks */
    static cilist io___6 = { 0, 6, 0, fmt_99, 0 };
    static cilist io___7 = { 0, 6, 0, fmt_98, 0 };


/* ***************************************************************C */
/* 								C */
/*     This is a package for experimenting with E866 data        C */
/*       -- the q**2 range is  1.25d0 < q**2 < 1d7,              C */
/*          the x range is still 1d-5 < x < 1d0                  C */
/* 								C */
/*     As before, x times the parton distribution is returned,   C */
/*     q is the scale in GeV, MSbar factorization is assumed,    C */
/*     and Lambda(MSbar,nf=4) = 300 MeV for all modes            C */
/*                       ie alpha_s(M_Z)=0.1175			C */
/*                                                               C */
/*     Mode =1 gives a set which is like MRS(G) as far as        C */
/*     describing NA51 data and the prelim. E866 data            C */
/*     Mode =2 gives a set gives a better description of         C */
/*     E866 but gives lousy NA51 result.                         C */
/*     Mode =3 gives a set which is even better for E866         C */
/*     but is even worse for NA51                                C */
/*     Mode =4 gives a set which tries to be a compromise        C */
/*     between E866 and NA51 but really fails for both           C */
/*                                                               C */
/*             >>>>>>>>  CROSS CHECK  <<<<<<<<                   C */
/*     The first  datafile B02GRID.DAT (read in on stream 33     C */
/*     The second datafile B03GRID.DAT (read in on stream 34     C */
/*     The third  datafile B04GRID.DAT (read in on stream 31     C */
/*     The fourth datafile B05GRID.DAT (read in on stream 32     C */
/* ***************************************************************C */
    q2 = *q * *q;
    if (q2 < qsqmin || q2 > qsqmax) {
	s_wsfe(&io___6);
	e_wsfe();
    }
    if (*x < xmin || *x > xmax) {
	s_wsfe(&io___7);
	e_wsfe();
    }
    if (*mode == 1) {
	mrsr1_(x, &q2, upv, dnv, usea, dsea, str, chm, bot, glu);
    } else if (*mode == 2) {
	mrsr2_(x, &q2, upv, dnv, usea, dsea, str, chm, bot, glu);
    } else if (*mode == 3) {
	mrsr3_(x, &q2, upv, dnv, usea, dsea, str, chm, bot, glu);
    } else if (*mode == 4) {
	mrsr4_(x, &q2, upv, dnv, usea, dsea, str, chm, bot, glu);
    }
    return 0;
} /* mrs97x_ */

/* Subroutine */ int mrsr1_(x, qsq, upv, dnv, usea, dsea, str, chm, bot, glu)
doublereal *x, *qsq, *upv, *dnv, *usea, *dsea, *str, *chm, *bot, *glu;
{
    /* Initialized data */

    static doublereal xmin = 1e-5;
    static doublereal xmax = 1.;
    static doublereal qsqmin = 1.25;
    static doublereal qsqmax = 1e7;
    static doublereal n0[8] = { 3.,4.,5.,9.,9.,9.,9.,9. };
    static integer init = 0;
    static doublereal xx[49] = { 1e-5,2e-5,4e-5,6e-5,8e-5,1e-4,2e-4,4e-4,6e-4,
	    8e-4,.001,.002,.004,.006,.008,.01,.014,.02,.03,.04,.06,.08,.1,
	    .125,.15,.175,.2,.225,.25,.275,.3,.325,.35,.375,.4,.425,.45,.475,
	    .5,.525,.55,.575,.6,.65,.7,.75,.8,.9,1. };
    static doublereal qq[37] = { 1.25,1.5,2.,2.5,3.2,4.,5.,6.4,8.,10.,12.,18.,
	    26.,40.,64.,100.,160.,240.,400.,640.,1e3,1800.,3200.,5600.,1e4,
	    1.8e4,3.2e4,5.6e4,1e5,1.8e5,3.2e5,5.6e5,1e6,1.8e6,3.2e6,5.6e6,1e7 
	    };

    /* Format strings */
    static char fmt_50[] = "(8f10.5)";

    /* System generated locals */
    doublereal d__1;
    olist o__1;

    /* Builtin functions */
    integer f_open(), s_rsfe(), do_fio(), e_rsfe();
    double pow_dd(), d_lg10();

    /* Local variables */
    static doublereal a, b, f[14896]	/* was [8][49][38] */, g[8];
    static integer i__, j, k, m, n;
    static doublereal xsave, q2save, fac, xxx;

    /* Fortran I/O blocks */
    static cilist io___20 = { 0, 33, 0, fmt_50, 0 };


    xsave = *x;
    q2save = *qsq;
    if (init != 0) {
	goto L10;
    }
    o__1.oerr = 0;
    o__1.ounit = 33;
    o__1.ofnmlen = 11;
    o__1.ofnm = "b02grid.dat";
    o__1.orl = 0;
    o__1.osta = "old";
    o__1.oacc = 0;
    o__1.ofm = 0;
    o__1.oblnk = 0;
    f_open(&o__1);
    for (n = 1; n <= 48; ++n) {
	for (m = 1; m <= 37; ++m) {
	  s_rsfe(&io___20);
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 400], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 399], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 398], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 397], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 396], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 394], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 395], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 393], (ftnlen)sizeof(
		    doublereal));
	    e_rsfe();
/* notation: 1=uval 2=val 3=glue 4=usea 5=chm 6=str 7=btm 8=dsea */
	    for (i__ = 1; i__ <= 8; ++i__) {
/* L25: */
		d__1 = 1. - xx[n - 1];
		f[i__ + (n + m * 49 << 3) - 401] /= pow_dd(&d__1, &n0[i__ - 1]
			);
	    }
/* L20: */
	}
    }
    for (j = 1; j <= 22; ++j) {
	d__1 = xx[j - 1] / xx[22];
	xx[j - 1] = d_lg10(&d__1) + xx[22];
	for (i__ = 1; i__ <= 8; ++i__) {
	    if (i__ == 5 || i__ == 7) {
		goto L31;
	    }
	    for (k = 1; k <= 37; ++k) {
/* L30: */
		d__1 = f[i__ + (j + k * 49 << 3) - 401] / f[i__ + (k * 49 + 
			23 << 3) - 401];
		f[i__ + (j + k * 49 << 3) - 401] = d_lg10(&d__1) + f[i__ + (k 
			* 49 + 23 << 3) - 401];
	    }
L31:
	    ;
	}
    }
    for (i__ = 1; i__ <= 8; ++i__) {
	for (m = 1; m <= 37; ++m) {
/* L40: */
	    f[i__ + (m * 49 + 49 << 3) - 401] = 0.;
	}
    }
    init = 1;
L10:
    if (*x < xmin) {
	*x = xmin;
    }
    if (*x > xmax) {
	*x = xmax;
    }
    if (*qsq < qsqmin) {
	*qsq = qsqmin;
    }
    if (*qsq > qsqmax) {
	*qsq = qsqmax;
    }
    xxx = *x;
    if (*x < xx[22]) {
	d__1 = *x / xx[22];
	xxx = d_lg10(&d__1) + xx[22];
    }
    n = 0;
L70:
    ++n;
    if (xxx > xx[n]) {
	goto L70;
    }
    a = (xxx - xx[n - 1]) / (xx[n] - xx[n - 1]);
    m = 0;
L80:
    ++m;
    if (*qsq > qq[m]) {
	goto L80;
    }
    b = (*qsq - qq[m - 1]) / (qq[m] - qq[m - 1]);
    for (i__ = 1; i__ <= 8; ++i__) {
	g[i__ - 1] = (1. - a) * (1. - b) * f[i__ + (n + m * 49 << 3) - 401] + 
		(1. - a) * b * f[i__ + (n + (m + 1) * 49 << 3) - 401] + a * (
		1. - b) * f[i__ + (n + 1 + m * 49 << 3) - 401] + a * b * f[
		i__ + (n + 1 + (m + 1) * 49 << 3) - 401];
	if (n >= 23) {
	    goto L65;
	}
	if (i__ == 5 || i__ == 7) {
	    goto L65;
	}
	fac = (1. - b) * f[i__ + (m * 49 + 23 << 3) - 401] + b * f[i__ + ((m 
		+ 1) * 49 + 23 << 3) - 401];
	d__1 = g[i__ - 1] - fac;
	g[i__ - 1] = fac * pow_dd(&c_b29, &d__1);
L65:
	d__1 = 1. - *x;
	g[i__ - 1] *= pow_dd(&d__1, &n0[i__ - 1]);
/* L60: */
    }
    *upv = g[0];
    *dnv = g[1];
    *usea = g[3];
    *dsea = g[7];
    *str = g[5];
    *chm = g[4];
    *glu = g[2];
    *bot = g[6];
    *x = xsave;
    *qsq = q2save;
    return 0;
} /* mrsr1_ */

/* Subroutine */ int mrsr2_(x, qsq, upv, dnv, usea, dsea, str, chm, bot, glu)
doublereal *x, *qsq, *upv, *dnv, *usea, *dsea, *str, *chm, *bot, *glu;
{
    /* Initialized data */

    static doublereal xx[49] = { 1e-5,2e-5,4e-5,6e-5,8e-5,1e-4,2e-4,4e-4,6e-4,
	    8e-4,.001,.002,.004,.006,.008,.01,.014,.02,.03,.04,.06,.08,.1,
	    .125,.15,.175,.2,.225,.25,.275,.3,.325,.35,.375,.4,.425,.45,.475,
	    .5,.525,.55,.575,.6,.65,.7,.75,.8,.9,1. };
    static doublereal qq[37] = { 1.25,1.5,2.,2.5,3.2,4.,5.,6.4,8.,10.,12.,18.,
	    26.,40.,64.,100.,160.,240.,400.,640.,1e3,1800.,3200.,5600.,1e4,
	    1.8e4,3.2e4,5.6e4,1e5,1.8e5,3.2e5,5.6e5,1e6,1.8e6,3.2e6,5.6e6,1e7 
	    };
    static doublereal xmin = 1e-5;
    static doublereal xmax = 1.;
    static doublereal qsqmin = 1.25;
    static doublereal qsqmax = 1e7;
    static doublereal n0[8] = { 3.,4.,5.,9.,9.,9.,9.,9. };
    static integer init = 0;

    /* Format strings */
    static char fmt_50[] = "(8f10.5)";

    /* System generated locals */
    doublereal d__1;
    olist o__1;

    /* Builtin functions */
    integer f_open(), s_rsfe(), do_fio(), e_rsfe();
    double pow_dd(), d_lg10();

    /* Local variables */
    static doublereal a, b, f[14896]	/* was [8][49][38] */, g[8];
    static integer i__, j, k, m, n;
    static doublereal xsave, q2save, fac, xxx;

    /* Fortran I/O blocks */
    static cilist io___42 = { 0, 34, 0, fmt_50, 0 };


    xsave = *x;
    q2save = *qsq;
    if (init != 0) {
	goto L10;
    }
    o__1.oerr = 0;
    o__1.ounit = 34;
    o__1.ofnmlen = 11;
    o__1.ofnm = "b03grid.dat";
    o__1.orl = 0;
    o__1.osta = "old";
    o__1.oacc = 0;
    o__1.ofm = 0;
    o__1.oblnk = 0;
    f_open(&o__1);
    for (n = 1; n <= 48; ++n) {
	for (m = 1; m <= 37; ++m) {
	    s_rsfe(&io___42);
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 400], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 399], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 398], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 397], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 396], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 394], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 395], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 393], (ftnlen)sizeof(
		    doublereal));
	    e_rsfe();
/* notation: 1=uval 2=val 3=glue 4=usea 5=chm 6=str 7=btm 8=dsea */
	    for (i__ = 1; i__ <= 8; ++i__) {
/* L25: */
		d__1 = 1. - xx[n - 1];
		f[i__ + (n + m * 49 << 3) - 401] /= pow_dd(&d__1, &n0[i__ - 1]
			);
	    }
/* L20: */
	}
    }
    for (j = 1; j <= 22; ++j) {
	d__1 = xx[j - 1] / xx[22];
	xx[j - 1] = d_lg10(&d__1) + xx[22];
	for (i__ = 1; i__ <= 8; ++i__) {
	    if (i__ == 5 || i__ == 7) {
		goto L31;
	    }
	    for (k = 1; k <= 37; ++k) {
/* L30: */
		d__1 = f[i__ + (j + k * 49 << 3) - 401] / f[i__ + (k * 49 + 
			23 << 3) - 401];
		f[i__ + (j + k * 49 << 3) - 401] = d_lg10(&d__1) + f[i__ + (k 
			* 49 + 23 << 3) - 401];
	    }
L31:
	    ;
	}
    }
    for (i__ = 1; i__ <= 8; ++i__) {
	for (m = 1; m <= 37; ++m) {
/* L40: */
	    f[i__ + (m * 49 + 49 << 3) - 401] = 0.;
	}
    }
    init = 1;
L10:
    if (*x < xmin) {
	*x = xmin;
    }
    if (*x > xmax) {
	*x = xmax;
    }
    if (*qsq < qsqmin) {
	*qsq = qsqmin;
    }
    if (*qsq > qsqmax) {
	*qsq = qsqmax;
    }
    xxx = *x;
    if (*x < xx[22]) {
	d__1 = *x / xx[22];
	xxx = d_lg10(&d__1) + xx[22];
    }
    n = 0;
L70:
    ++n;
    if (xxx > xx[n]) {
	goto L70;
    }
    a = (xxx - xx[n - 1]) / (xx[n] - xx[n - 1]);
    m = 0;
L80:
    ++m;
    if (*qsq > qq[m]) {
	goto L80;
    }
    b = (*qsq - qq[m - 1]) / (qq[m] - qq[m - 1]);
    for (i__ = 1; i__ <= 8; ++i__) {
	g[i__ - 1] = (1. - a) * (1. - b) * f[i__ + (n + m * 49 << 3) - 401] + 
		(1. - a) * b * f[i__ + (n + (m + 1) * 49 << 3) - 401] + a * (
		1. - b) * f[i__ + (n + 1 + m * 49 << 3) - 401] + a * b * f[
		i__ + (n + 1 + (m + 1) * 49 << 3) - 401];
	if (n >= 23) {
	    goto L65;
	}
	if (i__ == 5 || i__ == 7) {
	    goto L65;
	}
	fac = (1. - b) * f[i__ + (m * 49 + 23 << 3) - 401] + b * f[i__ + ((m 
		+ 1) * 49 + 23 << 3) - 401];
	d__1 = g[i__ - 1] - fac;
	g[i__ - 1] = fac * pow_dd(&c_b29, &d__1);
L65:
	d__1 = 1. - *x;
	g[i__ - 1] *= pow_dd(&d__1, &n0[i__ - 1]);
/* L60: */
    }
    *upv = g[0];
    *dnv = g[1];
    *usea = g[3];
    *dsea = g[7];
    *str = g[5];
    *chm = g[4];
    *glu = g[2];
    *bot = g[6];
    *x = xsave;
    *qsq = q2save;
    return 0;
} /* mrsr2_ */

/* Subroutine */ int mrsr3_(x, qsq, upv, dnv, usea, dsea, str, chm, bot, glu)
doublereal *x, *qsq, *upv, *dnv, *usea, *dsea, *str, *chm, *bot, *glu;
{
    /* Initialized data */

    static doublereal xx[49] = { 1e-5,2e-5,4e-5,6e-5,8e-5,1e-4,2e-4,4e-4,6e-4,
	    8e-4,.001,.002,.004,.006,.008,.01,.014,.02,.03,.04,.06,.08,.1,
	    .125,.15,.175,.2,.225,.25,.275,.3,.325,.35,.375,.4,.425,.45,.475,
	    .5,.525,.55,.575,.6,.65,.7,.75,.8,.9,1. };
    static doublereal qq[37] = { 1.25,1.5,2.,2.5,3.2,4.,5.,6.4,8.,10.,12.,18.,
	    26.,40.,64.,100.,160.,240.,400.,640.,1e3,1800.,3200.,5600.,1e4,
	    1.8e4,3.2e4,5.6e4,1e5,1.8e5,3.2e5,5.6e5,1e6,1.8e6,3.2e6,5.6e6,1e7 
	    };
    static doublereal xmin = 1e-5;
    static doublereal xmax = 1.;
    static doublereal qsqmin = 1.25;
    static doublereal qsqmax = 1e7;
    static doublereal n0[8] = { 3.,4.,5.,9.,9.,9.,9.,9. };
    static integer init = 0;

    /* Format strings */
    static char fmt_50[] = "(8f10.5)";

    /* System generated locals */
    doublereal d__1;
    olist o__1;

    /* Builtin functions */
    integer f_open(), s_rsfe(), do_fio(), e_rsfe();
    double pow_dd(), d_lg10();

    /* Local variables */
    static doublereal a, b, f[14896]	/* was [8][49][38] */, g[8];
    static integer i__, j, k, m, n;
    static doublereal xsave, q2save, fac, xxx;

    /* Fortran I/O blocks */
    static cilist io___64 = { 0, 31, 0, fmt_50, 0 };


    xsave = *x;
    q2save = *qsq;
    if (init != 0) {
	goto L10;
    }
    o__1.oerr = 0;
    o__1.ounit = 31;
    o__1.ofnmlen = 11;
    o__1.ofnm = "b04grid.dat";
    o__1.orl = 0;
    o__1.osta = "old";
    o__1.oacc = 0;
    o__1.ofm = 0;
    o__1.oblnk = 0;
    f_open(&o__1);
    for (n = 1; n <= 48; ++n) {
	for (m = 1; m <= 37; ++m) {
	    s_rsfe(&io___64);
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 400], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 399], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 398], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 397], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 396], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 394], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 395], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 393], (ftnlen)sizeof(
		    doublereal));
	    e_rsfe();
/* notation: 1=uval 2=val 3=glue 4=usea 5=chm 6=str 7=btm 8=dsea */
	    for (i__ = 1; i__ <= 8; ++i__) {
/* L25: */
		d__1 = 1. - xx[n - 1];
		f[i__ + (n + m * 49 << 3) - 401] /= pow_dd(&d__1, &n0[i__ - 1]
			);
	    }
/* L20: */
	}
    }
    for (j = 1; j <= 22; ++j) {
	d__1 = xx[j - 1] / xx[22];
	xx[j - 1] = d_lg10(&d__1) + xx[22];
	for (i__ = 1; i__ <= 8; ++i__) {
	    if (i__ == 5 || i__ == 7) {
		goto L31;
	    }
	    for (k = 1; k <= 37; ++k) {
/* L30: */
		d__1 = f[i__ + (j + k * 49 << 3) - 401] / f[i__ + (k * 49 + 
			23 << 3) - 401];
		f[i__ + (j + k * 49 << 3) - 401] = d_lg10(&d__1) + f[i__ + (k 
			* 49 + 23 << 3) - 401];
	    }
L31:
	    ;
	}
    }
    for (i__ = 1; i__ <= 8; ++i__) {
	for (m = 1; m <= 37; ++m) {
/* L40: */
	    f[i__ + (m * 49 + 49 << 3) - 401] = 0.;
	}
    }
    init = 1;
L10:
    if (*x < xmin) {
	*x = xmin;
    }
    if (*x > xmax) {
	*x = xmax;
    }
    if (*qsq < qsqmin) {
	*qsq = qsqmin;
    }
    if (*qsq > qsqmax) {
	*qsq = qsqmax;
    }
    xxx = *x;
    if (*x < xx[22]) {
	d__1 = *x / xx[22];
	xxx = d_lg10(&d__1) + xx[22];
    }
    n = 0;
L70:
    ++n;
    if (xxx > xx[n]) {
	goto L70;
    }
    a = (xxx - xx[n - 1]) / (xx[n] - xx[n - 1]);
    m = 0;
L80:
    ++m;
    if (*qsq > qq[m]) {
	goto L80;
    }
    b = (*qsq - qq[m - 1]) / (qq[m] - qq[m - 1]);
    for (i__ = 1; i__ <= 8; ++i__) {
	g[i__ - 1] = (1. - a) * (1. - b) * f[i__ + (n + m * 49 << 3) - 401] + 
		(1. - a) * b * f[i__ + (n + (m + 1) * 49 << 3) - 401] + a * (
		1. - b) * f[i__ + (n + 1 + m * 49 << 3) - 401] + a * b * f[
		i__ + (n + 1 + (m + 1) * 49 << 3) - 401];
	if (n >= 23) {
	    goto L65;
	}
	if (i__ == 5 || i__ == 7) {
	    goto L65;
	}
	fac = (1. - b) * f[i__ + (m * 49 + 23 << 3) - 401] + b * f[i__ + ((m 
		+ 1) * 49 + 23 << 3) - 401];
	d__1 = g[i__ - 1] - fac;
	g[i__ - 1] = fac * pow_dd(&c_b29, &d__1);
L65:
	d__1 = 1. - *x;
	g[i__ - 1] *= pow_dd(&d__1, &n0[i__ - 1]);
/* L60: */
    }
    *upv = g[0];
    *dnv = g[1];
    *usea = g[3];
    *dsea = g[7];
    *str = g[5];
    *chm = g[4];
    *glu = g[2];
    *bot = g[6];
    *x = xsave;
    *qsq = q2save;
    return 0;
} /* mrsr3_ */

/* Subroutine */ int mrsr4_(x, qsq, upv, dnv, usea, dsea, str, chm, bot, glu)
doublereal *x, *qsq, *upv, *dnv, *usea, *dsea, *str, *chm, *bot, *glu;
{
    /* Initialized data */

    static doublereal xx[49] = { 1e-5,2e-5,4e-5,6e-5,8e-5,1e-4,2e-4,4e-4,6e-4,
	    8e-4,.001,.002,.004,.006,.008,.01,.014,.02,.03,.04,.06,.08,.1,
	    .125,.15,.175,.2,.225,.25,.275,.3,.325,.35,.375,.4,.425,.45,.475,
	    .5,.525,.55,.575,.6,.65,.7,.75,.8,.9,1. };
    static doublereal qq[37] = { 1.25,1.5,2.,2.5,3.2,4.,5.,6.4,8.,10.,12.,18.,
	    26.,40.,64.,100.,160.,240.,400.,640.,1e3,1800.,3200.,5600.,1e4,
	    1.8e4,3.2e4,5.6e4,1e5,1.8e5,3.2e5,5.6e5,1e6,1.8e6,3.2e6,5.6e6,1e7 
	    };
    static doublereal xmin = 1e-5;
    static doublereal xmax = 1.;
    static doublereal qsqmin = 1.25;
    static doublereal qsqmax = 1e7;
    static doublereal n0[8] = { 3.,4.,5.,9.,9.,9.,9.,9. };
    static integer init = 0;

    /* Format strings */
    static char fmt_50[] = "(8f10.5)";

    /* System generated locals */
    doublereal d__1;
    olist o__1;

    /* Builtin functions */
    integer f_open(), s_rsfe(), do_fio(), e_rsfe();
    double pow_dd(), d_lg10();

    /* Local variables */
    static doublereal a, b, f[14896]	/* was [8][49][38] */, g[8];
    static integer i__, j, k, m, n;
    static doublereal xsave, q2save, fac, xxx;

    /* Fortran I/O blocks */
    static cilist io___86 = { 0, 32, 0, fmt_50, 0 };


    xsave = *x;
    q2save = *qsq;
    if (init != 0) {
	goto L10;
    }
    o__1.oerr = 0;
    o__1.ounit = 32;
    o__1.ofnmlen = 11;
    o__1.ofnm = "b05grid.dat";
    o__1.orl = 0;
    o__1.osta = "old";
    o__1.oacc = 0;
    o__1.ofm = 0;
    o__1.oblnk = 0;
    f_open(&o__1);
    for (n = 1; n <= 48; ++n) {
	for (m = 1; m <= 37; ++m) {
	    s_rsfe(&io___86);
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 400], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 399], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 398], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 397], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 396], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 394], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 395], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 49 << 3) - 393], (ftnlen)sizeof(
		    doublereal));
	    e_rsfe();
/* notation: 1=uval 2=val 3=glue 4=usea 5=chm 6=str 7=btm 8=dsea */
	    for (i__ = 1; i__ <= 8; ++i__) {
/* L25: */
		d__1 = 1. - xx[n - 1];
		f[i__ + (n + m * 49 << 3) - 401] /= pow_dd(&d__1, &n0[i__ - 1]
			);
	    }
/* L20: */
	}
    }
    for (j = 1; j <= 22; ++j) {
	d__1 = xx[j - 1] / xx[22];
	xx[j - 1] = d_lg10(&d__1) + xx[22];
	for (i__ = 1; i__ <= 8; ++i__) {
	    if (i__ == 5 || i__ == 7) {
		goto L31;
	    }
	    for (k = 1; k <= 37; ++k) {
/* L30: */
		d__1 = f[i__ + (j + k * 49 << 3) - 401] / f[i__ + (k * 49 + 
			23 << 3) - 401];
		f[i__ + (j + k * 49 << 3) - 401] = d_lg10(&d__1) + f[i__ + (k 
			* 49 + 23 << 3) - 401];
	    }
L31:
	    ;
	}
    }
    for (i__ = 1; i__ <= 8; ++i__) {
	for (m = 1; m <= 37; ++m) {
/* L40: */
	    f[i__ + (m * 49 + 49 << 3) - 401] = 0.;
	}
    }
    init = 1;
L10:
    if (*x < xmin) {
	*x = xmin;
    }
    if (*x > xmax) {
	*x = xmax;
    }
    if (*qsq < qsqmin) {
	*qsq = qsqmin;
    }
    if (*qsq > qsqmax) {
	*qsq = qsqmax;
    }
    xxx = *x;
    if (*x < xx[22]) {
	d__1 = *x / xx[22];
	xxx = d_lg10(&d__1) + xx[22];
    }
    n = 0;
L70:
    ++n;
    if (xxx > xx[n]) {
	goto L70;
    }
    a = (xxx - xx[n - 1]) / (xx[n] - xx[n - 1]);
    m = 0;
L80:
    ++m;
    if (*qsq > qq[m]) {
	goto L80;
    }
    b = (*qsq - qq[m - 1]) / (qq[m] - qq[m - 1]);
    for (i__ = 1; i__ <= 8; ++i__) {
	g[i__ - 1] = (1. - a) * (1. - b) * f[i__ + (n + m * 49 << 3) - 401] + 
		(1. - a) * b * f[i__ + (n + (m + 1) * 49 << 3) - 401] + a * (
		1. - b) * f[i__ + (n + 1 + m * 49 << 3) - 401] + a * b * f[
		i__ + (n + 1 + (m + 1) * 49 << 3) - 401];
	if (n >= 23) {
	    goto L65;
	}
	if (i__ == 5 || i__ == 7) {
	    goto L65;
	}
	fac = (1. - b) * f[i__ + (m * 49 + 23 << 3) - 401] + b * f[i__ + ((m 
		+ 1) * 49 + 23 << 3) - 401];
	d__1 = g[i__ - 1] - fac;
	g[i__ - 1] = fac * pow_dd(&c_b29, &d__1);
L65:
	d__1 = 1. - *x;
	g[i__ - 1] *= pow_dd(&d__1, &n0[i__ - 1]);
/* L60: */
    }
    *upv = g[0];
    *dnv = g[1];
    *usea = g[3];
    *dsea = g[7];
    *str = g[5];
    *chm = g[4];
    *glu = g[2];
    *bot = g[6];
    *x = xsave;
    *qsq = q2save;
    return 0;
} /* mrsr4_ */