<< . .

. 44
( : 51)



. . >>

fclose(fp2);
}


A1.5 Numerical integration of Fourier“Bessel terms

The following code accepts a list of the ¬rst 300 zeros of the Bessel function J 0 (r ) from the local
directory and performs numerical integration to calculate the value of the terms in Eq. (2.57) for the
case of a volcanic cone with r f = 200 m, rr = 500 m, and rc = 1500 m.

int nsum;
float *zeros;

float func(x)
float x;
{
return x*x*bessj0(zeros[nsum]*x);
}

main()
{ FILE *fp0,*fp1;
double a,zero=0.0,s1,s2,s3,rf,rr,rc;

fp0=fopen("besselzeros","r");
fp1=fopen("besselzerosintegral","w");
zeros=vector(1,300);
a=3000.0;
rf=200.0/a;
rr=500.0/a;
rc=1500.0/a;
for (nsum=1;nsum<=300;nsum++)
{fscanf(fp0,"%f",&zeros[nsum]);
s1=qsimp(func,zero,rf);
s2=qsimp(func,zero,rr);
s3=qsimp(func,zero,rc);
A1.5 NUMERICAL INTEGRATION OF FOURIER“BESSEL TERMS 233


fprintf(fp1,"%f %12.10f %12.10f %12.10f\n",zeros[nsum],s1,s2,s3);}
fclose(fp1);
fclose(fp2);
}

The following code accepts output from the code above, computes the results of the Bessel series
in Eq. (2.57), and prints the results to an output ¬le.
main()
{
FILE *fp1,*fp2;
float rf,rr,rc,*f,kappaT,deltax,a,*prefact,*topo,*J0rr,*J0rc,*J1rr,*J1rc;
float *J0rf,*J1rf,*Intrf,*J0r0,*J1r0,*Intrr,*Intrc;
int lattice_size_x,i,nsum;

fp1=fopen("besselzerosintegral","r");
fp2=fopen("volcanicconeoutput","w");
kappaT=0.0001; /* dimensionless scaled to a^2 */
lattice_size_x=300;
deltax=10.0;
a=lattice_size_x*deltax;
rf=200.0/a;
rr=500.0/a;
rc=1500.0/a;
topo=vector(1,lattice_size_x);
f=vector(1,lattice_size_x);
zeros=vector(1,300);
J0rc=vector(1,300);
J1rc=vector(1,300);
Intrc=vector(1,300);
J0rr=vector(1,300);
J1rr=vector(1,300);
Intrr=vector(1,300);
J0rf=vector(1,300);
J1rf=vector(1,300);
Intrf=vector(1,300);
prefact=vector(1,300);
for (nsum=1;nsum<=300;nsum++)
{fscanf(fp1,"%f %f %f %f",&zeros[nsum],&Intrf[nsum],&Intrr[nsum],&Intrc[nsum]);
J0rc[nsum]=bessj0(rc*zeros[nsum]);
J1rc[nsum]=bessj1(rc*zeros[nsum]);
J0rr[nsum]=bessj0(rr*zeros[nsum]);
J1rr[nsum]=bessj1(rr*zeros[nsum]);
J0rf[nsum]=bessj0(rf*zeros[nsum]);
J1rf[nsum]=bessj1(rf*zeros[nsum]);}
for (i=1;i<=lattice_size_x;i++)
{if (i*deltax/a<rf)
f[i]=1+(rf-rr)/(rc-rr);
else if (i*deltax/a<rr)
f[i]=1+(i*deltax/a-rr)/(rc-rr);
234 CODES FOR SOLVING THE DIFFUSION EQUATION


else if (i*deltax/a<rc)
f[i]=1+(rr-i*deltax/a)/(rc-rr);
else f[i]=0;
/* print initial condition first */
fprintf(fp2,"%f %12.10f\n",i*deltax,f[i]);}
for (i=1;i<=lattice_size_x;i++)
topo[i]=0;
for (nsum=1;nsum<=300;nsum++)
prefact[nsum]=2.0/zeros[nsum]*exp(-kappaT*zeros[nsum]*zeros[nsum])/
(bessj1(zeros[nsum])*bessj1(zeros[nsum]));
for (i=1;i<=lattice_size_x;i++)
{for (nsum=1;nsum<=300;nsum++)
topo[i]+=prefact[nsum]*bessj0(zeros[nsum]*i*deltax/a)*(rf*rf/(rc-rr)*
J1rf[nsum]-2*rr*rr/(rc-rr)*J1rr[nsum]+(1+rr/(rc-rr))*
rc*J1rc[nsum]+2*zeros[nsum]/(rc-rr)*Intrr[nsum]-zeros[nsum]/
(rc-rr)*(Intrf[nsum]+Intrc[nsum]));
fprintf(fp2,"%f %12.10f\n",i*deltax,topo[i]);}
fclose(fp1);
fclose(fp2);
}
Appendix 2




