<< . .

. 47
( : 51)



. . >>

w[2*i*lattice_size_y+2*j+1]*=fact;
w[2*i*lattice_size_y+2*j+2]*=fact;
w[2*i*lattice_size_y+2*lattice_size_y-2*j+1]*=fact;
w[2*i*lattice_size_y+2*lattice_size_y-2*j+2]*=fact;
w[2*lattice_size_x*lattice_size_y-2*(i-1)*lattice_size_y-
2*(lattice_size_y-j)+1]*=fact;
w[2*lattice_size_x*lattice_size_y-2*(i-1)*lattice_size_y-
2*(lattice_size_y-j)+2]*=fact;
w[2*lattice_size_x*lattice_size_y-2*lattice_size_y-
2*(i-1)*lattice_size_y+2*(lattice_size_y-j)+1]*=fact;
w[2*lattice_size_x*lattice_size_y-2*lattice_size_y-
2*(i-1)*lattice_size_y+2*(lattice_size_y-j)+2]*=fact;}
fourn(w,nn,2,-1);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fprintf(fp2,"%f\n",w[2*(i-1)*lattice_size_y+2*j-1]
/(lattice_size_x*lattice_size_y));
fclose(fp1);
fclose(fp2);
}


A4.4 ADI technique applied to glacial isostatic response modeling

The following code uses the ADI technique to solve the ¬‚exure equation with spatially variable
rigidity. The code accepts two ¬les: westtelcomb3km (a 3 km/pixel grid of elastic thickness values) and
westelamask3km (a mask of equilibrium lines recti¬ed to the same coordinate system as the elastic
thickness grid). These ¬les can be downloaded from the CUP website. The ADI method is repeated
5000 times in this example (the method converges slowly). In practice, the number of iterations
should be determined by the required level of accuracy.
#define PI 3.141592653589793
#define adiiterations 5000

int idum,*iup,*idown,*jup,*jdown,lattice_size_y,lattice_size_x;
float **D,delrho,deltax,deltay,**load,**deflect,**deflect2,te;
260 CODES FOR SOLVING THE FLEXURE EQUATION


void solveimpl()
{
int i,j,count;
float d,**ax,**ay,*xx,*xy,*rx,*ry,**alx,**aly,tot;
unsigned long *indx,*indy;

indx=lvector(1,lattice_size_x);
indy=lvector(1,lattice_size_y);
ax=matrix(1,lattice_size_x,1,5);
ay=matrix(1,lattice_size_y,1,5);
xx=vector(1,lattice_size_x);
xy=vector(1,lattice_size_y);
rx=vector(1,lattice_size_x);
ry=vector(1,lattice_size_y);
alx=matrix(1,lattice_size_x,1,2);
aly=matrix(1,lattice_size_y,1,2);
for (count=1;count<=adiiterations;count++)
{for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
deflect2[i][j]=deflect[i][j];
tot=0;
for (j=1;j<=lattice_size_y;j++) {
for (i=1;i<=lattice_size_x;i++)
if ((i>2)&&(i<lattice_size_x-1))
{rx[i]=-load[i][j];
ax[i][1]=D[i][j];
ax[i][2]=-4*D[i][j];
ax[i][3]=12*D[i][j]+delrho;
rx[i]+=-D[i][j]*(deflect2[i][jup[jup[j]]]+deflect2[i][jdown[jdown[j]]]-
4*(deflect2[i][jup[j]]+deflect2[i][jdown[j]]));
ax[i][4]=-4*D[i][j];
ax[i][5]=D[i][j];}
else
{rx[i]=0;
ax[i][1]=0;
ax[i][2]=0;
ax[i][3]=1;
ax[i][4]=0;
ax[i][5]=0;}
bandec(ax,lattice_size_x,2,2,alx,indx,&d);
banbks(ax,lattice_size_x,2,2,alx,indx,rx);
for (i=1;i<=lattice_size_x;i++)
deflect[i][j]=rx[i];}
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
deflect2[i][j]=deflect[i][j];
for (i=1;i<=lattice_size_x;i++) {
for (j=1;j<=lattice_size_y;j++)
if ((j>2)&&(j<lattice_size_y-1))
A4.4 ADI TECHNIQUE APPLIED TO GLACIAL ISOSTATIC RESPONSE MODELING 261


{ry[j]=-load[i][j];
ay[j][1]=D[i][j];
ay[j][2]=-4*D[i][j];
ay[j][3]=12*D[i][j]+delrho;
ry[j]+=-D[i][j]*(deflect2[iup[iup[i]]][j]+deflect2[idown[idown[i]]][j]-
4*(deflect2[iup[i]][j]+deflect2[idown[i]][j]));
ay[j][4]=-4*D[i][j];
ay[j][5]=D[i][j];}
else
{ry[j]=0;
ay[j][1]=0;
ay[j][2]=0;
ay[j][3]=1;
ay[j][4]=0;
ay[j][5]=0;}
bandec(ay,lattice_size_y,2,2,aly,indy,&d);
banbks(ay,lattice_size_y,2,2,aly,indy,ry);
for (j=1;j<=lattice_size_y;j++)
deflect[i][j]=ry[j];}
}
}

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

