在线等待啊。。。

#include <stdio.h>  
#include <math.h>
#include "utilities.h"
#define NVAR 4
#define NSTEP 1000
#define PI 3.14159
void derivs(x,y,dydx,ro,zo,alpha,lam1,lam2,Bflag,rflag) 
float x, y[], dydx[];
double ro,zo,alpha, lam1,lam2;
int *Bflag, *rflag;
{
double g, Mrel, gamma, dens, densdot, denssw, dpwdz, deltdens, visc, D;
double Cs, Cdm, Cm, Z, ub, K, Bpos, Bneg, qz, propvec[7];
static double rD, qD, qrel, densrel, densref, massflux;

g = 981.0; /*gravitational accereration cm/s2 */
Mrel = 133.0; /*release C02 mass flux kg/s */


/*begin calculations*/
if(x == zo) 
{
properties(zo, propvec); /* determine initial values*/
densrel=propvec[2]; /* C02 density kg/m3 */
densref=propvec[0]; /* seawater density kg/m3 */
qrel=Mrel/densrel; /* release flow m3/s */
massflux=Mrel; /* release mass fluz kg/s */
}
properties(x, propvec);  /*determine co2 and ocean properties at x */
denssw=propvec[0];     /* seawater density kg/m3 */
dpwdz=propvec[1];      /* seawater density gradient kg/m4 */
dens= propvec[2];      /* co2 density kg/m3 */
densdot = propvec[3];   /* co2 density gradient kg/m9 */
Cs   = propvec[4];     /* co2 solubility kg/m3 */
D    = propvec[5];     /* molecular diffusion coefficient cm2/s */
visc = propvec[6];     /* water viscosity cm2/s */


/* z=x, r=y[1], delta rho w = y[2], b=y[3], Um=y[4] */
if(y[1]<=0.01) {       /*If r is near zero, indicate absence of bubbles*/
y[1] = 0.00001;
*rflag = 1;
}
if(y[2] > 1000000000.0) /* Numerical limit to delta rho w */
  y[2] = 1000000000.0;

if(y[3]>1000000000.0) /* Numerical limit to b */
  y[3] = 1000000000.0;

if(y[4]<0.0)         /* Numerical limit to Um */
  y[4]=0.00001;

if(y[3]<0.0)         /* Numerical limit to b*/
  y[3]=0.00001;


/* Solve for local properties */
  Z=0.217*2*y[1]*cqrt(g/(visc*visc);
if(zo > 500.0) {         /* if starts as liquid */
 if(x > 500.0) {         /* liquid phase */
   qz = qrel*(y[1]/ro)*(y[1]/ro)*(y[1]/ro);
   /*Clift Ub*/
   ub=sqrt(g*2.0*y[1]*(denssw-dens)/denssw)*0.00711;
   /*Teyball Ub*/
   /*ub=O.2392*pow(((denssw-dens)/ 1000.0),0.28)*pow((dens/1000.0),-0.45);*/
    rD = y[1];
  qD = qz;
}
if(x <= 500.0 ) {        /* change to vapor */ 
  qrel = qD*5.156;
  ro = rD*1.728;
  qz = qrel*(y[1]/ro)*(y[1]/ro)*(y[1]/ro);
  ub = (1.0/100.0)*(2.355)*(108.4/Z+sqrt(Z/0.5479));/*m/s*/
}
}


if(zo <= 500.0){       /* if starts as vapor */
qz = qrel*(y[1]/ro)*(y[1]/ro)*(y[1]/ro);
ub = (1.0/100.0)*(2.355)*(108.4/Z+sqrt(Z/0.5479));/*m/s*/
}
/* WCK */
/*K = (4.917e-3)*sqrt(100.0*ub)/(sqrt(2.O*y[1])); */
/*Clift K for Cap Bubbles*/
K = 1.25*sqrt(sqrt(g*(denssw-dens)/denssw))*sqrt(D)/(sqrt(sqrt(2.0*y[1])));
 /*Cussler K*/
/*K = 0.30919*cbrt((denssw- dens)/ denssw);*/
Cm = ((qz/(PI*y[3]*y[3]*lam1*lam1))/((y[4]/(1+lam1*lam1))+ub)); 
deltdens = (densref-dens)/densref;
gamma = 1.0;


/*PEELING MODEL */

Bpos = Cm*deltdens*lam1*lam1;
Bneg = y[2]*lam2*lam2/densref;
if((Bpos<Bneg)&&(y[1]>0.01))
{
y[2] =0.5*y[2];
y[3] =0.707*y[3]; /*b*/
/*y[4] = [4l];  /*U*/
*Bflag =1;
massflux =0.5*massflux+0.5*dens*qz;
}


/* C02 Accumulation */
/*Cdm = 2* (masflu - qz*dens)/(y[4]*PI*y[3]*y[3]);
printf("%10.6lf %10.6lf \n", x, Cdm);*/


/* solve differential equations */
/* dydx[1] = dr/dz, 2->dens, 3->b, 4->um */

dydx[1] = (K*Cs*0.85/(dens*(y[4]+ub)))-(y[1]*densdot/(3*dens)); /*cm/s*/
dydx[2] = -((1+lam2*lam2/(lam2*lam2))*dpwdz-(2*alpha*y[2]/y[3]));
dydx[3] = -((2*alpha)-(g*y[3]/(gamma*100.0*y[4]*y[4]))*(Cm*deltdens*laml*laml-(y[2]*lam2*lam2/densref)));
dydx[4] = -((2*g/(gamma*100.0*y[4]))*(Cm*deltdens*laml*laml-(y[2]*lam2*lam2/densref))-(2*alpha*y[4]/y[3]));
} /* end of derivs()*/


/*PROPERTY FUNCTION*/
properties(z, densptr)
double z, *densptr;
{
double press, temp;
int i;

/* PROPERTIES FOR HIGH GRADIENT PROFILE */
/* z, densa, dens, Ca, D, visc, P, T */
static double densA[17][8] = {
{  -1.0, 1025.15, 1.82, 1.74, 1.9e-5, 0.01, 1.00, 19.0 },
{ 100.0, 1025.15, 19.2, 16.2, 1.9e-5, 0.01, 10., 19.0 },
{ 200.0, 1025.41, 41.2, 30.4, 1.9e-5, 0.01, 20., 18.0 },
{ 300.0, 1025.68, 67.5, 42.1, 1.9e-5, 0.01, 30., 17.0 },
{ 400.0, 1025.94, 101.3, 50.6, 1.9e-5, 0.01, 40., 16.0 },
{ 500.0, 1026.2 , 160.0, 54.4, 1.9e-5, 0.01, 50.0, 15.0 },
{ 500.0, 1026.2 , 825.0, 55.4, 1.9e-5, 0.01, 50.0, 15.0 }, 
{ 600.0, 1026.46, 853.1, 60.0, 1.9e-5, 0.01, 60., 12.8 },
{ 700.0, 1026.72, 885.2, 60.0, 1.9e-5, 0.01, 70., 10.6 },
{ 800.0, 1026.98, 912.0, 60.0, 1.9e-5, 0.01, 80.0, 8.4 },
{ 900.0, 1027.24, 934.6, 60.0, 1.9e-5, 0.01, 90, 6.2 },
{1000.0, 1027.5 , 953.5, 60.0, 1.9e-5, 0.01, 100., 4.0},
{1200.0, 1027.52, 966.7, 60.0, 1.9e-5, 0.01, 120., 3.6},
{1400.0, 1027.54, 981.0, 60.0, 1.9e-5, 0.01, 140., 3.2},
{1600.0, 1027.56, 992.6, 60.0, 1.9e-5, 0.01, 160., 2.8},
{2000.0, 1027.6, 1013.2, 60.0, 1.9e-5, 0.01, 200., 2.0},
{3000.0, 1027.7, 1050.3, 60.0, 1.9e-5, 0.01, 300., 1.5} 
};

i=0;
while (z > densA[i][0])
i++;
/*denssw*/
*densptr = ((densA[i][1]-densA[i-1][1])*(z-densA[i-1][0])/(densA[i][0]-densA[i-1][0]))+densA[i-1][1];
/*dpwdz*/ 
*(densptr+1) = (densA[i][1]-densA[i-1][1])/(densA[i][0]-densA[i-1][0]);
/*dens*/
*(densptr+2) = ((densA[i][2]-densA[i-1][2])*(z-densA[i-1][0])/(densA[i][0]-densA[i-1][0]))+densA[i-1][2];
/*densdot */
*(densptr+3) = (densA[i][2]-densA[i-1][2])/(densA[i][0]-densA[i-1][0]);
/*Cs */
*(densptr+4) = ((densA[i][3]-densA[i-1][3])*(z-densA[i-1][0])/(densA[i][0]-densA[i-1][0]))+densA[i-1][3];
/*D*/
*(densptr+5) = ((densA[i][4]-densA[i-1][4])*(z-densA[i-1][0])/(densA[i][0]-densA[i-1][0]))+densA[i-1][4];
/*visc*/
*(densptr+6) = ((densA[i][5]-densA[i-1][5])*(z-densA[i-1][0])/(densA[i][0]-densA[i-1][0]))+densA[i-1][5];
/*pressure*/
press = ((densA[i][6]-densA[i-1][6])*(z-densA[i-1][0])/(densA[i][0]-densA[i-1][0])) + densA[i-1][6];
/*temperature*/
temp = ((densA[i][7]-densA[i-1][7])*(z-densA[i-1][0])/(densA[i][0]-densA[i-1][0])) + densA[i-1][7];
} /*END OF PROPERTY FUNCTION*/

extern float **y, *x;
float **y=0, *xx=0;

/********* integrator driver *******************/


void rkdrive(vstart, nvar, xl, x2, nstep, derivs, alpha, laml, lam2)
int nvar, nstep;
float vstart[], xl, x2, alpha, laml, lam2;
void (*derivs)();
{
int i, k;
float x,h;
double ro,zo;
float *v, *vout, *dv, *vector();
void rk4(), error(), free_vector();

v=vector(1,nvar);
vout = vector(1, nvar);
dv = vector(1, nvar);


ro = vstart[1];
zo = xl;


for (i=1; i<=nvar; i++) {
v[i]=vstart[i];
y[i][1]=v[i];
}
xx[1]=xl;
x=xl;
h=(x2-xl)/nstep;
for (k=1; k<=nstep; k++){
(*derivs)(x,v,dv,ro,zo, alpha, lam1, lam2);
rk4(v, dv, nvar, x, h, vout, derivs, ro, zo, alpha, laml, lam2);
if (x+h==x) nrerror("Step size too small in routine RKDRIVE");
x+=h;
xx[k+1]=x;
for (i=1; i<=nvar; i++) {
v[i] = vout[i];
y[i][k+1] = v[i];
}
}
free_vector(dv,1,nvar);
free_vector(vout,1,nvar);
free_vector(v,1,nvar);
} /** end of RKDRIVE **/




/* Fourth Order Runge-Kutta Numerical Integrator*/
void rk4(y,dydx,n,x,h,yout,derivs,ro,zo, alpha, laml, lam2)
float y[], dydx[], x, h, yout[];
void (*derivs)();
double ro, zo, alpha, laml, lam2;
{
int i, Bflag, rflag;
float xh, hh, h6, *dym, *dyt, *yt, *vector();


dym = vector(1,n);
dyt = vector(1,n);
yt = vector(1,n);
hh = h*0.5;
h6 = h/6.0;
xh = x + hh;
for (i=1;i<=n;i++) yt[i]=y[i]+hh*dydx[i];
(*derivs)(xh,yt,dyt,ro,zo,alpha,laml,lam2,&Bflag,&rflag);
for (i=1;i<=n;i++) yt[i]=y[i]+hh*dyt[i];
(*derivs)(xh, yt, dym, ro, zo, alpha, laml, lam2, &Bflag, &rflag);
for (i=1; i<=n;i++) {
yt[i]=y[i]+h*dym[i];
dym[i] += dyt[i];
}
(*derivs)(x+h, yt, dyt, ro, zo, alpha, laml, lam2, &Bflag, &rflag);
for (i=1; i<=n; i++)
yout[i] = y[i]+h6*(dydx[i]+dyt[i]+2.0*dym[i]);

/* limiters and modifiers */
if(rflag == 1)
yout[1]=0.0;

/* peeling model */
if((Bflag==1)&&(rflag!=1)) {
yout[2]=0.5*yout[2];
yout[3]=0.707*yout[3]; /*b*/
/*yout[4] = yout[4]; /*U*/
}

if((500.0-h >= x)&&(x>500.0)) { /*** phase change to vapor ***/
y[1]=1.728*y[1];
yout[1] = y[1];
}
free_vector(yt,1 ,n);
free_vector(dyt,1,n);
free_vector(dym,l,n);
} /** end of RK4 **/


main() /* Beginning of main */
{
FILE *fpwz, *fopen();
int j;
float x1,x2,*vstart,alpha,lam1,lam2;
vstart=vector(l, NVAR);
xx = vector(l, NSTEP+1);
y =matrix(l,NVAR,1,NSTEP+1);
    printf("Input initial depth [m ]:");
 scanf("%f",&xl);
x2=0.0;
printf("\nInput ro [cm] : ");
scanf("%f",&vstart[1]);
vstart[2] =0.0; /*delta pw*/
printf("\nInput bo [m] : ");
scanf("%f",&vstart[3]);
printf("\nlnput umo [m/s] : ");
scanf("%f",&vstart[4]);
printf("\nInput entrainment coefficient, alpha: ");
scanf("%f",&alpha);
printf("\nlnput C02 volume spreading ratio, laml: ");
scanf("%f",&lam1);
printf("\nlnput water density spreading ratio, lam2: ");
scanf("%f",&lam2);

/* call integrator */
rkdrive(vstart, NVAR, xl, x2, NSTEP, derivs, alpha, lam1, lam2);
/* output */
fpwz = fopen("plumeout", "w");
for(j=1;j<=NSTEP+1; j++)
{
/*y[1][i]=r y[2][il=deltap y[3][i]=b y[4][i]= Um*/
if(y[4][j] > 0) { 
     printf(fpwz, "%f %f %f %f %f \n",xx[j], y[1][j], y[2][j], y[3][j], y[4][j]);
   }
}
printf("\n Output in file: plumeout. \n");

free_matrix(y,1,NVAR,1,NSTEP+1);
free_vector(xx,1,NSTEP+1);
free_vector(vstart, 1,NVAR);
fclose(fpwz);
} /*end of main*/