Codes for ¬‚ow routing


A2.1 Filling in pits and ¬‚ats in a DEM

The following code is used to eliminate pits and ¬‚ats from a DEM. The output from this type of
processing is sometimes called a hydrologically corrected DEM. The code is currently constructed to
accept a DEM from the local directory named inputdem that has 300 columns and 300 rows. In this
routine, all of the grid points are tested to determine whether they have a downslope neighbor. If
they do not have a downslope neighbor then they must be a ¬‚at or the bottom of a pit. In such
cases, the local elevation value is incremented by a small value fillincrement and the neighbors
are searched, recursively, for additional pits and ¬‚ats. Any additional pits and ¬‚ats that are found
also have their elevation values incremented. This procedure stops when all of the pixels in the pit
or ¬‚at have a downslope neighbor (i.e. the areas drains). For very large pits or ¬‚ats, it is possible
for this routine to run out of memory by calling itself too many times. In such cases, the value
of fillincrement can be raised, and/or a counter can be set up fillinpitsandflats that limits the
number of times that the routine can call itself before returning to the main program.
#define fillincrement 0.01

float **topo,
int lattice_size_x,lattice_size_y,*iup,*idown,*jup,*jdown;

void fillinpitsandflats(i,j)
int i,j;
{ float min;

min=topo[i][j];
if (topo[iup[i]][j]<min) min=topo[iup[i]][j];
if (topo[idown[i]][j]<min) min=topo[idown[i]][j];
if (topo[i][jup[j]]<min) min=topo[i][jup[j]];
if (topo[i][jdown[j]]<min) min=topo[i][jdown[j]];
if (topo[iup[i]][jup[j]]<min) min=topo[iup[i]][jup[j]];
if (topo[idown[i]][jup[j]]<min) min=topo[idown[i]][jup[j]];
if (topo[idown[i]][jdown[j]]<min) min=topo[idown[i]][jdown[j]];
if (topo[iup[i]][jdown[j]]<min) min=topo[iup[i]][jdown[j]];
if ((topo[i][j]<=min)&&(i>1)&&(j>1)&&(i<lattice_size_x)&&(j<lattice_size_y))
236 CODES FOR FLOW ROUTING


{topo[i][j]=min+fillincrement;
fillinpitsandflats(i,j);
fillinpitsandflats(iup[i],j);
fillinpitsandflats(idown[i],j);
fillinpitsandflats(i,jup[j]);
fillinpitsandflats(i,jdown[j]);
fillinpitsandflats(iup[i],jup[j]);
fillinpitsandflats(idown[i],jup[j]);
fillinpitsandflats(idown[i],jdown[j]);
fillinpitsandflats(iup[i],jdown[j]);}
}

