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

#include "f2c.h"

/* Common Block Declarations */

struct {
    doublereal x[16], y[16];
} coor_;

#define coor_1 coor_

struct {
    doublereal fx[16], fy[16];
} force_;

#define force_1 force_

struct {
    doublereal side, sideh, rcutsq;
} para_;

#define para_1 para_

struct {
    doublereal vx[16], vy[16];
} vel_;

#define vel_1 vel_

struct {
    doublereal pv[100];
} meas1_;

#define meas1_1 meas1_

struct {
    doublereal gr[100];
} meas2_;

#define meas2_1 meas2_

/* Table of constant values */

static integer c__3 = 3;
static integer c__1 = 1;
static integer c__5 = 5;
static integer c__0 = 0;
static integer c__9 = 9;

/*     This program runs velocity verlet in 2d */
/*     Constant Energy Program */
/* Main program */ MAIN__()
{
    /* Format strings */
    static char fmt_1000[] = "(4(1x,f12.6))";
    static char fmt_2000[] = "(6(1x,f12.6))";

    /* System generated locals */
    integer i__1, i__2;
    doublereal d__1;
    olist o__1;

    /* Builtin functions */
    integer f_open(), s_rsle(), do_lio(), e_rsle(), s_wsle(), e_wsle(), 
	    s_wsfe(), do_fio(), e_wsfe();

    /* Local variables */
    extern /* Subroutine */ int fcal_();
    static doublereal ekin, ener, time, grav[100], dist, rmin, temp;
    extern /* Subroutine */ int corr_();
    static doublereal vmax, pvav[100], rcut;
    extern doublereal rsub_();
    static doublereal etot, rmin2;
    static integer i__, j, k, iseed, nside, nmeas, istep;
    static doublereal d1, d2, virav, press;
    static integer nevry;
    static doublereal vxsum, vysum, ep, dt;
    static integer ir;
    static doublereal rr, ekinav, xx, yy, epotav;
    static integer ijk, jkl, neq;
    static doublereal rho, vir;
    extern /* Subroutine */ int kinetic_();

    /* Fortran I/O blocks */
    static cilist io___1 = { 0, 7, 0, 0, 0 };
    static cilist io___25 = { 0, 6, 0, 0, 0 };
    static cilist io___26 = { 0, 9, 0, fmt_1000, 0 };
    static cilist io___33 = { 0, 8, 0, 0, 0 };
    static cilist io___34 = { 0, 8, 0, fmt_1000, 0 };
    static cilist io___36 = { 0, 8, 0, fmt_1000, 0 };
    static cilist io___46 = { 0, 11, 0, 0, 0 };
    static cilist io___47 = { 0, 11, 0, fmt_2000, 0 };
    static cilist io___48 = { 0, 12, 0, fmt_1000, 0 };
    static cilist io___49 = { 0, 13, 0, fmt_1000, 0 };
    static cilist io___50 = { 0, 10, 0, fmt_1000, 0 };


    o__1.oerr = 0;
    o__1.ounit = 7;
    o__1.ofnmlen = 8;
    o__1.ofnm = "md_e.inp";
    o__1.orl = 0;
    o__1.osta = 0;
    o__1.oacc = 0;
    o__1.ofm = 0;
    o__1.oblnk = 0;
    f_open(&o__1);
    o__1.oerr = 0;
    o__1.ounit = 8;
    o__1.ofnmlen = 8;
    o__1.ofnm = "md_e.out";
    o__1.orl = 0;
    o__1.osta = 0;
    o__1.oacc = 0;
    o__1.ofm = 0;
    o__1.oblnk = 0;
    f_open(&o__1);
    o__1.oerr = 0;
    o__1.ounit = 9;
    o__1.ofnmlen = 8;
    o__1.ofnm = "md_e.ini";
    o__1.orl = 0;
    o__1.osta = 0;
    o__1.oacc = 0;
    o__1.ofm = 0;
    o__1.oblnk = 0;
    f_open(&o__1);
    o__1.oerr = 0;
    o__1.ounit = 10;
    o__1.ofnmlen = 8;
    o__1.ofnm = "md_e.fin";
    o__1.orl = 0;
    o__1.osta = 0;
    o__1.oacc = 0;
    o__1.ofm = 0;
    o__1.oblnk = 0;
    f_open(&o__1);
    o__1.oerr = 0;
    o__1.ounit = 11;
    o__1.ofnmlen = 8;
    o__1.ofnm = "md_e.res";
    o__1.orl = 0;
    o__1.osta = 0;
    o__1.oacc = 0;
    o__1.ofm = 0;
    o__1.oblnk = 0;
    f_open(&o__1);
    o__1.oerr = 0;
    o__1.ounit = 12;
    o__1.ofnmlen = 7;
    o__1.ofnm = "md_e.gr";
    o__1.orl = 0;
    o__1.osta = 0;
    o__1.oacc = 0;
    o__1.ofm = 0;
    o__1.oblnk = 0;
    f_open(&o__1);
    o__1.oerr = 0;
    o__1.ounit = 13;
    o__1.ofnmlen = 7;
    o__1.ofnm = "md_e.pv";
    o__1.orl = 0;
    o__1.osta = 0;
    o__1.oacc = 0;
    o__1.ofm = 0;
    o__1.oblnk = 0;
    f_open(&o__1);
    s_rsle(&io___1);
    do_lio(&c__3, &c__1, (char *)&iseed, (ftnlen)sizeof(integer));
    do_lio(&c__3, &c__1, (char *)&ir, (ftnlen)sizeof(integer));
    do_lio(&c__5, &c__1, (char *)&para_1.side, (ftnlen)sizeof(doublereal));
    do_lio(&c__5, &c__1, (char *)&rcut, (ftnlen)sizeof(doublereal));
    do_lio(&c__5, &c__1, (char *)&rmin, (ftnlen)sizeof(doublereal));
    do_lio(&c__5, &c__1, (char *)&temp, (ftnlen)sizeof(doublereal));
    do_lio(&c__5, &c__1, (char *)&dt, (ftnlen)sizeof(doublereal));
    do_lio(&c__5, &c__1, (char *)&vmax, (ftnlen)sizeof(doublereal));
    do_lio(&c__3, &c__1, (char *)&neq, (ftnlen)sizeof(integer));
    do_lio(&c__3, &c__1, (char *)&nmeas, (ftnlen)sizeof(integer));
    do_lio(&c__3, &c__1, (char *)&nevry, (ftnlen)sizeof(integer));
    e_rsle();
    para_1.sideh = para_1.side / 2.;
    para_1.rcutsq = rcut * rcut;
    rmin2 = rmin * rmin;
/* Computing 2nd power */
    d__1 = para_1.side;
    rho = 16. / (d__1 * d__1);
    d1 = dt * (float)12. * dt;
    d2 = dt * (float)12.;
/*      Initilization */
    time = 0.;
    i__1 = -iseed;
    rr = rsub_(&i__1);
    if (ir == 1) {
/*       PLACE THEM RANDOMLY */
	coor_1.x[0] = 0.;
	coor_1.y[0] = 0.;
	for (i__ = 2; i__ <= 16; ++i__) {
L200:
	    coor_1.x[i__ - 1] = para_1.sideh * (rsub_(&c__0) * 2. - 1.);
	    coor_1.y[i__ - 1] = para_1.sideh * (rsub_(&c__0) * 2. - 1.);
	    i__1 = i__ - 1;
	    for (j = 1; j <= i__1; ++j) {
		xx = coor_1.x[i__ - 1] - coor_1.x[j - 1];
		if (xx < -para_1.sideh) {
		    xx += para_1.side;
		}
		if (xx > para_1.sideh) {
		    xx -= para_1.side;
		}
		yy = coor_1.y[i__ - 1] - coor_1.y[j - 1];
		if (yy < -para_1.sideh) {
		    yy += para_1.side;
		}
		if (yy > para_1.sideh) {
		    yy -= para_1.side;
		}
		dist = xx * xx + yy * yy;
		if (dist < rmin2) {
		    goto L200;
		}
	    }
	}
    } else {
/*       PLACE THEM ON A Lattice */
	nside = (integer) para_1.side;
	i__1 = nside;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    i__2 = nside;
	    for (j = 1; j <= i__2; ++j) {
		k = i__ + (j - 1) * nside;
		if (k > 16) {
		    goto L300;
		}
		coor_1.x[k - 1] = (doublereal) i__ - para_1.sideh;
		coor_1.y[k - 1] = (doublereal) j - para_1.sideh;
	    }
	}
L300:
	;
    }
    s_wsle(&io___25);
    do_lio(&c__9, &c__1, "I have placed the particles", 27L);
    e_wsle();
    for (i__ = 1; i__ <= 16; ++i__) {
	s_wsfe(&io___26);
	do_fio(&c__1, (char *)&coor_1.x[i__ - 1], (ftnlen)sizeof(doublereal));
	do_fio(&c__1, (char *)&coor_1.y[i__ - 1], (ftnlen)sizeof(doublereal));
	e_wsfe();
    }
    for (i__ = 1; i__ <= 16; ++i__) {
	vel_1.vx[i__ - 1] = (rsub_(&c__0) * 2. - 1.) * vmax;
	vel_1.vy[i__ - 1] = (rsub_(&c__0) * 2. - 1.) * vmax;
    }
    vxsum = 0.;
    vysum = 0.;
    for (i__ = 1; i__ <= 16; ++i__) {
	vxsum += vel_1.vx[i__ - 1];
	vysum += vel_1.vy[i__ - 1];
    }
    vxsum /= 16.;
    vysum /= 16.;
    for (i__ = 1; i__ <= 16; ++i__) {
	vel_1.vx[i__ - 1] -= vxsum;
	vel_1.vy[i__ - 1] -= vysum;
    }
