sinuous.pas

PROGRAM SINUOUS(INFILE,INDATA,OUTFILE,AGEDAT,PARAMS,LIST,DETAIL,OUTDATA,INSEDIMENT,DEBUG);
{SINUOUS - River Meander Simulation
Copyright (C) 2021 Alan D. Howard
Developer can be contacted by ahoward@psi.edu
This program is free software; you can redistribute it and/or modify it under the terms of the
 GNU General Public License as published by the Free Software Foundation; either version 3
 of the License, or (at your option) any later version. 
This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
 without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 See the GNU General Public License for more details. 
You should have received a copy of the GNU General Public License along with this program;
 if not, write to the Free Software Foundation, Inc.,
 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. }

CONST
  NUM_CASES = 20;
  STAT_INTERVAL = 0.5;
  MIN_STAT = 1.25;
  MAX_STAT = 10.0;
  INFINSMALL = 1.0E-09;
  PI = 3.141593;
  ACONST = 0.077;
  NINTYDEGREE = 1.5707963;
  HALF_LIFE = 0.693147;
  ANGLE_CONVERT = 0.17453E-01;
  RADIAN = 0.0174532;
  MAX_FIT = 21;
  MAX_LAG = 30;
  MIN_LAG = 3;
  HALFPI = 1.5708;
  DVPOWER = 0.3333;
  YPOINTS = 200;
  XPOINTS = 700;
  JCENTER = 100;
  ISTART = 5;
  LEASTFRAC = 0.05;
  INITIALAGE = -100;
  ELEVATIONSCALE = 1000.0;
  DISTCATAGORIES = 20;
  IDISTCATAGORIES = 40;
  MAX_TERMS = 15;
  MAXRATECURVE=1.5;
  MAXPOINTS=5000;
  GRAVITY=9.81;
  SUBMERGED_SPECIFIC_GRAVITY = 1.65;
  WATER_DENSITY = 1000.0;
  SEDIMENT_POROSITY=0.5;
  SECONDS_PER_YEAR = 31.536E+06;
  SEDIMENT_GRADIENT_EXPONENT = 0.7;
  {SEDIMENT_TRANSPORT_EXPONENT = 1.5;}
  {TRANSPORT_FACTOR=8.0;}
  {FLOW_FRACTION=0.05;}
  SEDIMENT_1_EXPONENT=0.5;
  SEDIMENT_2_EXPONENT=0.3;
  {TRANSPORT_CRITICAL_DIM_SHEAR = 0.047;}
  MAX_LOCATIONS=10000;
  MANNING_STRICKLER_COEFFICIENT = 8.1;
  SKIN_FRACTION = 1.0;
  {UPWIND = 1.0;}
  GRAIN_ROUGHNESS_FACTOR = 5.0;
  {BEDMATERIAL_RELATIVE_DEPTH_FRACTION = 0.25;}
  {MAXIMUM_FAULTS = 2;}
  {RESET_NUMBER = 5;}
  INITIAL_ELEVATION = 20000;
TYPE
  COEFTYPE = ARRAY[1..MAX_TERMS,1..MAX_TERMS] OF REAL;
  VECTYPE = ARRAY[1..MAX_TERMS] OF REAL;
  IVECTYPE = ARRAY[1..MAX_TERMS] OF LONGINT;
  DISTTYPE = ARRAY[1..DISTCATAGORIES] OF REAL;
  AGETYPE = ARRAY [1..XPOINTS,1..YPOINTS] OF LONGINT;
  IDISTTYPE = ARRAY[0..IDISTCATAGORIES,0..IDISTCATAGORIES] OF LONGINT;
  DERIVTYPE = ARRAY[0..MAX_TERMS] OF REAL;
  BACKWATERTYPE = ARRAY[1..MAX_LOCATIONS] OF DOUBLE;
  BACKWATERBAD = ARRAY[1..MAX_LOCATIONS] OF BOOLEAN;
  DATALINK = ^DATAPOINT;
    DATAPOINT = RECORD
      XVALUE : REAL;
      YVALUE : REAL;
      ELAPSE : REAL;
      RESISTANCE : REAL;
      BANKDEPTH : REAL;
		BANKVELOCITY : REAL;
		CURVSENSE : LONGINT;
		LOCALGRADIENT : DOUBLE;
		LOCALDEPTH : DOUBLE;
		LOCALVELOCITY : DOUBLE;
		SEDIMENTFLUX : DOUBLE;
		BEDELEVATION : DOUBLE;
		BEDCHANGE : DOUBLE;
		LOCALCFRICTION : DOUBLE;
		UPSTREAM_DISTANCE : REAL ;
		MIGRATION_RATE : REAL;
		CURVATURE : REAL;
		CSCURVATURE : REAL;
		LOCALWIDTH : REAL;
		FROUDE : REAL;
      NEXT   : DATALINK;
    END;
  DATAFORM = ARRAY[0..MAX_LAG] OF REAL;
  STATARRAY = ARRAY[1..NUM_CASES] OF REAL;
  NSTATARRAY = ARRAY[1..NUM_CASES] OF LONGINT;
  BOUNDTYPE = ARRAY[1..11,1..2] OF REAL;
  POINTTYPE = ARRAY[1..11,1..2] of DATALINK;
  AREATYPE = ARRAY[1..10] OF REAL;
  SEPTYPE = ARRAY[1..MAXPOINTS] OF REAL;
  REGRESS_TYPE = ARRAY[1..3] OF REAL;
VAR
	 DAVG,DSD : REAL;
  SEPARATION : SEPTYPE;
  ENDSEGBOUND,NOWSEGBOUND : BOUNDTYPE;
  ENDPOINT,NOWPOINT : POINTTYPE;
  ENDAREA,NOWAREA : AREATYPE;
  NSEPARATION,FIRSTWRITE,LASTWRITE : LONGINT;
  RATEMULTUSE, NDIM, IDEPOSIT : LONGINT;
  USERATEMULT, USEDEPOSIT : BOOLEAN;
  MAXINTERVAL, DEPTHWEIGHT, VELOCITYWEIGHT : REAL;
  D1TERM, D2TERM : DERIVTYPE;
  DISTCOEFS,DUMMYCOEFS :COEFTYPE;
  DUMMYVEC, FIRSTCOEFS, SECONDCOEFS :VECTYPE;
  DETERMINANT : REAL;
  PIVOT : IVECTYPE;
  TDEPTH, CSCURVATURE,D1,D2,D3,CDEPTH, FDEPTH : REAL;
  DISTWEIGHTS : DISTTYPE;
  IRDIST : IDISTTYPE;
  BANK,KVAL,IDIST, AGE, ELEVATION : AGETYPE;
  ELEVRATE, FLOODRATE, BANKRATE, BANKDISTFACTOR : REAL;
  DEPTHFRACTION,MAXWIDTHCHECK, MAXELEV, ELEVCHANNEL : REAL;
  IRANGE,ELEVMAX, CHANNELELEV, WATERELEV, MEANDEPTH : LONGINT;
  DISCH_FACTOR, DISCHARGE : REAL;
  LOGMEDIAN, LOGMEAN, STDEV, LOGSTDEV : REAL;
  SEED, ISEED1 : REAL;
  STATSUM : STATARRAY;
  STATNUM : NSTATARRAY;
  FIRSTPRINT : LONGINT;
  TOTLAG : LONGINT;
  PLOTINT, CALCINT : LONGINT;
  PLOT_T,CALC_T : BOOLEAN;
  TWOWGT : REAL;
  FIRSTONE, SECOND_ONE : BOOLEAN;
  MAXNEWRATE,MAXRATE,RATEMULT,NOM_TIME,NOM_RATE,SLOPE_FACTOR : REAL;
  STREAM,ENDSTREAM , BASE,ENDBASE : DATALINK;
  OUTDATA,AGEDAT,PARAMS, OUTTIME,OUTFILE,INDATA , INFILE ,INEND,
     DISTLIST,LIST, DETAIL,OUTRAW1,OUTRAW2, OUTCAPT1, OUTCAPT2 , INSEDIMENT,VECDEBUG, AVG_VALS, JP_PARAMETERS: TEXT;
  ITERATION , UPDATEINT , PRINTINT , LOOPCHKINT : LONGINT;
  LOOPCHECKFACTOR , SCALEFACTOR, INTERVAL , RUNNUMBER : REAL;
  NRANGE, ITERNUM : LONGINT;
  MAXPLOT : LONGINT;
  RIVERSIZE,STARTRIVERSIZE,ENDRIVERSIZE : LONGINT;
  WIDTH , YAVERAGE : REAL;
  RESULT : LONGINT;
  ESCAPE , FORMFEED : CHAR;
  AVG_LENGTH, BIAS : REAL;
  MIN_ADJ_FACTOR : REAL;
  MAX_ADJ_FACTOR : REAL;
  LIMITWIDTH, BIASFACTOR, PRINTTIME : BOOLEAN;
  ABSCURVATURE , CURVATURE : REAL;
  RW , X1 , X2 , X3 , Y1 , Y2 , Y3 , X4, Y4 : REAL;
  XCHANGE , YCHANGE , XCENTER , YCENTER , RADIUS : REAL;
  UPLENGTH, DOWNLENGTH , NEWRATE  :  REAL;
  SENSE , NEWSENSE : LONGINT;
  START_TIME : REAL;
  STRAIGHT, FOLLOW : REAL;
  ANGLE, XS, YS : REAL;
  MIGRATION : REAL;
  XSENSE  : LONGINT;
  INPUT_PRINT : LONGINT;
  TOT_TIME : REAL;
  VALLEY_WIDTH : REAL;
  YUPPER, YLOWER : REAL;
  IPLACE : LONGINT;
  XC1,YC1,XC2,YC2,ANGLE1, ANGLE2 : REAL;
  PASTCURVATURE,  DISTANCE, PASTCS, PASTCVEL,PASTFVEL : DATAFORM;
  DETAILPRINT : LONGINT;
  DETAILPAR : BOOLEAN;
  INDEX1,INDEX2 : LONGINT;
  DIGITS : PACKED ARRAY [1..10] OF CHAR;
  NAMEOFFILE : PACKED ARRAY [1..12] OF CHAR;
  DESCRIPTION : PACKED ARRAY [1..70] OF CHAR;
{  DIGITS : STRING [10] ;
  NAMEOFFILE : STRING [12] ;
  DESCRIPTION : STRING [70] ;}
  MINITERATION : LONGINT;
  MEANLENGTH, MAX_DISTANCE : REAL;
  ERROR : REAL;
  DUMMY, IREADDEPOSIT : LONGINT;
  MAXTOT_RATE, TOT_RATE, OLDTOT_RATE : REAL;
  SIMILPAR : REAL;
  RES_TYPE : LONGINT;
  RES_AVERAGE : BOOLEAN;
  CFRICTION, WDRATIO,ACROSS,FROUDE,ASCROSS,DELTA,EXS20,EXS1,EXS : REAL;
  NOMAEXPONENT,AEXPONENT,NOMCSEXPONENT,
          CSEXPONENT,CSWEIGHT,CWEIGHT : REAL;
  MCONST,BETACONST,FNOMWGHT,FWEIGHT,GAMMA, F_RATE,C_RATE : REAL;
  FCONST1,FCONST2,FCONST3,FCONST4,FCONST5 : REAL;
  SKIPPOINTS, USENEWPARKER, FVELO_USE : LONGINT;
  DOSKIPPOINTS,NEWPARKER,USEFVELOCITY : BOOLEAN;
  VELOCITYFACTORUSE,MAXIMUMELEVATIONUSE,AVERAGEELEVATIONUSE,
	ANGLEFACTORUSE,PROBCUTOFFUSE,BANKELEVUSE,LINEARRATEUSE : LONGINT;
  USEVELOCITYFACTOR,USEMAXIMUMELEVATION,USEAVERAGEELEVATION,
	USEANGLEFACTOR,USEPROBCUTOFF,USEBANKELEV,USELINEARRATE : BOOLEAN;
	ANGLEFACTOR,VELOCITYFACTOR,ELEVATIONFACTOR,DISTANCEFACTOR : REAL;
	PROBMULT,MINPERPDIST, AGGRADERATE, DEPOSITDECAY : REAL;
	BANKERODE : ARRAY[1..XPOINTS,1..YPOINTS] OF REAL;
	PAST_SEDLOADS : ARRAY[1..50] OF DOUBLE;
	SEDLOADS_INDEX, NUMBER_SEDLOADS : LONGINT;
	NOMINAL_SEDFLUX : DOUBLE;
	ERODIBILITY,ERODENOMINAL,ERODEPLUG,ERODEVALLEY,ELEVPLUGMAX : REAL;
	PLUGMAXELEV,PLUGDISTANCE,READERODIBILITY : LONGINT;
	OUTBANKTYPE,INERODEFILE,OUTERODEFILE,DEBUG, TRANSFER,BACKREPORT,SEDIMENT_FLAG : TEXT;
	USEVALLEY,DOPLUGERODE,USESTICKY : BOOLEAN;
	{USEENDSTREAM : BOOLEAN;	}   {NEW 2007}
	ENDSTREAMINDEX : LONGINT; {NEW 2007}
  SUMABSRATE,SUMNUM : REAL;
  SEDIMENT_ROUTE, DO_INITIAL_GRADIENT , TRANSPORT_THRESHOLD, SEDIMENT_DIFFUSION, FIRST_USE : BOOLEAN;
  TRANSPORT_ITERATIONS, III : LONGINT;
  INITIAL_GRADIENT,CHANNEL_WIDTH : REAL;
  UPSTREAM_SEDIMENT_DISCHARGE, UPSTREAM_SEDIMENT_FLUX, SEDIMENT_SIZE : REAL;
  SEDIMENT_FACTOR, DEPTH_CONSTANT, VELOCITY_CONSTANT, WD_CONSTANT : REAL;
  FR_CONSTANT, AS_CONSTANT, DIM_SHEAR_CONST, SED_C1, SED_C2, EQUILIBRIUM_GRADIENT,TRANSPORT_INCREMENT : REAL;
  FLOWDEPTH, FLOWVELOCITY,E1,E2,E3,E4, THRESHOLD_GRADIENT, THRESHOLD_DEPTH : REAL;
  AVERAGE_DEPTH,NUM_DEPTH, INITIAL_DEPTH, YEARLY_SEDIMENT_FLUX : REAL;
  REG_DISTANCES,REG_ELEVATIONS : REGRESS_TYPE;
  SUMX,SUMY,SUMXY,SUMN,SUMX2,SUMY2, REGRESS_GRADIENT, DEFAULT_DEPTH, TRANSPORT_RATIO, QPERW, ROUGHNESS : DOUBLE;
  NETBEDCHANGE, ABSBEDCHANGE, NETBEDNUMBER, SHEAR_EXCEED_FACTOR, DOWNSTREAM_DEPTH, SEDIMENT_SPECIFIC_GRAVITY, THRESHOLD_SEDIMENT_SIZE : DOUBLE;
  H_DEPTH,H_PRED,DELX_PRED,UPDIST_PRED,DH_DXPRED,U_PRED,DH_DXCORR,TAU_PRED,EL_PRED,GRAD_PRED,FLUX_PRED,ELDIF_PRED, LCFRICTION, H_ELEV : BACKWATERTYPE;
  EL_PROVISIONAL : BACKWATERTYPE;
  VEC_POINTS  : LONGINT;
  NUMBER_FAULTS, NUMBER_OK : LONGINT;
  SUM_FLUX, NUMBER_SUM_FLUX, DIMENSIONLESS_FLUX_RATIO, SUM_DISCH_FACTOR,NUM_DISCH_FACTOR,AVG_DISCH_FACTOR, SUM_DISCHARGE : REAL ;
  MAX_CFWEIGHT, SUM_CFWEIGHT,NUM_CFWEIGHT,SUMSQ_CFWEIGHT, SD_CFWEIGHT : REAL;
  FAULT_STATE, END_SIMULATION, USE_PROVISIONAL, NORMAL_APPROXIMATION, CURVATURE_WIDTH_VARIATION, BEDLOAD_WIDTH_VARIATION, BEDLOAD_INDIRECT_VARIATION : BOOLEAN;
  NONLINEAR_BANK_EROSION : BOOLEAN;
  BANK_EROSION_EXPONENT : DOUBLE;
  BADPOINT : BACKWATERBAD;
  SEDIMENT_TRANSPORT_EXPONENT,TRANSPORT_FACTOR,FLOW_FRACTION,TRANSPORT_CRITICAL_DIM_SHEAR, UPWIND, DEFAULT_WDRATIO,DEFAULT_FROUDE, ITERATION_TIME : DOUBLE;
  BEDMATERIAL_RELATIVE_DEPTH_FRACTION, DISCHARGE_FACTOR, CRITICAL_USTAR, CURVATURE_WIDTH_SLOPE, BEDLOAD_WIDTH_CONSTANT, ABSCSCURVATURE : DOUBLE;
  MAXIMUM_FAULTS, RESET_NUMBER, SMOOTH_METHOD, SMOOTH_LENGTH, PROVISIONAL_USE, MAX_FAULT_ITERATIONS : LONGINT;
{*********************************************************************************}
{*********************************************************************************}
FUNCTION RAN3(DUMMY : LONGINT) : REAL;

VAR        TEMP : REAL;

BEGIN
        RAN3 := RANDOM;
       { TEMP := ISEED1 * 16807.0;
        ISEED1 := ABS(TEMP - 2147483647.0*TRUNC(TEMP/2147483647.0));
        RAN3 := ISEED1 * 4.6566128752458E-10; }

END;
{*********************************************************************************}
{*********************************************************************************}
FUNCTION NORMAL(DUMMY : LONGINT) : REAL;

VAR V1, V2, SS : REAL;
    J : LONGINT;

BEGIN

     REPEAT
  V1 := 2.0*RAN3(J)-1.0;
        V2 := 2.0*RAN3(J)-1.0;
        SS := V1*V1 + V2*V2;
     UNTIL SS<1.0;
     SS := SQRT(-2.0*LN(SS)/SS);
     NORMAL := V1*SS;
END;
{*********************************************************************************}
{*********************************************************************************}
FUNCTION LOGNORMAL(LOGMEDIAN, STDEV : REAL) : REAL;

VAR J : LONGINT;
    Y : REAL;
BEGIN
      J := 1;
      Y := NORMAL(J);
      LOGNORMAL := EXP(Y * STDEV) * LOGMEDIAN;
END;

{*********************************************************************************}
{*********************************************************************************}
PROCEDURE READ1PARAMETERS;
BEGIN
{First of three procedures to read the parameter values in the meander_parameters.prm file}

  RANDOMIZE;
  RIVERSIZE := 0;
  READ (INFILE, DESCRIPTION);
  READLN (INFILE);
  WRITELN (LIST,DESCRIPTION);
  READ (INFILE, ITERATION);
  READLN (INFILE);
  WRITELN (LIST,'ITERATIONS = ' ,ITERATION);
  READ (INFILE, UPDATEINT);
  READLN (INFILE);
  WRITELN (LIST,'UPDATE INTERVAL = ', UPDATEINT);
  READ (INFILE, PRINTINT);
  READLN (INFILE);
  WRITELN (LIST,'PRINT INTERVAL = ', PRINTINT);
  READ (INFILE, LOOPCHKINT);
  READLN (INFILE);
  WRITELN (LIST,'LOOP-CHECK INTERVAL = ', LOOPCHKINT:1);
  READ (INFILE, PLOTINT);
  READLN (INFILE);
  WRITELN (LIST,'DATA OUTPUT INTERVAL = ', PLOTINT:1);
   WRITELN ('DATA OUTPUT INTERVAL = ', PLOTINT:1);
  READ (INFILE, MINITERATION);
  READLN (INFILE);
        WRITELN('READ1PARAMETERS');
  WRITELN (LIST,'NO. OF ITERATIONS UNTIL DATA OUTPUT = ',MINITERATION);
  READ (INFILE, NOM_TIME);
  READLN (INFILE);
  WRITELN (LIST,'NOMINAL TIME INCREMENT IN YEARS = ', NOM_TIME:10);
  READ (INFILE, NOM_RATE);
  READLN (INFILE);
  WRITELN (LIST,'NOMINAL MAXIMUM BANK EROSION RATE IN METERS/YR= ', NOM_RATE:10);
  READ (INFILE, MAXINTERVAL);
  READLN (INFILE);
  WRITELN (LIST,'MAXIMUM TIME INTERVAL= ', MAXINTERVAL:10);
  READ(INFILE,RATEMULTUSE);
  READLN(INFILE);
  WRITELN(LIST,'SELECTION PARAMETER FOR SINUOSITY CORRECTION= ',RATEMULTUSE);
  IF RATEMULTUSE>0 THEN USERATEMULT := TRUE ELSE USERATEMULT := FALSE;
  READ (INFILE, SLOPE_FACTOR);
  READLN (INFILE);
  WRITELN (LIST,'SINUOSITY-GRADIENT EXPONENT = ', SLOPE_FACTOR:10);
  READ (JP_PARAMETERS, CHANNEL_WIDTH);
  READLN (JP_PARAMETERS);
  WRITELN('READ1PARAMETERS');
  WRITELN (LIST,'CHANNEL_WIDTH = ', CHANNEL_WIDTH:10);
  WRITELN ('CHANNEL_WIDTH = ', CHANNEL_WIDTH:10);
  READ (INFILE, MIN_ADJ_FACTOR);
  READLN (INFILE);
  WRITELN (LIST,'MIN_ADJ_FACTOR = ', MIN_ADJ_FACTOR:10);
  READ (INFILE, MAX_ADJ_FACTOR);
  READLN (INFILE);
  WRITELN (LIST,'MAX_ADJ_FACTOR = ', MAX_ADJ_FACTOR:10);
  MEANLENGTH := (MIN_ADJ_FACTOR + MAX_ADJ_FACTOR)/2.0;
  READ (INFILE, MAX_DISTANCE);
  READLN (INFILE);
  WRITELN (LIST,'MAX_DISTANCE = ', MAX_DISTANCE:10);
  MAX_DISTANCE := MAX_DISTANCE * CHANNEL_WIDTH;
        WRITELN('READ1PARAMETERS');
END; (*READ1PARAMETERS*)
{*********************************************************************************}
{*********************************************************************************}

PROCEDURE READ2PARAMETERS;
VAR USESEDIMENTROUTE : LONGINT;
BEGIN
{Second procedure to read parameter values}
 WRITELN('READ2PARAMETERS');
  READ (INFILE, NDIM);
  READLN (INFILE);
  WRITELN (LIST,'NO. TERMS IN DERIVATIVES= ', NDIM);
  READ (INFILE, USENEWPARKER);
  READLN (INFILE);
  WRITELN (LIST,'1 IF USE NEW PARKER= ', USENEWPARKER:10);
  IF USENEWPARKER>0 THEN NEWPARKER := TRUE ELSE NEWPARKER := FALSE;
  READ (INFILE, FVELO_USE);
  READLN (INFILE);
  WRITELN (LIST,'1 IF USE F VELOCITY= ', FVELO_USE:10);
  IF FVELO_USE>0 THEN USEFVELOCITY := TRUE ELSE USEFVELOCITY := FALSE;
  READ (JP_PARAMETERS, CFRICTION);
  READLN (JP_PARAMETERS);
  READ(SEDIMENT_FLAG, USESEDIMENTROUTE);
  READLN(SEDIMENT_FLAG);
  IF USESEDIMENTROUTE >0 THEN SEDIMENT_ROUTE := TRUE ELSE SEDIMENT_ROUTE :=FALSE;
  WRITELN(LIST,'USESEDIMENTROUTE (0=NO, 1=YES) =',USESEDIMENTROUTE:10);
  WRITELN (LIST,'COEFF. FRICTION= ', CFRICTION:10);
  READ (JP_PARAMETERS, WDRATIO);
  READLN (JP_PARAMETERS);
  WRITELN (LIST,'WIDTH-DEPTH RATIO = ', WDRATIO:10);
  READ (JP_PARAMETERS, FROUDE);
  READLN (JP_PARAMETERS);
  WRITELN (LIST,'FROUDE # = ', FROUDE:10);
  DEFAULT_DEPTH := WDRATIO/CHANNEL_WIDTH;
  DEFAULT_FROUDE := FROUDE;
  DEFAULT_WDRATIO := WDRATIO;
  READ (INFILE, MCONST);
  READLN (INFILE);
  WRITELN (LIST,'MCONST = ', MCONST:10);
  READ (INFILE, BETACONST);
  READLN (INFILE);
  WRITELN (LIST,'BETA = ', BETACONST:10);
  READ (INFILE, DEPTHWEIGHT);
  READLN (INFILE);
  WRITELN (LIST,'DEPTH PERTURBATION WEIGHT= ', DEPTHWEIGHT:10);
  READ (INFILE, VELOCITYWEIGHT);
  READLN (INFILE);
  WRITELN (LIST,'VELOCITY PERTURBATION WEIGHT= ', VELOCITYWEIGHT:10);
 { VELOCITYWEIGHT := VELOCITYWEIGHT + 1.0;}
  READ (INFILE, LOOPCHECKFACTOR);
  READLN (INFILE);
  WRITELN (LIST,'LOOPCHECKFACTOR = ', LOOPCHECKFACTOR:10);
  READ (INFILE, BIAS);
  READLN (INFILE);
  WRITELN (LIST,'DOWNSTREAM BIAS = ', BIAS:10);
  READ (INFILE, START_TIME);
  READLN (INFILE);
  WRITELN (LIST,'STARTING TIME = ', START_TIME:10);
  READ (INFILE, INPUT_PRINT);
  READLN (INFILE);
  WRITELN (LIST,'INPUT_PRINT = ', INPUT_PRINT);
  READ (INFILE, VALLEY_WIDTH);
  READLN (INFILE);
  WRITELN (LIST,'VALLEY_WIDTH = ', VALLEY_WIDTH:10);
  READ (INFILE, TWOWGT);
  READLN (INFILE);
  WRITELN (LIST,'2ND POINT WEIGHT = ', TWOWGT:10);
  READ (INFILE, DETAILPRINT);
  READLN (INFILE);
  WRITELN (LIST,'DETAILPRINT INDEX = ', DETAILPRINT);
  READ (INFILE, FIRSTWRITE);
  READLN (INFILE);
  WRITELN (LIST,'FIRST POINT FOR DETAILED OUTPUT = ', FIRSTWRITE);
  READ (INFILE, LASTWRITE);
  READLN (INFILE);
  WRITELN (LIST,'LAST POINT FOR DETAILED OUTPUT = ', LASTWRITE);
  READ (INFILE, FIRSTPRINT);
  READLN (INFILE);
  WRITELN (LIST,'FIRSTPRINT INDEX = ', FIRSTPRINT);
