/* shockt.f -- translated by f2c (version 20031025).
   You must link the resulting object file with libf2c:
	on Microsoft Windows system, link with libf2c.lib;
	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
	or, if you install libf2c.a in a standard place, with -lf2c -lm
	-- in that order, at the end of the command line, as in
		cc *.o -lf2c -lm
	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

		http://www.netlib.org/f2c/libf2c.zip
*/

#include "f2c.h"

/* Common Block Declarations */

struct {
    real g, r1, p1, e1, u1, r4, p4, e4, u4, xd0;
} setups_;

#define setups_1 setups_

struct {
    integer mx;
    real dx, x[201], cfl, dt;
    integer nlast;
    real delta;
} cmpcnd_;

#define cmpcnd_1 cmpcnd_

struct {
    real r__[201], p[201], e[201], u[201], q[603]	/* was [3][201] */, 
	    dq[603]	/* was [3][201] */;
} flwdat_;

#define flwdat_1 flwdat_

/* Table of constant values */

static integer c__9 = 9;
static integer c__1 = 1;
static integer c__4 = 4;
static integer c__3 = 3;
static real c_b61 = 1.f;

/* Main program */ int MAIN__(void)
{
    /* Format strings */
    static char fmt_501[] = "(5f15.4)";

    /* System generated locals */
    integer i__1, i__2;
    real r__1, r__2, r__3, r__4;
    olist o__1;
    cllist cl__1;

    /* Builtin functions */
    integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), 
	    e_wsle(void), s_rsle(cilist *), e_rsle(void);
    double sqrt(doublereal), r_sign(real *, real *);
    integer f_open(olist *), s_wsfe(cilist *), do_fio(integer *, char *, 
	    ftnlen), e_wsfe(void), f_clos(cllist *);
    /* Subroutine */ int s_stop(char *, ftnlen);

    /* Local variables */
    static integer i__, n;
    static real s, b1, b2, a11, a12, a13, a21, a22, a23, a31, a32, a33, cm[
	    201], hm[201], x10, x40, um[201], qq, ff1, al1[201], al2[201], 
	    al3[201], ff2, cm2, dq1, dq2, dq3, ph1[201], ph2[201], ph3[201], 
	    dx2, tv1[201], tv2[201], tv3[201], dxh, xdo, tvd1, tvd2, tvd3, 
	    dash, time, vref, ramda;

    /* Fortran I/O blocks */
    static cilist io___1 = { 0, 6, 0, 0, 0 };
    static cilist io___4 = { 0, 6, 0, 0, 0 };
    static cilist io___5 = { 0, 6, 0, 0, 0 };
    static cilist io___6 = { 0, 6, 0, 0, 0 };
    static cilist io___7 = { 0, 6, 0, 0, 0 };
    static cilist io___8 = { 0, 6, 0, 0, 0 };
    static cilist io___9 = { 0, 6, 0, 0, 0 };
    static cilist io___13 = { 0, 6, 0, 0, 0 };
    static cilist io___14 = { 0, 5, 0, 0, 0 };
    static cilist io___16 = { 0, 6, 0, 0, 0 };
    static cilist io___17 = { 0, 6, 0, 0, 0 };
    static cilist io___18 = { 0, 5, 0, 0, 0 };
    static cilist io___19 = { 0, 6, 0, 0, 0 };
    static cilist io___20 = { 0, 6, 0, 0, 0 };
    static cilist io___24 = { 0, 6, 0, 0, 0 };
    static cilist io___61 = { 0, 6, 0, 0, 0 };
    static cilist io___62 = { 0, 6, 0, 0, 0 };
    static cilist io___63 = { 0, 60, 0, 0, 0 };
    static cilist io___64 = { 0, 60, 0, 0, 0 };
    static cilist io___65 = { 0, 60, 0, fmt_501, 0 };
    static cilist io___66 = { 0, 6, 0, 0, 0 };
    static cilist io___67 = { 0, 6, 0, 0, 0 };


    s_wsle(&io___1);
    do_lio(&c__9, &c__1, ">> SOLVE SHOCK-TUBE PROBLEM USING TVD METHOD", (
	    ftnlen)44);
    e_wsle();