/*      Initial Kinetic Energy */
    kinetic_(&ekin);
/*      Initial Potential Energy */
    fcal_(&ep, &vir);
    etot = ekin + ep;
    s_wsle(&io___33);
    do_lio(&c__9, &c__1, "time,ekin,ep,etot", 17L);
    e_wsle();
    s_wsfe(&io___34);
    do_fio(&c__1, (char *)&time, (ftnlen)sizeof(doublereal));
    do_fio(&c__1, (char *)&ekin, (ftnlen)sizeof(doublereal));
    do_fio(&c__1, (char *)&ep, (ftnlen)sizeof(doublereal));
    do_fio(&c__1, (char *)&etot, (ftnlen)sizeof(doublereal));
    e_wsfe();
/*        EQUILIBRATION STARTS! */
    i__1 = neq;
    for (istep = 1; istep <= i__1; ++istep) {
	for (i__ = 1; i__ <= 16; ++i__) {
	    coor_1.x[i__ - 1] = coor_1.x[i__ - 1] + dt * vel_1.vx[i__ - 1] + 
		    d1 * force_1.fx[i__ - 1];
	    if (coor_1.x[i__ - 1] < -para_1.sideh) {
		coor_1.x[i__ - 1] += para_1.side;
	    }
	    if (coor_1.x[i__ - 1] > para_1.sideh) {
		coor_1.x[i__ - 1] -= para_1.side;
	    }
	    coor_1.y[i__ - 1] = coor_1.y[i__ - 1] + dt * vel_1.vy[i__ - 1] + 
		    d1 * force_1.fy[i__ - 1];
	    if (coor_1.y[i__ - 1] < -para_1.sideh) {
		coor_1.y[i__ - 1] += para_1.side;
	    }
	    if (coor_1.y[i__ - 1] > para_1.sideh) {
		coor_1.y[i__ - 1] -= para_1.side;
	    }
	    vel_1.vx[i__ - 1] += d2 * force_1.fx[i__ - 1];
	    vel_1.vy[i__ - 1] += d2 * force_1.fy[i__ - 1];
	}
	fcal_(&ep, &vir);
	for (i__ = 1; i__ <= 16; ++i__) {
	    vel_1.vx[i__ - 1] += d2 * force_1.fx[i__ - 1];
	    vel_1.vy[i__ - 1] += d2 * force_1.fy[i__ - 1];
	}
	time += dt;
/*       Check Energies during Equilibration */
	kinetic_(&ekin);
	etot = ekin + ep;
	s_wsfe(&io___36);
	do_fio(&c__1, (char *)&time, (ftnlen)sizeof(doublereal));
	do_fio(&c__1, (char *)&ekin, (ftnlen)sizeof(doublereal));
	do_fio(&c__1, (char *)&ep, (ftnlen)sizeof(doublereal));
	do_fio(&c__1, (char *)&etot, (ftnlen)sizeof(doublereal));
	e_wsfe();
    }