END;(*READ2PARAMETERS*)
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE READ3PARAMETERS;
VAR STICKYUSE, WIDTHVARIATION,USE_NONLINEAR_BANK_EROSION : LONGINT;
BEGIN
{Third procedure to read parameter values}
WRITELN('READ3PARAMETERS');
  READ(INFILE,IDEPOSIT);
  READLN(INFILE);
  WRITELN(LIST,'1 IF USE FLOODPLAIN DEPOSITION ELSE 0=',IDEPOSIT);
  IF IDEPOSIT = 1 THEN USEDEPOSIT := TRUE ELSE USEDEPOSIT := FALSE;
  READ(INFILE,IREADDEPOSIT);
  READLN(INFILE);
  IF IREADDEPOSIT>0 THEN
    WRITELN(LIST,'FLOODPLAIN DATA READ FROM FILE')
  ELSE
    WRITELN(LIST,'FLOODPLAIN DATA RESET');
  READ(INFILE,ELEVRATE);
  READLN(INFILE);
  WRITELN(LIST,'RATE OF FLOODPLAIN DEPOSITION=',ELEVRATE:10);
  ELEVRATE := ELEVRATE*ELEVATIONSCALE;
  READ(INFILE,MAXELEV);
  READLN(INFILE);
  WRITELN(LIST,'MAXIMUM FLOODPLAIN ELEVATION=',MAXELEV:10);
  ELEVMAX := ROUND(ELEVATIONSCALE * MAXELEV);
  READ(INFILE,ELEVCHANNEL);
  READLN(INFILE);
  CHANNELELEV := ROUND(ELEVATIONSCALE * ELEVCHANNEL);
{  WATERELEV := CHANNELELEV+ROUND((-CHANNELELEV+ELEVMAX)/2.0);
  MEANDEPTH := WATERELEV - CHANNELELEV;}
  WRITELN(LIST,'ELEVATION OF CHANNEL BOTTOM=',ELEVCHANNEL:10);
  READ(INFILE,DEPTHFRACTION);
  READLN(INFILE);
  WRITELN(LIST,'MEAN WATER LEVEL AS FRACTION OF ELEV. DIFF. ');
  WRITELN(LIST,'  FROM CHANN. BOTTOM TO MAX. FLOODPL. ELEV.'
        ,DEPTHFRACTION:10);
  WRITELN(LIST,'MEAN WATER LEVEL=',(ELEVCHANNEL+
       (-ELEVCHANNEL+MAXELEV)*DEPTHFRACTION):10);
  WATERELEV := CHANNELELEV+ROUND((-CHANNELELEV+ELEVMAX)*DEPTHFRACTION);
  MEANDEPTH := WATERELEV - CHANNELELEV;
  READ(INFILE,FLOODRATE);
  READLN(INFILE);
  WRITELN(LIST,'RELATIVE FLOODPLAIN DEPOSITION RATE=',FLOODRATE:10);
  READ(INFILE,BANKRATE);
  READLN(INFILE);
  WRITELN(LIST,'RELATIVE OVERBANK DEPOSITION RATE=',BANKRATE:10);
  READ(INFILE,BANKDISTFACTOR);
  READLN(INFILE);
  WRITELN(LIST,'DISTANCE DECAY FACTOR FOR OVERBANK DEPOSITION=',
     BANKDISTFACTOR:10);
  BANKDISTFACTOR := -LN(0.5) / BANKDISTFACTOR;
  MAXWIDTHCHECK := -LN(0.05) / BANKDISTFACTOR;
  WRITELN(LIST,'MAXIMUM CHECK DISTANCE=',MAXWIDTHCHECK:10);
	READ(INFILE,PROBCUTOFFUSE);
	READLN(INFILE);
	IF PROBCUTOFFUSE>0 THEN USEPROBCUTOFF:= TRUE ELSE USEPROBCUTOFF :=FALSE;
	WRITELN(LIST,'USE PROBABILITY CUTOFF=',PROBCUTOFFUSE);
	READ(INFILE,VELOCITYFACTORUSE);
	IF VELOCITYFACTORUSE>0 THEN USEVELOCITYFACTOR := TRUE ELSE
		USEVELOCITYFACTOR := FALSE;
	READLN(INFILE);
	WRITELN(LIST,'VELOCITYFACTORUSE=',VELOCITYFACTORUSE);
	READ(INFILE,ANGLEFACTORUSE);
	IF ANGLEFACTORUSE>0 THEN USEANGLEFACTOR := TRUE ELSE
		USEANGLEFACTOR := FALSE;
	READLN(INFILE);
	WRITELN(LIST,'ANGLEFACTORUSE=',ANGLEFACTORUSE);

	READ(INFILE,MAXIMUMELEVATIONUSE);
	IF MAXIMUMELEVATIONUSE>0 THEN USEMAXIMUMELEVATION := TRUE ELSE
		USEMAXIMUMELEVATION := FALSE;
	READLN(INFILE);
	WRITELN(LIST,'MAXIMUMELEVATIONUSE=',MAXIMUMELEVATIONUSE);
	READ(INFILE,AVERAGEELEVATIONUSE);
	IF AVERAGEELEVATIONUSE>0 THEN USEAVERAGEELEVATION := TRUE ELSE
		USEAVERAGEELEVATION := FALSE;
	READLN(INFILE);
	WRITELN(LIST,'AVERAGEELEVATIONUSE=',AVERAGEELEVATIONUSE);
	READ(INFILE,ANGLEFACTOR);
	READLN(INFILE);
	WRITELN(LIST,'ANGLEFACTOR=',ANGLEFACTOR:10);
	READ(INFILE,VELOCITYFACTOR);
	READLN(INFILE);
	WRITELN(LIST,'VELOCITYFACTOR=',VELOCITYFACTOR:10);
	READ(INFILE,ELEVATIONFACTOR);
	READLN(INFILE);
	WRITELN(LIST,'ELEVATIONFACTOR=',ELEVATIONFACTOR:10);
	ELEVATIONFACTOR := ELEVATIONFACTOR / ELEVATIONSCALE;
	READ(INFILE,DISTANCEFACTOR);
	READLN(INFILE);
	WRITELN(LIST,'DISTANCEFACTOR=',DISTANCEFACTOR:10);
        DISTANCEFACTOR := -LN(0.5) / DISTANCEFACTOR;
	READ(INFILE,MINPERPDIST);
	READLN(INFILE);
	WRITELN(LIST,'MINPERPDIST=',MINPERPDIST:10);
	READ(INFILE,AGGRADERATE);
	READLN(INFILE);
	WRITELN(LIST,'AGGRADATION RATE=',AGGRADERATE:10);
	AGGRADERATE:= AGGRADERATE*ELEVATIONSCALE;
	READ(INFILE,DEPOSITDECAY);
	READLN(INFILE);
	WRITELN(LIST,'DEPOSITDECAY=',DEPOSITDECAY:10);
	DEPOSITDECAY := DEPOSITDECAY / SQR(ELEVATIONSCALE);
{	DEPOSITDECAY := (-LN(0.5) /DEPOSITDECAY)/ELEVATIONSCALE;
	WRITELN(LIST,'SCALED DEPOSITION FACTOR=',DEPOSITDECAY:10);}
	READ(INFILE,BANKELEVUSE);
	READLN(INFILE);
	WRITELN(LIST,'INDEX FOR INITIAL BANK ELEV (0=MEAN BED, 1=BANKELEV): ',BANKELEVUSE);
	IF BANKELEVUSE > 0 THEN USEBANKELEV := TRUE ELSE USEBANKELEV := FALSE;
	READ(INFILE,LINEARRATEUSE);
	READLN(INFILE);
	WRITELN(LIST,'INDEX FOR DEPOSITION RATE LAW (1=LINEAR, 0=ESQUARE): ',LINEARRATEUSE);
	IF LINEARRATEUSE >0 THEN USELINEARRATE := TRUE ELSE USELINEARRATE := FALSE;
	IF NOT USELINEARRATE THEN ELEVRATE := ELEVRATE*ELEVATIONSCALE;
	READ(INFILE,SKIPPOINTS);
	READLN(INFILE);
	WRITELN(LIST,'INDEX FOR SKIPPING POINTS DURING CUTOFF CHECK (0=NO, 1=YES)');
	IF SKIPPOINTS>0 THEN DOSKIPPOINTS := TRUE ELSE DOSKIPPOINTS := FALSE;
	READ(INFILE,PROBMULT);
	READLN(INFILE);
	WRITELN(LIST,'MULTIPLICATIVE CUTOFF PROBABILITY FACTOR=',PROBMULT:10:4);
	READ(INFILE,ERODENOMINAL);
	READLN(INFILE);
	WRITELN(LIST,'NOMINAL BANK ERODIBILITY=',ERODENOMINAL:12:5);
	READ(INFILE,ERODEPLUG);
	READLN(INFILE);
	WRITELN(LIST,'BANK ERODIBILITY OF CLAY PLUGS=',ERODEPLUG:12:5);
	READ(INFILE,ERODEVALLEY);
	READLN(INFILE);
	WRITELN(LIST,'BANK ERODIBILITY OF VALLEY WALLS=',ERODEVALLEY:12:5);
	IF ERODEVALLEY<>ERODENOMINAL THEN USEVALLEY := TRUE ELSE USEVALLEY := FALSE;
	READ(INFILE,PLUGDISTANCE);
	READLN(INFILE);
	WRITELN(LIST,'DISTANCE FROM CHANNEL TO BEGINNING OF CLAY PLUG=',PLUGDISTANCE);
	READ(INFILE,READERODIBILITY);
	READLN(INFILE);
	WRITELN(LIST,'INDEX FOR READING ERODIBILITIES FROM FILE=',READERODIBILITY);
		READ(INFILE,STICKYUSE);
	READLN(INFILE);
	WRITELN(LIST,'INDEX FOR STICKYNESS OF VALLEY WALLS AND PLUGS (0=LESS,1=MORE)',STICKYUSE);
	IF STICKYUSE>0 THEN USESTICKY := TRUE ELSE USESTICKY := FALSE;
	read(infile,critical_ustar);
	readln(infile);
	writeln(list,'Critical value of velocity perturbation for bank erosion =',critical_ustar:12);
	READ(INFILE,WIDTHVARIATION); {0=CONSTANT WIDTH, 1= INDIRECT INFLUENCE, 2=DIRECT INFLUENCE}
	READLN(INFILE);
	IF (WIDTHVARIATION=1) THEN
	BEGIN
	    CURVATURE_WIDTH_VARIATION := TRUE;
		WRITELN(LIST,' CHANNEL WIDTH VARIES WITH CURVATURE');
	END ELSE
	BEGIN
	    CURVATURE_WIDTH_VARIATION := FALSE;
		WRITELN(LIST,'CHANNEL WIDTH IS CONSTANT');
	END;
	READ(INFILE,CURVATURE_WIDTH_SLOPE);
	READLN(INFILE);
	IF CURVATURE_WIDTH_VARIATION THEN WRITELN(LIST,'WIDTH-CURVATURE CONSTANT =',CURVATURE_WIDTH_SLOPE:15);
	READ(INFILE,USE_NONLINEAR_BANK_EROSION);
	READLN(INFILE);
	READ(INFILE,BANK_EROSION_EXPONENT);
	READLN(INFILE);
	IF USE_NONLINEAR_BANK_EROSION>0 THEN
	    NONLINEAR_BANK_EROSION := TRUE
	ELSE
	    NONLINEAR_BANK_EROSION := FALSE;
	IF NONLINEAR_BANK_EROSION THEN
	   WRITELN(LIST,' NON-LINEAR BANK EROSION WITH EXPONENT =',BANK_EROSION_EXPONENT);

END;
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE READ_ROUTE_PARAMETERS;
VAR USE_INITIAL_EQUILIBRIUM_GRADIENT, USE_TRANSPORT_THRESHOLD, USE_SEDIMENT_DIFFUSION,USE_GRAIN_ROUGHNESS, USE_NORMAL_APPROXIMATION : LONGINT;
BEDLOAD_WIDTH_USE : LONGINT;
BEGIN
{Reads the sediment routing parameters from the sediment_parameters.prm file}
   READ(INSEDIMENT,CHANNEL_WIDTH);
   READLN(INSEDIMENT);
   CHANNEL_WIDTH := CHANNEL_WIDTH;
   READ(INSEDIMENT,INITIAL_GRADIENT);
   READLN(INSEDIMENT);
   READ(INSEDIMENT, DISCHARGE);
   READLN(INSEDIMENT);
   READ(INSEDIMENT, SHEAR_EXCEED_FACTOR);
   READLN(INSEDIMENT);
   READ(INSEDIMENT, TRANSPORT_INCREMENT);
   READLN(INSEDIMENT);

	SEDIMENT_SPECIFIC_GRAVITY :=  SUBMERGED_SPECIFIC_GRAVITY + 1.0;
   TRANSPORT_ITERATIONS := ROUND(NOM_TIME/TRANSPORT_INCREMENT);
   WRITELN(LIST,'INITIAL TRANSPORT ITERATIONS =',TRANSPORT_ITERATIONS:10);
   QPERW := DISCHARGE/CHANNEL_WIDTH;
   INTERVAL := NOM_TIME;
   READ(INSEDIMENT,USE_GRAIN_ROUGHNESS);
	READLN(INSEDIMENT);
	READ(INSEDIMENT,ROUGHNESS);
	READLN(INSEDIMENT);
	READ(INSEDIMENT,SEDIMENT_TRANSPORT_EXPONENT);
	READLN(INSEDIMENT);
	READ(INSEDIMENT,TRANSPORT_FACTOR);
	READLN(INSEDIMENT);
	READ(INSEDIMENT,FLOW_FRACTION);
	READLN(INSEDIMENT);
	READ(INSEDIMENT,TRANSPORT_CRITICAL_DIM_SHEAR);
	READLN(INSEDIMENT);
	READ(INSEDIMENT,UPWIND);
	READLN(INSEDIMENT);
	READ(INSEDIMENT,BEDMATERIAL_RELATIVE_DEPTH_FRACTION);
	READLN(INSEDIMENT);
	READ(INSEDIMENT,MAXIMUM_FAULTS);
	READLN(INSEDIMENT);
	READ(INSEDIMENT,RESET_NUMBER);
	READLN(INSEDIMENT);
	READ(INSEDIMENT,SMOOTH_METHOD);
	READLN(INSEDIMENT);
	READ(INSEDIMENT,SMOOTH_LENGTH);
	READLN(INSEDIMENT);
	READ(INSEDIMENT,PROVISIONAL_USE);
	READLN(INSEDIMENT);
	READ(INSEDIMENT,MAX_FAULT_ITERATIONS);
	READLN(INSEDIMENT);
	READ(INSEDIMENT,USE_NORMAL_APPROXIMATION);
	READLN(INSEDIMENT);
	ROUGHNESS := ROUGHNESS/100.0;
	IF USE_NORMAL_APPROXIMATION>0 THEN
	    NORMAL_APPROXIMATION := TRUE
	ELSE
	    NORMAL_APPROXIMATION :=FALSE;
	IF USE_GRAIN_ROUGHNESS>0 THEN
	BEGIN
	     ROUGHNESS := SEDIMENT_SIZE*GRAIN_ROUGHNESS_FACTOR;
		 WRITELN(LIST,'USING GRAIN SIZE ROUGHNESS, COMPUTED ROUGHNESS = ',ROUGHNESS:15);
	END ELSE
	     WRITELN(LIST,'USING INPUT ROUGHNESS VALUE = ',ROUGHNESS:15);
	READ(INSEDIMENT,BEDLOAD_WIDTH_USE);
	READLN(INSEDIMENT);
	IF BEDLOAD_WIDTH_USE =1 THEN
	BEGIN
	    BEDLOAD_INDIRECT_VARIATION := TRUE;
		BEDLOAD_WIDTH_VARIATION := FALSE;
		WRITELN(LIST,'MIGRATION VARIES WITH BEDLOAD FLUX');
	END ELSE
	    IF BEDLOAD_WIDTH_USE = 2 THEN
		BEGIN
		BEDLOAD_WIDTH_VARIATION := TRUE;
		BEDLOAD_INDIRECT_VARIATION := FALSE;
		WRITELN(LIST,' WIDTH DEPENDS UPON SEDIMENT FLUX');
		END ELSE
	       BEGIN
	          BEDLOAD_WIDTH_VARIATION := FALSE;
			  BEDLOAD_INDIRECT_VARIATION := FALSE;
		      WRITELN(LIST,'CHANNEL WIDTH IS CONSTANT AND MIGRATION DOES NOT DEPEND UPON SEDIMENT FLUX');
	       END;
    READ(INSEDIMENT,BEDLOAD_WIDTH_CONSTANT);
	READLN(INSEDIMENT);
	IF BEDLOAD_WIDTH_VARIATION THEN WRITELN(LIST,'BEDLOAD WIDTH CONSTANT =',BEDLOAD_WIDTH_CONSTANT:15);
   IF TRANSPORT_ITERATIONS<1 THEN TRANSPORT_ITERATIONS :=1;
   DOWNSTREAM_DEPTH := EXP(0.3*LN((EXP(0.3333*LN(ROUGHNESS))*SQR(QPERW))/(SQR(MANNING_STRICKLER_COEFFICIENT)*GRAVITY*INITIAL_GRADIENT)));
   CFRICTION := EXP(0.3333*LN(ROUGHNESS))/(SQR(MANNING_STRICKLER_COEFFICIENT)*EXP(0.3333*LN(DOWNSTREAM_DEPTH)));
   SEDIMENT_SIZE := CFRICTION*SQR(QPERW)/(SQR(DOWNSTREAM_DEPTH)*SUBMERGED_SPECIFIC_GRAVITY*GRAVITY*TRANSPORT_CRITICAL_DIM_SHEAR*SHEAR_EXCEED_FACTOR);
   THRESHOLD_SEDIMENT_SIZE :=ROUGHNESS/GRAIN_ROUGHNESS_FACTOR;
   UPSTREAM_SEDIMENT_FLUX := SQRT(SUBMERGED_SPECIFIC_GRAVITY*GRAVITY*SEDIMENT_SIZE)*SEDIMENT_SIZE*
         TRANSPORT_FACTOR*EXP(SEDIMENT_TRANSPORT_EXPONENT*LN(TRANSPORT_CRITICAL_DIM_SHEAR*(SKIN_FRACTION*SHEAR_EXCEED_FACTOR-1.0)));
   YEARLY_SEDIMENT_FLUX := SECONDS_PER_YEAR * UPSTREAM_SEDIMENT_FLUX*FLOW_FRACTION;
   UPSTREAM_SEDIMENT_DISCHARGE := UPSTREAM_SEDIMENT_FLUX * CHANNEL_WIDTH;
   	WRITELN(LIST,'**** SEDIMENT TRANSPORT PARAMETERS ****');
	WRITELN(LIST,'CHANNEL WIDTH =',CHANNEL_WIDTH:15);
	WRITELN(LIST,'INITIAL GRADIENT =',INITIAL_GRADIENT:15);
	WRITELN(LIST,'SHEAR EXCEED FACTOR (TAU/TAUC) =',SHEAR_EXCEED_FACTOR:15);
	WRITELN(LIST,'DOWNSTREAM FLOW DEPTH =',DOWNSTREAM_DEPTH:15);
	WRITELN(LIST,'FRICTION COEFFICIENT (USTAR/U)^2 =',CFRICTION:15);
	WRITELN(LIST,'SEDIMENT SIZE =',SEDIMENT_SIZE:15);
	WRITELN(LIST, 'SEDIMENT SIZE FOR THRESHOLD CONDITIONS =',THRESHOLD_SEDIMENT_SIZE:15);
	WRITELN(LIST,'UPSTREAM SEDIMENT DISCHARGE =',UPSTREAM_SEDIMENT_DISCHARGE:15);
	WRITELN(LIST,'UPSTREAM SEDIMENT FLUX =',UPSTREAM_SEDIMENT_FLUX:15);
	WRITELN(LIST,'TRANSPORT INCREMENT, YEARS =',TRANSPORT_INCREMENT:15);
	{WRITELN(LIST,'USE INITIAL EQUILIBRIUM GRADIENT =',USE_INITIAL_EQUILIBRIUM_GRADIENT:15);}
	WRITELN(LIST,'INITIAL ITERATIONS PER SIMULATION TIME STEP =',TRANSPORT_ITERATIONS:15);
	WRITELN(LIST,'SEDIMENT POROSITY =',SEDIMENT_POROSITY:15);
	WRITELN(LIST,'GRAVITY=',GRAVITY:15);
	WRITELN(LIST,'SUBMERGED SEDIMENT SPECIFIC GRAVITY =',SUBMERGED_SPECIFIC_GRAVITY:15);
	WRITELN(LIST,'WATER DENSITY =',WATER_DENSITY:15);
	WRITELN(LIST,'SEDIMENT TRANSPORT EXPONENT =',SEDIMENT_TRANSPORT_EXPONENT:15);
	WRITELN(LIST,'TRANSPORT FACTOR =',TRANSPORT_FACTOR:15);
	WRITELN(LIST,'CRITICAL DIMENSIONLESS SHEAR =',TRANSPORT_CRITICAL_DIM_SHEAR:15);
	DEPTH_CONSTANT :=EXP(0.33333*LN(CFRICTION/GRAVITY))*EXP(0.66667*LN(DISCHARGE/CHANNEL_WIDTH));
 	THRESHOLD_GRADIENT := EXP(1.5*LN(TRANSPORT_CRITICAL_DIM_SHEAR*(SEDIMENT_SPECIFIC_GRAVITY-1.0)*SEDIMENT_SIZE/DEPTH_CONSTANT));
	THRESHOLD_DEPTH := DEPTH_CONSTANT/EXP(0.3333*LN(THRESHOLD_GRADIENT));
	WRITELN(LIST,'THRESHOLD CHANNEL GRADIENT =',THRESHOLD_GRADIENT:15);
	WRITELN(LIST,' THRESHOLD CHANNEL DEPTH=',THRESHOLD_DEPTH:15);
	WRITELN(LIST,'WATER DISCHARGE PER UNIT WIDTH = ',QPERW:15);
	WRITELN(LIST,'EXPONENT IN SEDIMENT TRANSPORT RELATIONSHIP =',SEDIMENT_TRANSPORT_EXPONENT:15);
	WRITELN(LIST,'MULTIPLICATIVE FACTOR IN SEDIMENT TRANSPORT RELATIONSHIP =',TRANSPORT_FACTOR:15);
	WRITELN(LIST,'FRACTION OF YEAR WITH EFFECTIVE SEDIMENT TRANSPORT = ',FLOW_FRACTION:15);
	WRITELN(LIST,'CRITICAL DIMENSIONLESS SHEAR =',TRANSPORT_CRITICAL_DIM_SHEAR:15);
	WRITELN(LIST,'UPWINDING COEFFICIENT IN SEDIMENT TRANSPORT (0.5 -- 1.0) =',UPWIND:15);
	WRITELN(LIST,'RELATIVE DEPTH OF BED SEDIMENT IN BANKS/POINT BARS =',BEDMATERIAL_RELATIVE_DEPTH_FRACTION:15);
	WRITELN(LIST,'NUMBER OF SUCCESSIVE NEGATIVE PREDICTED DEPTHS TO TERMINATE PROGRAM =',MAXIMUM_FAULTS:8);
	WRITELN(LIST,'NUMBER OF ITERATIONS WITH NO NEGATIVE PREDICTED DEPTHS TO CLEAR FAULTS =',RESET_NUMBER:8);
	WRITELN(LIST,'PROFILE SMOOTHING METHOD (0=NONE,1=FLATTEN,2=REGRESS) =',SMOOTH_METHOD:6);
	WRITELN(LIST,'INITIAL RANGE OF STREAM POINTS FOR SMOOTHING, N: RANGE=+N TO -N =',SMOOTH_LENGTH:6);
	IF PROVISIONAL_USE>0 THEN USE_PROVISIONAL := TRUE ELSE USE_PROVISIONAL := FALSE;
	WRITELN(LIST,'USE PROVISIONAL ELEVATIONS IN BACKWATER CALCULATION, (0=NO, 1=YES) =',PROVISIONAL_USE:5);
	WRITELN(LIST,'MAXIMUM NUMBER OF SMOOTHING ITERATIONS TO USE =',MAX_FAULT_ITERATIONS:8);
	WRITELN(LIST,'USE THE NORMAL APPROXIMATION FOR FLOW IN CHANNEL (0=NO,1=YES) = ',USE_NORMAL_APPROXIMATION:8);

	    USERATEMULT :=FALSE;
        DO_INITIAL_GRADIENT := TRUE;
        {writeln('end of sediment parameter readin');}


 END; (*READROUTEPARAMETERS*)
{*********************************************************************************}
{*********************************************************************************}

PROCEDURE FINDDISTANCE (VAR DISTANCE:REAL; X1, Y1, X2, Y2 :
            REAL);

BEGIN
{calculates the straight-line distance between two nodes from their X and Y coordinates}
  DISTANCE := SQRT ((SQR(X2 - X1)) + (SQR(Y2 - Y1)));

END;  (* FINDDISTANCESTANCE *)
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE WRITE_DEBUG;
VAR XD : REAL;
BEGIN
{Writes values of bed-elevation, depth, channel gradient, sediment flux, and friction of all active nodes to the file debug.dat.}
   ASSIGN(DEBUG,'DEBUG.DAT');
   REWRITE(DEBUG);
   STREAM := BASE;
   XD :=0.0;
   X1 := STREAM^.XVALUE;
   Y1 := STREAM^.YVALUE;
     WRITELN(DEBUG,STREAM^.BEDELEVATION:15,',',STREAM^.LOCALDEPTH:15,',',STREAM^.LOCALGRADIENT:15,',',STREAM^.SEDIMENTFLUX:15,',',
   STREAM^.LOCALCFRICTION:15,',',XD:15);
  WHILE (STREAM^.NEXT <> NIL) DO
   BEGIN
   STREAM := STREAM^.NEXT;
   X2 := STREAM^.XVALUE;
   Y2 := STREAM^.YVALUE;
   FINDDISTANCE(XD,X1,Y1,X2,Y2);
   WRITELN(DEBUG,STREAM^.BEDELEVATION:15,',',STREAM^.LOCALDEPTH:15,',',STREAM^.LOCALGRADIENT:15,',',STREAM^.SEDIMENTFLUX:15,',',
   STREAM^.LOCALCFRICTION:15,',',XD:15);
   X1 :=X2;
   Y1 := Y2;
   END;
   CLOSE(DEBUG);
END; (*WRITE_DEBUG*)
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE WRITE_VECTOR;
VAR I : LONGINT;
BEGIN
{Writes debugging values of a vector of points when smoothing the stream profile to the file vector.dat.
See procedure for a list of vector variables written.}
ASSIGN(VECDEBUG,'VECTOR.DAT');
REWRITE(VECDEBUG);
FOR I := 1 TO VEC_POINTS DO
   WRITELN(VECDEBUG, EL_PRED[I]:15,',',H_ELEV[I]:15,',',H_DEPTH[I]:15,',',U_PRED[I]:15,',',
      LCFRICTION[I]:15,',',FLUX_PRED[I]:15,',',ELDIF_PRED[I]:15,',',DELX_PRED[I]:15,',',
	  GRAD_PRED[I]:15,',',UPDIST_PRED[I]);
CLOSE(VECDEBUG);
END; (*WRITEVECTOR*)
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE SMOOTH_PROFILE( ILOC, N, ncall : LONGINT);
VAR I, IUP, IDOWN : LONGINT;
    CUM_DISTANCE, TOTAL_DISTANCE, EUP, EDOWN, ESLOPE,EINTERCEPT, EORIG : DOUBLE;
BEGIN
{This replaces the elevation of a node by fitting a constant-gradient profile starting
 from the elevation of a node N nodes upstream to the elevation of a node N nodes downstream
 and replacing the elevation of the central node by the value given by the constant-gradient profile.
 Used in sediment routing}
   IUP := ILOC - N;
   IF IUP < 1 THEN IUP :=1;
   IDOWN := ILOC + N;
   IF IDOWN > VEC_POINTS THEN IDOWN := VEC_POINTS;
   TOTAL_DISTANCE := 0.0;
   FOR I := IUP TO IDOWN -1 DO
      TOTAL_DISTANCE := TOTAL_DISTANCE + DELX_PRED[I];
   IF USE_PROVISIONAL THEN
   BEGIN
      EUP :=EL_PROVISIONAL[IUP];
      EDOWN := EL_PROVISIONAL[IDOWN];
   END ELSE
   BEGIN
      EUP :=EL_PRED[IUP];
      EDOWN := EL_PRED[IDOWN];
   END;
   EINTERCEPT := EUP;
   ESLOPE :=(EDOWN-EUP)/(TOTAL_DISTANCE);
   CUM_DISTANCE := DELX_PRED[IUP];
   WRITELN(LIST,'iteration=',iternum:10,' I=',ILOC:6,' NSMOOTH=',N:8,' ncall=',ncall:8);
    WRITELN('iteration=',iternum:10,' I=',ILOC:6,' NSMOOTH=',N:8,' ncall=',ncall:8);
   FOR I := IUP+1 TO IDOWN-1 DO
   BEGIN
      EORIG:=EL_PROVISIONAL[I];
      EL_PROVISIONAL[I] := EINTERCEPT + ESLOPE * CUM_DISTANCE;
	  WRITELN(LIST,' IUP=',IUP:8,' EUP=',EUP:15,' IDOWN=',IDOWN:8,' EDOWN=',EDOWN:15,' EORIG=',EORIG:15,' ENEW=',EL_PROVISIONAL[I]:15);
	   WRITELN(' IUP=',IUP:8,' EUP=',EUP:15,' IDOWN=',IDOWN:8,' EDOWN=',EDOWN:15,' EORIG=',EORIG:15,' ENEW=',EL_PROVISIONAL[I]:15);
	  CUM_DISTANCE := CUM_DISTANCE+DELX_PRED[I];
   END;
END; (*SMOOTH_PROFILE*)
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE REGRESS_PROFILE( ILOC, N, NCALL : LONGINT);
VAR I, IUP, IDOWN : LONGINT;
    TOTAL_DISTANCE, SUMN,SUMX,SUMY,SUMXY,SUMX2,SUMY2, ESLOPE,EINTERCEPT, EORIG : DOUBLE;
BEGIN
{This replaces the elevation of a node by fitting a straight line regression of the elevation of the N nodes upstream and N nodes downstream
 of the location and replacing the elevation of the central node by the value
 given by the regression. Used in sediment routing}
   IUP := ILOC - N;
   IF IUP < 1 THEN IUP :=1;
   IDOWN := ILOC + N;
   IF IDOWN > VEC_POINTS THEN IDOWN := VEC_POINTS;
   TOTAL_DISTANCE := 0.0;
   SUMX :=0.0;
   SUMY := 0.0;
   SUMX2 := 0.0;
   SUMY2 := 0.0;
   SUMXY :=0.0;
   SUMN := 0.0;
   FOR I := IUP TO IDOWN DO
   BEGIN
      SUMX := SUMX+ TOTAL_DISTANCE;
	  SUMX2 := SUMX2 + SQR(TOTAL_DISTANCE);
	  IF USE_PROVISIONAL THEN
	  BEGIN
	      SUMY := SUMY + EL_PROVISIONAL[I];
	      SUMY2 := SUMY2+ SQR(EL_PROVISIONAL[I]);
	      SUMXY := SUMXY + TOTAL_DISTANCE*EL_PROVISIONAL[I];
	  END ELSE
	  BEGIN
	      SUMY := SUMY + EL_PRED[I];
	      SUMY2 := SUMY2+ SQR(EL_PRED[I]);
	      SUMXY := SUMXY + TOTAL_DISTANCE*EL_PRED[I];
	  END;
	  SUMN := SUMN + 1.0;
      TOTAL_DISTANCE := TOTAL_DISTANCE + DELX_PRED[I];
   END;

   ESLOPE :=(SUMXY-SUMX*SUMY/SUMN)/(SUMX2-SQR(SUMX)/SUMN);
   EINTERCEPT :=SUMY/SUMN -ESLOPE*SUMX/SUMN;
   TOTAL_DISTANCE := 0.0;
   WRITELN(LIST,'iteration=',iternum:10,' I=',ILOC:6,' NSMOOTH=',N:8,' ncall =',ncall:8);
   WRITELN('iteration=',iternum:10,' I=',ILOC:6,' NSMOOTH=',N:8,' ncall =',ncall:8);
   FOR I := IUP TO IDOWN DO
   BEGIN
      EORIG := EL_PROVISIONAL[I];
      EL_PROVISIONAL[I] := EINTERCEPT + ESLOPE * TOTAL_DISTANCE;
	  WRITELN(LIST,' I=',I:6,' X=',TOTAL_DISTANCE:15,' EORIG =',EORIG:15,' ENEW=',EL_PROVISIONAL[I]:15);
	  WRITELN(' I=',I:6,' X=',TOTAL_DISTANCE:15,' EORIG =',EORIG:15,' ENEW=',EL_PROVISIONAL[I]:15);
	  TOTAL_DISTANCE := TOTAL_DISTANCE + DELX_PRED[I];
   END;
END; (*REGRESS_PROFILE*)
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE MAKE_VECTORS;
VAR I : LONGINT;
    DISTANCE : DOUBLE;
BEGIN
{This constructs vectors of distance downstream, bed elevation, and channel gradient for the river.
Used in sediment routing}
{WRITELN('MAKE_VECTORS');}
   STREAM := BASE;
   I :=0;
   WHILE STREAM^.NEXT^.NEXT <> NIL DO
   BEGIN
       I :=I+1;
       X1 := STREAM^.XVALUE;
       X2 := STREAM^.NEXT^.XVALUE;
       Y1 := STREAM^.YVALUE;
       Y2 := STREAM^.NEXT^.YVALUE;
	   FINDDISTANCE(DISTANCE,X1,Y1,X2,Y2);
       DELX_PRED[I] := DISTANCE;
       EL_PRED[I] := STREAM^.BEDELEVATION;
       STREAM := STREAM^.NEXT;
   END;
   I :=I+1;
       X1 := STREAM^.XVALUE;
       X2 := STREAM^.NEXT^.XVALUE;
       Y1 := STREAM^.YVALUE;
       Y2 := STREAM^.NEXT^.YVALUE;
	   FINDDISTANCE(DISTANCE,X1,Y1,X2,Y2);
       DELX_PRED[I] := DISTANCE;
       EL_PRED[I] := STREAM^.BEDELEVATION;
   STREAM := STREAM^.NEXT;
   I := I+1;
   VEC_POINTS := I;
   DELX_PRED[I] := DELX_PRED[I-1];
   UPDIST_PRED[I] := 0.0;
   EL_PRED[I] := STREAM^.BEDELEVATION;
   FOR I := (VEC_POINTS-1) DOWNTO 1 DO
   BEGIN
      UPDIST_PRED[I] := UPDIST_PRED[I+1] + DELX_PRED[I];
      GRAD_PRED[I] := (EL_PRED[I]-EL_PRED[I+1])/DELX_PRED[I];
   END;
   {WRITELN('END MAKE_VECTORS');}