main ()
{ FILE *fp1,*fp2;
int i,j;

fp1=fopen("inputdem","r");
fp2=fopen("filleddem","w");
setupgridneighbors();
lattice_size_x=300;
lattice_size_y=300;
topo=matrix(1,lattice_size_x,1,lattice_size_y);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fscanf(fp1,"%f ",&topo[i][j]);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fillinpitsandflats(i,j);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fprintf(fp2,"%f\n",topo[i][j]);
fclose(fp1);
fclose(fp2);
}




A2.2 MFD ¬‚ow routing method

The following code calculates the contributing area of a 30 m/pixel DEM using MFD ¬‚ow routing
assuming a local input ¬le named inputdem.


#define oneoversqrt2 0.707106781187

float **topo,**flow,**flow1,**flow2,**flow3,**flow4,**flow5,**flow6,**flow7,**flow8;
int lattice_size_x,lattice_size_y,*iup,*idown,*jup,*jdown;
A2.2 MFD FLOW ROUTING METHOD 237


void mfdflowroute(i,j)
int i,j;
{ float tot;

tot=0;
if (topo[i][j]>topo[iup[i]][j])
tot+=pow(topo[i][j]-topo[iup[i]][j],1.1);
if (topo[i][j]>topo[idown[i]][j])
tot+=pow(topo[i][j]-topo[idown[i]][j],1.1);
if (topo[i][j]>topo[i][jup[j]])
tot+=pow(topo[i][j]-topo[i][jup[j]],1.1);
if (topo[i][j]>topo[i][jdown[j]])
tot+=pow(topo[i][j]-topo[i][jdown[j]],1.1);
if (topo[i][j]>topo[iup[i]][jup[j]])
tot+=pow((topo[i][j]-topo[iup[i]][jup[j]])*oneoversqrt2,1.1);
if (topo[i][j]>topo[iup[i]][jdown[j]])
tot+=pow((topo[i][j]-topo[iup[i]][jdown[j]])*oneoversqrt2,1.1);
if (topo[i][j]>topo[idown[i]][jup[j]])
tot+=pow((topo[i][j]-topo[idown[i]][jup[j]])*oneoversqrt2,1.1);
if (topo[i][j]>topo[idown[i]][jdown[j]])
tot+=pow((topo[i][j]-topo[idown[i]][jdown[j]])*oneoversqrt2,1.1);
if (topo[i][j]>topo[iup[i]][j])
flow1[i][j]=pow(topo[i][j]-topo[iup[i]][j],1.1)/tot;
else flow1[i][j]=0;
if (topo[i][j]>topo[idown[i]][j])
flow2[i][j]=pow(topo[i][j]-topo[idown[i]][j],1.1)/tot;
else flow2[i][j]=0;
if (topo[i][j]>topo[i][jup[j]])
flow3[i][j]=pow(topo[i][j]-topo[i][jup[j]],1.1)/tot;
else flow3[i][j]=0;
if (topo[i][j]>topo[i][jdown[j]])
flow4[i][j]=pow(topo[i][j]-topo[i][jdown[j]],1.1)/tot;
else flow4[i][j]=0;
if (topo[i][j]>topo[iup[i]][jup[j]])
flow5[i][j]=pow((topo[i][j]-topo[iup[i]][jup[j]])*oneoversqrt2,1.1)/tot;
else flow5[i][j]=0;
if (topo[i][j]>topo[iup[i]][jdown[j]])
flow6[i][j]=pow((topo[i][j]-topo[iup[i]][jdown[j]])*oneoversqrt2,1.1)/tot;
else flow6[i][j]=0;
if (topo[i][j]>topo[idown[i]][jup[j]])
flow7[i][j]=pow((topo[i][j]-topo[idown[i]][jup[j]])*oneoversqrt2,1.1)/tot;
else flow7[i][j]=0;
if (topo[i][j]>topo[idown[i]][jdown[j]])
flow8[i][j]=pow((topo[i][j]-topo[idown[i]][jdown[j]])*oneoversqrt2,1.1)/tot;
else flow8[i][j]=0;
flow[iup[i]][j]+=flow[i][j]*flow1[i][j];
flow[idown[i]][j]+=flow[i][j]*flow2[i][j];
flow[i][jup[j]]+=flow[i][j]*flow3[i][j];
flow[i][jdown[j]]+=flow[i][j]*flow4[i][j];
238 CODES FOR FLOW ROUTING


flow[iup[i]][jup[j]]+=flow[i][j]*flow5[i][j];
flow[iup[i]][jdown[j]]+=flow[i][j]*flow6[i][j];
flow[idown[i]][jup[j]]+=flow[i][j]*flow7[i][j];
flow[idown[i]][jdown[j]]+=flow[i][j]*flow8[i][j];
}