/*      MEASUREMENT STARTS */
    ekinav = 0.;
    epotav = 0.;
    virav = 0.;
    for (i__ = 1; i__ <= 100; ++i__) {
	pvav[i__ - 1] = 0.;
	grav[i__ - 1] = 0.;
    }
    i__1 = nmeas;
    for (ijk = 1; ijk <= i__1; ++ijk) {
	i__2 = nevry;
	for (jkl = 1; jkl <= i__2; ++jkl) {
	    for (i__ = 1; i__ <= 16; ++i__) {
		coor_1.x[i__ - 1] = coor_1.x[i__ - 1] + dt * vel_1.vx[i__ - 1]
			 + d1 * force_1.fx[i__ - 1];
		if (coor_1.x[i__ - 1] < -para_1.sideh) {
		    coor_1.x[i__ - 1] += para_1.side;
		}
		if (coor_1.x[i__ - 1] > para_1.sideh) {
		    coor_1.x[i__ - 1] -= para_1.side;
		}
		coor_1.y[i__ - 1] = coor_1.y[i__ - 1] + dt * vel_1.vy[i__ - 1]
			 + d1 * force_1.fy[i__ - 1];
		if (coor_1.y[i__ - 1] < -para_1.sideh) {
		    coor_1.y[i__ - 1] += para_1.side;
		}
		if (coor_1.y[i__ - 1] > para_1.sideh) {
		    coor_1.y[i__ - 1] -= para_1.side;
		}
		vel_1.vx[i__ - 1] += d2 * force_1.fx[i__ - 1];
		vel_1.vy[i__ - 1] += d2 * force_1.fy[i__ - 1];
	    }
	    fcal_(&ep, &vir);
	    for (i__ = 1; i__ <= 16; ++i__) {
		vel_1.vx[i__ - 1] += d2 * force_1.fx[i__ - 1];
		vel_1.vy[i__ - 1] += d2 * force_1.fy[i__ - 1];
	    }
	    time += dt;
	}
	epotav += ep;
	virav += vir;
	corr_();
	for (i__ = 1; i__ <= 100; ++i__) {
	    grav[i__ - 1] += meas2_1.gr[i__ - 1];
	}
	kinetic_(&ekin);
	ekinav += ekin;
	for (i__ = 1; i__ <= 100; ++i__) {
	    pvav[i__ - 1] += meas1_1.pv[i__ - 1];
	}
    }
    epotav /= (doublereal) nmeas;
    ekinav /= (doublereal) nmeas;
    ener = epotav + ekinav;
    temp = ekinav / 16.;
    virav /= (doublereal) nmeas;