fp1=fopen("westtelcomb3km","r");
fp2=fopen("westsnowlinemask3km","r");
fp3=fopen("westdeflectcombnew","w");
lattice_size_x=626;
lattice_size_y=642;
deltax=3; /* km */
deltay=3; /* km */
delrho=6000; /* kg/m^3 */
deflect=matrix(1,lattice_size_x,1,lattice_size_y);
deflect2=matrix(1,lattice_size_x,1,lattice_size_y);
load=matrix(1,lattice_size_x,1,lattice_size_y);
D=matrix(1,lattice_size_x,1,lattice_size_y);
setupgridneighborsperiodic();
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
{fscanf(fp1,"%f",&te);
if ((i<5)||(i>lattice_size_x-5)||(j<5)||(j>lattice_size_y-5)) te=0;
D[i][j]=te*te*te*70000000/(deltax*deltax*deltax*deltax)/
(12*(1-0.25*0.25));
fscanf(fp2,"%f",&load[i][j]);
load[i][j]*=0.82; /* rho_c/rho_m */
if ((i<5)||(i>lattice_size_x-5)||(j<5)||(j>lattice_size_y-5)) load[i][j]=0;
deflect[i][j]=load[i][j];}
262 CODES FOR SOLVING THE FLEXURE EQUATION


solveimpl();
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fprintf(fp3,"%f\n",deflect[i][j]);
fclose(fp1);
fclose(fp2);
fclose(fp3);
}
Appendix 5




Codes for modeling non-Newtonian ¬‚ows

A5.1 2D radially symmetric lava ¬‚ow simulation

The following code implements the 2D radially symmetric gravity ¬‚ow model with temperature-
dependent viscosity described in Section 6.3.