main ()
{ FILE *fp1,*fp2;
int i,j,t,*topovecind;
float *topovec,delta;

fp1=fopen("inputdem","r");
fp2=fopen("MFDcontribarea","w");
setupgridneighbors();
lattice_size_x=300;
lattice_size_y=300;
delta=30.0; /* m */
topo=matrix(1,lattice_size_x,1,lattice_size_y);
topovec=vector(1,lattice_size_x*lattice_size_y);
topovecind=ivector(1,lattice_size_x*lattice_size_y);
flow=matrix(1,lattice_size_x,1,lattice_size_y);
flow1=matrix(1,lattice_size_x,1,lattice_size_y);
flow2=matrix(1,lattice_size_x,1,lattice_size_y);
flow3=matrix(1,lattice_size_x,1,lattice_size_y);
flow4=matrix(1,lattice_size_x,1,lattice_size_y);
flow5=matrix(1,lattice_size_x,1,lattice_size_y);
flow6=matrix(1,lattice_size_x,1,lattice_size_y);
flow7=matrix(1,lattice_size_x,1,lattice_size_y);
flow8=matrix(1,lattice_size_x,1,lattice_size_y);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
{fscanf(fp1,"%f ",&topo[i][j]);
flow[i][j]=delta*delta;}
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fillinpitsandflats(i,j);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
topovec[(j-1)*lattice_size_x+i]=topo[i][j];
indexx(lattice_size_x*lattice_size_y,topovec,topovecind);
t=lattice_size_x*lattice_size_y+1;
while (t>1)
{t--;
i=(topovecind[t])%lattice_size_x;
if (i==0) i=lattice_size_x;
j=(topovecind[t])/lattice_size_x+1;
if (i==lattice_size_x) j--;
mfdflowroute(i,j);}
for (j=1;j<=lattice_size_y;j++)
A2.3 SUCCESSIVE FLOW ROUTING WITH MFD METHOD 239


for (i=1;i<=lattice_size_x;i++)
fprintf(fp2,"%f\n",flow[i][j]);
fclose(fp1);
fclose(fp2);
}


A2.3 Successive ¬‚ow routing with MFD method

The following code performs a successive ¬‚ow routing analysis on an input DEM. The user is prompted
for the names of the input and output DEM ¬les, the size of the DEM, and the prescribed runoff
depth. The routine uses Manning™s equation to convert a prescribed discharge into an equivalent
depth, assuming steady ¬‚ow. As such, the local ¬‚ow depths depend sensitively on the local bed
slope. For some applications of this code it may be necessary to ¬rst smooth the input DEM in order
to remove topographic irregularities. In the example below it is assumed that the smallest bed slope
was 1% (any smaller value is made equal to 1%).

#define niterations 100
#define oneoversqrt2 0.707106781186

float deltax,**topo,**toposave,*topovec,**depth,**slope;
float **flow,**flow1,**flow2,**flow3,**flow4,**flow5,**flow6,**flow7,**flow8;
int lattice_size_x,lattice_size_y,*iup,*idown,*jup,*jdown;