/* Computing 2nd power */
    d__1 = para_1.side;
    press = (16. * temp + virav * (float).5) / (d__1 * d__1);
    s_wsle(&io___46);
    do_lio(&c__9, &c__1, "PE, KE, ENER, T,P,RHO", 21L);
    e_wsle();
    s_wsfe(&io___47);
    do_fio(&c__1, (char *)&epotav, (ftnlen)sizeof(doublereal));
    do_fio(&c__1, (char *)&ekinav, (ftnlen)sizeof(doublereal));
    do_fio(&c__1, (char *)&ener, (ftnlen)sizeof(doublereal));
    do_fio(&c__1, (char *)&temp, (ftnlen)sizeof(doublereal));
    do_fio(&c__1, (char *)&press, (ftnlen)sizeof(doublereal));
    do_fio(&c__1, (char *)&rho, (ftnlen)sizeof(doublereal));
    e_wsfe();
    for (i__ = 1; i__ <= 100; ++i__) {
	pvav[i__ - 1] /= (doublereal) nmeas;
	grav[i__ - 1] /= (doublereal) nmeas;
    }
    for (i__ = 1; i__ <= 100; ++i__) {
	s_wsfe(&io___48);
	d__1 = (doublereal) i__ * (float).1;
	do_fio(&c__1, (char *)&d__1, (ftnlen)sizeof(doublereal));
	do_fio(&c__1, (char *)&grav[i__ - 1], (ftnlen)sizeof(doublereal));
	e_wsfe();
	s_wsfe(&io___49);
	d__1 = (doublereal) i__ * (float).1;
	do_fio(&c__1, (char *)&d__1, (ftnlen)sizeof(doublereal));
	do_fio(&c__1, (char *)&pvav[i__ - 1], (ftnlen)sizeof(doublereal));
	e_wsfe();
    }
    for (i__ = 1; i__ <= 16; ++i__) {
	s_wsfe(&io___50);
	do_fio(&c__1, (char *)&coor_1.x[i__ - 1], (ftnlen)sizeof(doublereal));
	do_fio(&c__1, (char *)&coor_1.y[i__ - 1], (ftnlen)sizeof(doublereal));
	e_wsfe();
    }
} /* MAIN__ */