END; (*MAKE_VECTORS*)
{**************************************************************************************}
{*********************************************************************************}
PROCEDURE CHECKPROFILE;
VAR I, II, IL, IH, IPOINT , NSMOOTH, IICALL, NCALL : LONGINT;
    MUSTPROCESS : BOOLEAN;
    SQUARE_FROUDE,LOCALFRICTION,FNP,HPR,H_PREDICTED,FN,FR2 : DOUBLE;
BEGIN
{This examines predicted values of channel depth to see if negative or large flow depths
 are predicted and if so, calls SMOOTH_PROFILE or REGRESS_PROFILE}
NSMOOTH := SMOOTH_LENGTH;
MUSTPROCESS := FALSE;
BADPOINT[1] := FALSE;
FOR I := 1 TO VEC_POINTS DO EL_PROVISIONAL[I] := EL_PRED[I];
REPEAT
BEGIN
NUMBER_FAULTS:=0;
H_DEPTH[VEC_POINTS] := DOWNSTREAM_DEPTH;
IF (NOT USE_PROVISIONAL) THEN FOR I := 1 TO VEC_POINTS DO EL_PROVISIONAL[I] := EL_PRED[I];
FOR I:= 1 TO VEC_POINTS-1 DO
BEGIN
    if ((h_depth[vec_points+1-i]<=0.0) OR (h_depth[vec_points+1-i]>(3.0*DOWNSTREAM_DEPTH)))then
	begin
	        ncall :=1;
           	IF (SMOOTH_METHOD = 1) THEN
		  	     SMOOTH_PROFILE(VEC_POINTS+1-I,NSMOOTH,ncall)
		    ELSE
			     REGRESS_PROFILE(VEC_POINTS+1-I,NSMOOTH,ncall);
		   NUMBER_OK := 0;
		   NUMBER_FAULTS := NUMBER_FAULTS+1;
		   IICALL:= vec_points+1-i;
		   BREAK;
	end;
    SQUARE_FROUDE := SQR(QPERW)/(GRAVITY*EXP(3.0*LN(H_DEPTH[VEC_POINTS+1-I])));
	IF SQUARE_FROUDE >0.8 THEN
	BEGIN
	        ncall :=2;
	      	IF (SMOOTH_METHOD = 1) THEN
		  	     SMOOTH_PROFILE(VEC_POINTS+1-I,NSMOOTH,ncall)
		    ELSE
			     REGRESS_PROFILE(VEC_POINTS+1-I,NSMOOTH,ncall);
		   NUMBER_OK := 0;
		   NUMBER_FAULTS := NUMBER_FAULTS+1;
		   IICALL := vec_points+1-i;
		   BREAK;
    END;
	LCFRICTION[I] := EXP(-0.3333*LN(H_DEPTH[VEC_POINTS+1-I]/ROUGHNESS))/(SQR(MANNING_STRICKLER_COEFFICIENT));
	FNP := (EL_PROVISIONAL[VEC_POINTS-I] - EL_PROVISIONAL[VEC_POINTS+1-I]-LCFRICTION[I]*SQUARE_FROUDE*DELX_PRED[VEC_POINTS-I])/(1.0-SQUARE_FROUDE);
	HPR :=H_DEPTH[VEC_POINTS+1-I]-FNP;
	if hpr<= 0.0 then
        BEGIN
		    ncall :=3;
            IF (SMOOTH_METHOD = 1) THEN
		  	     SMOOTH_PROFILE(VEC_POINTS+1-I,NSMOOTH,ncall)
		    ELSE
			     REGRESS_PROFILE(VEC_POINTS+1-I,NSMOOTH,ncall);
		    NUMBER_OK := 0;
		     NUMBER_FAULTS := NUMBER_FAULTS+1;
			 IICALL:= vec_points+1-i;
		     BREAK;
         END;
    FR2 := SQR(QPERW)/(GRAVITY*EXP(3.0*LN(HPR)));
	FN := (EL_PROVISIONAL[VEC_POINTS-I]-EL_PROVISIONAL[VEC_POINTS+1-I] - LCFRICTION[I]* FR2* DELX_PRED[VEC_POINTS-I])/(1.0-FR2);
	H_DEPTH[VEC_POINTS-I] := H_DEPTH[VEC_POINTS+1-I] -0.5 *(FNP+FN);
	IF ((H_DEPTH[VEC_POINTS-I]<0.0) OR (H_DEPTH[VEC_POINTS-I]>(3.0*DOWNSTREAM_DEPTH)) OR (FR2>0.8)) THEN
	BEGIN
	     ncall :=4;
	     IF (SMOOTH_METHOD = 1) THEN
		  	     SMOOTH_PROFILE(VEC_POINTS+1-I,NSMOOTH,ncall)
		    ELSE
			     REGRESS_PROFILE(VEC_POINTS+1-I,NSMOOTH,ncall);
		   NUMBER_OK := 0;
		   NUMBER_FAULTS := NUMBER_FAULTS+1;
		   IICALL:= vec_points+1-i;
		   BREAK;
	END;
    END;
END;
NSMOOTH := NSMOOTH+1;
UNTIL ((NSMOOTH>(SMOOTH_LENGTH+MAX_FAULT_ITERATIONS)) OR (NUMBER_FAULTS =0));
IF (NSMOOTH > SMOOTH_LENGTH+MAX_FAULT_ITERATIONS) THEN
BEGIN
   FAULT_STATE := TRUE;
   II := IICALL;
   WRITELN('FAULT STATE, I= ',(II):5,' NSMOOTH=',(NSMOOTH-1):8,' H_DEPTH =',H_DEPTH[II]:15,' EL =',EL_PROVISIONAL[II]:15,
          ' EL ORIG=',EL_PRED[II]:15,' SQFR=',SQUARE_FROUDE:15,' FNP=',FNP:15,' HPR=',HPR:15,' FR2=',FR2:15,' FN=',FN:15);
   WRITELN(LIST,'FAULT STATE, I= ',(II):5,' NSMOOTH=',(NSMOOTH-1):8,' H_DEPTH =',H_DEPTH[II]:15,' EL =',EL_PROVISIONAL[II]:15,
          ' EL ORIG=',EL_PRED[I]:15,' SQFR=',SQUARE_FROUDE:15,' FNP=',FNP:15,' HPR=',HPR:15,' FR2=',FR2:15,' FN=',FN:15);
END
ELSE
  FOR I:= 1 TO VEC_POINTS-1 DO  EL_PRED[I] := EL_PROVISIONAL[I];



END; (*CHECKPROFILE*)
{**************************************************************************************}
{*********************************************************************************}
PROCEDURE PARKER_BACKWATER_TRANSPORT;
VAR I ,ii, NSMOOTH: LONGINT;
    SQUARE_FROUDE, FNP,  HPR, FR2, FN, DIMENSIONLESS_SHEAR, DIMENSIONLESS_TRANSPORT_FLUX, QTBACK,
	QTIT, QTFRNT, QTDIF,ELDIF, LOCALGRADIENT, LOCALDEPTH,DELX_USE: DOUBLE;
BEGIN
{This routes bed-sediment downstream by first constructing the backwater water elevation profile,
 calculating flow depth and velocity, and determining elevation changes by calculating sediment flux
 divergence  to deposit or erode the bed. This uses the Parker backwater profile
 sediment aggradation and degradation routine}
IF NORMAL_APPROXIMATION THEN
BEGIN
    FAULT_STATE :=FALSE;
    FOR I :=1 TO VEC_POINTS-1 DO
	BEGIN
	   LOCALGRADIENT := (EL_PRED[I]-EL_PRED[I+1])/DELX_PRED[I];
           IF LOCALGRADIENT > 0.0 THEN
	      H_DEPTH[I] := EXP(0.3*LN((EXP(0.3333*LN(ROUGHNESS))*SQR(QPERW))/(SQR(MANNING_STRICKLER_COEFFICIENT)*GRAVITY*LOCALGRADIENT)))
           ELSE
              H_DEPTH[I] := 1.5 * DOWNSTREAM_DEPTH;
	END;
	H_DEPTH[VEC_POINTS] := DOWNSTREAM_DEPTH;
END
ELSE
BEGIN
   IF (SMOOTH_METHOD <>0) THEN CHECKPROFILE;
   NSMOOTH := 3;
END;
IF NOT FAULT_STATE THEN
BEGIN

FOR I :=1 TO VEC_POINTS DO
BEGIN
   H_ELEV[I] := EL_PRED[I] + H_DEPTH[I];
   U_PRED[I] := QPERW/H_DEPTH[I];
END;
{WRITELN('H[1]=',H_DEPTH[1]:12,'  H[END]=',H_DEPTH[VEC_POINTS]:12,
  ' U[1]=',U_PRED[1]:12,' U[END]=',U_PRED[VEC_POINTS]:12);}
{STEP 2: FIND SHIELDS STRESS AND LOAD}
FOR I :=1 TO VEC_POINTS DO
BEGIN
IF NORMAL_APPROXIMATION THEN
    LOCALDEPTH :=H_DEPTH[I]
ELSE
BEGIN
   if h_depth[i] <=0.0 then
       LOCALDEPTH := DOWNSTREAM_DEPTH
   ELSE
      LOCALDEPTH :=H_DEPTH[I]; {writeln('oops, negative depth');}
END;
LCFRICTION[I] := (1.0/SQR(MANNING_STRICKLER_COEFFICIENT))/EXP(0.3333*LN(LOCALDEPTH/ROUGHNESS));
DIMENSIONLESS_SHEAR := LCFRICTION[I] *SQR(QPERW/LOCALDEPTH)/(SUBMERGED_SPECIFIC_GRAVITY* GRAVITY * SEDIMENT_SIZE);
{WRITELN('T=',DIMENSIONLESS_SHEAR:12);}
IF (DIMENSIONLESS_SHEAR*SKIN_FRACTION)>TRANSPORT_CRITICAL_DIM_SHEAR THEN
   DIMENSIONLESS_TRANSPORT_FLUX := TRANSPORT_FACTOR*EXP(SEDIMENT_TRANSPORT_EXPONENT*LN(SKIN_FRACTION*DIMENSIONLESS_SHEAR-TRANSPORT_CRITICAL_DIM_SHEAR))
ELSE
   DIMENSIONLESS_TRANSPORT_FLUX := 0.0;
FLUX_PRED[I] := SQRT(SUBMERGED_SPECIFIC_GRAVITY * GRAVITY*SEDIMENT_SIZE)* SEDIMENT_SIZE * DIMENSIONLESS_TRANSPORT_FLUX;
END;
{WRITELN('CF[1]=',LCFRICTION[1]:12,' CF[END]=',LCFRICTION[VEC_POINTS]:12,
' QS[1]=',FLUX_PRED[1]:12,' QS[END]=',FLUX_PRED[VEC_POINTS]:12);}
{STEP 3: EVOLVE THE BED PROFILE}
FOR I:=1 TO VEC_POINTS-1 DO
BEGIN
IF I=1 THEN
BEGIN
   QTBACK := UPSTREAM_SEDIMENT_FLUX;
   DELX_USE:= DELX_PRED[I];
END ELSE
BEGIN
   QTBACK := FLUX_PRED[I-1];
   DELX_USE:=0.5*(DELX_PRED[I-1]+DELX_PRED[I]);
END;
QTIT := FLUX_PRED[I];
QTFRNT :=FLUX_PRED[I+1];
QTDIF := UPWIND * (QTBACK - QTIT) + (1.0 - UPWIND) *(QTIT - QTFRNT);
ELDIF_PRED[I] :=  QTDIF * ITERATION_TIME/((1.0-SEDIMENT_POROSITY)*DELX_USE);
{WRITELN(' CHECK RATE =',QTDIF*ITERATION_TIME/DELX_USE:15);}
EL_PRED[I] := EL_PRED[I] + ELDIF_PRED[I];
  NETBEDCHANGE := NETBEDCHANGE + ELDIF_PRED[I];
      ABSBEDCHANGE := ABSBEDCHANGE + ABS(ELDIF_PRED[I]);
      NETBEDNUMBER := NETBEDNUMBER+1.0;
END;
QTDIF := FLUX_PRED[VEC_POINTS - 1] -FLUX_PRED[VEC_POINTS];
IF NUMBER_FAULTS>MAXIMUM_FAULTS THEN
BEGIN
 FAULT_STATE := TRUE;
 WRITELN(LIST,'EXITING PROGRAM BECAUSE OF LARGE SUPERCRITICAL FLOW REACH');
 WRITELN('EXITING PROGRAM BECAUSE OF LARGE SUPERCRITICAL FLOW REACH');
END;

IF NUMBER_OK >= RESET_NUMBER THEN NUMBER_FAULTS :=0;
NUMBER_OK := NUMBER_OK+1;
END;
{COMMENT OUT NEXT STATMENT BECAUSE ENDPOINT IS FIXED}
{EL_PRED[VEC_POINTS] := EL[VEC_POINTS] + TRANSPORT_INCREMENT*QTDIF * FLOW_FRACTION /((1.0-SEDIMENT_POROSITY)*DELX_PRED[VEC_POINTS]);}
END;
{*********************************************************************************}
{*********************************************************************************}

{*********************************************************************************}
{*********************************************************************************}
PROCEDURE PUT_VECTORS;
VAR I: LONGINT;
GRAD_CARRY : DOUBLE;
BEGIN
{This determines vector values of depth, velocity, bed elevation, and sediment
 flux based upon changes predicted by PARKER_BACKWATER_TRANSPORT}
{WRITELN('START PUT_VECTORS');}
    STREAM := BASE;
	IF (VEC_POINTS<>RIVERSIZE) THEN
	    WRITE('OOPS, VECPOINTS =',VEC_POINTS:7, '  RIVERSIZE =',RIVERSIZE:7);
    FOR I := 1 TO VEC_POINTS-1 DO
	BEGIN
	   STREAM^.LOCALDEPTH := H_DEPTH[I];
	   STREAM^.LOCALVELOCITY := U_PRED[I];
	   STREAM^.BEDELEVATION := EL_PRED[I];
	   STREAM^.SEDIMENTFLUX := FLUX_PRED[I];
	   SUM_FLUX := SUM_FLUX + FLUX_PRED[I];
	   NUMBER_SUM_FLUX := NUMBER_SUM_FLUX + 1.0;
	   STREAM^.LOCALCFRICTION := LCFRICTION[I];
           STREAM^.LOCALGRADIENT := (EL_PRED[I]-EL_PRED[I+1])/DELX_PRED[I];
		   STREAM^.BEDCHANGE := ELDIF_PRED[I];
           GRAD_CARRY := STREAM^.LOCALGRADIENT;
	   STREAM^.UPSTREAM_DISTANCE := UPDIST_PRED[I];
	   STREAM := STREAM^.NEXT;
	END;
	   STREAM^.LOCALDEPTH := H_DEPTH[VEC_POINTS];
	   STREAM^.LOCALVELOCITY := U_PRED[VEC_POINTS];
	   STREAM^.BEDELEVATION := EL_PRED[VEC_POINTS];
	   STREAM^.SEDIMENTFLUX := FLUX_PRED[VEC_POINTS];
	   SUM_FLUX := SUM_FLUX + FLUX_PRED[I];
	   NUMBER_SUM_FLUX := NUMBER_SUM_FLUX + 1.0;
           STREAM^.LOCALGRADIENT := GRAD_CARRY;
		   STREAM^.BEDCHANGE := ELDIF_PRED[VEC_POINTS];
	   STREAM^.LOCALCFRICTION :=LCFRICTION[VEC_POINTS];
	   STREAM^.UPSTREAM_DISTANCE := UPDIST_PRED[VEC_POINTS];
END; (*PUT_VECTORS*)
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE EVOLVE_STREAM_PROFILE;
VAR I : LONGINT;
AVG_SEDFLUX, SUM_SEDFLUX : DOUBLE;

BEGIN
{This procedure evolves the channel bed profile and sediment flux by calling
 PARKER_BACKWATER_TRANSPORT several iterations per master channel meander iterations.
 The number of subiterations is given by the parameter TRANSPORT_ITERATIONS}
{WRITELN('START EVOLVE_STREAM_PROFILE');}
   {TRANSPORT_ITERATIONS := ROUND(INTERVAL/TRANSPORT_INCREMENT);
   IF TRANSPORT_ITERATIONS<1 THEN TRANSPORT_ITERATIONS :=1;}
 {WRITELN('TRANSPORT ITERATIONS = ',TRANSPORT_ITERATIONS:10,'INTERVAL =',INTERVAL:10);}
   MAKE_VECTORS;

   FOR I :=1 TO TRANSPORT_ITERATIONS DO
   BEGIN
      PARKER_BACKWATER_TRANSPORT;
	  IF FAULT_STATE THEN EXIT;
   END;
   PUT_VECTORS;
   AVG_SEDFLUX := 0.0;
   FOR I := 1 TO VEC_POINTS DO
       AVG_SEDFLUX := AVG_SEDFLUX + FLUX_PRED[I];
   AVG_SEDFLUX := AVG_SEDFLUX/VEC_POINTS;
   SEDLOADS_INDEX := SEDLOADS_INDEX +1;
   IF SEDLOADS_INDEX > 50 THEN SEDLOADS_INDEX := 1;
   PAST_SEDLOADS[SEDLOADS_INDEX] := AVG_SEDFLUX;
   IF NUMBER_SEDLOADS <50 THEN NUMBER_SEDLOADS := NUMBER_SEDLOADS +1;
   FOR I := 1 TO NUMBER_SEDLOADS DO
       SUM_SEDFLUX := SUM_SEDFLUX + PAST_SEDLOADS[I];
   NOMINAL_SEDFLUX := SUM_SEDFLUX /NUMBER_SEDLOADS;
   {WRITELN(' NOMINAL SEDFLUX = ',NOMINAL_SEDFLUX:15);}


   {IF (ITERNUM MOD PLOTINT = 0) AND (ITERNUM > MINITERATION) THEN
   BEGIN
      WRITE_DEBUG;
	  WRITE_VECTOR;
   END;}
{WRITELN('END NEW EVOLVE STREAM PROFILE');}
END; (*EVOLVE_STREAM_PROFILE*)

{*********************************************************************************}
{*********************************************************************************}

PROCEDURE SET_INITIAL_ELEVATIONS;
VAR LASTELEV : REAL;
I : LONGINT;
BEGIN
{This creates the initial stream profile for simulations using sediment routing
 using the INITIAL_GRADIENT parameter and sets initial values of velocity, flow depth,
 channel width, and bed friction using initial values that are read in}
WRITELN('INITIAL ELEVS');
FOR I := 1 TO 50 DO
    PAST_SEDLOADS[I] := 0.0;
SEDLOADS_INDEX := 0;
NUMBER_SEDLOADS := 0;
NOMINAL_SEDFLUX := 0.0;
STREAM := BASE;
STREAM^.BEDELEVATION := 0.0;
X1 := STREAM^.XVALUE;
X2 := STREAM^.NEXT^.XVALUE;
Y1 := STREAM^.YVALUE;
Y2 := STREAM^.NEXT^.YVALUE;
FINDDISTANCE(UPLENGTH,X1,Y1,X2,Y2);
{IF DO_INITIAL_GRADIENT THEN}
      STREAM^.NEXT^.BEDELEVATION := STREAM^.BEDELEVATION-UPLENGTH*INITIAL_GRADIENT;
	  {ELSE
	  STREAM^.NEXT^.BEDELEVATION := STREAM^.BEDELEVATION-UPLENGTH*EQUILIBRIUM_GRADIENT;}


WHILE (STREAM^.NEXT^.NEXT <> NIL) DO
BEGIN
  STREAM := STREAM^.NEXT;
  X1 := X2;
  Y1 := Y2;
  X2 := STREAM^.NEXT^.XVALUE;
  Y2 := STREAM^.NEXT^.YVALUE;
  FINDDISTANCE(UPLENGTH,X1,Y1,X2,Y2);
  {IF DO_INITIAL_GRADIENT THEN}
      STREAM^.NEXT^.BEDELEVATION := STREAM^.BEDELEVATION-UPLENGTH*INITIAL_GRADIENT;
	  {ELSE
	  STREAM^.NEXT^.BEDELEVATION := STREAM^.BEDELEVATION-UPLENGTH*EQUILIBRIUM_GRADIENT;}

  LASTELEV := STREAM^.NEXT^.BEDELEVATION;
END;
STREAM := BASE;
STREAM^.BEDELEVATION := STREAM^.BEDELEVATION-LASTELEV;
WRITELN(LIST,'UPSTREAM BED ELEVATION =',STREAM^.BEDELEVATION:15);
STREAM^.LOCALDEPTH:= DOWNSTREAM_DEPTH;
STREAM^.LOCALVELOCITY := QPERW/DOWNSTREAM_DEPTH;
STREAM^.LOCALCFRICTION := CFRICTION;
STREAM^.LOCALWIDTH := CHANNEL_WIDTH;
WHILE (STREAM^.NEXT <> NIL) DO
BEGIN
  STREAM := STREAM^.NEXT;
  STREAM^.LOCALDEPTH:= DOWNSTREAM_DEPTH;
  STREAM^.LOCALVELOCITY := QPERW/DOWNSTREAM_DEPTH;
  STREAM^.LOCALCFRICTION := CFRICTION;
  STREAM^.LOCALWIDTH := CHANNEL_WIDTH;
  STREAM^.BEDELEVATION := STREAM^.BEDELEVATION-LASTELEV;
END;
WRITELN(LIST,' DOWNSTREAM BED ELEVATION =',STREAM^.BEDELEVATION:15);
END; (*SET_INTITIAL_ELEVATIONS*)
{********************************************************************}
{*********************************************************************************}



PROCEDURE WRITEBANKDATA;
VAR I,J : LONGINT;
BEGIN
{This writes an output file of LONGINT values depending upon the whether the floodplain cell has normal or plug erodibility}
	FOR I := 1 TO XPOINTS DO
	BEGIN
	  FOR J := 1 TO YPOINTS DO
	    IF ((BANKERODE[I,J] = ERODEVALLEY) OR (BANKERODE[I,J] = -ERODEVALLEY)) AND (ERODEVALLEY<>ERODENOMINAL) THEN
	        WRITE(OUTBANKTYPE,'1')
	    ELSE IF ((BANKERODE[I,J] = ERODEPLUG) OR (BANKERODE[I,J] = -ERODEPLUG)) AND (ERODEPLUG<>ERODENOMINAL) THEN
	        WRITE(OUTBANKTYPE,'2')
	    ELSE
	        WRITE(OUTBANKTYPE,'0');
	    WRITELN(OUTBANKTYPE);
	END;
END;
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE FINDAGE;
VAR I,J,IBEGIN,JBEGIN,IFINISH,JFINISH : LONGINT;
    DX,DY,IX,IY,IDX,IDY,INCR,XX,YY,IPX,IPY : LONGINT;
    PLOT : BOOLEAN;
    meanage : real;
BEGIN
{This increments the age of each floodplain cell from the start of the simulation
 or the last time the channel has migrated through that cell, whichever is smaller.
 It also determines the elevation of the channel bank if the cell has been migrated through}
  FOR I := 1 TO XPOINTS DO
  FOR J := 1 TO YPOINTS DO
  IF AGE[I,J] >= INITIALAGE THEN
    AGE[I,J] := AGE[I,J] + UPDATEINT;
  STREAM := BASE;
  IBEGIN := ROUND(STREAM^.XVALUE/CHANNEL_WIDTH) + ISTART;
  JBEGIN := ROUND(STREAM^.YVALUE/CHANNEL_WIDTH) + JCENTER;
  STREAM := STREAM^.NEXT;
  WHILE (STREAM^.NEXT <> NIL) DO
  BEGIN
      IFINISH := ROUND(STREAM^.XVALUE/CHANNEL_WIDTH) + ISTART;
      JFINISH := ROUND(STREAM^.YVALUE/CHANNEL_WIDTH) + JCENTER;
      IF (IBEGIN >= 1) AND (IBEGIN <= XPOINTS) AND
         (IFINISH >= 1) AND (IFINISH <= XPOINTS) AND
         (JBEGIN >= 1) AND (JBEGIN <= YPOINTS) AND
         (JFINISH >= 1) AND (JFINISH <= YPOINTS) THEN
      BEGIN
    DX := IFINISH - IBEGIN;
    DY := JFINISH - JBEGIN;
    IF DX >= 0 THEN IDX := 1 ELSE IDX := -1;
    IF DY >= 0 THEN IDY := 1 ELSE IDY := -1;
          IX := ABS(DX);
    IY := ABS(DY);
          IF IX > IY THEN INCR := IX ELSE INCR := IY;
    IPX := IBEGIN;
    IPY := JBEGIN;
    XX := 0;
    YY := 0;
    AGE[IPX,IPY] := 0;
    ELEVATION[IPX,IPY]:=CHANNELELEV;
    IF USEBANKELEV THEN
                BANK[IPX,IPY] :=
                   ROUND(STREAM^.BANKDEPTH * MEANDEPTH)
    ELSE
    		BANK[IPX,IPY] := 0;
    FOR I := 0 TO INCR DO
    BEGIN
        XX := XX + IX;
        YY := YY + IY;
        PLOT := FALSE;
        IF XX > INCR THEN
        BEGIN
      PLOT := TRUE;
      XX := XX - INCR;
      IPX := IPX + IDX;
        END;
        IF YY > INCR THEN
        BEGIN
      PLOT := TRUE;
      YY := YY - INCR;
      IPY := IPY + IDY;
        END;
        IF PLOT THEN
                    BEGIN
      AGE[IPX,IPY] := 0;
      ELEVATION[IPX,IPY] := CHANNELELEV;
      IF USEBANKELEV THEN
                        BANK[IPX,IPY] :=
                           ROUND(STREAM^.BANKDEPTH * MEANDEPTH)
      ELSE
      			BANK[IPX,IPY] := 0;
        END;
    END;
      END;
      STREAM := STREAM^.NEXT;
      IBEGIN := IFINISH;
      JBEGIN := JFINISH;
  END;
  meanage := 0.0;
  FOR I := 1 TO XPOINTS DO
  FOR J := 1 TO YPOINTS DO
  BEGIN
	IF (AGE[I,J]>0) AND (BANKERODE[I,J]<0.0) THEN BANKERODE[I,J] := ERODENOMINAL;
   IF (AGE[I,J]=0) AND (BANKERODE[I,J]>0.0) THEN BANKERODE[I,J] := -BANKERODE[I,J];
   meanage := meanage + age[i,j];
  END;
  meanage :=meanage/(xpoints*ypoints);
  writeln('mean age=',meanage:10);
END;
{*********************************************************************************}
{*********************************************************************************}

PROCEDURE FINDCHANNELDISTANCE(I,J : LONGINT; VAR DISTANCE : REAL);
VAR IPOINT, JPOINT : LONGINT;
VAR XSTREAM, YSTREAM  : REAL;
VAR K : LONGINT;
BEGIN
{Given a cell I,J location this calculates the real distance from the channel to the cell}
  DISTANCE := MAXWIDTHCHECK;
  IPOINT := (I-ISTART);
  JPOINT := (J-JCENTER);
    IF KVAL[I,J] > 0 THEN
    BEGIN
  K := 1;
  STREAM := BASE;
  WHILE (STREAM^.NEXT <> NIL) DO
  BEGIN
     IF (K = KVAL[I,J])  THEN
     BEGIN
    XSTREAM := STREAM^.XVALUE/CHANNEL_WIDTH;
    YSTREAM := STREAM^.YVALUE/CHANNEL_WIDTH;
    DISTANCE := SQRT(SQR(IPOINT-XSTREAM)
                    +SQR(JPOINT-YSTREAM));
     END;
     STREAM := STREAM^.NEXT;
     K := K+1;
  END;
  IF DISTANCE > MAXWIDTHCHECK THEN
           WRITELN('IP=',IPOINT,' JP=',JPOINT,
    ' D=',DISTANCE:10);

    END;
END;
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE DEPOSIT;
VAR I,J,ICAT,IMAX,ENOW : LONGINT;
     AVGVAL,NSUM,EXPARG,DISTANCE, DEPOSIT_RATE : REAL;