main()
{ FILE *fp1;
int lattice_size_x,*iup,*idown,i;
float deltar,nu,ro,rn,*height,*heightold,*temp,*tempold,*flux;
float simulationlength,outputinterval,timestep,elapsedtime,maxdelt,maxdelh,del;

fp1=fopen("viscousflownu100","w");
lattice_size_x=260;
deltar=0.02;
nu=100.0;
ro=1.0;
rn=1.1;
iup=ivector(1,lattice_size_x);
idown=ivector(1,lattice_size_x);
height=vector(1,lattice_size_x);
heightold=vector(1,lattice_size_x);
temp=vector(1,lattice_size_x);
tempold=vector(1,lattice_size_x);
flux=vector(1,lattice_size_x);
simulationlength=10.0;
outputinterval=0.1;
for (i=1;i<=lattice_size_x;i++)
{temp[i]=exp(-pow(i*deltar/rn,20));
tempold[i]=temp[i];
if (i*deltar<=ro) height[i]=pow(4.0/(1+nu)*log(rn/ro),0.25);
else if (i*deltar<=rn) height[i]=
pow(4.0/(1+nu)*log(rn/(i*deltar)),0.25); else height[i]=0.001;
heightold[i]=height[i];
flux[i]=0;
264 CODES FOR MODELING NON-NEWTONIAN FLOWS


iup[i]=i+1;idown[i]=i-1;}
iup[lattice_size_x]=lattice_size_x;idown[1]=1;
timestep=0.00001;
elapsedtime=0.0;
while (elapsedtime<simulationlength)
{maxdelt=0;maxdelh=0;
flux[1]=1;
for (i=1;i<=lattice_size_x;i++)
if (i*deltar<ro) flux[i]=1;
else flux[i]=i*(1+nu*tempold[i])*pow(0.5*(heightold[idown[i]]+heightold[i]),3)*
(heightold[idown[i]]-heightold[i]);
for (i=1;i<=lattice_size_x-1;i++)
{del=timestep/(i*deltar*deltar)*(flux[i]-flux[iup[i]]);
height[i]+=del;
if (fabs(del)>maxdelh) maxdelh=del;
del=0;
if (heightold[i]>0)
{del=timestep*(flux[i]*(tempold[idown[i]]-tempold[i])/
(heightold[i]*i*deltar*deltar)-
tempold[i]/(heightold[i]*heightold[i]));
if (i*deltar<ro) del=0;
temp[i]+=del;} else del=0;
if (fabs(del)>maxdelt) maxdelt=fabs(del);
if (temp[i]<0) temp[i]=0;
if (temp[i]>1) temp[i]=1;}
for (i=1;i<=lattice_size_x;i++)
if (i*deltar<ro) height[i]=height[(int)(ro/deltar)];
elapsedtime+=timestep;
if ((maxdelt>0.1)||(maxdelh>0.1))
{elapsedtime-=timestep;
timestep/=2.0;
for (i=1;i<=lattice_size_x;i++)
{height[i]=heightold[i];
temp[i]=tempold[i];}}
else
if ((maxdelt<0.01)&&(maxdelh<0.01)) timestep*=1.2;
for (i=1;i<=lattice_size_x;i++)
{heightold[i]=height[i];
tempold[i]=temp[i];}}
for (i=1;i<=lattice_size_x;i++)
fprintf(fp1,"%d %f %f\n",i,height[i],temp[i]);
fclose(fp1);
}



A5.2 Sandpile method for ice-sheet and glacier reconstruction

The following code implements the sandpile method of ice-sheet or glacier reconstruction using the
Greenland ice sheet as an example. The code accepts the input ¬les greenlandsmallbed (a 10 km/pixel
A5.2 SANDPILE METHOD FOR ICE-SHEET AND GLACIER RECONSTRUCTION 265


DEM of the bed topography) and greenlandsmallmask (a 10 km/pixel mask grid de¬ning locations of
ice coverage) from the local directory and outputs the ¬le greenlandsmallsurf. The routine checkmin
is just a simple version of fillinpitsandflats that only applies to the masked area (i.e. the area
with ice cover).

#define fillincrement 0.1

float taugam,**topo,block,*heightvec,**mask,**height,delta;
int *heightvecind,*iup,*jup,*idown,*jdown,lattice_size_x,lattice_size_y,change;

void addiflessthanthreshold(i,j)
int i,j;
{ float slope,slopex,slopey;

height[i][j]+=block;
slopex=height[i][j]-height[iup[i]][j];
if (height[i][j]-height[idown[i]][j]>slopex) slopex=height[i][j]-height[idown[i]][j];
slopey=height[i][j]-height[i][jup[j]];
if (height[i][j]-height[i][jdown[j]]>slopey) slopey=height[i][j]-height[i][jdown[j]];
slope=sqrt(slopex*slopex+slopey*slopey)/delta;
/* this version implements slope-dependent threshold
tau(S)=15S^0.55 */
if (pow(slope,0.45)*(height[i][j]-topo[i][j])>taugam*15)
height[i][j]-=block;
else change=1;
}

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

min=10000;
if (mask[i][j]>0.1) {
if (height[iup[i]][j]<min) min=height[iup[i]][j];
if (height[idown[i]][j]<min) min=height[idown[i]][j];
if (height[i][jup[j]]<min) min=height[i][jup[j]];
if (height[i][jdown[j]]<min) min=height[i][jdown[j]];
if (height[i][j]<min) {height[i][j]=min+0.1;
checkmin(iup[i],j);
checkmin(idown[i],j);
checkmin(i,jup[j]);
checkmin(i,jdown[j]);} }
}

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

