MODEL:
 ! Model of flow in an uncovered stream, 
e.g., small river, sewer, irrigation canal, etc.;
! Given a desired flow in volume/second, and
  a slope in length/length, find the width and
  height of a rectangular channel to achieve
  this flow at minimum cost of materials;
! Keywords: Hydraulics, Channel flow, Water flow,
  Gauckler-Manning-Strickler, River flow,
  Irrigation channel, Sewer flow;
DATA:
 S=0.001; ! Slope of the channel;
 Q=100;   ! Required flow in meters^3 or feet^3/ second;
 k = 1;   !units conversion factor, = 1 for SI units (meters)
         = 1.48592 = (3.28083499 ft/m)^(1/3), for English units (feet);
! Some bounds on b and h;
 blo = .5; bhi = 100;
 hlo = .5; hhi = 100;
ENDDATA

SUBMODEL STREAMFLO:
 ! Parameters:
    Q = discharge rate in meters^3 or feet^3/ second,
    S = slope of water surface, or head loss in length/length,
    n = Gauckler-Manning-Strickler roughness coefficient,
         higher means a slower flow,
         concrete has a coefficient of about .014,
    k = units conversion factor, = 1 for SI units (meters)
         = 1.4859 for English units (feet).
 Variables:
    b = channel width,
    P = wetted perimeter of rectangular cross section,
        e.g., amount of concrete needed,
    A = area of cross section,
    R = hydraulic radius,
    h = height of fluid in channel,
 ;

 min=P; ! Min sum of sidewalls + base;

 P=b+2*h;
 A=b*h; ! Area;
 P*R=A; ! Compute hydraulic radius;
! The Gauckler-Manning-Strickler formula for discharge rate gives
 flow velocity = (R^(2/3))*(S^0.5)/n, so discharge rate
           Q = A*(R^(2/3))*(S^0.5)/n, or ;
 k*A*(R^(2/3))=n*Q/(S^0.5);

! This model is hard to solve, so be sure to turn on
 either global or multi-start;
! Some reasonable bounds on b and h;
 @BND(blo, b, bhi);
 @BND(hlo, h, hhi);
ENDSUBMODEL


CALC:
 @SET("TERSEO",2); ! Turn off default output;
 @WRITE('Channel height, h, and width b, for various coefficients, n:',@NEWLINE(1));
 @WRITE(' to achieve a target flow volume of ',Q,' cubic ');
 @IFC( k #EQ# 1:
   @WRITE('meters/sec.',@NEWLINE(1));
  @ELSE
   @WRITE('feet/sec.',@NEWLINE(1));
     );
 @WRITE('   n         b       h       P ',@NEWLINE(1));
 FIRST = .01;  ! Smallest value of n to try;
 LAST = .03;   ! Largest value of n to try;
 STEP = .01;   ! Step size over n;
 n = FIRST;
 @WHILE( n #LE# LAST:
   @SOLVE( STREAMFLO);
   @WRITE( @FORMAT(n,'6.3f'),'  ',@FORMAT(b,'7.3f'),' ',@FORMAT(h,'7.3f'),'   ',@FORMAT(p,'7.3f'), @NEWLINE(1));
   n = n + STEP;
       );
ENDCALC