BEGIN
{This deposits sediment on the floodplain as a function of the distance from each
 cell to the nearest channel and resets the elevation to the channel bed elevation
 when the channel migrates through the cell}
        IMAX := TRUNC(MAXWIDTHCHECK)+1;
  AVGVAL := 0.0;
  NSUM := 0.0;
  FOR I := 1 TO XPOINTS DO
  FOR J := 1 TO YPOINTS DO
  BEGIN
   IF (ELEVATION[I,J]< INITIAL_ELEVATION) THEN
   BEGIN
	IF (AGE[I,J] > 0) THEN
	BEGIN
		ELEVATION[I,J] := ELEVATION[I,J] - ROUND(AGGRADERATE*UPDATEINT);
		DISTANCE := IDIST[I,J];
		IF DISTANCE>IMAX THEN DISTANCE := IMAX;
		ELEVATION[I,J] := ELEVATION[I,J]+BANK[I,J];
		IF ABS(BANK[I,J])>INITIAL_ELEVATION THEN WRITELN(I,' ',J,' BANK=',BANK[I,J]);
		BANK[I,J] := 0;
		ENOW := ELEVATION[I,J];
		IF ENOW>INITIAL_ELEVATION THEN WRITELN(I,' ',J,' ENOW=',ENOW,' BANK=',BANK[I,J]);
		ICAT := ROUND(DISTANCE) + 1;
	IF USELINEARRATE THEN
		DEPOSIT_RATE := ELEVRATE * (ELEVMAX - ENOW)
			* (FLOODRATE + BANKRATE * DISTWEIGHTS[ICAT])
	ELSE BEGIN
		IF ENOW>CHANNELELEV THEN
		BEGIN
			EXPARG := SQR(ENOW-CHANNELELEV) * DEPOSITDECAY;
                        IF EXPARG<0.0 THEN EXPARG :=0.0;
			AVGVAL := AVGVAL + EXPARG;
			NSUM := NSUM+1.0;
			IF EXPARG<15.0 THEN
				DEPOSIT_RATE := ELEVRATE * EXP(-EXPARG)
				* (FLOODRATE + BANKRATE * DISTWEIGHTS[ICAT])
			ELSE
				DEPOSIT_RATE := 0.0;
		END ELSE
     			DEPOSIT_RATE := ELEVRATE * (FLOODRATE + BANKRATE * DISTWEIGHTS[ICAT]);
	END;
		IF DEPOSIT_RATE > 0.0 THEN
			ELEVATION[I,J] :=  ELEVATION[I,J] +
			ROUND(DEPOSIT_RATE * UPDATEINT);
		{IF (ITERNUM = ITERATION) THEN
		BEGIN
		   WRITELN('I=',I,' J=',J,' E=',ENOW:10,' ENEW=',ELEVATION[I,J]:10,' DEP=',DEPOSIT_RATE:10);
		   WRITELN('EXPARG=',EXPARG:10,' ICAT=',ICAT,' DWGT=',DISTWEIGHTS[ICAT]:10);
		END;
		IF (ELEVATION[I,J]<(-50000)) OR ((DEPOSIT_RATE * UPDATEINT)>20000)THEN
		BEGIN
		  WRITELN(I,' ',J,' E=',ELEVATION[I,J],' DR=',DEPOSIT_RATE * UPDATEINT:10,' B=',BANK[I,J]);
		  WRITELN('ENOW=',ENOW,' DIST=',DISTWEIGHTS[ICAT]:10,' ICAT=',ICAT);
		END;}
      END; {AGE>}

       {IF AGE[I,J] = 0 THEN
		ELEVATION[I,J] := CHANNELELEV;}
	END; {INITIAL_ELEVATION}
   END; {FOR I J}
END;
{*********************************************************************************}
{*********************************************************************************}

PROCEDURE FINDRANGE;
VAR K,I,J,ISTREAM,JSTREAM,II,JJ, IIDIST : LONGINT;

BEGIN
{This locates the channel node closest to each floodplain cell and the LONGINT
 distance from the channel to the cell}
{  WRITELN('FINDRANGE');  }
  FOR I := 1 TO XPOINTS DO
  FOR J := 1 TO YPOINTS DO
  BEGIN
     IDIST[I,J] := 20000;
     KVAL[I,J] := 0;
  END;
  STREAM := BASE;
  K := 1;
  WHILE (STREAM^.NEXT <> NIL) DO
  BEGIN
    ISTREAM := ROUND(STREAM^.XVALUE/CHANNEL_WIDTH)+ ISTART;
    JSTREAM := ROUND(STREAM^.YVALUE/CHANNEL_WIDTH)+ JCENTER;
    FOR II := -IRANGE TO IRANGE DO
    FOR JJ := -IRANGE TO IRANGE DO
    BEGIN
        I := ISTREAM + II;
        J := JSTREAM + JJ;
        IF (I>0) THEN
        IF (I<=XPOINTS) THEN
        IF (J>0) THEN
        IF (J<=YPOINTS) THEN
        BEGIN
      IIDIST := IRDIST[II+NRANGE,JJ+NRANGE];
      IF IIDIST < IDIST[I,J] THEN
      BEGIN
         IDIST[I,J] := IIDIST;
         KVAL[I,J] := K;
      END;
       END;
    END;
    STREAM := STREAM^.NEXT;
    K := K+1;
  END;
END;
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE SETUPFLOODPLAIN;
VAR I,J,II,JJ : LONGINT;
   VAR TEMP : REAL;
		 JLOWER,JUPPER : LONGINT;
BEGIN
{This initializes the variables defining the state of the floodplain, including
 the floodplain age, erodibility of the floodplain deposits}
  FOR I := 1 TO XPOINTS DO
  FOR J := 1 TO YPOINTS DO
  BEGIN
       AGE[I,J] := INITIALAGE;
             ELEVATION[I,J] := INITIAL_ELEVATION;
       BANK[I,J] := 0;
  END;
  IF VALLEY_WIDTH >0.0 THEN
  BEGIN
  	JUPPER := ROUND(VALLEY_WIDTH/CHANNEL_WIDTH) + JCENTER;
  	JLOWER := ROUND(-VALLEY_WIDTH/CHANNEL_WIDTH) + JCENTER;
   FOR I := 1 TO XPOINTS DO
   BEGIN
     FOR J := 1 TO JLOWER DO
     BEGIN
		BANKERODE[I,J] := ERODEVALLEY;
     END;
     FOR J := (JLOWER+1) TO JUPPER DO
     BEGIN
		BANKERODE[I,J] := ERODENOMINAL;
     END;
     FOR J := (JUPPER+1) TO YPOINTS DO
     BEGIN
		BANKERODE[I,J] := ERODEVALLEY;
     END;
	END;
  END ELSE IF READERODIBILITY < 1 THEN
  BEGIN
   FOR I := 1 TO XPOINTS DO
   FOR J := 1 TO YPOINTS DO
		BANKERODE[I,J] := ERODEVALLEY;
  END ELSE
  BEGIN
   FOR I := 1 TO XPOINTS DO
   FOR J := 1 TO YPOINTS DO
		BANKERODE[I,J] := ERODENOMINAL;
  END;
  IF IREADDEPOSIT > 0 THEN
BEGIN
	  ASSIGN(AGEDAT,'alev.dat');
     RESET(AGEDAT);
     WHILE NOT EOF(AGEDAT) DO
                BEGIN
      READ(AGEDAT,I,J);
      READLN(AGEDAT,AGE[I,J],ELEVATION[I,J],IDIST[I,J]);
                END;
     CLOSE(AGEDAT);
        END;
  FOR I := 1 TO DISTCATAGORIES DO
    DISTWEIGHTS[I] := EXP(-(I-1) * BANKDISTFACTOR);
  IRANGE := IDISTCATAGORIES DIV 2 - 1;
  NRANGE := IRANGE;
  FOR I := 0 TO (2*IRANGE) DO
  FOR J := 0 TO (2*IRANGE) DO
  BEGIN
      II := I - IRANGE;
      JJ := J - IRANGE;
      IRDIST[I,J] := ROUND(SQRT(II*II+JJ*JJ));
  END;
	IF READERODIBILITY>0 THEN
	BEGIN
		ASSIGN(INERODEFILE,'inebank.dat');
		RESET(INERODEFILE);
		REPEAT
			READLN(INERODEFILE,II,JJ,TEMP);
			BANKERODE[II,JJ]:= TEMP;
		UNTIL EOF(INERODEFILE);
	END;
END;
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE WRITEWIDTH;
VAR ISTART,IEND,IWIDTH,I,J : LONGINT;
    WIDTHSUM,WIDTHNUM : REAL;
BEGIN
{This writes out the vertical width of the floodplain has been reworked by
 channel migration.  Note that this assumes the channel initially
 has a vertical extent of zero (i.e., follows the X axis)}
   WIDTHNUM:=0.0;
   WIDTHSUM:=0.0;
  FOR I := 50 TO XPOINTS DO
  BEGIN
      ISTART := -1;
      IEND := ISTART;
      FOR J := 1 TO YPOINTS DO
      BEGIN
       IF (ISTART<0) AND (AGE[I,J]>0) THEN
           ISTART := J;
       IF AGE[I,J]>0 THEN
           IEND := J;
      END;
      IF ITERNUM<15 THEN
         WRITELN(' I=',I,' ISTART=',ISTART,' IEND=',IEND);
      IWIDTH := IEND - ISTART + 1;
      WIDTHSUM := WIDTHSUM + IWIDTH;
      WIDTHNUM := WIDTHNUM + 1.0;
  END;
  WRITELN(OUTTIME,' ',(WIDTHSUM/WIDTHNUM):13);
END;
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE AGE_ELEVATION;
BEGIN
{This is the master routine for modeling the floodplain, calling
 FINDAGE, FINDRANGE, DEPOSIT, and WRITEWIDTH sequentially}
  FINDAGE;
  FINDRANGE;
  DEPOSIT;
  WRITEWIDTH;
{  WRITELN('DONE');}
END;
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE WRITEAGEELEV;
CONST
	AGECATS = 20;
	ELEVCATS = 20;
	DISTCATS = 31;
	AGEINTERVAL  = 250.0;
	DISTINTERVAL  = 1.0;
	ELEVINTERVAL = 2000.0;
       {ELEVBASE = CHANNELELEV;}
	AGEBASE  = 0;
	DISTBASE = 0;
VAR I, J, IVAL, ELEVBASE : LONGINT;
    THERANGE,MAXVAL,MINVAL : REAL;
	AGESUM : ARRAY[1..AGECATS] OF LONGINT;
	ELEVSUM : ARRAY[0..ELEVCATS] OF LONGINT;
	DISTSUM : ARRAY[1..DISTCATS] OF LONGINT;
	AGEELEV : ARRAY[1..AGECATS] OF REAL;
	NAGEELEV : ARRAY[1..AGECATS] OF LONGINT;
	DISTAGE : ARRAY[1..DISTCATS] OF REAL;
	NDISTAGE : ARRAY[1..DISTCATS] OF LONGINT;
	DISTELEV : ARRAY[1..DISTCATS] OF REAL;
	NDISTELEV : ARRAY[1..DISTCATS] OF LONGINT;
	OUTAGELEV : TEXT;
	TCOUNT, IAGE,IELEV,IIDIST,AGELOC,DISTLOC,ELEVLOC : LONGINT;
	RCOUNT : REAL;
BEGIN
{This routine provides the primary output of the state of the floodplain,
 printing, for each I,J location the location, age, elevation, and
 LONGINT distance to the nearest channel.  It also writes binary image files
 of floodplain age and elevation, as well as several other summary data files}
  ASSIGN(AGEDAT,'aelev.dat');
  REWRITE(AGEDAT);
  ELEVBASE := CHANNELELEV;
  FOR I := 1 TO XPOINTS DO
  FOR J := 1 TO YPOINTS DO
        IF AGE[I,J]>=0 THEN
    WRITELN(AGEDAT,I:3,' ',J:3,' ',AGE[I,J]:5,' ',ELEVATION[I,J]:6
    ,' ',IDIST[I,J]:5);
  CLOSE(AGEDAT);
   ASSIGN(OUTRAW1,'elev.PGM');
   REWRITE(OUTRAW1);
   {TESTING PGM}
   WRITE(OUTRAW1,'P5 ');
   WRITE(OUTRAW1,XPOINTS:10);
   WRITE(OUTRAW1,' ');
   WRITE(OUTRAW1,YPOINTS:10);
   WRITE(OUTRAW1,' 255 ');
   MAXVAL := -1.0E+25;
   MINVAL := -MAXVAL;
   FOR I := 1 TO XPOINTS DO
   FOR J := 1 TO YPOINTS DO
   BEGIN
     IF AGE[I,J]>=0 THEN
     BEGIN
     	IF ELEVATION[I,J]< MINVAL THEN MINVAL := ELEVATION[I,J];
     	IF ELEVATION[I,J]> MAXVAL THEN MAXVAL := ELEVATION[I,J];
     END;
   END;
   THERANGE := MAXVAL-MINVAL;
   writeln(minval:10,' ',maxval:10,' ',therange:10);
   FOR J := 1 TO YPOINTS DO
   FOR I := 1 TO XPOINTS DO
   BEGIN
     IF AGE[I,J]>=0 THEN
         IVAL := TRUNC((ELEVATION[I,J]-MINVAL)*254.0/THERANGE)+1
     ELSE
         IVAL := 0;
     WRITE(OUTRAW1,CHR(IVAL));
   END;
   CLOSE(OUTRAW1);
	ASSIGN(OUTRAW2,'age.PGM');
   REWRITE(OUTRAW2);
   WRITE(OUTRAW2,'P5 ');
   WRITE(OUTRAW2,XPOINTS:10);
   WRITE(OUTRAW2,' ');
   WRITE(OUTRAW2,YPOINTS:10);
   WRITE(OUTRAW2,' 255 ');
   MAXVAL := -1.0E+25;
   MINVAL := -MAXVAL;
   FOR I := 1 TO XPOINTS DO
   FOR J := 1 TO YPOINTS DO
   BEGIN
     IF AGE[I,J]>=0 THEN
     BEGIN
     	IF AGE[I,J]< MINVAL THEN MINVAL := AGE[I,J];
     	IF AGE[I,J]> MAXVAL THEN MAXVAL := AGE[I,J];
     END;
   END;
   THERANGE := MAXVAL-MINVAL;
   FOR J := 1 TO YPOINTS DO
   FOR I := 1 TO XPOINTS DO
   BEGIN
     IF AGE[I,J]>=0 THEN
         IVAL := TRUNC((AGE[I,J]-MINVAL)*254.0/THERANGE)+1
     ELSE
         IVAL := 0;
     WRITE(OUTRAW2,CHR(IVAL));
   END;
   CLOSE(OUTRAW2);
	FOR I := 1 TO AGECATS DO
	BEGIN
		AGESUM[I] := 0;
		AGEELEV[I] := 0.0;
		NAGEELEV[I] := 0;
	END;
	FOR I := 0 TO ELEVCATS DO
	BEGIN
		ELEVSUM[I] := 0;
	END;
	FOR I := 1 TO DISTCATS DO
	BEGIN
		DISTSUM[I] := 0;
		DISTAGE[I] := 0.0;
		NDISTAGE[I] := 0;
		DISTELEV[I] := 0.0;
		NDISTELEV[I] := 0;
	END;
	ASSIGN(OUTAGELEV,'agesum.dat');
	REWRITE(OUTAGELEV);
	TCOUNT := 0;
	RCOUNT := 0.0;
  FOR I := 1 TO XPOINTS DO
  FOR J := 1 TO YPOINTS DO
        IF AGE[I,J]>=0 THEN
	   BEGIN
		IAGE := AGE[I,J];
		IELEV := ELEVATION[I,J];
		IIDIST := IDIST[I,J];
		AGELOC := TRUNC((IAGE-AGEBASE)/AGEINTERVAL)+1;
		IF AGELOC<1 THEN AGELOC := 1;
		IF AGELOC>AGECATS THEN AGELOC := AGECATS;
		ELEVLOC := TRUNC((IELEV-ELEVBASE)/ELEVINTERVAL)+1;
		IF ELEVLOC<0 THEN ELEVLOC := 0;
		IF ELEVLOC>ELEVCATS THEN ELEVLOC := ELEVCATS;
		DISTLOC := TRUNC((IIDIST-DISTBASE)/DISTINTERVAL)+1;
		IF DISTLOC<1 THEN DISTLOC := 1;
		IF DISTLOC>DISTCATS THEN DISTLOC := DISTCATS;
		AGESUM[AGELOC] := AGESUM[AGELOC]+1;
		ELEVSUM[ELEVLOC] := ELEVSUM[ELEVLOC]+1;
		DISTSUM[DISTLOC] := DISTSUM[DISTLOC]+1;
		AGEELEV[AGELOC] := AGEELEV[AGELOC] + IELEV/1000.0;
		NAGEELEV[AGELOC] := NAGEELEV[AGELOC] + 1;
		DISTAGE[DISTLOC] := DISTAGE[DISTLOC]+IAGE;
		NDISTAGE[DISTLOC] := NDISTAGE[DISTLOC]+1;
		DISTELEV[DISTLOC] := DISTELEV[DISTLOC]+IELEV/1000.0;
		NDISTELEV[DISTLOC] := NDISTELEV[DISTLOC]+1;
		TCOUNT := TCOUNT + 1;
		RCOUNT := RCOUNT + 1.0;
	END;
	WRITELN(OUTAGELEV,TCOUNT);
	WRITELN(OUTAGELEV);
	WRITELN(OUTAGELEV,'    0 ',ELEVBASE/1000.0:12:2,' ',ELEVSUM[0]:7);
	FOR I := 1 TO ELEVCATS DO
		WRITELN(OUTAGELEV, I:5,' ',(ELEVBASE+(I-1)*ELEVINTERVAL)/1000.0:12:2,' ',ELEVSUM[I]:7);
	WRITELN(OUTAGELEV);
	FOR I := 1 TO AGECATS DO
	BEGIN
		WRITE(OUTAGELEV,I:5,' ',(AGEBASE+(I-1)*AGEINTERVAL):12:2,' ',AGESUM[I]:7);
		IF NAGEELEV[I]>0 THEN AGEELEV[I] := AGEELEV[I]/NAGEELEV[I];
		WRITELN(OUTAGELEV,' ',AGEELEV[I]:12:3);
	END;
	WRITELN(OUTAGELEV);
	FOR I := 1 TO DISTCATS DO
	BEGIN
		WRITE(OUTAGELEV,I:5,' ',(DISTBASE+(I-1)*DISTINTERVAL):12:2,' ',DISTSUM[I]:7);
		IF NDISTAGE[I]>0 THEN DISTAGE[I] := DISTAGE[I]/NDISTAGE[I];
		IF NDISTELEV[I]>0 THEN DISTELEV[I] := DISTELEV[I]/NDISTELEV[I];
		WRITELN(OUTAGELEV,' ',DISTAGE[I]:12:3,' ',DISTELEV[I]:12:3);
	END;
	CLOSE(OUTAGELEV);
	ASSIGN(OUTERODEFILE,'outebank.dat');
	REWRITE(OUTERODEFILE);
  FOR I := 1 TO XPOINTS DO
  FOR J := 1 TO YPOINTS DO
  BEGIN
  	IF BANKERODE[I,J]<>ERODENOMINAL THEN
		WRITELN(OUTERODEFILE,I,' ',J,' ',BANKERODE[I,J]:12:5);
  END;
  CLOSE(OUTERODEFILE);
END;
{*********************************************************************************}
{*********************************************************************************}

PROCEDURE GETVALUES;

VAR DUMMY : REAL;
    IDUMMY : LONGINT;

BEGIN
{This reads several simulation parameters from the initial.prm file}

    READ(PARAMS,SEED);
    READLN(PARAMS);
    READ(PARAMS, SIMILPAR);
    READLN(PARAMS);
    READ(PARAMS,ERROR);
    READLN(PARAMS);
    READ(PARAMS,LOGSTDEV);
    READLN(PARAMS);
    READ(PARAMS,RES_TYPE);
    READLN(PARAMS);
    IF RES_TYPE = 1 THEN RES_AVERAGE := TRUE ELSE RES_AVERAGE := FALSE;
{    WRITELN('SEED=',SEED:12);
    WRITELN('SIMILPAR=',SIMILPAR:12);
    WRITELN('ERROR=',ERROR:12);
    WRITELN('DISCH ERROR=',LOGSTDEV:12);}
    DUMMY := RAN3(IDUMMY); {RANDOM(IDUMMY);}
END; {GETVALUES}
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE STATVAL(II : LONGINT;  RATE : REAL);

BEGIN
{Adds an observation to the statistical database}
       STATSUM [II] := STATSUM [II] + RATE;
       STATNUM [II] := STATNUM [II] + 1;
END; (*STATVAL*)
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE STATISTICS;


VAR II : LONGINT;
     RW, RATE : REAL;

BEGIN
{Adds observations about migration rate and radius/curvature of planform to database}
    RW := RADIUS / CHANNEL_WIDTH;
    RATE := NEWRATE * NEWSENSE/(INTERVAL * CHANNEL_WIDTH);
    IF (RW < 1.0) THEN
      BEGIN
       II := 1;
       STATVAL(II,RATE)
      END
     ELSE IF (RW < MAX_STAT) THEN
       BEGIN
     II := TRUNC(RW * 2.0);
     STATVAL(II,RATE)
       END (*ELSE IF*)
     ELSE BEGIN
       II := NUM_CASES;
       STATVAL(II,RATE)
     END; (*IF*)
  END; (* STATISTICS *)
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE STATINITIALIZE;

VAR JJ : LONGINT;

BEGIN
{Initializes the statistics variables}
      FOR JJ := 1 TO NUM_CASES DO
        BEGIN
          STATSUM [JJ] := 0.0;
          STATNUM [JJ] := 0;
        END; (*FOR JJ*)
END; (*STATINITIALIZE*)
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE STATPRINT;

VAR II, JJ : LONGINT;
    NTOT : REAL;
    DUMMY : REAL;

BEGIN
{Writes out the statistics about radius/width and migration rate to the list file}
     WRITELN(LIST,'RADIUS/WIDTH,  AVERAGE MIGRATION RATE,',
          ' RELATIVE FREQUENCY');
     WRITELN(LIST);
     NTOT := 0;
     FOR II := 1 TO NUM_CASES DO
         NTOT := NTOT + STATNUM[II];
     IF STATNUM[1] > 0 THEN DUMMY := STATSUM[1]/STATNUM[1]
    ELSE
       DUMMY := 0.0;
     WRITELN(LIST,'R/W= ','    <1.0    ',' RATE= ', DUMMY:10,
          ' F= ',STATNUM[1]/NTOT:10);

     FOR JJ := 2 TO (NUM_CASES-1) DO
      BEGIN
    IF STATNUM[JJ] > 0 THEN DUMMY := STATSUM[JJ]/STATNUM[JJ]
       ELSE DUMMY := 0.0;
    WRITELN(LIST,'R/W= ',(MIN_STAT+(JJ-2)*STAT_INTERVAL):10,' RATE= ',
             DUMMY:10,
               ' F= ',STATNUM[JJ]/NTOT:10);
      END; (*FOR*)
     IF STATNUM[NUM_CASES] > 0 THEN
    DUMMY := STATSUM[NUM_CASES]/STATNUM[NUM_CASES]
    ELSE DUMMY := 0.0;
     WRITELN(LIST,'R/W= ','    >10.0   ',' RATE= ', DUMMY:10,
          ' F= ',STATNUM[NUM_CASES]/NTOT:10);
END;(*STATPRINT*)
{*********************************************************************************}
{*********************************************************************************}



PROCEDURE FINDPERPDIS(VAR PERP: REAL;X2,Y2,XS,YS,XE,YE: REAL);
VAR SLOPE, INTERCEPT : REAL;
BEGIN
{Finds the perpendicular distance between a point and a line defined by two other points}
    IF ABS(XS - XE) < 1.0E-4 THEN
       PERP := ABS(X2-XE)
    ELSE
    BEGIN
       SLOPE := (YE-YS)/(XE-XS);
       INTERCEPT := YE - SLOPE * XE;
       PERP := ABS((X2*SLOPE-Y2+INTERCEPT)/SQRT(SQR(SLOPE)+1));
    END;
{WRITELN(X2:8:3,' ',Y2:8:3,' ',XS:8:3,' ',YS:8:3,' ',XE:8:3,
     ' ',YE:8:3,' ',PERP:10:3);}
END;
{*********************************************************************************}
{*********************************************************************************}

PROCEDURE CURVANGLE(X1,X2,X3,Y1,Y2,Y3 : REAL;VAR ANGLE : REAL);

VAR
   NUMERATOR , DIVISOR  : REAL;
   X21,Y21,X32,Y32 : REAL;

BEGIN
{Given three successive nodes on the centerline, calculates the acute angle
 of intersection of the lines joining the first two points
 and the line joining the second two points}
   X21 := X2 - X1;
   Y21 := Y2 - Y1;
   X32 := X3 - X2;
   Y32 := Y3 - Y2;
   DIVISOR := X21 * X32 + Y21 * Y32;
   NUMERATOR := Y21 * X32 - X21 * Y32;
   IF (DIVISOR = 0.0)
     THEN
    IF (NUMERATOR = 0.0)
      THEN
         ANGLE := 0.0
      ELSE
         IF (NUMERATOR > 0.0)
        THEN
           ANGLE := HALFPI
        ELSE
           ANGLE := -HALFPI
     ELSE
    IF (DIVISOR > 0.0)
       THEN
          ANGLE := ARCTAN(NUMERATOR / DIVISOR)
       ELSE
          ANGLE := -ARCTAN(NUMERATOR / DIVISOR);
   IF (ANGLE >= 0.0)
    THEN
      SENSE := 1
    ELSE
      SENSE := -1;
END; (*CURVANGLEGLE*)
{*********************************************************************************}
{*********************************************************************************}


PROCEDURE INCLUDEANGLE(X1,X2,X3,Y1,Y2,Y3 : REAL;VAR ANGLE : REAL);

VAR
   NUMERATOR , DIVISOR  : REAL;
   X21,Y21,X31,Y31 : REAL;

BEGIN
{Like CURVANGLE, but calculates the interior angle}
   X21 := X2 - X1;
   Y21 := Y2 - Y1;
   X31 := X3 - X1;
   Y31 := Y3 - Y1;
   DIVISOR := X21 * X31 + Y21 * Y31;
   NUMERATOR := Y21 * X31 - X21 * Y31;
{   WRITELN('NUM=', NUMERATOR:12:3,' ','DIV=',DIVISOR:12:3);}
   IF (DIVISOR = 0.0)
   THEN
      IF (NUMERATOR = 0.0)
      THEN
         ANGLE := 0.0
      ELSE
         IF (NUMERATOR > 0.0)
         THEN
           ANGLE := HALFPI
         ELSE
           ANGLE := -HALFPI
   ELSE
      IF (DIVISOR > 0.0)
      THEN
          ANGLE := ARCTAN(NUMERATOR / DIVISOR)
      ELSE
          IF (NUMERATOR>0.0)
          THEN
              ANGLE := PI -ARCTAN(-NUMERATOR / DIVISOR)
          ELSE
              ANGLE := (-PI -ARCTAN(-NUMERATOR / DIVISOR));
END; (*INCLUDEANGLE*)
{*********************************************************************************}
{*********************************************************************************}

PROCEDURE FINDANGLE (VAR ANGLE : REAL; X1, Y1, X2, Y2 : REAL);


BEGIN
{Given two points, calculates the angular orientation of the line extending
 from the first point to the second}
  IF (X1 = X2) OR (Y1 = Y2)
  THEN
    IF X1 = X2
    THEN
      IF Y1 < Y2
      THEN
        ANGLE :=  HALFPI
      ELSE
        ANGLE := -HALFPI
    ELSE
      IF X1 < X2
      THEN
        ANGLE := 0
      ELSE
        ANGLE := PI
  ELSE
    BEGIN
      ANGLE := ARCTAN ((Y2 - Y1) / (X2 - X1));
      IF ((X2 - X1) < 0) AND ((Y2 - Y1) < 0)
      THEN
        ANGLE := ANGLE + PI
      ELSE
        IF ((X2 - X1) < 0)
        THEN
          ANGLE := ANGLE + PI;
    END; (*IF*)
END; (* FINDANGLEGLE *)
{*********************************************************************************}
{*********************************************************************************}


PROCEDURE FINDPLUGS;
VAR I,J,IW,IE,JN,JS,DISTANCE,TCOUNT : LONGINT;
BEGIN
{Determines whether a location on the floodplain is underlain by sediment
 of normal erodibility or is the location of a clay plug deposited in a cutoff loop}
	FINDRANGE;
  FOR I := 1 TO XPOINTS DO
  FOR J := 1 TO YPOINTS DO
  BEGIN
	   IF AGE[I,J] = 0 THEN
		BEGIN
		  DISTANCE := IDIST[I,J];
		  IF DISTANCE>=PLUGDISTANCE THEN BANKERODE[I,J] := ERODEPLUG;
		END;
  END;
  FOR I := 1 TO XPOINTS DO
  FOR J := 1 TO YPOINTS DO
  BEGIN
    TCOUNT := 0;
    IW:=I-1;
    IE:=I+1;
    JN:=J-1;
    JS:=J+1;
    IF (BANKERODE[I,J]=ERODEVALLEY) THEN
    BEGIN
       IF I>1 THEN
         BEGIN
           IF BANKERODE[IW,J]=ERODENOMINAL THEN TCOUNT := TCOUNT+1;
	   IF J>1 THEN
	      IF BANKERODE[IW,JN]=ERODENOMINAL THEN TCOUNT := TCOUNT+1;
	   IF J<YPOINTS THEN
	      IF BANKERODE[IW,JS]=ERODENOMINAL THEN TCOUNT := TCOUNT+1;
	 END;
       IF J>1 THEN
         IF BANKERODE[I,JN]=ERODENOMINAL THEN TCOUNT := TCOUNT+1;
       IF J<YPOINTS THEN
	 IF BANKERODE[I,JS]=ERODENOMINAL THEN TCOUNT := TCOUNT+1;
       IF I< XPOINTS THEN
         BEGIN
           IF BANKERODE[IE,J]=ERODENOMINAL THEN TCOUNT := TCOUNT+1;
	   IF J>1 THEN
	      IF BANKERODE[IE,JN]=ERODENOMINAL THEN TCOUNT := TCOUNT+1;
	   IF J<YPOINTS THEN
	      IF BANKERODE[IE,JS]=ERODENOMINAL THEN TCOUNT := TCOUNT+1;
	 END;
    IF TCOUNT>5 THEN BANKERODE[I,J] := ERODENOMINAL;
    END;
  END;
END;
{*********************************************************************************}
{*********************************************************************************}

PROCEDURE FINDSTREAMLOC(VAR ILOC,JLOC : LONGINT; XLOC,YLOC : REAL);
BEGIN
{Given an X,Y location finds the I,J location of a floodplain cell}
  ILOC := ROUND(XLOC/CHANNEL_WIDTH) + ISTART;
  JLOC := ROUND(YLOC/CHANNEL_WIDTH) + JCENTER;
END;
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE FINDERODIBILITY(XLOC,YLOC : REAL; VAR XCHANGE,YCHANGE : REAL);
VAR ILOC,JLOC : LONGINT;
	 TEMP1,TEMP2 : REAL;
	 DOCHANGE : BOOLEAN;
BEGIN
{Determines the bank erodibility of an X,Y location of a stream node.
  Returns the default erodibility if floodplain modeling is not invoked}
IF USEDEPOSIT THEN
BEGIN
        DOCHANGE := FALSE;
	FINDSTREAMLOC(ILOC,JLOC,XLOC,YLOC);