void calculatealongchannelslope(i,j)
int i,j;
{ float down;

down=0;
if (topo[iup[i]][j]-topo[i][j]<down) down=topo[iup[i]][j]-topo[i][j];
if (topo[idown[i]][j]-topo[i][j]<down) down=topo[idown[i]][j]-topo[i][j];
if (topo[i][jup[j]]-topo[i][j]<down) down=topo[i][jup[j]]-topo[i][j];
if (topo[i][jdown[j]]-topo[i][j]<down) down=topo[i][jdown[j]]-topo[i][j];
if ((topo[iup[i]][jup[j]]-topo[i][j])*oneoversqrt2<down)
down=(topo[iup[i]][jup[j]]-topo[i][j])*oneoversqrt2;
if ((topo[idown[i]][jup[j]]-topo[i][j])*oneoversqrt2<down)
down=(topo[idown[i]][jup[j]]-topo[i][j])*oneoversqrt2;
if ((topo[iup[i]][jdown[j]]-topo[i][j])*oneoversqrt2<down)
down=(topo[iup[i]][jdown[j]]-topo[i][j])*oneoversqrt2;
if ((topo[idown[i]][jdown[j]]-topo[i][j])*oneoversqrt2<down)
down=(topo[idown[i]][jdown[j]]-topo[i][j])*oneoversqrt2;
slope[i][j]=fabs(down)/deltax;
}

main()
{ FILE *fp1,*fp2;
float runoff,manningsn,*topovec;
char inputfile[100],outputfile[100];
int i,j,k,t,*topovecind;
240 CODES FOR FLOW ROUTING


scanf("%s",&inputfile);
scanf("%s",&outputfile);
scanf("%d",&lattice_size_x);
scanf("%d",&lattice_size_y);
scanf("%f",&deltax); /* m */
scanf("%f",&runoff); /* m */
scanf("%f",&manningsn);
fp1=fopen(inputfile,"r");
fp2=fopen(outputfile,"w");
setupgridneighbors();
topo=matrix(1,lattice_size_x,1,lattice_size_y);
toposave=matrix(1,lattice_size_x,1,lattice_size_y);
topovec=vector(1,lattice_size_x*lattice_size_y);
topovecind=ivector(1,lattice_size_x*lattice_size_y);
depth=matrix(1,lattice_size_x,1,lattice_size_y);
slope=matrix(1,lattice_size_x,1,lattice_size_y);
flow=matrix(1,lattice_size_x,1,lattice_size_y);
flow1=matrix(1,lattice_size_x,1,lattice_size_y);
flow2=matrix(1,lattice_size_x,1,lattice_size_y);
flow3=matrix(1,lattice_size_x,1,lattice_size_y);
flow4=matrix(1,lattice_size_x,1,lattice_size_y);
flow5=matrix(1,lattice_size_x,1,lattice_size_y);
flow6=matrix(1,lattice_size_x,1,lattice_size_y);
flow7=matrix(1,lattice_size_x,1,lattice_size_y);
flow8=matrix(1,lattice_size_x,1,lattice_size_y);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fscanf(fp1,"%f",&topo[i][j]);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fillinpitsandflats(i,j);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
{toposave[i][j]=topo[i][j];
calculatechannelslope(i,j);
if (slope[i][j]<0.01) slope[i][j]=0.01;}
for (k=1;k<=niterations;k++)
{printf("iteration = %d/%d\n",k,niterations);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
{flow[i][j]=pow(runoff,1.66667)*sqrt(slope[i][j])*deltax/manningsn;
topovec[(j-1)*lattice_size_x+i]=topo[i][j];}
indexx(lattice_size_x*lattice_size_y,topovec,topovecind);
t=lattice_size_x*lattice_size_y+1;
while (t>1)
{t--;
i=(topovecind[t])%lattice_size_x;
if (i==0) i=lattice_size_x;
j=(topovecind[t])/lattice_size_x+1;

<< . .

. 44
( : 51)



. . >>