/* ......(1) SHOCK-TUBE SETUP */
    setups_1.g = 1.4f;
    setups_1.r1 = .25f;
    setups_1.p1 = .2154f;
    setups_1.e1 = setups_1.p1 / ((setups_1.g - 1.f) * setups_1.r1);
    setups_1.u1 = 0.f;
    setups_1.r4 = 1.f;
    setups_1.p4 = 1.f;
    setups_1.e4 = setups_1.p4 / ((setups_1.g - 1.f) * setups_1.r4);
    setups_1.u4 = 0.f;
    setups_1.xd0 = 0.f;
    x40 = -.6f;
    x10 = .6f;
    s_wsle(&io___4);
    do_lio(&c__9, &c__1, " ", (ftnlen)1);
    e_wsle();
    s_wsle(&io___5);
    do_lio(&c__9, &c__1, ">> SHOCK-TUBE SETUP", (ftnlen)19);
    e_wsle();
    s_wsle(&io___6);
    do_lio(&c__9, &c__1, "  HIGH PRESSURE SIDE : PRESSURE = ", (ftnlen)34);
    do_lio(&c__4, &c__1, (char *)&setups_1.p4, (ftnlen)sizeof(real));
    e_wsle();
    s_wsle(&io___7);
    do_lio(&c__9, &c__1, "                       DENSITY  = ", (ftnlen)34);
    do_lio(&c__4, &c__1, (char *)&setups_1.r4, (ftnlen)sizeof(real));
    e_wsle();
    s_wsle(&io___8);
    do_lio(&c__9, &c__1, "   LOW PRESSURE SIDE : PRESSURE = ", (ftnlen)34);
    do_lio(&c__4, &c__1, (char *)&setups_1.p1, (ftnlen)sizeof(real));
    e_wsle();
    s_wsle(&io___9);
    do_lio(&c__9, &c__1, "                       DENSITY  = ", (ftnlen)34);
    do_lio(&c__4, &c__1, (char *)&setups_1.p4, (ftnlen)sizeof(real));
    e_wsle();
/* ......(2) GRID GENERATION */
    cmpcnd_1.mx = 121;
    cmpcnd_1.dx = (x10 - x40) / (real) (cmpcnd_1.mx - 1);
    dx2 = cmpcnd_1.dx * 2.f;
    dxh = cmpcnd_1.dx / 2.f;
    i__1 = cmpcnd_1.mx + 1;
    for (i__ = 1; i__ <= i__1; ++i__) {
	cmpcnd_1.x[i__ - 1] = x40 + cmpcnd_1.dx * (real) (i__ - 1);
/* L200: */
    }
/* ......(3) TIME SETUP SIZE */
    s_wsle(&io___13);
    do_lio(&c__9, &c__1, ">> INPUT CFL NUMBER ===>", (ftnlen)24);
    e_wsle();
    s_rsle(&io___14);
    do_lio(&c__4, &c__1, (char *)&cmpcnd_1.cfl, (ftnlen)sizeof(real));
    e_rsle();
    vref = 1.f;
    cmpcnd_1.dt = cmpcnd_1.cfl * cmpcnd_1.dx / vref;
    s_wsle(&io___16);
    do_lio(&c__9, &c__1, "TIME STEP SIZE =", (ftnlen)16);
    do_lio(&c__4, &c__1, (char *)&cmpcnd_1.dt, (ftnlen)sizeof(real));
    e_wsle();
    s_wsle(&io___17);
    do_lio(&c__9, &c__1, ">> INPUT NO. OF TIME STEPS ===>", (ftnlen)31);
    e_wsle();
    s_rsle(&io___18);
    do_lio(&c__3, &c__1, (char *)&cmpcnd_1.nlast, (ftnlen)sizeof(integer));
    e_rsle();