IF (ILOC>0) AND (ILOC<XPOINTS) AND (JLOC>0) AND (JLOC<YPOINTS) THEN
BEGIN
	ERODIBILITY := ABS(BANKERODE[ILOC,JLOC]);
{	IF ERODIBILITY=ERODEPLUG THEN
	BEGIN
		 IF ELEVATION[ILOC,JLOC]<PLUGMAXELEV THEN
		 BEGIN
			 TEMP1 := ROUND(ELEVATION[ILOC,JLOC]-CHANNELELEV) /
			      ROUND(PLUGMAXELEV-CHANNELELEV);
			 IF TEMP1<0.0 THEN TEMP1 := 0.0;
			 ERODIBILITY := TEMP1*ERODEPLUG + (1.0-TEMP1)*ERODENOMINAL;
		 END;
	END;}
   IF ERODIBILITY<>ERODENOMINAL THEN
   BEGIN
     IF USESTICKY THEN
     BEGIN
        IF (XCHANGE>0.0) AND (ILOC<XPOINTS) THEN
	BEGIN
	   IF (YCHANGE<0.0) AND (JLOC>1) THEN
	   BEGIN
	       IF ABS(BANKERODE[ILOC+1,JLOC-1]) = ERODIBILITY THEN DOCHANGE := TRUE;
	       IF ABS(BANKERODE[ILOC,JLOC-1]) = ERODIBILITY THEN DOCHANGE := TRUE;
	   END;
	   IF ABS(BANKERODE[ILOC+1,JLOC]) = ERODIBILITY THEN DOCHANGE := TRUE;
	   IF (YCHANGE>0.0) AND (JLOC<YPOINTS) THEN
	   BEGIN
	       IF ABS(BANKERODE[ILOC+1,JLOC+1]) = ERODIBILITY THEN DOCHANGE := TRUE;
	       IF ABS(BANKERODE[ILOC,JLOC+1]) = ERODIBILITY THEN DOCHANGE := TRUE;
    	   END;
	END ELSE IF (XCHANGE<0.0) AND (ILOC>1) THEN
	BEGIN
	   IF (YCHANGE<0.0) AND (JLOC>1) THEN
	   BEGIN
	       IF ABS(BANKERODE[ILOC-1,JLOC-1]) = ERODIBILITY THEN DOCHANGE := TRUE;
    	       IF ABS(BANKERODE[ILOC,JLOC-1]) = ERODIBILITY THEN DOCHANGE := TRUE;
	   END;
	   IF ABS(BANKERODE[ILOC-1,JLOC]) = ERODIBILITY THEN DOCHANGE := TRUE;
	   IF (YCHANGE>0.0) AND (JLOC<YPOINTS) THEN
	   BEGIN
	       IF ABS(BANKERODE[ILOC-1,JLOC+1]) = ERODIBILITY THEN DOCHANGE := TRUE;
	       IF ABS(BANKERODE[ILOC,JLOC+1]) = ERODIBILITY THEN DOCHANGE := TRUE;
	   END;
	END ELSE
	BEGIN
	   IF (YCHANGE<0.0) AND (JLOC>1) THEN
	       IF ABS(BANKERODE[ILOC,JLOC-1]) = ERODIBILITY THEN DOCHANGE := TRUE;
	   IF (YCHANGE>0.0) AND (JLOC<YPOINTS) THEN
	       IF ABS(BANKERODE[ILOC,JLOC+1]) = ERODIBILITY THEN DOCHANGE := TRUE;
        END;
     END ELSE
     BEGIN
        IF (XCHANGE>0.0) AND (ILOC<XPOINTS) THEN
	BEGIN
	   IF (YCHANGE<0.0) AND (JLOC>1) THEN
	   BEGIN
	       IF (BANKERODE[ILOC+1,JLOC-1]) = ERODIBILITY THEN DOCHANGE := TRUE;
	       IF (BANKERODE[ILOC,JLOC-1]) = ERODIBILITY THEN DOCHANGE := TRUE;
	   END;
	   IF (BANKERODE[ILOC+1,JLOC]) = ERODIBILITY THEN DOCHANGE := TRUE;
	   IF (YCHANGE>0.0) AND (JLOC<YPOINTS) THEN
	   BEGIN
	       IF (BANKERODE[ILOC+1,JLOC+1]) = ERODIBILITY THEN DOCHANGE := TRUE;
	       IF (BANKERODE[ILOC,JLOC+1]) = ERODIBILITY THEN DOCHANGE := TRUE;
    	   END;
	END ELSE IF (XCHANGE<0.0) AND (ILOC>1) THEN
	BEGIN
	   IF (YCHANGE<0.0) AND (JLOC>1) THEN
	   BEGIN
	       IF (BANKERODE[ILOC-1,JLOC-1]) = ERODIBILITY THEN DOCHANGE := TRUE;
    	       IF (BANKERODE[ILOC,JLOC-1]) = ERODIBILITY THEN DOCHANGE := TRUE;
	   END;
	   IF (BANKERODE[ILOC-1,JLOC]) = ERODIBILITY THEN DOCHANGE := TRUE;
	   IF (YCHANGE>0.0) AND (JLOC<YPOINTS) THEN
	   BEGIN
	       IF (BANKERODE[ILOC-1,JLOC+1]) = ERODIBILITY THEN DOCHANGE := TRUE;
	       IF (BANKERODE[ILOC,JLOC+1]) = ERODIBILITY THEN DOCHANGE := TRUE;
	   END;
	END ELSE
	BEGIN
	   IF (YCHANGE<0.0) AND (JLOC>1) THEN
	       IF (BANKERODE[ILOC,JLOC-1]) = ERODIBILITY THEN DOCHANGE := TRUE;
	   IF (YCHANGE>0.0) AND (JLOC<YPOINTS) THEN
	       IF (BANKERODE[ILOC,JLOC+1]) = ERODIBILITY THEN DOCHANGE := TRUE;
        END;
     END;
	IF DOCHANGE THEN
	BEGIN
	    XCHANGE := XCHANGE * ERODIBILITY;
	    YCHANGE := YCHANGE * ERODIBILITY;
	END ELSE
	BEGIN
	    XCHANGE := XCHANGE * ERODENOMINAL;
	    YCHANGE := YCHANGE * ERODENOMINAL;
	END;
    END ELSE
    BEGIN
	    XCHANGE := XCHANGE * ERODENOMINAL;
	    YCHANGE := YCHANGE * ERODENOMINAL;
    END;
END ELSE BEGIN
	    XCHANGE := XCHANGE * ERODENOMINAL;
	    YCHANGE := YCHANGE * ERODENOMINAL;
END;
END;
END;
{*********************************************************************************}
{*********************************************************************************}

PROCEDURE FINDCENTER (X1, X2, X3, Y1, Y2, Y3 : REAL;
              VAR XCENTER, YCENTER : REAL);

BEGIN
{Finds the central X,Y location defined by three successive stream nodes.}
  IF (((X3 - X1) * (Y2 - Y1)) - ((X2 - X1) * (Y3 - Y1))) = 0.0
    THEN
    YCENTER := (((X2 - X1) * (SQR(X1) - SQR(X3) + SQR(Y1) -
           SQR(Y3)) * 0.5) - ((X3 - X1) * (SQR(X1) - SQR(X2)
    + SQR(Y1) - SQR(Y2)) * 0.5)) / INFINSMALL
  ELSE
    YCENTER := (((X2 - X1) * (SQR(X1) - SQR(X3) + SQR(Y1) -
           SQR(Y3)) * 0.5) - ((X3 - X1) * (SQR(X1) - SQR(X2)
    + SQR(Y1) - SQR(Y2)) * 0.5)) / (((X3 - X1) *
    (Y2 - Y1)) - ((X2 - X1) * (Y3 - Y1)));
  IF (X2 - X1) = 0.0
    THEN
    XCENTER := -(((SQR(X1) - SQR(X2) + SQR(Y1) - SQR(Y2)) * 0.5)
    + ((Y2 - Y1) * YCENTER)) / INFINSMALL
  ELSE
    XCENTER := -(((SQR(X1) - SQR(X2) + SQR(Y1) - SQR(Y2)) * 0.5)
    + ((Y2 - Y1) * YCENTER)) / (X2 - X1);

END; (* FINDCENTER *)
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE ADELETE (LASTONE : BOOLEAN);

BEGIN
{Deletes the next channel node downstream from a given node, and resets the
 X,Y location of the given to be midway between its original position
 and that of the new next downstream node}
  IF NOT LASTONE
    THEN
    BEGIN
      STREAM^.NEXT := STREAM^.NEXT^.NEXT;
      STREAM^.XVALUE := (STREAM^.XVALUE + STREAM^.NEXT^.XVALUE)/2.0;
      STREAM^.YVALUE := (STREAM^.YVALUE + STREAM^.NEXT^.YVALUE)/2.0;
      RIVERSIZE := RIVERSIZE - 1;
    END; (* THEN *)
END; (* ADELETE *)
{*********************************************************************************}
{*********************************************************************************}

PROCEDURE FINDNEWPOINT (VAR XXNEWPOINT, YYNEWPOINT : REAL; RADIUS,
            XCENTER, YCENTER, XMIDPOINT, YMIDPOINT : REAL);

VAR
    SLOPE2 : REAL;

BEGIN
{Determines the X,Y location of a new point inserted between two existing stream
 nodes based upon the central location and radius of curvature}
  IF XMIDPOINT = XCENTER
    THEN
    BEGIN
      XXNEWPOINT := XCENTER;
      IF YMIDPOINT < YCENTER
    THEN
    YYNEWPOINT := YCENTER - RADIUS
      ELSE
    YYNEWPOINT := YCENTER + RADIUS
    END  (* THEN *)
  ELSE
    BEGIN
      SLOPE2 := (YMIDPOINT - YCENTER) / (XMIDPOINT - XCENTER);

      IF XMIDPOINT > XCENTER
        THEN
        BEGIN
          YYNEWPOINT := YCENTER + (RADIUS * (SIN (ARCTAN
                        (SLOPE2))));
          XXNEWPOINT := XCENTER + (RADIUS * (COS (ARCTAN
                        (SLOPE2))));
        END (* THEN *)
      ELSE
        BEGIN
          YYNEWPOINT := YCENTER - (RADIUS * (SIN (ARCTAN
                        (SLOPE2))));
          XXNEWPOINT := XCENTER - (RADIUS * (COS (ARCTAN
                        (SLOPE2))));
        END; (* IF *)
    END; (* IF *)
END; (* FINDNEWPOINT *)
{*********************************************************************************}
{*********************************************************************************}




PROCEDURE INSERTS (FIRSTONE, LASTONE : BOOLEAN; XPRECEDING,
                  YPRECEDING : REAL);


VAR
    XCENTER , YCENTER , XMIDPOINT , YMIDPOINT , RADIUS , XNEWPOINT ,
    YNEWPOINT , X1NEWPOINT , Y1NEWPOINT , X2NEWPOINT , Y2NEWPOINT ,
    ELAPSAVG,RESAVG, FRICAVG, VELAVG, WIDAVG : REAL;
	GRADAVG,DEPTHAVG,FLUXAVG,ELEVAVG,BEDCHANGEAVG, CRVAVG,CSAVG : REAL;
    POINT : DATALINK;

BEGIN
{This is the main code for inserting a new channel node between two existing nodes
 when the inter-node distance has increased beyond a critical value.
 It also assigns values to the variables associated with the new node by averaging
 the values for the existing upstream and downstream nodes.
 Calls FINDCENTER, FINDDISTANCE, FINDNEWPOINT}
  XMIDPOINT := STREAM^.XVALUE + (STREAM^.NEXT^.XVALUE - STREAM^
  .XVALUE) / 2.0;
  YMIDPOINT := STREAM^.YVALUE + (STREAM^.NEXT^.YVALUE - STREAM^
  .YVALUE) /2.0;
  X1 := STREAM^.XVALUE;
  Y1 := STREAM^.YVALUE;
  X2 := STREAM^.NEXT^.XVALUE;
  Y2 := STREAM^.NEXT^.YVALUE;
  ELAPSAVG :=(STREAM^.ELAPSE+STREAM^.NEXT^.ELAPSE)/2.0;
  RESAVG :=(STREAM^.RESISTANCE+STREAM^.NEXT^.RESISTANCE)/2.0;
  GRADAVG :=(STREAM^.LOCALGRADIENT+STREAM^.NEXT^.LOCALGRADIENT)/2.0;
  DEPTHAVG :=(STREAM^.LOCALDEPTH+STREAM^.NEXT^.LOCALDEPTH)/2.0;
  FLUXAVG :=(STREAM^.SEDIMENTFLUX+STREAM^.NEXT^.SEDIMENTFLUX)/2.0;
  ELEVAVG :=(STREAM^.BEDELEVATION+STREAM^.NEXT^.BEDELEVATION)/2.0;
  FRICAVG :=(STREAM^.LOCALCFRICTION+STREAM^.NEXT^.LOCALCFRICTION)/2.0;
  VELAVG  :=(STREAM^.LOCALVELOCITY+STREAM^.NEXT^.LOCALVELOCITY)/2.0;
  WIDAVG  :=(STREAM^.LOCALWIDTH+STREAM^.NEXT^.LOCALWIDTH)/2.0;
  BEDCHANGEAVG  :=(STREAM^.BEDCHANGE+STREAM^.NEXT^.BEDCHANGE)/2.0;
  CRVAVG := (STREAM^.CURVATURE+STREAM^.NEXT^.CURVATURE)/2.0;
  CSAVG := (STREAM^.CSCURVATURE+STREAM^.NEXT^.CSCURVATURE)/2.0;
  IF (FIRSTONE) OR (LASTONE)
    THEN
    BEGIN
      IF FIRSTONE
    THEN
    FINDCENTER (STREAM^.XVALUE, STREAM^.NEXT^.XVALUE,
    STREAM^.NEXT^.NEXT^.XVALUE, STREAM^.YVALUE,
    STREAM^.NEXT^.YVALUE, STREAM^.NEXT^.NEXT^.YVALUE,
    XCENTER, YCENTER)
      ELSE
    FINDCENTER (XPRECEDING, STREAM^.XVALUE,
    STREAM^.NEXT^.XVALUE, YPRECEDING, STREAM^.YVALUE,
    STREAM^.NEXT^.YVALUE, XCENTER, YCENTER);
      FINDDISTANCE (RADIUS, XCENTER, YCENTER, STREAM^.XVALUE,
      STREAM^.YVALUE);
      FINDNEWPOINT (XNEWPOINT, YNEWPOINT, RADIUS, XCENTER,
            YCENTER, XMIDPOINT, YMIDPOINT);
    END (* THEN *)
  ELSE
    BEGIN
      FINDCENTER (XPRECEDING, STREAM^.XVALUE,
      STREAM^.NEXT^.XVALUE, YPRECEDING, STREAM^.YVALUE,
      STREAM^.NEXT^.YVALUE, XCENTER, YCENTER);
      FINDDISTANCE (RADIUS, XCENTER, YCENTER, STREAM^.XVALUE,
      STREAM^.YVALUE);
      FINDNEWPOINT (X1NEWPOINT, Y1NEWPOINT, RADIUS, XCENTER,
            YCENTER, XMIDPOINT, YMIDPOINT);
      FINDCENTER (STREAM^.XVALUE, STREAM^.NEXT^.XVALUE,
      STREAM^.NEXT^.NEXT^.XVALUE, STREAM^.YVALUE,
      STREAM^.NEXT^.YVALUE, STREAM^.NEXT^.NEXT^.YVALUE,
      XCENTER, YCENTER);
      FINDDISTANCE (RADIUS, XCENTER, YCENTER, STREAM^.XVALUE,
      STREAM^.YVALUE);
      FINDNEWPOINT (X2NEWPOINT, Y2NEWPOINT, RADIUS, XCENTER,
            YCENTER, XMIDPOINT, YMIDPOINT);
      XNEWPOINT := (X1NEWPOINT + X2NEWPOINT) / 2.0;
      YNEWPOINT := (Y1NEWPOINT + Y2NEWPOINT) / 2.0;
    END; (*IF*)
  NEW (POINT);
  POINT^.NEXT := STREAM^.NEXT;
  STREAM^.NEXT := POINT;
  POINT^.XVALUE := XNEWPOINT;
  POINT^.YVALUE := YNEWPOINT;
  POINT^.BANKDEPTH := 0.0;
  POINT^.BANKVELOCITY := 0.0;
  POINT^.CURVSENSE := 1;
  POINT^.LOCALGRADIENT :=GRADAVG;
  POINT^.LOCALDEPTH := DEPTHAVG;
  POINT^.SEDIMENTFLUX := FLUXAVG;
  POINT^.BEDELEVATION := ELEVAVG;
  POINT^.LOCALCFRICTION :=FRICAVG;
  POINT^.LOCALVELOCITY :=VELAVG;
  POINT^.LOCALWIDTH := WIDAVG;
  POINT^.BEDCHANGE := BEDCHANGEAVG;
  POINT^.CURVATURE := CRVAVG;
  POINT^.CSCURVATURE := CSAVG;
  {WRITELN(LIST,'NEW, X = ',XNEWPOINT:15,',',' Y =',YNEWPOINT:15,',',' XP =',XPRECEDING:15,',',' YP =',YPRECEDING:15);
  WRITELN(LIST,'    NEW, X1 =',X1:15,' Y1=',Y1:15,',',' X2 =',X2:15,',',' Y2=',Y2:15);}
  IF RES_AVERAGE THEN
  BEGIN
    POINT^.ELAPSE := ELAPSAVG;
    POINT^.RESISTANCE := RESAVG;
  END ELSE
  BEGIN
    POINT^.ELAPSE := 0.0;
    POINT^.RESISTANCE := EXP(ERROR * NORMAL(DUMMY));
  END;
  RIVERSIZE := RIVERSIZE + 1;
END; (* INSERTS *)
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE ADJUST;

VAR
    FIRSTONE , LASTONE : BOOLEAN;
    DISTANCE : REAL;
    XPRECEDING , YPRECEDING : REAL;
    PREVRES, RESAVERAGE, LOCALRES : REAL;

BEGIN
{This cycles through the nodes in the channel database and determines if a node
 needs to be inserted or deleted, and also adjusts the bank resistance if the
 channel has migrated more than one channel width}
  LASTONE := FALSE;
  FIRSTONE := TRUE;
  STREAM := BASE;
  PREVRES := 1.0;
  WHILE STREAM^.
    NEXT <> NIL DO
    BEGIN
      LOCALRES := STREAM^.RESISTANCE;
      IF STREAM^.ELAPSE > CHANNEL_WIDTH THEN
         BEGIN
            STREAM^.ELAPSE := 0.0;
            IF RES_AVERAGE THEN
            BEGIN
              RESAVERAGE := (PREVRES + STREAM^.NEXT^.RESISTANCE)/2.0;
              STREAM^.RESISTANCE := RESAVERAGE *
                   EXP(ERROR * NORMAL(DUMMY));
            END ELSE
              STREAM^.RESISTANCE := EXP(ERROR * NORMAL(DUMMY));
         END;
      DISTANCE := SQRT (SQR (STREAM^.NEXT^.XVALUE - STREAM^
      .XVALUE) + SQR (STREAM^.NEXT^.YVALUE -
      STREAM^.YVALUE));
      IF DISTANCE > (MAX_ADJ_FACTOR * CHANNEL_WIDTH)
    THEN
    INSERTS (FIRSTONE,
        LASTONE, XPRECEDING, YPRECEDING)
      ELSE
    IF DISTANCE < (MIN_ADJ_FACTOR * CHANNEL_WIDTH)
      THEN
      ADELETE (LASTONE);
      FIRSTONE := FALSE;
      XPRECEDING := STREAM^.XVALUE;
      YPRECEDING := STREAM^.YVALUE;
      PREVRES := LOCALRES;
      STREAM := STREAM^.NEXT;
      IF (STREAM^.NEXT <> NIL)
    THEN
    IF (STREAM^.NEXT^.NEXT = NIL)
                 THEN
                 LASTONE := TRUE;
    END; (*WHILE*)
END; (* ADJUST *)
{*********************************************************************************}
{*********************************************************************************}

PROCEDURE CUTOFF (VAR CHECKPOINT : DATALINK; LASTX,LASTY:REAL);

VAR
    TRAVERSER, MINCHECKPOINT : DATALINK;
    DELCOUNT : LONGINT;
    X1,X2,X3,Y1,Y2,Y3,DS,DU,DD,DCROSS,ANGLE,CURV,MINCURV : REAL;
    OLDCURV,XBASE,YBASE : REAL;
    OK1,OK2 : BOOLEAN;

BEGIN
{This procedure implements neck cutoffs by eliminating channel nodes of the
 abandoned loop.  The procedure is initiated with the X,Y location
 of the current node and of the downstream location (CHECKPOINT)
 that will become the new next downstream location}
  TRAVERSER := CHECKPOINT;
  X1 := STREAM^.XVALUE;
  X2 := CHECKPOINT^.XVALUE;
  Y1 := STREAM^.YVALUE;
  Y2 := CHECKPOINT^.YVALUE;
IF DETAILPAR THEN WRITELN(LIST,'ORIGINAL: X1=',X1:10,' Y1=',Y1:10,
     ' X2=',X2:10,' Y2=',Y2:10);
  XBASE := X1;
  YBASE := Y1;
  FINDDISTANCE(DS,XBASE,YBASE,X1,Y1);
  FINDDISTANCE(DU,X1,Y1,X2,Y2);
  MINCURV := 1.0E+15;
  OK1 := TRUE;
  OK2 := TRUE;
  REPEAT
    IF TRAVERSER^.NEXT <> NIL THEN
    BEGIN
       X3 := TRAVERSER^.NEXT^.XVALUE;
       Y3 := TRAVERSER^.NEXT^.YVALUE;
       FINDDISTANCE(DCROSS,X1,Y1,X3,Y3);
       CURVANGLE(XBASE,X1,X2,YBASE,Y2,Y3,ANGLE);
       OLDCURV := ABS(2.0*ANGLE/(DS+DU));
       FINDDISTANCE(DD,X2,Y2,X3,Y3);
       CURVANGLE(X1,X2,X3,Y1,Y2,Y3,ANGLE);
       CURV := ABS(2.0*ANGLE/(DU+DD));
       WRITELN(OUTCAPT1,ITERNUM:6,' ',X3:10:2,
       ' ',Y3:10:2,' ',OLDCURV:10:4,
       ' ',CURV:10:4,' ',MINCURV:10:4,' ',DCROSS:10:2);
       IF (OLDCURV + CURV) < MINCURV THEN
       BEGIN
         MINCURV := OLDCURV + CURV;
         CHECKPOINT := TRAVERSER;
     WRITELN(OUTCAPT1,CHECKPOINT^.XVALUE:10:2,' ',
       CHECKPOINT^.YVALUE:10:2);
       END;
       TRAVERSER := TRAVERSER^.NEXT;
       X2 := X3;
       Y2 := Y3;
       DU := DD;
       IF DCROSS > (CHANNEL_WIDTH * LOOPCHECKFACTOR) THEN OK1 := FALSE;
    END ELSE OK2 := FALSE;
  UNTIL (NOT OK1) OR (NOT OK2);
  IF OK2 THEN
  BEGIN
     X1 := STREAM^.XVALUE;
     Y1 := STREAM^.YVALUE;
     STREAM^.XVALUE := (LASTX+X1+CHECKPOINT^.XVALUE)/3.0;
     STREAM^.YVALUE := (LASTY+Y1+CHECKPOINT^.YVALUE)/3.0;
     IF CHECKPOINT^.NEXT<>NIL THEN
     BEGIN
       CHECKPOINT^.XVALUE :=(X1+CHECKPOINT^.XVALUE
                +CHECKPOINT^.NEXT^.XVALUE)
                            /3.0;
       CHECKPOINT^.YVALUE :=(Y1+CHECKPOINT^.YVALUE
                +CHECKPOINT^.NEXT^.YVALUE)
                            /3.0;
     END;
     WRITELN(OUTCAPT1,STREAM^.XVALUE:10:2,
     ' ',STREAM^.YVALUE:10:2,
     ' ',CHECKPOINT^.XVALUE:10:2,' ',CHECKPOINT^.YVALUE:10:2);
  END;
  DELCOUNT := -1;
  TRAVERSER := STREAM;
  REPEAT
    TRAVERSER := TRAVERSER^.NEXT;
    DELCOUNT := DELCOUNT + 1;
  UNTIL TRAVERSER = CHECKPOINT;
  STREAM^.NEXT := CHECKPOINT;
  RIVERSIZE := ABS(RIVERSIZE - DELCOUNT);
  WRITELN (OUTCAPT1,DELCOUNT:5, ' ',
       ' ',RIVERSIZE:5);
  WRITELN(OUTCAPT1);
END; (* CUTOFF *)
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE DOCUTOFF(STREAM,CHECKPOINT : DATALINK; LASTX,LASTY,DISTANCE : REAL);
VAR
	XS,YS,XE,YE,SLP,INTR,ANGL,XINC,YN,XN,YINC,LDIST,LCHECK : REAL;
	STILLPOINTS, STILLOK : BOOLEAN;
	TRAVERSER, LASTTRAVERSER, POINT : DATALINK;
	I,DELCOUNT,NPOINTS : LONGINT;
	CHECKLASTX,CHECKLASTY : REAL;

