(* 4/11/95  modified to compute average value of absolute
        error between flux and accumulation curves *) 
massglobal/:massglobal[file1_String]:=Block[{dat,time,accum,influx,outflux,netflux,x,z,v,gr1,gr2,gr4,gr5,gr7,graphlist={}},
(*	SetOptions[ListPlot,DisplayFunction->Identity,Background->g10,Frame->True]; *)
(*	Run["a.out infleet out1 out2"];  *)
(* Plot global mass balance *)
	dat=ReadList[file1,{Number,Number,Number,Number}];
	time=take2[dat];time=time/3600.0;
	accum=take3[dat];accum=accum/10^6;
	netflux=take4[dat];netflux=netflux/10^6;
        error=Apply[Plus,Abs[accum-netflux]]/Length[accum];
        Print["Average of |mass_error| for global domain = ",error];
        Print[];
(* comment out this Write if don't want to go to output file. 
        NOTE!  Must open the output file (and close it) from
        the Mathematica session and must call it outfl *)
(*        Write[outfl,error];   *)
	gr1=ListPlot[xyjoin[time,netflux],PlotStyle->{t1,d1},FrameLabel->{"Time (hours from start)","Cumulative Mass (10^6)","",""},PlotLabel->"Entire Pond (accum is solid)"];
	gr2=ListPlot[xyjoin[time,accum],PlotStyle->{t1}];
	gr5=Show[gr1,gr2,Axes->False]; disp[gr5];
	Return[gr5]
	]

(* (* Plot global momentum balance *)
	dat=ReadList[file2,{Number,Number,Number,Number}];
	time=take2[dat];time=time/3600.0;
	accum=take3[dat];accum=accum/10^6;
	netflux=take4[dat];netflux=netflux/10^6;
	gr1=ListPlot[xyjoin[time,netflux],PlotStyle->{t1,d1},FrameLabel->{"Time (hours)","Cumulative Momentum (10^6)","",""},PlotLabel->"Entire Pond (accum is solid)"];
	gr2=ListPlot[xyjoin[time,accum],PlotStyle->{t1}];
	gr7=Show[gr1,gr2,Axes->False]; disp[gr7];
	Return[{gr5,gr7}]
	] *)

(* (* Plot H-v curve for last time step *)
	dat=ReadList["~/newton/oneD/out1",{Number,Number,Number,Number}];
	x=take2[dat]; z=take3[dat]; v=take4[dat];
	x=x/3.28/1000; z=z/3.28; v=v/3.28;
	gr1=ListPlot[xyjoin[x,z],PlotStyle->{t1,d0},FrameLabel->{"x (km)","Elevation (m) or Velocity (m/sec)"},PlotLabel->"Vel/Elev Response (last time step; z is solid)"];
	gr2=ListPlot[xyjoin[x,v],PlotStyle->{t1,d12}];
	gr4=Show[gr1,gr2,Axes->False]; disp[gr4];
(*	gr4=Show[gr4,bg->g10]; gr5=Show[gr5,bg->g10];*)
(*	lpr[gr4]; lpr[gr5];  *)
	temp={gr4,gr5,gr7};
	SetOptions[ListPlot,DisplayFunction:>$DisplayFunction];
	Return[temp]
	] *)