/* ......(4) TVD SOLUTION */
/* ......4-1 INITIAL CONDITION */
    s_wsle(&io___19);
    do_lio(&c__9, &c__1, " ", (ftnlen)1);
    e_wsle();
    s_wsle(&io___20);
    do_lio(&c__9, &c__1, ">> NUMERICAL SIMULATION STARTS !!", (ftnlen)33);
    e_wsle();
    i__1 = cmpcnd_1.mx;
    for (i__ = 1; i__ <= i__1; ++i__) {
	if (cmpcnd_1.x[i__ - 1] < xdo) {
	    flwdat_1.r__[i__ - 1] = setups_1.r4;
	    flwdat_1.p[i__ - 1] = setups_1.p4;
	    flwdat_1.e[i__ - 1] = setups_1.e4;
	    flwdat_1.u[i__ - 1] = setups_1.u4;
	} else {
	    flwdat_1.r__[i__ - 1] = setups_1.r1;
	    flwdat_1.p[i__ - 1] = setups_1.p1;
	    flwdat_1.e[i__ - 1] = setups_1.e1;
	    flwdat_1.u[i__ - 1] = setups_1.u1;
	}
	flwdat_1.q[i__ * 3 - 3] = flwdat_1.r__[i__ - 1];
	flwdat_1.q[i__ * 3 - 2] = flwdat_1.r__[i__ - 1] * flwdat_1.u[i__ - 1];
/* Computing 2nd power */
	r__1 = flwdat_1.u[i__ - 1];
	flwdat_1.q[i__ * 3 - 1] = flwdat_1.r__[i__ - 1] * flwdat_1.e[i__ - 1] 
		+ flwdat_1.r__[i__ - 1] * .5f * (r__1 * r__1);
/* L410: */
    }