/* Subroutine */ int fcal_(epot, virial)
doublereal *epot, *virial;
{
    static integer i__, j;
    static doublereal force, xcomp, ycomp, r2, xi, yi, xx, yy, sr2, sr6, eij, 
	    sr12, vij;

    *epot = 0.;
    *virial = 0.;
    for (i__ = 1; i__ <= 16; ++i__) {
	force_1.fx[i__ - 1] = 0.;
	force_1.fy[i__ - 1] = 0.;
    }
    for (i__ = 1; i__ <= 15; ++i__) {
	xi = coor_1.x[i__ - 1];
	yi = coor_1.y[i__ - 1];
	for (j = i__ + 1; j <= 16; ++j) {
	    xx = xi - coor_1.x[j - 1];
	    yy = yi - coor_1.y[j - 1];
	    if (xx < -para_1.sideh) {
		xx += para_1.side;
	    }
	    if (xx > para_1.sideh) {
		xx -= para_1.side;
	    }
	    if (yy < -para_1.sideh) {
		yy += para_1.side;
	    }
	    if (yy > para_1.sideh) {
		yy -= para_1.side;
	    }
	    r2 = xx * xx + yy * yy;
	    if (r2 > para_1.rcutsq) {
		goto L10;
	    }
	    sr2 = 1. / r2;
	    sr6 = sr2 * sr2 * sr2;
	    sr12 = sr6 * sr6;
	    eij = sr12 - sr6;
	    *epot += eij;
	    vij = eij + sr12;
	    *virial += vij;
	    force = vij * sr2;
	    xcomp = xx * force;
	    force_1.fx[i__ - 1] += xcomp;
	    force_1.fx[j - 1] -= xcomp;
	    ycomp = yy * force;
	    force_1.fy[i__ - 1] += ycomp;
	    force_1.fy[j - 1] -= ycomp;
L10:
	    ;
	}
    }
    *epot *= 4.;
    *virial *= 24.;
    return 0;
} /* fcal_ */

/* Subroutine */ int kinetic_(ekin)
doublereal *ekin;
{
    /* Builtin functions */
    double sqrt();

    /* Local variables */
    static integer i__, j;
    static doublereal v, v2;

    for (i__ = 1; i__ <= 100; ++i__) {
	meas1_1.pv[i__ - 1] = 0.;
    }
    *ekin = 0.;
    for (i__ = 1; i__ <= 16; ++i__) {
	v2 = vel_1.vx[i__ - 1] * vel_1.vx[i__ - 1] + vel_1.vy[i__ - 1] * 
		vel_1.vy[i__ - 1];
	*ekin += v2;
	v = sqrt(v2);
	if (v <= 10.) {
	    j = (integer) (v / (float).1);
	    meas1_1.pv[j - 1] += 1.;
	}
    }
    *ekin *= (float).5;
    return 0;
} /* kinetic_ */

/* Subroutine */ int corr_()
{
    /* Builtin functions */
    double atan(), sqrt();

    /* Local variables */
    static integer i__, j;
    static doublereal r__, r2, pi, xi, yi, rr, xx, yy;

    pi = atan(1.) * 4.;
    for (i__ = 1; i__ <= 100; ++i__) {
	meas2_1.gr[i__ - 1] = 0.;
    }
    for (i__ = 1; i__ <= 15; ++i__) {
	xi = coor_1.x[i__ - 1];
	yi = coor_1.y[i__ - 1];
	for (j = i__ + 1; j <= 16; ++j) {
	    xx = xi - coor_1.x[j - 1];
	    yy = yi - coor_1.y[j - 1];
	    if (xx < -para_1.sideh) {
		xx += para_1.side;
	    }
	    if (xx > para_1.sideh) {
		xx -= para_1.side;
	    }
	    if (yy < -para_1.sideh) {
		yy += para_1.side;
	    }
	    if (yy > para_1.sideh) {
		yy -= para_1.side;
	    }
	    r2 = xx * xx + yy * yy;
/*        compute g(r) */
	    rr = sqrt(r2);
	    if (rr <= 10.) {
		j = (integer) (rr / (float).1);
		meas2_1.gr[j - 1] += 1.;
	    }
/* L10: */
	}
    }
    for (i__ = 1; i__ <= 100; ++i__) {
	r__ = (doublereal) i__ * (float).1;
	meas2_1.gr[i__ - 1] /= pi * 2. * r__ * (float).1;
    }
    return 0;
} /* corr_ */