BEGIN
{This procedure implements chute cutoffs by eliminating the cutoff nodes
 and inserting new nodes along the cutoff path and assigns
 variable values associated with the new nodes}
	XS := STREAM^.XVALUE;
	YS := STREAM^.YVALUE;
	XE := CHECKPOINT^.XVALUE;
	YE := CHECKPOINT^.YVALUE;
	CHECKLASTX := XE;
	CHECKLASTY := YE;
	NPOINTS := 0;
	TRAVERSER :=STREAM;
	STILLOK := TRUE;
	LCHECK := MIN_ADJ_FACTOR*CHANNEL_WIDTH;
	LASTTRAVERSER := STREAM;
	IF TRAVERSER <> CHECKPOINT THEN STILLPOINTS := TRUE
		ELSE STILLPOINTS := FALSE;
	FINDDISTANCE(LDIST,XS,YS,XE,YE);
{if distance<2.0 then
begin
  writeln('i=1, x=',xs:10:3,' y=',ys:10:3);
  i := 1;
  while traverser <> checkpoint do
  begin
     traverser := traverser^.next;
     i := i+1;
     writeln('i=',i,' x=',traverser^.xvalue:10:3,' y=',traverser^.yvalue:10:3);
  end;
  traverser := stream;
end;}
	IF LDIST<LCHECK THEN STILLOK := FALSE;
      IF STILLOK THEN
      BEGIN
	IF ABS(XS-XE)>0.00001 THEN
	BEGIN
		SLP := (YE-YS)/(XE-XS);
		INTR := YS - SLP * XS;
		ANGL := ARCTAN(SLP);
		XINC := COS(ANGL) * MEANLENGTH;
{if distance<2.0 then writeln('slp=',slp:10:5,' xinc=',XINC:10:5);}
		IF ((XINC>0.0) AND ((XINC+XS) < XE)) OR ((XINC<0.0) AND ((XINC+XS) > XE))
		THEN BEGIN
			REPEAT
				IF STILLPOINTS THEN
				BEGIN
					LASTTRAVERSER := TRAVERSER;
					TRAVERSER :=TRAVERSER^.NEXT;
{if distance<2.0 then writeln('traverser advanced');}
				END;
				XN := XS + XINC;
				YN := XS * SLP + INTR;
{if distance<2.0 then writeln('xn=',xn:10:3,' yn=',yn:10:3);}
				FINDDISTANCE(LDIST,XN,YN,XE,YE);
				IF LDIST<LCHECK THEN
				BEGIN
					STILLOK := FALSE;
					TRAVERSER := LASTTRAVERSER;
{if distance<2.0 then writeln('not stillok - traverser decremented');}
				END;
			    IF STILLOK THEN
			    BEGIN
				NPOINTS := NPOINTS+1;
				IF TRAVERSER <> CHECKPOINT THEN
				BEGIN
{if distance<2.0 then writeln('reassigned x at',XN:10,' y',YN:10);}
					TRAVERSER^.XVALUE := XN;
					TRAVERSER^.YVALUE := YN;
					CHECKLASTX := XN;
					CHECKLASTY := YN;
				END ELSE BEGIN
					NEW(POINT);
{if distance<2.0 then writeln('new x at',XN:10,' y',YN:10);}
  					POINT^.NEXT := TRAVERSER;
  					LASTTRAVERSER^.NEXT := POINT;
  					POINT^.XVALUE := XN;
  					POINT^.YVALUE := YN;
					CHECKLASTX := XN;
					CHECKLASTY := YN;
					POINT^.BANKDEPTH := 0.0;
					POINT^.BANKVELOCITY := 0.0;
					POINT^.CURVSENSE := 1;
					POINT^.LOCALGRADIENT := 0.0;
					POINT^.LOCALDEPTH := TRAVERSER^.LOCALDEPTH;
					POINT^.LOCALCFRICTION := TRAVERSER^.LOCALCFRICTION;
					POINT^.LOCALVELOCITY := TRAVERSER^.LOCALVELOCITY;
					POINT^.SEDIMENTFLUX := 0.0;
					POINT^.BEDELEVATION := TRAVERSER^.BEDELEVATION;
					POINT^.LOCALWIDTH := CHANNEL_WIDTH;
					POINT^.BEDCHANGE := TRAVERSER^.BEDCHANGE;
					POINT^.CURVATURE := TRAVERSER^.CURVATURE;
					POINT^.CSCURVATURE := TRAVERSER^.CSCURVATURE;
  					IF RES_AVERAGE THEN
  					BEGIN
    					POINT^.ELAPSE := TRAVERSER^.ELAPSE;
    					POINT^.RESISTANCE := TRAVERSER^.RESISTANCE;
  					END ELSE
  					BEGIN
    					POINT^.ELAPSE := 0.0;
    					POINT^.RESISTANCE := EXP(ERROR * NORMAL(DUMMY));
  					END;
  					RIVERSIZE := RIVERSIZE + 1;
					STILLPOINTS := FALSE;
					LASTTRAVERSER := POINT;
				END;
				XS := XN;
			    END;
			UNTIL (((XINC>0.0) AND (XN>=XE)) OR ((XINC<0.0) AND (XN<=XE))) OR (NOT STILLOK);
			IF (TRAVERSER <> CHECKPOINT) AND STILLPOINTS THEN
			BEGIN
				DELCOUNT := -1;
				POINT := TRAVERSER;
				REPEAT
					POINT := POINT^.NEXT;
					DELCOUNT := DELCOUNT + 1;
				UNTIL POINT = CHECKPOINT;
{if distance<2.0 then writeln(delcount,' points deleted');}
				TRAVERSER^.NEXT := CHECKPOINT;
				CHECKLASTX := TRAVERSER^.XVALUE;
				CHECKLASTY := TRAVERSER^.YVALUE;
				RIVERSIZE := ABS(RIVERSIZE - DELCOUNT);
			END;
		END;
	END ELSE
	BEGIN
		YINC := MEANLENGTH * (YE-YS)/ABS(YE-YS);
		IF ((YINC>0.0) AND ((YINC+YS) < YE)) OR ((YINC<0.0) AND ((YINC+YS) > YE))
		THEN BEGIN
			REPEAT
				IF STILLPOINTS THEN
				BEGIN
					LASTTRAVERSER := TRAVERSER;
					TRAVERSER :=TRAVERSER^.NEXT;
				END;
				XN := XS;
				YN := YS+ YINC;
				FINDDISTANCE(LDIST,XN,YN,XE,YE);
				IF LDIST<LCHECK THEN
				BEGIN
					STILLOK := FALSE;
					TRAVERSER := LASTTRAVERSER;
				END;
				NPOINTS := NPOINTS+1;
			    IF STILLOK THEN
			    BEGIN
				IF TRAVERSER <> CHECKPOINT THEN
				BEGIN
					TRAVERSER^.XVALUE := XN;
					TRAVERSER^.YVALUE := YN;
					CHECKLASTX := XN;
					CHECKLASTY := YN;
				END ELSE BEGIN
					NEW(POINT);
  					POINT^.NEXT := TRAVERSER;
  					LASTTRAVERSER^.NEXT := POINT;
  					POINT^.XVALUE := XN;
  					POINT^.YVALUE := YN;
					CHECKLASTX := XN;
					CHECKLASTY := YN;
					POINT^.BANKDEPTH := 0.0;
					POINT^.BANKVELOCITY := 0.0;
					POINT^.CURVSENSE := 1;
		            POINT^.LOCALGRADIENT := 0.0;
					POINT^.LOCALDEPTH := TRAVERSER^.LOCALDEPTH;
					POINT^.LOCALCFRICTION := TRAVERSER^.LOCALCFRICTION;
					POINT^.LOCALVELOCITY := TRAVERSER^.LOCALVELOCITY;
					POINT^.SEDIMENTFLUX := 0.0;
					POINT^.LOCALWIDTH := CHANNEL_WIDTH;
					POINT^.BEDELEVATION := TRAVERSER^.BEDELEVATION;
					POINT^.CURVATURE := TRAVERSER^.CURVATURE;
					POINT^.CSCURVATURE := TRAVERSER^.CSCURVATURE;
					POINT^.BEDCHANGE := 0.0;
  					IF RES_AVERAGE THEN
  					BEGIN
    					POINT^.ELAPSE := TRAVERSER^.ELAPSE;
    					POINT^.RESISTANCE := TRAVERSER^.RESISTANCE;
  					END ELSE
  					BEGIN
    					POINT^.ELAPSE := 0.0;
    					POINT^.RESISTANCE := EXP(ERROR * NORMAL(DUMMY));
  					END;
  					RIVERSIZE := RIVERSIZE + 1;
					STILLPOINTS := FALSE;
					LASTTRAVERSER := POINT;
				END;
				YS := YN;
			    END;
			UNTIL (((YINC>0.0) AND (YN>=YE)) OR ((YINC<0.0) AND (YN<=YE))) OR (NOT STILLOK);
			IF (TRAVERSER <> CHECKPOINT) AND STILLPOINTS THEN
			BEGIN
				DELCOUNT := -1;
				POINT := TRAVERSER;
				REPEAT
					POINT := POINT^.NEXT;
					DELCOUNT := DELCOUNT + 1;
				UNTIL POINT = CHECKPOINT;
				TRAVERSER^.NEXT := CHECKPOINT;
				CHECKLASTX := TRAVERSER^.XVALUE;
				CHECKLASTY := TRAVERSER^.YVALUE;
				RIVERSIZE := ABS(RIVERSIZE - DELCOUNT);
			END;
		END;
	END;
   FINDDISTANCE(LDIST,LASTX,LASTY,STREAM^.XVALUE,STREAM^.YVALUE);
   IF LDIST<MAX_ADJ_FACTOR*CHANNEL_WIDTH THEN
   BEGIN
      STREAM^.XVALUE := (LASTX+STREAM^.XVALUE+STREAM^.NEXT^.XVALUE)/3.0;
      STREAM^.YVALUE := (LASTY+STREAM^.YVALUE+STREAM^.NEXT^.YVALUE)/3.0;
   END ELSE
      WRITELN('PROBLEM IN STREAM AVERAGING: DIST=',LDIST:10:4,'  XU=',LASTX:10:3,' YU=',LASTY:10:3,
          ' XSTR=',STREAM^.XVALUE:10:3,' YSTR=',STREAM^.YVALUE:10:3);

   IF CHECKPOINT^.NEXT <> NIL THEN
   BEGIN
      FINDDISTANCE(LDIST,CHECKLASTX,CHECKLASTY,CHECKPOINT^.XVALUE,CHECKPOINT^.YVALUE);
      IF LDIST<MAX_ADJ_FACTOR*CHANNEL_WIDTH THEN
      BEGIN
        TRAVERSER := CHECKPOINT;
        TRAVERSER^.XVALUE := (CHECKLASTX+TRAVERSER^.XVALUE+TRAVERSER^.NEXT^.XVALUE)/3.0;
        TRAVERSER^.YVALUE := (CHECKLASTY+TRAVERSER^.YVALUE+TRAVERSER^.NEXT^.YVALUE)/3.0;
      END ELSE
        WRITELN('PROBLEM IN CHECKPOINT AVERAGING: DIST=',LDIST:10:4,' XU=',CHECKLASTX:10:3,' YU=',CHECKLASTY:10:3,
          ' XCHK=',CHECKPOINT^.XVALUE:10:3,' YCHK=',CHECKPOINT^.YVALUE:10:3);
   END;
   END;
END;
{*********************************************************************************}
{*********************************************************************************}

PROCEDURE CHECKPATH(STREAM,CHECKPOINT: DATALINK;LASTX,LASTY,
			NEXTX,NEXTY,DISTANCE : REAL; VAR PROBFACT,
			PATHLENGTH, MAXPERPDIST : REAL; VAR MAXIMUMELEV : LONGINT;
			VAR AVERAGEELEV,
			THEANGLE, ANGLECORRECT, VELOCITYCORRECT, ELEVATIONCORRECT
			: REAL);
VAR
	LOCATION : DATALINK;
	PERP : REAL;
	NTERMS, IBEGIN, JBEGIN, IFINISH, JFINISH : LONGINT;
    DX,DY,IX,IY,IDX,IDY,INCR,XX,YY,IPX,IPY, I : LONGINT;
    XS,YS,XE,YE,NEWX,NEWY,XL,YL,XN,YN : REAL;
    FOUNDVALLEY,PLOT : BOOLEAN;
    NEWDISTANCE,EXPWEIGHT,
    	RANDCOMP,LOCDIST : REAL;
BEGIN
{This procedure checks the probability of a chute cutoff will occur between
 two nodes along the centerline}
        FOUNDVALLEY := FALSE;
IF USEDEPOSIT THEN
BEGIN
	IBEGIN := ROUND(STREAM^.XVALUE/CHANNEL_WIDTH) + ISTART;
	JBEGIN := ROUND(STREAM^.YVALUE/CHANNEL_WIDTH) + JCENTER;
	IFINISH := ROUND(CHECKPOINT^.XVALUE/CHANNEL_WIDTH) + ISTART;
	JFINISH := ROUND(CHECKPOINT^.YVALUE/CHANNEL_WIDTH) + JCENTER;
END;
	PROBFACT := -2.0;
	XS := STREAM^.XVALUE;
	YS := STREAM^.YVALUE;
	XE := CHECKPOINT^.XVALUE;
	YE := CHECKPOINT^.YVALUE;
	PATHLENGTH := 0.0;
	LOCATION := STREAM^.NEXT;
	XL := LOCATION^.XVALUE;
	YL := LOCATION^.YVALUE;
	FINDDISTANCE(LOCDIST,XS,YS,XL,YL);
	PATHLENGTH := PATHLENGTH+LOCDIST;
	LOCATION := STREAM^.NEXT;
	XN := LOCATION^.XVALUE;
	YN := LOCATION^.YVALUE;
	FINDDISTANCE(LOCDIST,XL,YL,XN,YN);
	PATHLENGTH := PATHLENGTH+LOCDIST;
	XL := XN;
	YL := YN;
	MAXPERPDIST := 0.0;
	REPEAT
		NEWX := LOCATION^.XVALUE;
		NEWY := LOCATION^.YVALUE;
		FINDPERPDIS(PERP,NEWX,NEWY,XS,YS,XE,YE);
		IF PERP>MAXPERPDIST THEN MAXPERPDIST := PERP;
		LOCATION := LOCATION^.NEXT;
		XN := LOCATION^.XVALUE;
		YN := LOCATION^.YVALUE;
		FINDDISTANCE(LOCDIST,XL,YL,XN,YN);
		PATHLENGTH := PATHLENGTH+LOCDIST;
		XL := XN;
		YL := YN;
	UNTIL LOCATION = CHECKPOINT;
IF MAXPERPDIST>(MINPERPDIST*CHANNEL_WIDTH) THEN
BEGIN
{writeln(' maxperpdist');}
	NTERMS := 0;
	MAXIMUMELEV := -20000;
	AVERAGEELEV := 0.0;
IF USEDEPOSIT THEN
BEGIN
	IF (IBEGIN >= 1) AND (IBEGIN <= XPOINTS) AND
		(IFINISH >= 1) AND (IFINISH <= XPOINTS) AND
		(JBEGIN >= 1) AND (JBEGIN <= YPOINTS) AND
		(JFINISH >= 1) AND (JFINISH <= YPOINTS) THEN
	BEGIN
		DX := IFINISH - IBEGIN;
		DY := JFINISH - JBEGIN;
		IF DX >= 0 THEN IDX := 1 ELSE IDX := -1;
		IF DY >= 0 THEN IDY := 1 ELSE IDY := -1;
		IX := ABS(DX);
		IY := ABS(DY);
		IF IX > IY THEN INCR := IX ELSE INCR := IY;
		IPX := IBEGIN;
		IPY := JBEGIN;
		XX := 0;
		YY := 0;
		FOR I := 0 TO INCR DO
		BEGIN
			XX := XX + IX;
			YY := YY + IY;
			PLOT := FALSE;
			IF XX > INCR THEN
			BEGIN
				PLOT := TRUE;
				XX := XX - INCR;
				IPX := IPX + IDX;
			END;
			IF YY > INCR THEN
			BEGIN
				PLOT := TRUE;
				YY := YY - INCR;
				IPY := IPY + IDY;
			END;
			IF PLOT THEN
			BEGIN
				IF USEVALLEY AND (BANKERODE[IPX,IPY] = ERODEVALLEY) THEN FOUNDVALLEY := TRUE;
				IF ELEVATION[IPX,IPY] > MAXIMUMELEV THEN MAXIMUMELEV :=
					ELEVATION[IPX,IPY];
				AVERAGEELEV := AVERAGEELEV + ELEVATION[IPX,IPY];
				NTERMS := NTERMS + 1;
			END;
		END;
	END;
	IF NTERMS > 0 THEN AVERAGEELEV := AVERAGEELEV / NTERMS;
END;
	INCLUDEANGLE(LASTX,NEXTX,CHECKPOINT^.XVALUE,LASTY,NEXTY,
		CHECKPOINT^.YVALUE, THEANGLE);
	IF USEANGLEFACTOR THEN ANGLECORRECT := -COS(ABS(THEANGLE))*ANGLEFACTOR
		ELSE ANGLECORRECT := 0.0;
	IF USEVELOCITYFACTOR THEN
	BEGIN
		IF THEANGLE <> 0.0 THEN
			VELOCITYCORRECT := STREAM^.CURVSENSE*(THEANGLE/ABS(THEANGLE))*
				VELOCITYFACTOR * STREAM^.BANKVELOCITY
		ELSE
			VELOCITYCORRECT := 0.0
	END ELSE
		VELOCITYCORRECT := 0.0;
	ELEVATIONCORRECT := 0.0;
	IF USEMAXIMUMELEVATION THEN
		IF MAXIMUMELEV>CHANNELELEV THEN
			ELEVATIONCORRECT := (MAXIMUMELEV-CHANNELELEV)	*
				ELEVATIONFACTOR
		ELSE
			ELEVATIONCORRECT := 0.0;
	IF USEAVERAGEELEVATION THEN
		IF AVERAGEELEV>CHANNELELEV THEN
			ELEVATIONCORRECT := (AVERAGEELEV- CHANNELELEV)	*
				ELEVATIONFACTOR
		ELSE
			ELEVATIONCORRECT := 0.0;
	NEWDISTANCE := PROBMULT*(PATHLENGTH-DISTANCE)/PATHLENGTH;
	IF NEWDISTANCE > 0.0 THEN
	BEGIN
		EXPWEIGHT := DISTANCE * DISTANCEFACTOR/CHANNEL_WIDTH + VELOCITYCORRECT + ANGLECORRECT
		+ ELEVATIONCORRECT;
		IF EXPWEIGHT >0.0 THEN
			IF EXPWEIGHT<15.0 THEN
				PROBFACT := NEWDISTANCE*EXP(-EXPWEIGHT)
			ELSE
				PROBFACT := 0.0
		ELSE
			PROBFACT := NEWDISTANCE;
	END;
	IF FOUNDVALLEY THEN PROBFACT := 0.0;
END;
END;
{*********************************************************************************}
{*********************************************************************************}

PROCEDURE CHECKLOOP;


VAR
    CHECKPOINT, MAXCHECKPOINT,MAXSTARTPOINT : DATALINK;
    DISTANCE,RANDCOMP : REAL;
    K, NCHECK,NCONSIDER : LONGINT;
    NECKCUTOFF,OK1,OK2 ,FIRSTCASE, LOOPSFOUND,MAKECUTOFF,ENDOFLOOP,CHUTESFOUND : BOOLEAN;
    NSKIP,KK : LONGINT;
    LASTX, LASTY,NEXTX,NEXTY,PROBFACT,MAXPROBFACT,MAXDISTANCE,MAXSTRX,MAXSTRY,
       MAXNEXTX,MAXNEXTY,LOOPDISTANCE ,STRX,STRY: REAL;
	 PTHL,MPDIST,AVGE,THEANG,ACORR,VCORR,ECORR : REAL;
	 MAXPTHL,MAXMPDIST,MAXAVGE,MAXTHEANG,MAXACORR,MAXVCORR,MAXECORR : REAL;
	 MAXSENSE,MAXMAXE,MAXE : LONGINT;
	 MAXVEL,MAXENDX,MAXENDY,MAXLASTX,MAXLASTY : REAL;
	 CLOSEST_APPROACH, CLOSEST_COMP : real;

BEGIN
{; This cycles downstream from node to node and determines if a neck
 or chute cutoff event will occur and calls the appropriate routine}
  NCHECK:= 0;
  LOOPSFOUND := FALSE;
  CHUTESFOUND := FALSE;
  FIRSTCASE := TRUE;
  CLOSEST_APPROACH := 1.0E25;
  OK1 := TRUE;
  STREAM := BASE;
  LASTX := STREAM^.XVALUE;
  LASTY := STREAM^.YVALUE;
  STRX := LASTX;
  STRY := LASTY;
  NEXTX := STREAM^.NEXT^.XVALUE;
  NEXTY := STREAM^.NEXT^.YVALUE;
  WHILE (STREAM^.NEXT <> NIL) AND OK1 DO
  BEGIN
        NCHECK := NCHECK+1;
		NCONSIDER := 0;
		CHECKPOINT := STREAM;
		NECKCUTOFF := FALSE;
		MAXPROBFACT := -1.0;
		MAXCHECKPOINT := CHECKPOINT;
		MAXDISTANCE := 1.0E+25;
		CLOSEST_COMP := CHANNEL_WIDTH*LOOPCHECKFACTOR;
		FOR K := 1 TO 3 DO
			IF CHECKPOINT^.NEXT <> NIL
			THEN
				CHECKPOINT := CHECKPOINT^.NEXT
			ELSE
				OK1 := FALSE;
		IF OK1 THEN
		BEGIN
			OK2 := TRUE;
			WHILE OK2 DO
			BEGIN
			    NCONSIDER := NCONSIDER+1;
				DISTANCE := SQRT ((SQR (CHECKPOINT^.XVALUE -
					STREAM^.XVALUE)) + (SQR (CHECKPOINT^.
					YVALUE - STREAM^.YVALUE)));
                                {IF (DISTANCE<(2.0*CHANNEL_WIDTH)) THEN WRITELN('D=',DISTANCE:15);}
				IF DISTANCE <CLOSEST_APPROACH THEN CLOSEST_APPROACH := DISTANCE;
				IF DISTANCE > MAX_DISTANCE  THEN OK2 := FALSE;
				IF DOSKIPPOINTS THEN
				BEGIN
			          NSKIP := TRUNC(DISTANCE/MEANLENGTH) -2;
				  IF NSKIP < 1 THEN NSKIP :=1;
				  IF NSKIP > 4 THEN NSKIP :=4;
				END ELSE
				  NSKIP := 1;
				MAKECUTOFF := FALSE;
				IF DISTANCE <= (CHANNEL_WIDTH * LOOPCHECKFACTOR) THEN
				BEGIN
					MAKECUTOFF := TRUE;
					WRITELN('CUTOFF');
				END ELSE BEGIN
					IF USEPROBCUTOFF THEN
					BEGIN
						CHECKPATH(STREAM,CHECKPOINT,STRX,STRY,
						NEXTX,NEXTY,DISTANCE,PROBFACT,PTHL,MPDIST,
						MAXE,AVGE,THEANG,ACORR,VCORR,ECORR);
						IF PROBFACT > MAXPROBFACT THEN
						BEGIN
      {  FINDDISTANCE(LOOPDISTANCE,STREAM^.XVALUE,STREAM^.YVALUE,CHECKPOINT^.XVALUE,CHECKPOINT^.YVALUE);
						IF LOOPDISTANCE> 2.5*CHANNEL_WIDTH THEN
						BEGIN}
							MAXPROBFACT := PROBFACT;
							MAXCHECKPOINT := CHECKPOINT;
							MAXDISTANCE := DISTANCE;
							MAXSTARTPOINT := STREAM;
							MAXLASTX := LASTX;
							MAXLASTY := LASTY;
							MAXSTRX := STRX;
							MAXSTRY := STRY;
							MAXNEXTX := NEXTX;
							MAXNEXTY := NEXTY;
							MAXENDX := CHECKPOINT^.XVALUE;
							MAXENDY := CHECKPOINT^.YVALUE;
							MAXMPDIST := MPDIST;
							MAXPTHL := PTHL;
							MAXMAXE := MAXE;
							MAXAVGE := AVGE;
							MAXTHEANG := THEANG;
							MAXACORR := ACORR;
							MAXVCORR := VCORR;
							MAXECORR := ECORR;
							MAXSENSE := STREAM^.CURVSENSE;
							MAXVEL := STREAM^.BANKVELOCITY;
					       {	END;  }
						END;
					END;
				END;
				IF MAKECUTOFF THEN
				BEGIN
					LOOPSFOUND := TRUE;
					IF FIRSTCASE THEN
					BEGIN
						FIRSTCASE := FALSE;
					END;(*IF*)
					CUTOFF (CHECKPOINT,LASTX,LASTY);
					NECKCUTOFF := TRUE;
					FOR K := 1 TO 3 DO
						IF CHECKPOINT^.NEXT <> NIL	THEN
							CHECKPOINT := CHECKPOINT^.NEXT
						ELSE
							OK2 := FALSE;
{writeln('neck cutoff, x=',stream^.xvalue:10:3,' y=',stream^.yvalue:10:3);
writeln('new checkpoint, x=',checkpoint^.xvalue:10:3,' y=',checkpoint^.yvalue:10:3);}
				END (* THEN*)
				ELSE
				BEGIN
					KK := 0;
					REPEAT
						IF CHECKPOINT^.NEXT <> NIL THEN
							CHECKPOINT := CHECKPOINT^.NEXT
						ELSE
							OK2 := FALSE;
							KK := KK + 1;
					UNTIL (KK = NSKIP) OR (NOT OK2);
				END; (*ELSE*)
			END; (* WHILE *)
			IF (NOT MAKECUTOFF) AND (NOT NECKCUTOFF) THEN
			BEGIN
				IF MAXPROBFACT > 0.0 THEN
				BEGIN
					RANDCOMP := RAN3(DUMMY) {RANDOM(0)};
					IF RANDCOMP < MAXPROBFACT THEN
					BEGIN
writeln('cutoff from x=',STREAM^.XVALUE:12,' y=',STREAM^.YVALUE:12,' to x=',
MAXCHECKPOINT^.XVALUE:12,' y=',MAXCHECKPOINT^.YVALUE:12);
writeln('checkpoint x=',CHECKPOINT^.XVALUE:12,' y=',CHECKPOINT^.YVALUE:12);
writeln('randcomp=',RANDCOMP:12,' maxprobfact=',MAXPROBFACT:12,' maxdistance=',
MAXDISTANCE:12);
if maxstartpoint <> stream then writeln('PROBLEM : MAXSTARTPOINT <> STREAM, MXS x=',
maxstartpoint^.xvalue:10:3,' y=',maxstartpoint^.yvalue:10:3);
						DOCUTOFF(STREAM,MAXCHECKPOINT,MAXLASTX,MAXLASTY,MAXDISTANCE);
						WRITELN(OUTCAPT2,ITERNUM:6,' ',MAXSTRX:10:4,' ',MAXSTRY:10:4,' ',
						  MAXNEXTX:10:4,' ',MAXNEXTY:10:4,' ',MAXENDX:10:4,' ',MAXENDY:10:4);
						WRITELN(OUTCAPT2,MAXPROBFACT:10:5,' ',MAXDISTANCE:10:4,' ',MAXPTHL:10:4,
						  ' ',MAXMPDIST:10:4);
						WRITELN(OUTCAPT2,MAXMAXE/1000.0:10:4,' ',MAXAVGE/1000.0:10:4,' ',MAXECORR:10:4);
						WRITELN(OUTCAPT2,MAXTHEANG/0.01745329:10:2,' ',MAXACORR:10:4);
						WRITELN(OUTCAPT2,MAXVEL:10:4,' ',MAXSENSE:5,' ',MAXVCORR:10:4);
						WRITELN(OUTCAPT2);
						STREAM := MAXCHECKPOINT;
						CHECKPOINT := MAXCHECKPOINT;
						STRX := STREAM^.XVALUE;
						STRY := STREAM^.YVALUE;
						IF STREAM^.NEXT <> NIL THEN
						BEGIN
							NEXTX := STREAM^.NEXT^.XVALUE;
							NEXTY := STREAM^.NEXT^.YVALUE;
						END ELSE
							OK1 := FALSE;
writeln('new stream x=',STREAM^.XVALUE:12,' y=',STREAM^.YVALUE:12);
						LOOPSFOUND := TRUE;
						IF FIRSTCASE THEN
						BEGIN
							FIRSTCASE := FALSE;
						END;(*IF*)
						FOR K := 1 TO 3 DO
							IF CHECKPOINT^.NEXT <> NIL	THEN
								CHECKPOINT := CHECKPOINT^.NEXT
							ELSE
								OK1 := FALSE;
writeln('new checkpoint x=',CHECKPOINT^.XVALUE:12,' y=',CHECKPOINT^.YVALUE:12);
					END;
				END;
			END;
			LASTX := STREAM^.XVALUE;
			LASTY := STREAM^.YVALUE;
			IF STREAM^.NEXT <> NIL THEN
			BEGIN
				STREAM := STREAM^.NEXT;
				STRX := STREAM^.XVALUE;
				STRY := STREAM^.YVALUE;
				IF STREAM^.NEXT <> NIL THEN
				BEGIN
					NEXTX := STREAM^.NEXT^.XVALUE;
					NEXTY := STREAM^.NEXT^.YVALUE;
				END ELSE
					OK1 := FALSE;
			END ELSE
				OK1 := FALSE;
		END;(*IF*)
		{WRITELN('NCONSIDER=',NCONSIDER:10);}
	END; (*WHILE*)
	{WRITELN('CLOSEST_APPROACH=',CLOSEST_APPROACH:15,' NEED TO BE CLOSER THAN =',CLOSEST_COMP:15,' NCHECK=',NCHECK:10);}
	IF LOOPSFOUND THEN
	BEGIN
	END; (*IF*)
END; (* CHECKLOOP *)
{*********************************************************************************}
{*********************************************************************************}

PROCEDURE INITIALIZE;



VAR
    X , Y,EL, RE, ELEV : REAL;{}
    I : LONGINT;


BEGIN
{This reads in the initial values of the X,Y locations of stream nodes,
 the elapsed time at the beginning of the simulation, and the initial bank resistance}
  YAVERAGE := 0.0;
  BASE := NIL;
  WHILE NOT EOF(INDATA) DO
    BEGIN
      READ (INDATA,X);
      IF NOT EOF (INDATA)
    THEN
    BEGIN
      READ (INDATA,Y);
      {READ (INDATA,EL);
      READ (INDATA,RE);
	  READ (INDATA,ELEV);}{}
      IF EOF (INDATA)
        THEN
        WRITELN (LIST,'WARNING : ODD NO. OF ENTRIES')
      ELSE
        BEGIN
          IF (INPUT_PRINT > 0) THEN
        WRITELN (LIST,'INDATA : ',X:10,'   ',Y:10);
          YAVERAGE := YAVERAGE + Y;
          NEW (STREAM);
          RIVERSIZE := RIVERSIZE + 1;
          STREAM^.XVALUE := X;
          STREAM^.YVALUE := Y;
          STREAM^.ELAPSE := 0.0;
          STREAM^.RESISTANCE := 1.0;
          STREAM^.NEXT := BASE;
	  STREAM^.BANKDEPTH := 0.0;
	  STREAM^.BANKVELOCITY := 0.0;
	  STREAM^.CURVSENSE := 1;
	  STREAM^.LOCALDEPTH := 0.0;
	  STREAM^.SEDIMENTFLUX := 0.0;
	  STREAM^.LOCALGRADIENT :=0.0;
	  STREAM^.LOCALCFRICTION := CFRICTION;
	  STREAM^.LOCALVELOCITY := 0.0;
	  STREAM^.BEDELEVATION := 0.0;{}
	  STREAM^.UPSTREAM_DISTANCE := 0.0;
	  STREAM^.MIGRATION_RATE := 0.0;
	  STREAM^.LOCALWIDTH := CHANNEL_WIDTH;
	  STREAM^.CURVATURE := 0.0;
	  STREAM^.CSCURVATURE := 0.0;
	  STREAM^.BEDCHANGE := 0.0;
          BASE := STREAM;
          READLN(INDATA);
        END;
    END;
    END;

  ENDBASE := NIL;
  WRITELN(' INITIAL RIVER SIZE=',RIVERSIZE: 5);
     YAVERAGE := YAVERAGE / RIVERSIZE;
    WRITELN (LIST);
    WRITELN (LIST,'YAVERAGE = ',YAVERAGE:10);
    WRITELN (LIST,'RIVERSIZE = ',RIVERSIZE);

    WRITELN(LIST,' ACROSS=',ACROSS:10);
    ASCROSS := (4.0*181.0) *(2.0*EXS*EXS+4.0*EXS/5.0+1.0/15.0)/
         (WDRATIO * WDRATIO *EXS1);
    WRITELN(LIST,' INITIAL ASCROSS=',ASCROSS:10);
   FOR I := 1 TO NDIM DO
   BEGIN
     FIRSTCOEFS[I] := 0.0;
     SECONDCOEFS[I] := 0.0;
   END;
   FIRSTCOEFS[2] := 1;
   SECONDCOEFS[3] := 2;
END;
{*********************************************************************************}
{*********************************************************************************}

PROCEDURE SETUP_ARRAY;

VAR  I: LONGINT;

BEGIN
{This sets up the upstream autoregression array of distances, C_velocities, and F_velocities}
   FOR I :=0 TO MAX_LAG DO
    BEGIN
        DISTANCE[I] := -1.0;
        PASTCVEL[I] := 0.0;
        PASTFVEL[I] := 0.0;
    END; (*FOR*)
END;(*SETUP_ARRAY*)
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE ERODE(YPOSITION, ANGLE : REAL;
                      VAR XCHANGE,YCHANGE,TOT_RATE : REAL);

VAR
   NINTYANGLE : REAL;

BEGIN
{This calculates the normalized erosion rate based on the
 curvature and calculated flow and depth perturbations}
   IF (SENSE*NEWSENSE = 1)
       THEN
     NINTYANGLE := ANGLE + HALFPI
       ELSE
     NINTYANGLE := ANGLE - HALFPI;
   IF BIASFACTOR
      THEN
    TOT_RATE := DISCH_FACTOR *  RATEMULT *
           NEWRATE * (1.0 - BIAS * COS(ANGLE)) / 2.0
      ELSE
    TOT_RATE := DISCH_FACTOR *  RATEMULT  * NEWRATE / 2.0;
{   IF BIASFACTOR
      THEN
    TOT_RATE := DISCH_FACTOR *  RATEMULT *
           NEWRATE * (1.0 - BIAS * COS(ANGLE)) / (2.0 * RESVALUE)
      ELSE
    TOT_RATE := DISCH_FACTOR *  RATEMULT * NEWRATE / (2.0 * RESVALUE);}
   XCHANGE := TOT_RATE * COS(NINTYANGLE);
   YCHANGE := TOT_RATE * SIN(NINTYANGLE);
   IF LIMITWIDTH THEN
      BEGIN
     IF (YPOSITION > YUPPER) AND (YCHANGE > 0.0)
         THEN BEGIN
           YCHANGE := 0.0;
           TOT_RATE := XCHANGE;
         END;
     IF (YPOSITION < YLOWER) AND (YCHANGE < 0.0)
         THEN BEGIN
           YCHANGE := 0.0;
           TOT_RATE := XCHANGE;
         END;
      END;
	  SUMABSRATE := SUMABSRATE+TOT_RATE;
      SUMNUM := SUMNUM+1.0;
    IF TOT_RATE > MAXTOT_RATE THEN MAXTOT_RATE := TOT_RATE;