fp1=fopen("greenlandsmallbed","r");
fp2=fopen("greenlandsmallmask","r");
266 CODES FOR MODELING NON-NEWTONIAN FLOWS


fp3=fopen("greenlandsmallsurf","w");
lattice_size_x=150;
lattice_size_y=280;
delta=10000.0; /* km */
taugam=100000.0/(920.0*9.8); /* tau = 10^5 Pa */
height=matrix(1,lattice_size_x,1,lattice_size_y);
mask=matrix(1,lattice_size_x,1,lattice_size_y);
topo=matrix(1,lattice_size_x,1,lattice_size_y);
setupgridneighbors();
heightvec=vector(1,lattice_size_x*lattice_size_y);
heightvecind=ivector(1,lattice_size_x*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]);
fscanf(fp2,"%f",&mask[i][j]);
height[i][j]=topo[i][j];}
block=100; /* start at 100 m block size */
while (block>0.9) /* stop when block size < 1 m */
{change=1;
while (change>0)
{change=0;
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
heightvec[(j-1)*lattice_size_x+i]=height[i][j];
indexx(lattice_size_x*lattice_size_y,heightvec,heightvecind);
t=0;
while (t<lattice_size_x*lattice_size_y)
{t++;
i=(heightvecind[t])%lattice_size_x;
if (i==0) i=lattice_size_x;
j=(heightvecind[t])/lattice_size_x+1;
if (mask[i][j]>0.1)
addiflessthanthreshold(i,j);}
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
checkmin(i,j);}
block=block/3.33;}
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
{if (height[i][j]>0)
fprintf(fp3,"%f\n",height[i][j]);
else fprintf(fp3,"0\n");}
fclose(fp1);
fclose(fp2);
fclose(fp3);
}
Appendix 6




Codes for modeling instabilities


A6.1 Werner™s eolian dune model

The following code implements Werner™s eolian dune model described in Section 7.4 with parameters
l = 5, ps = 0.6, pns = 0.4, and a shadow-zone angle equal to half of the angle of repose. The code was
written by the author based on the description in Werner (1995). The initial condition is a ¬‚at bed
with three sand slabs at each grid point.


int **height,*iup,*jup,*idown,*jdown,lattice_size_x,lattice_size_y;
float **mask,thresh;

void avalanchedown(i,j)
int i,j;
{ float hshadow;
int i2,violate,min,mini,minj;

violate=0;min=height[i][j];mini=i;minj=j;
if (height[i][j]-height[iup[i]][j]>thresh)
{violate=1;
if (height[iup[i]][j]<min)
{min=height[iup[i]][j];mini=iup[i];minj=j;}}
if (height[i][j]-height[idown[i]][j]>thresh)
{violate=1;
if (height[idown[i]][j]<min)
{min=height[idown[i]][j];mini=idown[i];minj=j;}}
if (height[i][j]-height[i][jup[j]]>thresh)
{violate=1;
if (height[i][jup[j]]<min)
{min=height[i][jup[j]];mini=i;minj=jup[j];}}
if (height[i][j]-height[i][jdown[j]]>thresh)
{violate=1;
if (height[i][jdown[j]]<min)
{min=height[i][jdown[j]];mini=i;minj=jdown[j];}}
if (violate==1)
{height[i][j]--;
268 CODES FOR MODELING INSTABILITIES


height[mini][minj]++;
avalanchedown(mini,minj);}
else
{hshadow=height[i][j]-0.5;
i2=iup[i];
while (height[i2][j]<hshadow)
{if (mask[i2][j]<hshadow) mask[i2][j]=hshadow;
i2=iup[i2];hshadow-=0.5;}}
}

void avalancheup(i,j)
int i,j;
{ float hshadow;
int i2,violate,min,mini,minj;

violate=0;min=height[i][j];mini=i;minj=j;
if (height[iup[i]][j]-height[i][j]>thresh)

<< . .

. 47
( : 51)



. . >>