doublereal rsub_(iseed)
integer *iseed;
{
    /* Initialized data */

    static logical ini = FALSE_;

    /* System generated locals */
    doublereal ret_val;

    /* Local variables */
    static integer i__, j, k, jrand, istack[55], idx;

/* 	subtraction method of random number generation (Box 7-2) */
/* 	Define three constants m10=10^9, i_s=21 and i_r=30. */
/* 	Set up a stack $S(k)$ of length $L_k$ for integers. */
/* ================================================================ */
/*        ADDED BY Amit Chakrabarti (12-13-94) to solve problem with */
/*        SILICON GRAPHICS MACHINES */
/* ================================================================ */
/* 	reinitializes if ISEED < 0 */
    if (! ini || *iseed < 0) {
/* 	initialization */
/* 	Store the the absolute value of the seed $s_d$ in $S(L_k)$ */
	*iseed = abs(*iseed);
	istack[54] = *iseed;
/* 		Let $k=1$ and $j=s_d$. */
	j = *iseed;
	k = 1;
/* 		carry out the following steps to store the stack: */
	for (i__ = 1; i__ <= 54; ++i__) {
/* 			Let $\ell_x=i_s*\ell \mod L_k$. */
	    idx = i__ * 21 % 55;
/* 			Store the value $k$ in $S(\ell_x)$. */
	    istack[idx - 1] = k;
/* 			Replace $k$ by $j-k$. If $j<0$ increase $j$ by 10^9 */
	    k = j - k;
	    if (k < 0) {
		k += 1000000000;
	    }
/* 			Let $j=S(\ell_x)$ */
	    j = istack[idx - 1];
/* L50: */
	}
/* 		Randomize the stack by carrying out the following steps */
/* 			three times: */
	for (j = 1; j <= 3; ++j) {
	    for (i__ = 1; i__ <= 55; ++i__) {
/* 			For $\ell=1$ to 24 ($\ell+i_r \leq L_k$), */
/* 				replace $S(\ell)$ by $S(\ell)-S(\ell+31)$. */
/* 				For $\ell=25$ to $L_k$, replace $S(\ell)$ */
/* 					by $S(\ell)-S(\ell-24)$. */
		istack[i__ - 1] -= istack[(i__ + 30) % 55];
/* 				If $S(\ell)$ becomes negative, */
/* 					increase $S(\ell)$ by $10^9$. */
		if (istack[i__ - 1] < 0) {
		    istack[i__ - 1] += 1000000000;
		}
/* L80: */
	    }
/* L90: */
	}
/* 		Set stack counter $j_s = 0$. */
	jrand = 0;
/* 		Mark the stack as initialized. */
	ini = TRUE_;
    }
/* 	normal entry point of rsub */

/* 	Increase the stack counter $j_s$ by 1. */
    ++jrand;
    if (jrand > 55) {
/* 	Replace the stack by repeating the do 80 loop */
	for (i__ = 1; i__ <= 55; ++i__) {
	    istack[i__ - 1] -= istack[(i__ + 30) % 55];
	    if (istack[i__ - 1] < 0) {
		istack[i__ - 1] += 1000000000;
	    }
/* L300: */
	}
/* 		Set stack counter $j_s = 1$. */
	jrand = 1;
    }
/* 	Return the value of (float) $S(j_s) \times 10^{-9}$ */
/* 		as the uniform random number in the interval $[0,1]$. */
    ret_val = (doublereal) istack[jrand - 1] / 1e9;
/* write(6,*) jrand, rsub */
    return ret_val;
} /* rsub_ */

/* Main program alias */ int md_e__ () { MAIN__ (); }