END; (* ERODE *)
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE LUDECOMP;
VAR
   NLESS1,I,PIVPLACE,IPLUS1,J,L,KCOL, JCOL,JROW,TEMPPIVOT : LONGINT;
   SAVEIT,RATIO,THEVALUE : REAL;
BEGIN
{This solves for the discriminant of the upstream distance coefficients}
  DETERMINANT := 1.0;
  NLESS1 := NDIM-1;
  FOR I := 1 TO NDIM DO PIVOT[I] := I;
  FOR I := 1 TO NLESS1 DO
  BEGIN
    IPLUS1 := I+1;
    PIVPLACE := I;
    FOR J := IPLUS1 TO NDIM DO
       IF ABS(DISTCOEFS[PIVPLACE,I])<ABS(DISTCOEFS[J,I])
            THEN PIVPLACE := J;
    IF (PIVPLACE <> I) THEN
    BEGIN
      TEMPPIVOT := PIVOT[I];
      PIVOT[I] := PIVOT[PIVPLACE];
      FOR JCOL := 1 TO NDIM DO
      BEGIN
        SAVEIT := DISTCOEFS[I,JCOL];
        DISTCOEFS[I,JCOL] := DISTCOEFS[PIVPLACE,JCOL];
        DISTCOEFS[PIVPLACE,JCOL] := SAVEIT;
      END;
      PIVOT[PIVPLACE] := TEMPPIVOT;
      DETERMINANT := - DETERMINANT;
    END;
    FOR JROW := IPLUS1 TO NDIM DO
    BEGIN
      IF DISTCOEFS[JROW,I] <> 0.0 THEN
      BEGIN
          DISTCOEFS[JROW,I] :=
               DISTCOEFS[JROW,I] / DISTCOEFS[I,I];
          FOR KCOL := IPLUS1 TO NDIM DO
             DISTCOEFS[JROW,KCOL] := DISTCOEFS[JROW,KCOL] -
                 DISTCOEFS[JROW,I]*DISTCOEFS[I,KCOL];
      END;
    END;
  END;
  FOR I := 1 TO NDIM DO
    DETERMINANT := DETERMINANT* DISTCOEFS[I,I];
END;
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE MATRIXSOLVE;

VAR
    IROW,JCOL,I : LONGINT;
    XTEMP : VECTYPE;
    SUM : REAL;
BEGIN
{This solves the matrices determining upstream weights for calculating
 C_velocities and F_velocities setting up the weighting vectors D1TERM
 (1st derivative) and D2TERM (2nd derivative)}
   FOR I := 1 TO NDIM DO
     XTEMP[I] := DUMMYVEC[PIVOT[I]];
   FOR IROW := 2 TO NDIM DO
   BEGIN
     SUM := XTEMP[IROW];
     FOR JCOL := 1 TO IROW-1 DO
         SUM := SUM - DUMMYCOEFS[IROW,JCOL] * XTEMP[JCOL];
     XTEMP[IROW] := SUM;
   END;
   DUMMYVEC[NDIM] := XTEMP[NDIM]/DUMMYCOEFS[NDIM,NDIM];
   FOR IROW := NDIM - 1 DOWNTO 1 DO
   BEGIN
     SUM := XTEMP[IROW];
     FOR JCOL := IROW+1 TO NDIM DO
       SUM := SUM - DUMMYCOEFS[IROW,JCOL] * DUMMYVEC[JCOL];
     DUMMYVEC[IROW] := SUM/DUMMYCOEFS[IROW,IROW];
   END;
END;
{*********************************************************************************}
{*********************************************************************************}

PROCEDURE DERIVWEIGHTS;
VAR POWER,TERM : REAL;
    I,J : LONGINT;
BEGIN
{This sets up coefficients for evaluating the upstream weightings for
 calculating the C_velocities, and F_velocities, allowing for
 uneven distances between upstream nodes}
   IF DISTANCE[NDIM-2]>0.0 THEN
   BEGIN
     FOR J := 1 TO NDIM DO
     BEGIN
       DISTCOEFS[1,J] := 1.0;
       IF J = 1 THEN TERM := 0.0 ELSE TERM := TERM - DISTANCE[J-2];
       POWER := TERM;
       DISTCOEFS[2,J] := TERM;
       FOR I := 3 TO NDIM DO
       BEGIN
         POWER := POWER * TERM;
         DISTCOEFS[I,J] := POWER;
       END;
     END;
     LUDECOMP;
     FOR I := 1 TO NDIM DO
       FOR J := 1 TO NDIM DO
          DUMMYCOEFS[I,J] := DISTCOEFS[I,J];
     FOR I := 1 TO NDIM DO
       DUMMYVEC[I] := FIRSTCOEFS[I];
     MATRIXSOLVE;
     FOR I := 1 TO NDIM DO
       D1TERM[I-1] := DUMMYVEC[I];
     FOR I := 1 TO NDIM DO
       FOR J := 1 TO NDIM DO
          DUMMYCOEFS[I,J] := DISTCOEFS[I,J];
     FOR I := 1 TO NDIM DO
       DUMMYVEC[I] := SECONDCOEFS[I];
     MATRIXSOLVE;
     FOR I := 1 TO NDIM DO
       D2TERM[I-1] := DUMMYVEC[I];
   END;
END;
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE CALCCS;
VAR I : LONGINT;
    TEMP : REAL;
BEGIN
{This calculates the secondary circulation term}
  IF DISTANCE[NDIM-2] > 0.0 THEN
  BEGIN
    TEMP := 0.0;
    FOR I := 1 TO NDIM-1 DO TEMP := TEMP + PASTCS[I]*D1TERM[I];
    CSCURVATURE := (CURVATURE * CSEXPONENT - TEMP)
       / (D1TERM[0] +CSEXPONENT);
  END ELSE CSCURVATURE := CURVATURE;
  PASTCS[0] := CSCURVATURE;
  ABSCSCURVATURE := ABS(CSCURVATURE);
END;
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE FVELOCITY;
VAR
    TERM1,TERM2,TERM3,TERM4,TERM5 : REAL;
    I : LONGINT;
BEGIN
{This calculates the velocity and depth perturbations due to bedforms
 and sediment transport}
   IF DISTANCE[NDIM-2]>0.0 THEN
   BEGIN
     TERM1 := 0.0;
     TERM2 := 0.0;
     FOR I := 1 TO NDIM-1 DO
     BEGIN
        TERM1 := TERM1 + D2TERM[I] * PASTFVEL[I];
        TERM2 := TERM2 + D1TERM[I] * PASTFVEL[I];
     END;
     TERM3 := 0.0;
     TERM4 := 0.0;
     TERM5 := 0.0;
     FOR I := 0 TO NDIM-1 DO
     BEGIN
        TERM3 := TERM3 + D1TERM[I] * PASTCVEL[I];
        TERM4 := TERM4 + D1TERM[I] * PASTCURVATURE[I];
        TERM5 := TERM5 + D1TERM[I] * PASTCS[I];
     END;
     F_RATE := (-TERM1 - FCONST1*TERM2 + FCONST3*TERM3 + FCONST4*TERM4
              + FCONST5*TERM5)
              /(D2TERM[0]+FCONST1*D1TERM[0]+FCONST2);
     PASTFVEL[0] := F_RATE;
     TERM1 := 0.0;
     FOR I := 0 TO NDIM-1 DO
       TERM1 := TERM1 + PASTFVEL[I]*D1TERM[I];
     FDEPTH := {2.0*}(TERM1/FWEIGHT + 2.0 * F_RATE);
   END ELSE BEGIN
     F_RATE := 0.0;
     FDEPTH := 0.0;
     PASTFVEL[0] := 0.0;
   END;
END;
{*********************************************************************************}
{*********************************************************************************}

PROCEDURE RATESENSE;

VAR XSENSE : LONGINT;

BEGIN
{This determines whether the migration is to the left or right of the downstream direction}
IF NEWRATE = 0.0
    THEN NEWSENSE := 1
    ELSE IF NEWRATE > 0.0 THEN
        IF SENSE < 0 THEN NEWSENSE := -1
            ELSE NEWSENSE := 1
        ELSE IF SENSE < 0 THEN NEWSENSE := 1
            ELSE NEWSENSE := -1;

END; (*RATESENSENSE*)
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE PRINTDETAIL;


BEGIN
{This prints out a report of the major variables for a given stream location,
 including the X,Y coordinates, curvature, local gradient,
 velocity and depth perturbations}
    WRITELN(DETAIL,STREAM^.XVALUE:15,',',STREAM^.YVALUE:15,',',
        CURVATURE:15,',',C_RATE:15,',',
        F_RATE:15,',',CSCURVATURE:15,',',
		STREAM^.LOCALGRADIENT:15,',',STREAM^.LOCALDEPTH:15,',',STREAM^.SEDIMENTFLUX:15,',',STREAM^.BEDELEVATION:15,',',
        CDEPTH:15,',',FDEPTH:15,',',SENSE,',',NEWSENSE);
END;(*PRINTDETAIL*)
{*********************************************************************************}
{*********************************************************************************}

PROCEDURE LASTSTUFF;
VAR I : LONGINT;

BEGIN
{This migrates the stream and prints out details if it is time to do so}
      IF DETAILPAR AND (ITERNUM MOD PLOTINT=0) THEN
        BEGIN
           IF (IPLACE >= FIRSTWRITE) AND (IPLACE<=LASTWRITE) THEN
               PRINTDETAIL;
        END;
      IF FIRSTONE OR SECOND_ONE THEN BEGIN
          STREAM^.XVALUE := X2 + XC1 * TWOWGT;
          STREAM^.YVALUE := Y2 + YC1 * TWOWGT;
          STREAM^.ELAPSE:=STREAM^.ELAPSE+TOT_RATE*TWOWGT;
          IF NOT FIRSTONE THEN SECOND_ONE := FALSE;
          FIRSTONE := FALSE;
      END ELSE BEGIN
          STREAM^.XVALUE := X2 + 2.0*XC1;
          STREAM^.YVALUE := Y2 + 2.0*YC1;
          STREAM^.ELAPSE := STREAM^.ELAPSE + TOT_RATE + OLDTOT_RATE;
      END;
      UPLENGTH := DOWNLENGTH;
      X1 := X2;
      Y1 := Y2;
	  E1 := E2;{}
      X2 := X3;
      Y2 := Y3;
	  E2 := E3;{}
      X3 := X4;
      Y3 := Y4;
	  E3 := E4;{}
      IPLACE := IPLACE + 1;
      ANGLE1 := ANGLE2;
      IF PRINTTIME THEN STATISTICS;
END; (*LASTSTUFFUFF*)
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE CVELOCITY;
VAR I : LONGINT;
    TERM1, TERM2 : REAL;
BEGIN
{This calculates the velocity and depth perturbations due to planform curvature}
  IF DISTANCE[NDIM-2]>0.0 THEN
  BEGIN
    TERM1 := CURVATURE * D1TERM[0];
    TERM2 := 0.0;
    FOR I := 1 TO NDIM-1 DO
    BEGIN
      TERM2 := TERM2 + PASTCVEL[I]*D1TERM[I];
      TERM1 := TERM1 + PASTCURVATURE[I]*D1TERM[I];
    END;
    IF NEWPARKER THEN
       NEWRATE := (-TERM2+0.5*(-EXS20*TERM1+CWEIGHT*CURVATURE
                   +CSWEIGHT*CSCURVATURE))
                   /(D1TERM[0]+AEXPONENT)
    ELSE
       NEWRATE := (-TERM2+0.5*(TERM1+CWEIGHT*CURVATURE))
                  /(D1TERM[0]+AEXPONENT);
  END ELSE NEWRATE := 0.0;
END;
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE RATEAVERAGE;

VAR
       I  : LONGINT;
       DEPTHVAL,ABS_CFWEIGHT : REAL;
BEGIN
{This calculates the nominal migration rate based upon the calulated velocity
 and depth perturbations and other simulation parameters}
    IF NEWPARKER THEN
    BEGIN
      CALCCS;
	  STREAM^.CSCURVATURE := CSCURVATURE;
      CDEPTH := ACROSS * CSCURVATURE/2.0;
      CVELOCITY;
    END ELSE BEGIN
      CVELOCITY;
      CDEPTH := ACROSS * CURVATURE;
    END;
   IF USEFVELOCITY THEN
   BEGIN
      PASTCVEL[0] := NEWRATE;
      FVELOCITY;
      PASTFVEL[0] := F_RATE;
      C_RATE := NEWRATE;
      TDEPTH := CDEPTH+FDEPTH;
		STREAM^.BANKVELOCITY := C_RATE+F_RATE;
      IF TDEPTH<-0.9 THEN TDEPTH := -0.9;
      IF TDEPTH>0.9 THEN TDEPTH := 0.9;
      {DEPTHVAL := EXP(DEPTHWEIGHT*LN(TDEPTH));
      NEWRATE := (NEWRATE+F_RATE)*INTERVAL*0.01;
      IF NEWRATE > 0.0 THEN NEWRATE := NEWRATE*DEPTHVAL
          ELSE NEWRATE := NEWRATE/DEPTHVAL;}
      NEWRATE := (VELOCITYWEIGHT * (NEWRATE+F_RATE) +
                    DEPTHWEIGHT * (TDEPTH));
	  IF NEWRATE >= 0.0 THEN
	      IF NEWRATE> CRITICAL_USTAR THEN
		      NEWRATE :=NEWRATE-CRITICAL_USTAR
		  ELSE
		      NEWRATE := 0.0
	  ELSE
	      IF NEWRATE <= -CRITICAL_USTAR THEN
		      NEWRATE := NEWRATE+CRITICAL_USTAR
		  ELSE
		      NEWRATE := 0.0;

      {IF NEWRATE>MAXNEWRATE THEN
          NEWRATE := MAXNEWRATE*INTERVAL*0.01
      ELSE NEWRATE := NEWRATE*INTERVAL*0.01;}
	  NEWRATE := NEWRATE*NOM_TIME*NOM_RATE;
   END ELSE
   BEGIN
      PASTCVEL[0] := NEWRATE;
    {  FVELOCITY;
      PASTFVEL[0] := F_RATE; }
      C_RATE := NEWRATE;
	STREAM^.BANKVELOCITY := C_RATE;
      TDEPTH := CDEPTH;
      IF TDEPTH<-0.9 THEN TDEPTH := -0.9;
      IF TDEPTH>0.9 THEN TDEPTH := 0.9;
{      DEPTHVAL := EXP(DEPTHWEIGHT*LN(TDEPTH));
		STREAM^.BANKVELOCITY := NEWRATE;
    NEWRATE := NEWRATE*INTERVAL*0.01;
    IF NEWRATE > 0.0 THEN NEWRATE := NEWRATE*DEPTHVAL
          ELSE NEWRATE := NEWRATE/DEPTHVAL;}
      NEWRATE := (VELOCITYWEIGHT * NEWRATE +
                    DEPTHWEIGHT * (TDEPTH));
	   IF NEWRATE >= 0.0 THEN
	      IF NEWRATE> CRITICAL_USTAR THEN
		      NEWRATE :=NEWRATE-CRITICAL_USTAR
		  ELSE
		      NEWRATE := 0.0
	  ELSE
	      IF NEWRATE <= -CRITICAL_USTAR THEN
		      NEWRATE := NEWRATE+CRITICAL_USTAR
		  ELSE
		      NEWRATE := 0.0;


      {IF NEWRATE>MAXNEWRATE THEN
          NEWRATE := MAXNEWRATE*INTERVAL*0.01
      ELSE NEWRATE := NEWRATE*INTERVAL*0.01;}
	  NEWRATE := NEWRATE*NOM_TIME*NOM_RATE;
   END;
   IF NONLINEAR_BANK_EROSION THEN
      IF NEWRATE>0.0 THEN
	     NEWRATE := EXP(BANK_EROSION_EXPONENT*LN(NEWRATE))
	  ELSE IF NEWRATE<0.0 THEN
	     NEWRATE := -EXP(BANK_EROSION_EXPONENT*(LN(ABS(NEWRATE))))
	  ELSE NEWRATE :=0.0;
    IF ((SEDIMENT_ROUTE) AND (BEDLOAD_INDIRECT_VARIATION)) THEN
	  IF ((ITERNUM>10) {AND(STREAM^.SEDIMENTFLUX > NOMINAL_SEDFLUX)}) THEN
	   NEWRATE := NEWRATE*(1.0+(STREAM^.SEDIMENTFLUX/NOMINAL_SEDFLUX)*CSCURVATURE*BEDLOAD_WIDTH_CONSTANT);
	  ABS_CFWEIGHT := ABS(NEWRATE);
	  IF (ABS_CFWEIGHT>MAX_CFWEIGHT) THEN MAX_CFWEIGHT := ABS_CFWEIGHT;
	  SUM_CFWEIGHT := SUM_CFWEIGHT +ABS_CFWEIGHT;
	  NUM_CFWEIGHT := NUM_CFWEIGHT + 1.0;
	  SUMSQ_CFWEIGHT :=SUMSQ_CFWEIGHT+SQR(ABS_CFWEIGHT);
    STREAM^.MIGRATION_RATE := NEWRATE;
    MIGRATION := ABS(NEWRATE);
{    IF MIGRATION > 0.0 THEN
       NEWRATE := EXP(VELOCITYWEIGHT*LN(MIGRATION))/NEWRATE;}
    IF NEWRATE= 0.0 THEN STREAM^.BANKDEPTH := TDEPTH ELSE
      STREAM^.BANKDEPTH := TDEPTH*ABS(NEWRATE)/NEWRATE;
	 STREAM^.CURVSENSE := SENSE;
    RATESENSE;
   FOR I := (NDIM-1) DOWNTO 0 DO
      BEGIN
     DISTANCE [ I+1 ] := DISTANCE [ I ];
     PASTCURVATURE [ I+1 ] := PASTCURVATURE [I] ;
     PASTCS[I+1] := PASTCS[I];
     PASTCVEL[I+1] := PASTCVEL[I];
     PASTFVEL[I+1] := PASTFVEL[I];
      END; (*FOR*)
   NEWRATE := ABS(NEWRATE);

END; {RATEAVERAGE}
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE FINDMAXNEWRATE;
VAR CVELOCITY,CDEPTH : REAL;
BEGIN
{Finds the maximum possible values of velocity and depth perturbations and erosion rate}
    CVELOCITY := ((CWEIGHT+CSWEIGHT)*MAXRATECURVE)/(2.0*AEXPONENT);
    CDEPTH := ACROSS * MAXRATECURVE/2.0;
    MAXNEWRATE := VELOCITYWEIGHT * CVELOCITY +
                    DEPTHWEIGHT * CDEPTH;
    WRITELN(LIST,' MAXIMUM CVELOCITY=',CVELOCITY:15:7,' MAX DEPTH=',CDEPTH:15:7,
          ' MAXIMUM EROSION RATE=',MAXNEWRATE:15:7);
    {CLOSE(LIST);}
END;

{*********************************************************************************}
{*********************************************************************************}

PROCEDURE LOCAL_PARAMETERS(PARAMETER_INDEX : BOOLEAN);
VAR TEMPDEPTH : DOUBLE;
BEGIN
{This calculates local values of several flow parameters as well as
 unique values if channel width is specified to vary downstream}
IF SEDIMENT_ROUTE THEN
BEGIN
  IF PARAMETER_INDEX THEN
  BEGIN
     IF BEDLOAD_WIDTH_VARIATION THEN
	 BEGIN
	      STREAM^.LOCALWIDTH := CHANNEL_WIDTH*(1.0-(STREAM^.SEDIMENTFLUX/NOMINAL_SEDFLUX)*ABSCSCURVATURE*BEDLOAD_WIDTH_CONSTANT);
		  IF STREAM^.LOCALWIDTH < CHANNEL_WIDTH*0.6 THEN STREAM^.LOCALWIDTH := CHANNEL_WIDTH*0.6;
		  WDRATIO := DEFAULT_WDRATIO*EXP((5.0/3.0)*LN(STREAM^.LOCALWIDTH/CHANNEL_WIDTH));
          FROUDE := DEFAULT_FROUDE/EXP(0.33333*LN(STREAM^.LOCALWIDTH/CHANNEL_WIDTH));
	 END ELSE
	 BEGIN
         FLOWDEPTH :=STREAM^.LOCALDEPTH;
         FLOWVELOCITY := STREAM^.LOCALVELOCITY;
         IF FLOWDEPTH>0.0 THEN
            WDRATIO := CHANNEL_WIDTH/FLOWDEPTH
         ELSE
             WDRATIO :=20.0;
	     FROUDE := FLOWVELOCITY/SQRT(GRAVITY*FLOWDEPTH);
     END;
	 CFRICTION := STREAM^.LOCALCFRICTION;
  END;
 { CFRICTION := 1/(SQR(MANNING_STRICKLER_COEFFICIENT)*EXP(0.3333*LN(FLOWDEPTH/ROUGHNESS)));}
  IF (FLOWDEPTH<=0.0) OR( FLOWVELOCITY<=0.0) OR (WDRATIO<=0.0) OR (CFRICTION<=0.0) THEN
   BEGIN
      WRITELN('PROBLEM WITH FLOW PARAMETERS, BED ELEVATION =',STREAM^.BEDELEVATION:15);
      WRITELN('DEPTH =',STREAM^.LOCALDEPTH:15,' VELOCITY =',STREAM^.LOCALVELOCITY:15);
      WRITELN('CFRICTION =',STREAM^.LOCALCFRICTION:15);
    END;
END ELSE
BEGIN
  STREAM^.LOCALWIDTH := CHANNEL_WIDTH*(1.0+ CURVATURE_WIDTH_SLOPE*ABSCURVATURE);
  WDRATIO := DEFAULT_WDRATIO*EXP((5.0/3.0)*LN(STREAM^.LOCALWIDTH/CHANNEL_WIDTH));
  FROUDE := DEFAULT_FROUDE/EXP(0.33333*LN(STREAM^.LOCALWIDTH/CHANNEL_WIDTH));
 {WRITELN(' H=',FLOWDEPTH:12,'WD=',WDRATIO:12,
    ' U=',FLOWVELOCITY:12,'CF=',CFRICTION:12);}
END;
  STREAM^.FROUDE := FROUDE;
  ASCROSS := (4.0*181.0) *(2.0*EXS*EXS+4.0*EXS/5.0+1.0/15.0)/
         (WDRATIO * WDRATIO *EXS1);
  CSEXPONENT := DELTA * CFRICTION *WDRATIO;
  CSWEIGHT := WDRATIO*CFRICTION*(ACROSS+ASCROSS);
  CWEIGHT := WDRATIO*CFRICTION*(SQR(FROUDE)
            *EXS20-1.0);
  AEXPONENT :=  2.0 * CFRICTION * WDRATIO;
  FWEIGHT := CFRICTION*WDRATIO;
  GAMMA := BETACONST * 4.0/(CFRICTION*WDRATIO*WDRATIO);
  FCONST1 := FWEIGHT *(3.0-MCONST+NINTYDEGREE*NINTYDEGREE*GAMMA);
  FCONST2 := FWEIGHT*FWEIGHT*2.0*NINTYDEGREE*NINTYDEGREE*GAMMA;
  FCONST3 := FWEIGHT*(MCONST-1.0);
  FCONST4 := -FWEIGHT*FROUDE*FROUDE*EXS20/2.0;
  FCONST5 := -FWEIGHT*ACROSS/2.0;

END; (*LOCAL_PARAMETERS*)
{*********************************************************************************}
{*********************************************************************************}

PROCEDURE CYCLE1;
VAR COUNT : LONGINT;

BEGIN
{This is the main procedure that at each location calls routines to
 calculate the velocity and depth perturbations, to determine the rate
 and directional sense of node migration, and to actually migrate the channel.
 It then moves to the next point downstream until the end is reached}
    {WRITELN('CYCLE1');}
    COUNT :=0;
    WHILE (STREAM^.NEXT^.NEXT^.NEXT <> NIL) DO
    BEGIN
        STREAM := STREAM^.NEXT;
        COUNT :=COUNT+1;
        {WRITELN('COUNT=',COUNT:10);}
        FINDDISTANCE(DOWNLENGTH,X2,Y2,X3,Y3);
       {WRITELN('^^');}
	{	STREAM^.LOCALGRADIENT := (E1-E3)/(UPLENGTH+DOWNLENGTH);}
      { IF STREAM^.LOCALGRADIENT<0.0 THEN STREAM^.NEXT^.LOCALGRADIENT :=0.0;}
        CURVANGLE(X1,X2,X3,Y1,Y2,Y3,ANGLE);
        {WRITELN('##');}
        if (uplength + downlength)=0.0 then
        begin
          curvature:=0.0;
             writeln('c=',curvature:15:7);
        end
        else
        begin
           CURVATURE := 2.0*ANGLE * CHANNEL_WIDTH / (UPLENGTH + DOWNLENGTH);
        end;
      PASTCURVATURE[0] := CURVATURE;
      DISTANCE[0] := UPLENGTH / CHANNEL_WIDTH;
      PASTCS[0] := CURVATURE;
      IF (SENSE=0) THEN WRITELN('NONSENSE');
      ABSCURVATURE := CURVATURE / SENSE;
	  STREAM^.CURVATURE :=CURVATURE;
      IF (ABSCURVATURE > 0.0)
      THEN
        RADIUS := CHANNEL_WIDTH / ABSCURVATURE
      ELSE
        RADIUS := 1.0E+10;
		{if stream^.localgradient<=0.0 then writeln('#');}
	  {IF SEDIMENT_ROUTE THEN LOCAL_PARAMETERS(STREAM^.LOCALGRADIENT);}
         { WRITELN('**');}
	  IF ((SEDIMENT_ROUTE) OR (CURVATURE_WIDTH_VARIATION) OR (BEDLOAD_WIDTH_VARIATION)) THEN LOCAL_PARAMETERS(FIRST_USE);
          {WRITELN('%%');}
      DERIVWEIGHTS;
      RATEAVERAGE;
      IF (NEWRATE > MAXRATE) THEN
    MAXRATE := NEWRATE;
      FOLLOW := FOLLOW + UPLENGTH;
      ERODE(Y2,ANGLE1,XC1,YC1,TOT_RATE);
      FINDERODIBILITY(X2,Y2,XC1,YC1);
{      ERODE(Y2,ANGLE1,STREAM^.RESISTANCE,XC1,YC1,TOT_RATE);}
        X4 := STREAM^.NEXT^.NEXT^.XVALUE;
        Y4 := STREAM^.NEXT^.NEXT^.YVALUE;
		E4 := STREAM^.NEXT^.NEXT^.BEDELEVATION;{}
        FINDANGLE(ANGLE2,X2,Y2,X4,Y4);
        LASTSTUFF;
    END;
    STREAM := STREAM^.NEXT;
END; (*CYCLE1*)
{*********************************************************************************}
{*********************************************************************************}

PROCEDURE ITERATE;