/* ......4-2 TIME MARCHING */
    i__1 = cmpcnd_1.nlast;
    for (n = 1; n <= i__1; ++n) {
	time = cmpcnd_1.dt * (real) n;
	if (n / 10 * 10 == n) {
	    s_wsle(&io___24);
	    do_lio(&c__9, &c__1, " ... STEP = ", (ftnlen)12);
	    do_lio(&c__3, &c__1, (char *)&n, (ftnlen)sizeof(integer));
	    e_wsle();
	}
/*          <1> I+1/2 VALUES */
	i__2 = cmpcnd_1.mx - 1;
	for (i__ = 1; i__ <= i__2; ++i__) {
/*             <== ROE'S AVERAGING */
	    dash = sqrt(flwdat_1.r__[i__] / flwdat_1.r__[i__ - 1]);
	    um[i__ - 1] = (dash * flwdat_1.u[i__] + flwdat_1.u[i__ - 1]) / (
		    dash + 1.f);
	    hm[i__ - 1] = (dash * (flwdat_1.q[(i__ + 1) * 3 - 1] + flwdat_1.p[
		    i__]) / flwdat_1.r__[i__] + (flwdat_1.q[i__ * 3 - 1] + 
		    flwdat_1.p[i__ - 1]) / flwdat_1.r__[i__ - 1]) / (dash + 
		    1.f);
/* Computing 2nd power */
	    r__1 = um[i__ - 1];
	    cm2 = (setups_1.g - 1.f) * (hm[i__ - 1] - r__1 * r__1 * .5f);
/* Computing MAX */
/* Computing MIN */
	    r__3 = setups_1.g * (setups_1.g - 1.f) * flwdat_1.e[i__ - 1], 
		    r__4 = setups_1.g * (setups_1.g - 1.f) * flwdat_1.e[i__];
	    r__1 = cm2, r__2 = dmin(r__3,r__4);
	    cm[i__ - 1] = sqrt((dmax(r__1,r__2)));
	    dq1 = flwdat_1.q[(i__ + 1) * 3 - 3] - flwdat_1.q[i__ * 3 - 3];
	    dq2 = flwdat_1.q[(i__ + 1) * 3 - 2] - flwdat_1.q[i__ * 3 - 2];
	    dq3 = flwdat_1.q[(i__ + 1) * 3 - 1] - flwdat_1.q[i__ * 3 - 1];
/*             <== COMPUTR R^-1 */
/* Computing 2nd power */
	    r__1 = cm[i__ - 1];
	    b2 = (setups_1.g - 1.f) / (r__1 * r__1);
/* Computing 2nd power */
	    r__1 = um[i__ - 1];
	    b1 = b2 * (r__1 * r__1) / 2.f;
	    a11 = (b1 + um[i__ - 1] / cm[i__ - 1]) * .5f;
	    a12 = (-b2 * um[i__ - 1] - 1.f / cm[i__ - 1]) * .5f;
	    a13 = b2 * .5f;
	    a21 = 1.f - b1;
	    a22 = b2 * um[i__ - 1];
	    a23 = -b2;
	    a31 = (b1 - um[i__ - 1] / cm[i__ - 1]) * .5f;
	    a32 = (-b2 * um[i__ - 1] + 1.f / cm[i__ - 1]) * .5f;
	    a33 = b2 * .5f;
/*             <== COMPUTE A=R^-1*DQ */
	    al1[i__ - 1] = a11 * dq1 + a12 * dq2 + a13 * dq3;
	    al2[i__ - 1] = a21 * dq1 + a22 * dq2 + a23 * dq3;
	    al3[i__ - 1] = a31 * dq1 + a32 * dq2 + a33 * dq3;
/* L4210: */
	}
/*          <2> LIMITER FUNCTION */
/*             MINMOD(X,Y,Z)=S*MAX[0,MIN(S*X,S*Y,S*Z)], S=SGN(X) */
/*          <== ENTROPY CORRECTION PARAMETER */
	cmpcnd_1.delta = .15f;
	i__2 = cmpcnd_1.mx - 2;
	for (i__ = 2; i__ <= i__2; ++i__) {
	    ramda = um[i__ - 1] - cm[i__ - 1];
/* Computing 2nd power */
	    r__1 = ramda;
	    ff1 = cmpcnd_1.dt / cmpcnd_1.dx * (r__1 * r__1);
	    ff2 = dabs(ramda);
	    if (ff2 < cmpcnd_1.delta) {
/* Computing 2nd power */
		r__1 = ramda;
/* Computing 2nd power */
		r__2 = cmpcnd_1.delta;
		ff2 = (r__1 * r__1 + r__2 * r__2) * .5f / cmpcnd_1.delta;
	    }
	    s = r_sign(&c_b61, &al1[1 - i__ - 1]);
/* Computing MAX */
/* Computing MIN */
	    r__3 = s * 2.f * al1[i__ - 2], r__4 = s * 2.f * al1[i__ - 1], 
		    r__3 = min(r__3,r__4), r__4 = s * 2.f * al1[i__], r__3 = 
		    min(r__3,r__4), r__4 = s * .5f * (al1[i__ - 2] + al1[i__])
		    ;
	    r__1 = 0.f, r__2 = dmin(r__3,r__4);
	    qq = s * dmax(r__1,r__2);
	    ph1[i__ - 1] = -ff1 * qq - ff2 * (al1[i__ - 1] - qq);
	    ramda = um[i__ - 1];
/* Computing 2nd power */
	    r__1 = ramda;
	    ff1 = cmpcnd_1.dt / cmpcnd_1.dx * (r__1 * r__1);
	    ff2 = dabs(ramda);
	    if (ff2 < cmpcnd_1.delta) {
/* Computing 2nd power */
		r__1 = ramda;
/* Computing 2nd power */
		r__2 = cmpcnd_1.delta;
		ff2 = (r__1 * r__1 + r__2 * r__2) * .5f / cmpcnd_1.delta;
	    }
	    s = r_sign(&c_b61, &al2[i__ - 2]);
/* Computing MAX */
/* Computing MIN */
	    r__3 = s * 2.f * al2[i__ - 2], r__4 = s * 2.f * al2[i__ - 1], 
		    r__3 = min(r__3,r__4), r__4 = s * 2.f * al2[i__], r__3 = 
		    min(r__3,r__4), r__4 = s * .5f * (al2[i__ - 2] + al2[i__])
		    ;
	    r__1 = 0.f, r__2 = dmin(r__3,r__4);
	    qq = s * dmax(r__1,r__2);
	    ph2[i__ - 1] = -ff1 * qq - ff2 * (al2[i__ - 1] - qq);
	    ramda = um[i__ - 1] + cm[i__ - 1];
/* Computing 2nd power */
	    r__1 = ramda;
	    ff1 = cmpcnd_1.dt / cmpcnd_1.dx * (r__1 * r__1);
	    ff2 = dabs(ramda);
	    if (ff2 < cmpcnd_1.delta) {
/* Computing 2nd power */
		r__1 = ramda;
/* Computing 2nd power */
		r__2 = cmpcnd_1.delta;
		ff2 = (r__1 * r__1 + r__2 * r__2) * .5f / cmpcnd_1.delta;
	    }
	    s = r_sign(&c_b61, &al3[i__ - 2]);
/* Computing MAX */
/* Computing MIN */
	    r__3 = s * 2.f * al3[i__ - 2], r__4 = s * 2.f * al3[i__ - 1], 
		    r__3 = min(r__3,r__4), r__4 = s * 2.f * al3[i__], r__3 = 
		    min(r__3,r__4), r__4 = s * .5f * (al3[i__ - 2] + al2[i__])
		    ;
	    r__1 = 0.f, r__2 = dmin(r__3,r__4);
	    qq = s * dmax(r__1,r__2);
	    ph3[i__ - 1] = -ff1 * qq - ff2 * (al3[i__ - 1] - qq);
/* L4220: */
	}
/*          <3> NUMERICAL FLUX */
	i__2 = cmpcnd_1.mx - 2;
	for (i__ = 2; i__ <= i__2; ++i__) {
	    tvd1 = ph1[i__ - 1] + ph2[i__ - 1] + ph3[i__ - 1];
	    tvd2 = um[i__ - 1] * (ph1[i__ - 1] + ph2[i__ - 1] + ph3[i__ - 1]) 
		    - cm[i__ - 1] * (ph1[i__ - 1] - ph3[i__ - 1]);
/* Computing 2nd power */
	    r__1 = um[i__ - 1];
	    tvd3 = (hm[i__ - 1] - um[i__ - 1] * cm[i__ - 1]) * ph1[i__ - 1] + 
		    r__1 * r__1 * .5f * ph2[i__ - 1] + (hm[i__ - 1] + um[i__ 
		    - 1] * cm[i__ - 1]) * ph3[i__ - 1];
	    tv1[i__ - 1] = (flwdat_1.q[i__ * 3 - 2] + flwdat_1.q[(i__ + 1) * 
		    3 - 2] + tvd1) / 2.f;
	    tv2[i__ - 1] = (flwdat_1.q[i__ * 3 - 2] * flwdat_1.u[i__ - 1] + 
		    flwdat_1.p[i__ - 1] + flwdat_1.q[(i__ + 1) * 3 - 2] * 
		    flwdat_1.u[i__] + flwdat_1.p[i__] + tvd2) / 2.f;
	    tv3[i__ - 1] = ((flwdat_1.q[i__ * 3 - 1] + flwdat_1.p[i__ - 1]) * 
		    flwdat_1.u[i__ - 1] + (flwdat_1.q[(i__ + 1) * 3 - 1] + 
		    flwdat_1.p[i__]) * flwdat_1.u[i__] + tvd3) / 2.f;
/* L4230: */
	}
/*          <4> DIFFERENCING */
	i__2 = cmpcnd_1.mx - 2;
	for (i__ = 3; i__ <= i__2; ++i__) {
	    flwdat_1.dq[i__ * 3 - 3] = -(cmpcnd_1.dt / cmpcnd_1.dx) * (tv1[
		    i__ - 1] - tv1[i__ - 2]);
	    flwdat_1.dq[i__ * 3 - 2] = -(cmpcnd_1.dt / cmpcnd_1.dx) * (tv2[
		    i__ - 1] - tv2[i__ - 2]);
	    flwdat_1.dq[i__ * 3 - 1] = -(cmpcnd_1.dt / cmpcnd_1.dx) * (tv3[
		    i__ - 1] - tv3[i__ - 2]);
/* L4240: */
	}
/*          <5> UP-DATE */
	i__2 = cmpcnd_1.mx - 2;
	for (i__ = 3; i__ <= i__2; ++i__) {
	    flwdat_1.q[i__ * 3 - 3] += flwdat_1.dq[i__ * 3 - 3];
	    flwdat_1.q[i__ * 3 - 2] += flwdat_1.dq[i__ * 3 - 2];
	    flwdat_1.q[i__ * 3 - 1] += flwdat_1.dq[i__ * 3 - 1];
/* L4250: */
	}
/*          <6> SET BOUNDARY CONDITIONS */
	flwdat_1.q[0] = setups_1.r4;
	flwdat_1.q[1] = setups_1.r4 * setups_1.u4;
/* Computing 2nd power */
	r__1 = setups_1.u4;
	flwdat_1.q[2] = setups_1.r4 * setups_1.e4 + setups_1.r4 * .5f * (r__1 
		* r__1);
	flwdat_1.q[3] = setups_1.r4;
	flwdat_1.q[4] = setups_1.r4 * setups_1.u4;
/* Computing 2nd power */
	r__1 = setups_1.u4;
	flwdat_1.q[5] = setups_1.r4 * setups_1.e4 + setups_1.r4 * .5f * (r__1 
		* r__1);
	flwdat_1.q[(cmpcnd_1.mx - 1) * 3 - 3] = setups_1.r1;
	flwdat_1.q[(cmpcnd_1.mx - 1) * 3 - 2] = setups_1.r1 * setups_1.u1;
/* Computing 2nd power */
	r__1 = setups_1.u1;
	flwdat_1.q[(cmpcnd_1.mx - 1) * 3 - 1] = setups_1.r1 * setups_1.e1 + 
		setups_1.r1 * .5f * (r__1 * r__1);
	flwdat_1.q[cmpcnd_1.mx * 3 - 3] = setups_1.r1;
	flwdat_1.q[cmpcnd_1.mx * 3 - 2] = setups_1.r1 * setups_1.u1;
/* Computing 2nd power */
	r__1 = setups_1.u1;
	flwdat_1.q[cmpcnd_1.mx * 3 - 1] = setups_1.r1 * setups_1.e1 + 
		setups_1.r1 * .5f * (r__1 * r__1);
/*          <7>PRIMITIVE VARIABLES */
	i__2 = cmpcnd_1.mx;
	for (i__ = 1; i__ <= i__2; ++i__) {
	    flwdat_1.r__[i__ - 1] = flwdat_1.q[i__ * 3 - 3];
	    flwdat_1.u[i__ - 1] = flwdat_1.q[i__ * 3 - 2] / flwdat_1.r__[i__ 
		    - 1];
/* Computing 2nd power */
	    r__1 = flwdat_1.u[i__ - 1];
	    flwdat_1.e[i__ - 1] = flwdat_1.q[i__ * 3 - 1] / flwdat_1.r__[i__ 
		    - 1] - r__1 * r__1 * .5f;
	    flwdat_1.p[i__ - 1] = (setups_1.g - 1.f) * flwdat_1.r__[i__ - 1] *
		     flwdat_1.e[i__ - 1];
/* L4270: */
	}

/* L4200: */
    }