BEGIN
{This sets up variables for the given iteration, calls CYCLE1 to calculate
 the migration for each node, calculates migration for the last node,
 and prints out summary details for the iteration}
{ IF DETAILPAR AND (ITERNUM MOD PLOTINT=0) THEN
    WRITELN(DETAIL,ITERNUM);}
  MAXTOT_RATE := 0.0;
  REPEAT
    DISCHARGE_FACTOR := LOGNORMAL(LOGMEDIAN,STDEV);
    DISCH_FACTOR := EXP(DVPOWER * LN(DISCHARGE_FACTOR));
  UNTIL DISCH_FACTOR > 0.1;
  IF DISCH_FACTOR > 3.0 THEN DISCH_FACTOR := 3.0;
  SUM_DISCHARGE := SUM_DISCHARGE + DISCHARGE_FACTOR;
  SUM_DISCH_FACTOR := SUM_DISCH_FACTOR +DISCH_FACTOR;
  NUM_DISCH_FACTOR := NUM_DISCH_FACTOR + 1.0;
  AEXPONENT := AEXPONENT / DISCH_FACTOR;
  CSEXPONENT := CSEXPONENT /DISCH_FACTOR;
  TOT_TIME := TOT_TIME + INTERVAL;
  MAXRATE := 0.0;
  FIRSTONE := TRUE;
  SECOND_ONE := TRUE;
  SETUP_ARRAY;
  FOLLOW := 0.0;
  STREAM := BASE;
  OLDTOT_RATE := 0.0;
  STREAM^.BANKDEPTH := 1.0;
  X1 := STREAM^.XVALUE;
  Y1 := STREAM^.YVALUE;
  E1 := STREAM^.BEDELEVATION;
  XS := X1;
  YS := Y1;
  X2 := STREAM^.NEXT^.XVALUE;
  Y2 := STREAM^.NEXT^.YVALUE;
  E2 := STREAM^.NEXT^.BEDELEVATION;
  X3 := STREAM^.NEXT^.NEXT^.XVALUE;
  Y3 := STREAM^.NEXT^.NEXT^.YVALUE;
  E3 := STREAM^.NEXT^.NEXT^.BEDELEVATION;
  FINDDISTANCE(UPLENGTH,X1,Y1,X2,Y2);
  FINDANGLE(ANGLE1,X1,Y1,X3,Y3);
  {STREAM^.LOCALGRADIENT :=(E1-E2)/UPLENGTH;}
{ IF STREAM^.LOCALGRADIENT<0.0 THEN STREAM^.LOCALGRADIENT := 0.0;}
  IPLACE := 2;
  CYCLE1; {CYXLW1 progresses doenstream and calculates flow, erosion rate, and migrates all stream nodes}
  FINDDISTANCE(DOWNLENGTH,X2,Y2,X3,Y3);
 { STREAM^.NEXT^.LOCALGRADIENT := (E2-E3)/DOWNLENGTH;}
  {IF STREAM^.NEXT^.LOCALGRADIENT<0.0 THEN STREAM^.NEXT^.LOCALGRADIENT :=0.0;}


  CURVANGLE(X1,X2,X3,Y1,Y2,Y3,ANGLE);
      CURVATURE := 2.0*ANGLE * CHANNEL_WIDTH / (UPLENGTH + DOWNLENGTH);
      PASTCURVATURE[0] := CURVATURE;
      DISTANCE[0] := UPLENGTH / CHANNEL_WIDTH;
      PASTCS[0] := CURVATURE;
      ABSCURVATURE := CURVATURE / SENSE;
	  STREAM^.NEXT^.CURVATURE :=CURVATURE;
      IF (ABSCURVATURE > 0.0)
    THEN
    RADIUS := CHANNEL_WIDTH / ABSCURVATURE
      ELSE
    RADIUS := 1.0E+10;
	{if stream^.localgradient<=0.0 then writeln('@');}
	IF ((SEDIMENT_ROUTE) OR (CURVATURE_WIDTH_VARIATION) OR(BEDLOAD_WIDTH_VARIATION)) THEN LOCAL_PARAMETERS(FIRST_USE);
	FIRST_USE := TRUE;
      RATEAVERAGE;
      IF (NEWRATE > MAXRATE) THEN
    MAXRATE := NEWRATE;
      FOLLOW := FOLLOW + UPLENGTH;
  FINDANGLE(ANGLE2,X2,Y2,X3,Y3);
  ERODE(Y3,ANGLE2,XC2,YC2,TOT_RATE);
  FINDERODIBILITY(X3,Y3,XC2,YC2);
  STREAM^.NEXT^.XVALUE := X3 + XC2* 2.0;
  STREAM^.NEXT^.YVALUE := Y3 + YC2* 2.0;
  STREAM^.NEXT^.ELAPSE := STREAM^.NEXT^.ELAPSE + TOT_RATE * 2.0;
  LASTSTUFF;
  FOLLOW := FOLLOW + UPLENGTH;
  FINDDISTANCE(STRAIGHT,X3,Y3,XS,YS);
{    WRITELN('I=',ITERNUM,' QF=',
       DISCH_FACTOR:10,' MT=',MAXTOT_RATE:10,' IV=',INTERVAL:10);}
    {WRITELN('I=',ITERNUM);}
  IF USERATEMULT THEN
      RATEMULT := EXP(SLOPE_FACTOR * LN(STRAIGHT / FOLLOW))
  ELSE RATEMULT := 1.0;
  IF ((NOT SEDIMENT_ROUTE) AND (NOT CURVATURE_WIDTH_VARIATION) AND (NOT BEDLOAD_WIDTH_VARIATION)) THEN
  BEGIN
  AEXPONENT := NOMAEXPONENT * RATEMULT;
  CSEXPONENT := NOMCSEXPONENT * RATEMULT;
  CSWEIGHT := WDRATIO*RATEMULT*CFRICTION*
            (ACROSS+ASCROSS/SQR(RATEMULT));
  CWEIGHT := WDRATIO*RATEMULT*CFRICTION*(SQR(FROUDE)
            *RATEMULT*RATEMULT*RATEMULT
            *EXS20-1.0);
  FWEIGHT := FNOMWGHT*RATEMULT;
  GAMMA := BETACONST * 4.0/(CFRICTION*WDRATIO
          *WDRATIO*RATEMULT*RATEMULT);
  FCONST1 := FWEIGHT *(3.0-MCONST+NINTYDEGREE*NINTYDEGREE*GAMMA);
  FCONST2 := FWEIGHT*FWEIGHT*2.0*NINTYDEGREE*NINTYDEGREE*GAMMA;
  FCONST3 := FWEIGHT*(MCONST-1.0);
  FCONST4 := -FWEIGHT*FROUDE*FROUDE*EXS20
           *RATEMULT*RATEMULT*RATEMULT/2.0;
  FCONST5 := -FWEIGHT*ACROSS/2.0;
  END;
  {INTERVAL := NOM_TIME * INTERVAL * CHANNEL_WIDTH * NOM_RATE /
          (MAXRATE * RATEMULT);
  IF INTERVAL > MAXINTERVAL * CHANNEL_WIDTH THEN
             INTERVAL := MAXINTERVAL * CHANNEL_WIDTH;}
  INTERVAL := NOM_TIME;

  IF ITERNUM MOD PRINTINT=0 THEN
    BEGIN
      WRITELN(LIST);
      WRITELN(LIST,'ITERATION = ',ITERNUM,
             '   RIVER SIZE = ',RIVERSIZE);
      AVG_LENGTH := FOLLOW/(RIVERSIZE - 1);
      WRITELN(LIST,'CURRENT TIME=',TOT_TIME:10);
      WRITELN(LIST,'STRAIGHT-LINE, CURVE-LENGTH, AND SINUOSITY ');
      WRITELN(LIST,STRAIGHT:10,'  ',FOLLOW:10,
             '  ',(FOLLOW / STRAIGHT):10);
      WRITELN(LIST,'AVERAGE STREAM LENGTH=',AVG_LENGTH:10);
    END;
{ IF DETAILPAR AND (ITERNUM MOD PLOTINT=0) THEN
    WRITELN(DETAIL);}
  IF ITERNUM MOD UPDATEINT=0 THEN
  BEGIN
      WRITELN(OUTTIME,ITERNUM:6,' ',TOT_TIME:13,' ',(FOLLOW/STRAIGHT):13,
        ' ',(SUMABSRATE/SUMNUM):13);
  END;
END;(*ITERATE*)
{*********************************************************************************}
{*********************************************************************************}
PROCEDURE WRITELONGFILE;

VAR I : LONGINT;
AVG_BANKPER,AVG_VELPER,AVG_GRAD,AVG_DEPTH,AVG_FLUX,AVG_VELOCITY,AVG_FRICTION,AVG_MIGRATE,AVG_CURVATURE,AVG_FROUDE, AVG_CHANGE,AVG_WIDTH, AVG_CSCURV : DOUBLE;

BEGIN
{This cycles downstream and outputs the values of local variables for each node,
 then outputs several average values of variables}
  AVG_BANKPER := 0.0;
  AVG_VELPER := 0.0;
  AVG_GRAD :=0.0;
  AVG_DEPTH := 0.0;
  AVG_FLUX := 0.0;
  AVG_VELOCITY := 0.0;
  AVG_FRICTION := 0.0;
  AVG_MIGRATE := 0.0;
  AVG_CURVATURE := 0.0;
  AVG_FROUDE := 0.0;
  AVG_CHANGE :=0.0;
  AVG_WIDTH := 0.0;
  AVG_CSCURV :=0.0;
  STREAM := BASE;
  FOR I := 1 TO RIVERSIZE DO
    BEGIN
      WRITELN(OUTFILE,STREAM^.XVALUE:14,',',STREAM^.YVALUE:14,',',STREAM^.UPSTREAM_DISTANCE:14
              ,',',STREAM^.ELAPSE:14,',',STREAM^.RESISTANCE:14
			  ,',',STREAM^.BANKDEPTH:14,',',STREAM^.BANKVELOCITY:14,',',STREAM^.CURVSENSE
			  ,',',STREAM^.LOCALGRADIENT:15,',',STREAM^.LOCALDEPTH:15,',',STREAM^.SEDIMENTFLUX:15,',',STREAM^.BEDELEVATION:15
			  ,',',STREAM^.LOCALVELOCITY:15,',',STREAM^.LOCALCFRICTION:15,',',STREAM^.MIGRATION_RATE:15
			  ,',',STREAM^.CURVATURE:15, ',',STREAM^.FROUDE:15,',',STREAM^.BEDCHANGE:15,',',STREAM^.LOCALWIDTH:15,',',STREAM^.CSCURVATURE:15);
	   AVG_BANKPER := AVG_BANKPER+ABS(STREAM^.BANKDEPTH);
	   AVG_VELPER := AVG_VELPER+ABS(STREAM^.BANKVELOCITY);
	   AVG_GRAD := AVG_GRAD+ STREAM^.LOCALGRADIENT;
	   AVG_DEPTH := AVG_DEPTH+STREAM^.LOCALDEPTH;
	   AVG_FLUX := AVG_FLUX +STREAM^.SEDIMENTFLUX;
	   AVG_VELOCITY := AVG_VELOCITY+STREAM^.LOCALVELOCITY;
	   AVG_FRICTION := AVG_FRICTION+STREAM^.LOCALCFRICTION;
	   AVG_MIGRATE := AVG_MIGRATE+ABS(STREAM^.MIGRATION_RATE);
	   AVG_CURVATURE := AVG_CURVATURE+ABS(STREAM^.CURVATURE);
	   AVG_FROUDE := AVG_FROUDE+STREAM^.FROUDE;
	   AVG_CHANGE := AVG_CHANGE+STREAM^.BEDCHANGE;
	   AVG_WIDTH := AVG_WIDTH + STREAM^.LOCALWIDTH;
	   AVG_CSCURV := AVG_CSCURV+ABS(STREAM^.CSCURVATURE);

      STREAM := STREAM^.NEXT;
    END;
	WRITELN(AVG_VALS, AVG_BANKPER/RIVERSIZE:15,',',AVG_VELPER/RIVERSIZE:15,',',AVG_GRAD/RIVERSIZE:15,',',
	AVG_DEPTH/RIVERSIZE:15,',',AVG_FLUX/RIVERSIZE:15,',',AVG_VELOCITY/RIVERSIZE:15,',',AVG_FRICTION/RIVERSIZE:15,',',AVG_MIGRATE/RIVERSIZE:15,',',AVG_CURVATURE/RIVERSIZE:15,',',AVG_FROUDE/RIVERSIZE:15
	,',',AVG_CHANGE/RIVERSIZE:15,',',DISCHARGE:15,',',CHANNEL_WIDTH:15,',',SHEAR_EXCEED_FACTOR:15,',',ROUGHNESS:15,',',SEDIMENT_TRANSPORT_EXPONENT:15,',',
	TRANSPORT_FACTOR:15,',',FLOW_FRACTION:15,',',TRANSPORT_CRITICAL_DIM_SHEAR:15,',',UPWIND:15,',',BEDMATERIAL_RELATIVE_DEPTH_FRACTION:15
	,',',SEDIMENT_SIZE:15,',',UPSTREAM_SEDIMENT_FLUX:15,',',AVG_WIDTH/RIVERSIZE:15,',',AVG_CSCURV/RIVERSIZE:15);

{writeln('done writelongfile');}
END;(*WRITELONGFILE*)
{*********************************************************************************}
{*********************************************************************************}

PROCEDURE WRITESHORTFILE;

VAR I,ICOUNT : LONGINT;

BEGIN
{This just writes out the X,Y values of all stream nodes}
 {WRITELN(OUTDATA, RIVERSIZE:10);}
  STREAM := BASE;
  FOR I := 1 TO RIVERSIZE DO
    BEGIN
      WRITELN(OUTDATA,STREAM^.XVALUE:14,' ',STREAM^.YVALUE:14);
      STREAM := STREAM^.NEXT;
    END;
END;(*WRITESHORTFILE*)

{*********************************************************************************}
{*********************************************************************************}

BEGIN (*MAIN PROGRAM*)
{This opens the major input and output files, reads simulation parameters
 and input data, calculates several additional fixed and initial parameters,
 sets up the floodplain if used, sets up the stream profile if sediment transport
 is used and then iterates over the specified number of time increments
 (calling ITERATE), outputs data at specific time increments, and finalizes
 the simulation and closes files.}
  ASSIGN(OUTTIME,'timestat.dat');
  REWRITE (OUTTIME);
  ASSIGN(INFILE,'meander_parameters.prm');
  RESET (INFILE);
  ASSIGN(INDATA,'indata.dat');
  RESET (INDATA);
   ASSIGN(PARAMS,'initial.prm');
  RESET(PARAMS);
  ASSIGN(INSEDIMENT,'sediment_parameters.prm');
  RESET(INSEDIMENT);
  ASSIGN(LIST,'list.dat');
  REWRITE (LIST);
  ASSIGN(DETAIL,'detail.dat');
  REWRITE (DETAIL);
  ASSIGN(OUTCAPT1,'neckcapt.dat');
  REWRITE (OUTCAPT1);
  ASSIGN(OUTCAPT2,'chutecapt.dat');
  REWRITE (OUTCAPT2);
  ASSIGN(OUTBANKTYPE,'banktype.dat');
  REWRITE (OUTBANKTYPE);
  ASSIGN(OUTDATA,'outdata.dat');
  REWRITE(OUTDATA);
  {ASSIGN(DISTLIST,'distance.dat');
  REWRITE(DISTLIST);}
  ASSIGN(OUTFILE,'plotdata.dat');
  REWRITE(OUTFILE);
  ASSIGN(AVG_VALS,'average_values.dat');
  REWRITE(AVG_VALS);
  ASSIGN(JP_PARAMETERS,'jp_parameters.dat');
  RESET(JP_PARAMETERS);
  ASSIGN(SEDIMENT_FLAG,'sediment_flag.prm');
  RESET(SEDIMENT_FLAG);
  GETVALUES;
  writeln('getvalues');
  ISEED1 := SEED/4.6566128752458E-10;
  DUMMY := 1;
  WRITELN(LIST,'STANDARD ERROR OF RESISTANCE=',ERROR:12);
  WRITELN(LIST,'STANDARD ERROR OF DISCHARGE=',LOGSTDEV:12);
  WRITELN(LIST,'SIMILARITY PARAMETER FOR RESISTANCE=',SIMILPAR:12);
  WRITELN(LIST,'RESISTANCE AVERAGING PARAMETER (0 FOR NONE)='
         ,RES_TYPE);
  LOGMEAN := 1.0;
  STDEV := SQRT(LN(SQR(LOGSTDEV/LOGMEAN)+1.0));
  LOGMEDIAN := LOGMEAN/EXP(SQR(STDEV)/2.0);
  WRITELN(LIST,'MEDIAN DISCHARGE RELATIVE TO MEAN=',LOGMEDIAN:12);
  WRITELN(LIST,'DISCHARGE EQUIVALENT NORMAL ST. ERR.=',STDEV:12);
  WRITELN(LIST,'RANDOM NUMBER SEED=',SEED:12);
{  SETUPFILE;}
  READ1PARAMETERS;
  READ2PARAMETERS;
  READ3PARAMETERS;
  CLOSE(JP_PARAMETERS);
  EXS1 := ACONST / SQRT(CFRICTION);
    EXS := EXS1 - 1.0/3.0;
    EXS20 := (EXS*EXS*EXS+EXS*EXS+2.0*EXS/5.0+2.0/35.0)/
  (EXS1*EXS1*EXS1);
    DELTA := ((EXS1*EXS1)*(EXS+0.25))/
        (EXS*EXS/12.0+EXS*11.0/360.0+1.0/504.0);
    WRITELN(LIST,'FIXED CONSTANTS:');
    WRITELN(LIST,'EXS1=',EXS1:10,' EXS=',EXS:10,' EXS20=',EXS20:10,
         ' DELTA=',DELTA:10);
    ACROSS := 2.0*(EXS+2.0/7.0)/(45.0*ACONST*ACONST*
            BETACONST*(EXS+1.0/3.0));

  writeln('readparameters');
  INITIALIZE;
  writeln('initialize');
   NOMAEXPONENT :=  2.0 * CFRICTION * WDRATIO;
  AEXPONENT := NOMAEXPONENT;
  NOMCSEXPONENT := DELTA * CFRICTION *WDRATIO;
  CSEXPONENT := NOMCSEXPONENT;
  CSWEIGHT := WDRATIO*CFRICTION*(ACROSS+ASCROSS);
  CWEIGHT := WDRATIO*CFRICTION*(SQR(FROUDE)
            *EXS20-1.0);
  FNOMWGHT := CFRICTION*WDRATIO;
  FWEIGHT := FNOMWGHT;
  GAMMA := BETACONST * 4.0/(CFRICTION*WDRATIO*WDRATIO);
  FCONST1 := FWEIGHT *(3.0-MCONST+NINTYDEGREE*NINTYDEGREE*GAMMA);
  FCONST2 := FWEIGHT*FWEIGHT*2.0*NINTYDEGREE*NINTYDEGREE*GAMMA;
  FCONST3 := FWEIGHT*(MCONST-1.0);
  FCONST4 := -FWEIGHT*FROUDE*FROUDE*EXS20/2.0;
  FCONST5 := -FWEIGHT*ACROSS/2.0;
  IF NOT SEDIMENT_ROUTE THEN
  BEGIN
  WRITELN(LIST);
  WRITELN(LIST,'SCALED PARAMETER VALUES');
  WRITELN(LIST,'AEXPONENT= ',AEXPONENT:10,
           ' CSEXPONENT=',CSEXPONENT:10);
  WRITELN(LIST,'CWEIGHT=',CWEIGHT:10,' CSWEIGHT=',CSWEIGHT:10);
  WRITELN(LIST,'FWEIGHT=',FWEIGHT:10,' FCONSTANTS:');
  WRITELN(LIST,FCONST1:10,' ',FCONST2:10,' ',FCONST3:10,' ',
           FCONST4:10,' ',FCONST5:10);
  END; {not_sedimentroute}
  FINDMAXNEWRATE;
  {WRITELN('FINDMAXNEWRATE');}
  INTERVAL := NOM_TIME;
  RATEMULT := 1.0;
  TOT_TIME := START_TIME;
  ESCAPE := CHR(27);
  FORMFEED := CHR(12);
  IF DETAILPRINT = 1 THEN DETAILPAR := TRUE ELSE DETAILPAR := FALSE;
  LIMITWIDTH := TRUE;
  IF VALLEY_WIDTH < 0.0
    THEN
       LIMITWIDTH := FALSE
    ELSE
      BEGIN
    YUPPER := YAVERAGE + VALLEY_WIDTH;
    YLOWER := YAVERAGE - VALLEY_WIDTH;
      END; {valley_width}
  IF BIAS > 0.0 THEN BIASFACTOR := TRUE ELSE BIASFACTOR := FALSE;
  WRITELN(LIST,FORMFEED);
  IF USEDEPOSIT THEN SETUPFLOODPLAIN;
  ADJUST;
    IF SEDIMENT_ROUTE THEN
  BEGIN
    READ_ROUTE_PARAMETERS;
	ASSIGN(JP_PARAMETERS,'jp_parameters.dat');
	REWRITE(JP_PARAMETERS);
	WRITELN(JP_PARAMETERS,CHANNEL_WIDTH:15,'  - stream width');
	WRITELN(JP_PARAMETERS,CFRICTION:15,' - coefficient of friction');
	WRITELN(JP_PARAMETERS,WDRATIO:15,'   - width-depth ratio');
	WRITELN(JP_PARAMETERS,FROUDE:15,'   - froude number');
	CLOSE(JP_PARAMETERS);
	IF DO_INITIAL_GRADIENT THEN SET_INITIAL_ELEVATIONS;
	WDRATIO := CHANNEL_WIDTH/DOWNSTREAM_DEPTH;
	FROUDE := QPERW/(SQRT(GRAVITY)*EXP(1.5*LN(DOWNSTREAM_DEPTH)));
        {WRITELN('SEDIMENT INITIALIZED, WDRATIO=',WDRATIO:10,' FROUDE=',FROUDE:10);}
  END; {sediment_route}
  PRINTTIME := FALSE;
  IF INPUT_PRINT >0 THEN
  BEGIN
        WRITELN(OUTDATA,TOT_TIME,' ',RIVERSIZE);
        WRITESHORTFILE;
  END; {input_print}
  {WRITELN('TOSEDIMENTROUTE');}
  IF SEDIMENT_ROUTE THEN
  BEGIN
     SET_INITIAL_ELEVATIONS;
	 AVERAGE_DEPTH := DOWNSTREAM_DEPTH;
	 FIRST_USE := FALSE;
	 WRITELN('AVERAGE CHANNEL DEPTH =',DOWNSTREAM_DEPTH:15);
	 FLOWDEPTH := DOWNSTREAM_DEPTH;
	 FLOWVELOCITY := QPERW/DOWNSTREAM_DEPTH;
	 WDRATIO := CHANNEL_WIDTH/DOWNSTREAM_DEPTH;
	 LOCAL_PARAMETERS(FIRST_USE);
	 WRITELN(LIST);

     WRITELN('FLOW DEPTH =',FLOWDEPTH:15);
     WRITELN(LIST,'FLOW VELOCITY =',FLOWVELOCITY:15);
     WRITELN(LIST,'WIDTH DEPTH RATIO =',WDRATIO:15);
     WRITELN(LIST,'FROUDE NO =',FROUDE:15);
     WRITELN(LIST,'ASCROSS =',ASCROSS:15);
  WRITELN(LIST,'SCALED PARAMETER VALUES');
  WRITELN(LIST,'AEXPONENT= ',AEXPONENT:10,
           ' CSEXPONENT=',CSEXPONENT:10);
  WRITELN(LIST,'CWEIGHT=',CWEIGHT:10,' CSWEIGHT=',CSWEIGHT:10);
  WRITELN(LIST,'FWEIGHT=',FWEIGHT:10,' FCONSTANTS:');
  WRITELN(LIST,FCONST1:10,' ',FCONST2:10,' ',FCONST3:10,' ',
           FCONST4:10,' ',FCONST5:10);
  END; {sediment_route}
 writeln('begin iterate');
   SUMABSRATE := 0.0;
   SUMNUM := 0.0;
   SUM_FLUX := 0.0;
   NUMBER_SUM_FLUX := 0.0;
   SUM_DISCH_FACTOR := 0.0;
   NUM_DISCH_FACTOR := 0.0;
   SUM_DISCHARGE := 0.0;
   MAX_CFWEIGHT := 0.0;
   SUM_CFWEIGHT := 0.0;
   NUM_CFWEIGHT := 0.0;
   SUMSQ_CFWEIGHT := 0.0;
   FAULT_STATE := FALSE;
   END_SIMULATION := FALSE;
  FOR ITERNUM := 1 TO ITERATION DO
    BEGIN
	IF NOT FAULT_STATE THEN
	BEGIN
	WRITELN('ITERATION =',ITERNUM:10,' TOTAL TIME=',TOT_TIME:15);
      IF ITERNUM MOD PRINTINT = 0 THEN
        BEGIN
          PRINTTIME := TRUE;
          STATINITIALIZE;
        END; {iternum mod printint}

      ITERATE;
      {WRITELN('##');}
      if SEDIMENT_ROUTE then
      BEGIN
         { WRITELN('INTERVAL =',INTERVAL:15,' TRANSPORT_INCREMENT =',TRANSPORT_INCREMENT:15);}
          TRANSPORT_RATIO := INTERVAL/TRANSPORT_INCREMENT;
          IF TRANSPORT_RATIO >2000.0 THEN
              TRANSPORT_ITERATIONS := 2000 ELSE
	  TRANSPORT_ITERATIONS := ROUND(INTERVAL/TRANSPORT_INCREMENT);
	  ITERATION_TIME := SECONDS_PER_YEAR*TRANSPORT_INCREMENT*FLOW_FRACTION/TRANSPORT_ITERATIONS;
	  {WRITELN('TRANSPORT_ITERATIONS=',TRANSPORT_ITERATIONS:15,' ITERATION TIME = ', ITERATION_TIME:15);}
        IF TRANSPORT_ITERATIONS<1 THEN TRANSPORT_ITERATIONS :=1;
      { FOR III := 1 TO 100 DO }{WRITELN('TRANSPORT_ITERATIONS =',TRANSPORT_ITERATIONS:15);}
      END; {sediment_route}
       {WRITELN('$$');}
      ADJUST;
      {IF (ITERNUM = 2) THEN
      BEGIN
      ASSIGN(OUTFILE,'tempout.out');
      REWRITE(OUTFILE);
      WRITELN(OUTFILE,ITERNUM,'  ITERATION NUMBER');
      WRITELONGFILE;
      CLOSE(OUTFILE);
	  END;}
       {WRITELN('@@');}
	  NETBEDCHANGE :=0.0;
	  NETBEDNUMBER :=0.0;
	  ABSBEDCHANGE := 0.0;
	  IF SEDIMENT_ROUTE THEN
	  BEGIN
	     EVOLVE_STREAM_PROFILE;
		 IF FAULT_STATE THEN BREAK;
	     NETBEDCHANGE := NETBEDCHANGE/NETBEDNUMBER;
	     ABSBEDCHANGE := ABSBEDCHANGE/NETBEDNUMBER;
	     {WRITELN('I =',ITERNUM:10,'  NET BED CHANGE =',NETBEDCHANGE:15,'  ABSOLUTE BED CHANGE=',ABSBEDCHANGE:15,'  NET BED COUNT =',NETBEDNUMBER:15);}
	  END; {sediment route}
	  IF ITERNUM = 1 THEN
	  BEGIN
	  ASSIGN(TRANSFER,'bedprofile.dat');
	  REWRITE(TRANSFER);
	  STREAM := BASE;
	  WHILE STREAM^.NEXT <> NIL DO
	  BEGIN
	    WRITELN(TRANSFER, STREAM^.XVALUE:15,',',STREAM^.YVALUE:15,',',STREAM^.BEDELEVATION:15);
          STREAM :=STREAM^.NEXT;
	  END; {steeam_nest}
	  CLOSE(TRANSFER);
	  END; {iternum}

       IF (ITERNUM MOD UPDATEINT = 0)
    THEN
    BEGIN
        ADJUST;
      IF USEDEPOSIT THEN
       BEGIN
         IF (ITERNUM MOD PLOTINT = 0) AND (ITERNUM > MINITERATION) THEN
    IRANGE := IDISTCATAGORIES DIV 2 - 1
   ELSE
    IRANGE := TRUNC(MAXWIDTHCHECK)+1;
         AGE_ELEVATION;
       END; {usedeposit}
    END; {uternum mod updateint}
    END {NOT FAULT STATE } ELSE
	   WRITELN(LIST,'FAULT STATE, CLOSING DOWN');

    IF ((ITERNUM MOD PLOTINT = 0) AND (ITERNUM > MINITERATION)) OR FAULT_STATE THEN
      BEGIN
       { FILENAME;}
	    WRITELN(LIST,'ITERATION =',ITERNUM:10,' TOTAL TIME=',TOT_TIME:15,' RIVER SIZE =',RIVERSIZE:10);
        WRITELN(OUTDATA,TOT_TIME,' ',RIVERSIZE);
        WRITESHORTFILE;
	    WRITEBANKDATA;
		WRITELN(OUTFILE,TOT_TIME,' ',RIVERSIZE);
		WRITELONGFILE;
		IF NOT FAULT_STATE THEN
		BEGIN
		IF (SUMABSRATE*SUMNUM) >0.0 THEN
		BEGIN
                if sediment_route then
                begin
		DIMENSIONLESS_FLUX_RATIO := (SUM_FLUX*NUMBER_SUM_FLUX*CHANNEL_WIDTH)/(BEDMATERIAL_RELATIVE_DEPTH_FRACTION*SUMABSRATE*SUMNUM);
		WRITELN(LIST,'RATIO OF FLUXES OF DOWNSTREAM SEDIMENT TRANSPORT TO MIGRATION BED REWORKING =',DIMENSIONLESS_FLUX_RATIO:15);
		SUM_FLUX :=0.0;
		NUMBER_SUM_FLUX :=0.0;
		SUMABSRATE :=0.0;
		SUMNUM :=0.0;
		AVG_DISCH_FACTOR := SUM_DISCH_FACTOR/NUM_DISCH_FACTOR;
		SUM_DISCHARGE := SUM_DISCHARGE /NUM_DISCH_FACTOR;
		WRITELN(LIST,'AVERAGE DISCHARGE = ',SUM_DISCHARGE:15,' AVERAGE DISCHARGE FACTOR=',AVG_DISCH_FACTOR:15);
		{WRITELN('SSQ =',SUMSQ_CFWEIGHT:13,' SUM =',SUM_CFWEIGHT:15,' NIM = ',NUM_CFWEIGHT:15);}
		IF NUM_CFWEIGHT>0.0 THEN
		SUMSQ_CFWEIGHT := SQRT((SUMSQ_CFWEIGHT-SQR(SUM_CFWEIGHT)/NUM_CFWEIGHT)/NUM_CFWEIGHT)
		ELSE
		SUMSQ_CFWEIGHT := 0.0;
		SUM_CFWEIGHT := SUM_CFWEIGHT/NUM_CFWEIGHT;
		WRITELN(LIST,'AVERAGE ABSOLUTE VELOCITY DEVIATION = ',SUM_CFWEIGHT:15,' MAXIMUM DEVIATION =',MAX_CFWEIGHT:15);
		WRITELN(LIST,'STANDARD DEVIATION OF ABSOLUTE VELOCITY DEVIATIONS =',SUMSQ_CFWEIGHT:15);
		SUM_DISCH_FACTOR := 0.0;
		NUM_DISCH_FACTOR := 0.0;
		SUM_DISCHARGE := 0.0;
		MAX_CFWEIGHT := 0.0;
        SUM_CFWEIGHT := 0.0;
        NUM_CFWEIGHT := 0.0;
        SUMSQ_CFWEIGHT := 0.0;
		END; {sediment_route}
        end;{sumabsrate+sumnum}
		END;{not fault_state}
        {INCREMENT;}
      END; {iternum mod plotint}

      IF PRINTTIME OR FAULT_STATE
    THEN
    BEGIN
		{ASSIGN(OUTFILE,'tempout.out');
      REWRITE(OUTFILE);
      WRITELN(OUTFILE,ITERNUM,'  ITERATION NUMBER');
      WRITELONGFILE;
      CLOSE(OUTFILE);}
      WRITELN(LIST);
      WRITELN(LIST,'UPDATE AT ITERATION ',ITERNUM);
      WRITELN(LIST,' AEXPONENT=',AEXPONENT*DISCH_FACTOR:10,
            ' CSEXPONENT=',CSEXPONENT*DISCH_FACTOR:10,
            ' CSWEIGHT=',CSWEIGHT:10,
            ' CWEIGHT=',CWEIGHT:10);
      WRITELN(LIST,'RATEMULT=',RATEMULT:10,
         ' INTERVAL=',INTERVAL:10);
      STATPRINT;
      WRITELN(LIST);
      PRINTTIME := FALSE;
	WRITE_DEBUG;
	WRITE_VECTOR;
         END; {printtime pr fault_state}
      IF (ITERNUM MOD LOOPCHKINT = 0) AND (ITERNUM <> ITERATION) AND (NOT FAULT_STATE)
    THEN
    BEGIN
      CHECKLOOP;
      ADJUST;
      if usedeposit then FINDPLUGS;
    END; {iternum mod loopcheckint}
	IF FAULT_STATE THEN BREAK;
    END;(*ITERATIONS*)
        IF USEDEPOSIT THEN WRITEAGEELEV;
   CLOSE(LIST);
   CLOSE(DETAIL);
	CLOSE(OUTCAPT1);
	CLOSE(OUTCAPT2);
        CLOSE(OUTDATA);
	CLOSE(OUTBANKTYPE);
   CLOSE(OUTTIME);
	{CLOSE(DISTLIST);}
	CLOSE(OUTFILE);
	CLOSE(AVG_VALS);
	CLOSE(SEDIMENT_FLAG);
END. {begin}