/* ......(5) WRITE TVD SOLUTION IN FILE */
    s_wsle(&io___61);
    do_lio(&c__9, &c__1, " ", (ftnlen)1);
    e_wsle();
    s_wsle(&io___62);
    do_lio(&c__9, &c__1, ">> WRITE CFD SOLUTION ON FILE.", (ftnlen)30);
    e_wsle();
    o__1.oerr = 0;
    o__1.ounit = 60;
    o__1.ofnmlen = 9;
    o__1.ofnm = "SOCK.DATA";
    o__1.orl = 0;
    o__1.osta = 0;
    o__1.oacc = 0;
    o__1.ofm = "FORMATTED";
    o__1.oblnk = 0;
    f_open(&o__1);
    s_wsle(&io___63);
    do_lio(&c__3, &c__1, (char *)&cmpcnd_1.mx, (ftnlen)sizeof(integer));
    do_lio(&c__4, &c__1, (char *)&cmpcnd_1.cfl, (ftnlen)sizeof(real));
    do_lio(&c__4, &c__1, (char *)&time, (ftnlen)sizeof(real));
    do_lio(&c__3, &c__1, (char *)&cmpcnd_1.nlast, (ftnlen)sizeof(integer));
    e_wsle();
    s_wsle(&io___64);
    do_lio(&c__9, &c__1, "   x   Density   Pressure   e   Velocity", (ftnlen)
	    40);
    e_wsle();
    i__1 = cmpcnd_1.mx;
    for (i__ = 1; i__ <= i__1; ++i__) {
	s_wsfe(&io___65);
	do_fio(&c__1, (char *)&cmpcnd_1.x[i__ - 1], (ftnlen)sizeof(real));
	do_fio(&c__1, (char *)&flwdat_1.r__[i__ - 1], (ftnlen)sizeof(real));
	do_fio(&c__1, (char *)&flwdat_1.p[i__ - 1], (ftnlen)sizeof(real));
	do_fio(&c__1, (char *)&flwdat_1.e[i__ - 1], (ftnlen)sizeof(real));
	do_fio(&c__1, (char *)&flwdat_1.u[i__ - 1], (ftnlen)sizeof(real));
	e_wsfe();
/* L500: */
    }
    cl__1.cerr = 0;
    cl__1.cunit = 60;
    cl__1.csta = 0;
    f_clos(&cl__1);

    s_wsle(&io___66);
    do_lio(&c__9, &c__1, " ", (ftnlen)1);
    e_wsle();
    s_wsle(&io___67);
    do_lio(&c__9, &c__1, ">> COMPLETED.", (ftnlen)13);
    e_wsle();


    s_stop("", (ftnlen)0);
    return 0;
} /* MAIN__ */