printf("(sound_velocity) INVALID iphi value (=%16.8e) at (%d,%d)\n",iphi(phi[i][j]),i,j);
		exit(-1);
	}
	c2=(1.-iphi(phi[i][j]))*cair2 + iphi(phi[i][j])*cwater2;
	if(c2<=0.){
		printf("Negative Sound Speed =%16.8e at (%2d,%2d), iphi=%16.8e, cair2=%16.8e, cwater2=%16.8e\n",c2, i,j,iphi(phi[i][j]),cair2, cwater2);
		exit(-1);
	}
	return(c2);
}
int bc_p(double pn[IMAX][JMAX]){
	int i,j;
	/********** B.C. for Pressure & Sound speed *****************************/
	/* LEFT side */
	/*... outlet ...*/
	for(j=2; j<=OUTLET+1; j++){
	pn[0 ][j]=pn[2 ][j];
	pn[1 ][j]=pn[2 ][j];
	}
	/*... wall ...*/
	for(j=OUTLET+2; j<=JMAX-3; j++){
	pn[0 ][j]=pn[3 ][j];
	pn[1 ][j]=pn[2 ][j];
	}
	/* RIGHT side */
	/*... inlet&open boundary ...*/
	for(j=2; j<=JMAX-3; j++){
	pn[IMAX-1][j]=pn[IMAX-3][j];
	pn[IMAX-2][j]=pn[IMAX-3][j];
	}
	/* BOTTM */
	for(i=2; i<=IMAX-3; i++){
	/* bottom */
	pn[i][0]=pn[i][3];
	pn[i][1]=pn[i][2];
	}
	/* TOP */
	for(i=0; i<=IMAX-1; i++){
	pn[i][JMAX-1]=Patm;
	pn[i][JMAX-2]=Patm;
	}
	/* corner */
	/* left & bottom */
	pn[0][0]=pn[3][3];
	pn[0][1]=pn[3][2];
	pn[1][0]=pn[2][3];
	pn[1][1]=pn[2][2];
	/* right & bottom */
	pn[IMAX-1][0]=pn[IMAX-4][3];
	pn[IMAX-1][1]=pn[IMAX-4][2];
	pn[IMAX-2][0]=pn[IMAX-3][3];
	pn[IMAX-2][1]=pn[IMAX-3][2];
	/********** B.C. for Pressure *****************************/
	return(0);
}


int psor(double rho[IMAX][JMAX], double u[IMAX][JMAX], double v[IMAX][JMAX], double phi[IMAX][JMAX], double p[IMAX][JMAX], double pn[IMAX][JMAX], int rtt){
	int i, j;
	int ite=0;
	double oldp[IMAX][JMAX], op;
	double res=10., rhs=0., bnm=0.,npt;
	double fdx, fdy, fdx2, fdy2, DT2, fDT, fDT2;
	double b0, b1, b2, b3, b4, b;
	double div;
	double rhol, rhor, rhou, rhod;
	double dx2, dy2, cs2, rhoav;
	for(i=0; i<IMAX; i++){
		for(j=0; j<JMAX; j++){
			oldp[i][j]=p[i][j];
		}
	}
	dx2=DX*DX;
	dy2=DY*DY;
	fdx=1./DX;
	if(DX<=0.){
		printf("invalid fdx(=%16.8e). aborted! in psor()\n",fdx);
		exit(-1);
	}
	fdy=1./DY;
	if(DY<=0.){
		printf("invalid fdy(=%16.8e). aborted! in psor()\n",fdy);
		exit(-1);
	}
	fdx2=1./dx2;
	if(pow(DX,2)<=0.){
		printf("invalid fdx2(=%16.8e). aborted! in psor()\n",fdx2);
		exit(-1);
	}
	fdy2=1./dy2;
	if(pow(DY,2)<=0.){
		printf("invalid fdy2(=%16.8e). aborted! in psor()\n",fdy2);
		exit(-1);
	}
	DT2=pow(DT,2);
	if(DT2<=0.){
		printf("invalid DT2(=%16.8e). aborted! in psor()\n",DT2);
		exit(-1);
	}
	fDT=1./DT;
	if(DT<=0.){
		printf("invalid fDT(=%16.8e). aborted! in psor()\n",fDT);
		exit(-1);
	}
	fDT2=1./DT2;
	if(fDT2<=0.){
		printf("invalid DT2(=%16.8e). aborted! in psor()\n",DT2);
		exit(-1);
	}
	while ( (res>EPS) && (ite<MITE)){
		for(i=0; i<IMAX; i++){
			for(j=0; j<JMAX; j++){
				p[i][j]=ALPHA*pn[i][j]+(1.-ALPHA)*p[i][j];
				/*
				sample
				pp p
				p pn
				yp oldp
				*/
			}
		}
		ite++;
		res=0.;
		for(j=2; j<=JMAX-3; j++){
			for(i=2; i<=IMAX-3; i++){
				cs2=sound_velocity(i,j,rho,oldp,phi);
				rhol=(rho[i][j]+rho[i-1][j ])*0.5;
				rhor=(rho[i][j]+rho[i+1][j ])*0.5;
				rhod=(rho[i][j]+rho[i ][j-1])*0.5;
				rhou=(rho[i][j]+rho[i ][j+1])*0.5;
				rhoav=0.25*(rhol+rhor+rhou+rhod);
				if(rhol<=0.){
					printf("Negative Density at [%d][%d] rhol=%30.22e\n",i,j,rhol);
					exit(-1);
				}
				if(rhor<=0.){
				printf("Negative Density at [%d][%d] rhor=%30.22e\n",i,j,rhor);
				exit(-1);
				}
				if(rhod<=0.){
				printf("Negative Density at [%d][%d] rhod=%30.22e\n",i,j,rhod);
				exit(-1);
				}
				if(rhou<=0.){
				printf("Negative Density at [%d][%d] rhou=%30.22e rho=%30.22e rho(j+1)=%30.22e\n",i,j,rhou,rho[i][j],rho[i][j+1]);
				exit(-1);
				}
				/*
				sample
				pp p
				p pn
				yp oldp
				*/
				/* original sample SOR ( not included sound speed ) */
				// div=1.0+GAMMA*oldp[i][j]*DT*DT*((1.0/rhor+1.0/rhol)/dx2+(1.0/rhou+1.0/rhod)/dy2);
				// pn[i][j]=(oldp[i][j]-GAMMA*oldp[i][j]*DT*((u[i+1][j]-u[i][j])/DX - (p[i+1][j]/rhor+p[i-1][j]/rhol)*DT/dx2
				// +(v[i][j+1]-v[i][j])/DY - (p[i][j+1]/rhou + p[i][j-1]/rhod)*DT/dy2))/div;
				/* (improved) sample SOR, including sound speed */
				div=1.0+cs2*rho[i][j]*DT*DT*((1.0/rhor+1.0/rhol)/dx2+(1.0/rhou+1.0/rhod)/dy2);
				pn[i][j]=(oldp[i][j]-cs2*rho[i][j]*DT*((u[i+1][j]-u[i][j])/DX - (p[i+1][j]/rhor+p[i-1][j]/rhol)*DT/dx2+(v[i][j+1]-v[i][j])/DY - (p[i][j+1]/rhou + p[i][j-1]/rhod)*DT/dy2))/div;
				if(pn[i][j]<=0.){
					printf("Negative Pressure at [%d][%d] p=%30.22e\n",i,j,pn[i][j]);
					return(-1);
				}
				res+=pow(p[i][j]-pn[i][j],2);
			}//for.i
		}//for.j
		if(ite%200==0)   printf("--(SOR)%4d st.ite., res.=%16.8e (%6.1f t.l.t. EPS)...(SOR)\n",ite,res,(res/EPS));
		/* B.C. */
		bc_p(pn);
	/* ... end of pressure iteration loop ... */
	}// while
	if(ite==MITE)   printf("--not converged within %4d steps iterations, residual=%16.8e....(SOR)\n",MITE,res);
	else            printf("--converged in %4d steps iterations. res=%16.8e...(SOR)\n",ite,res);
	for(i=0; i<IMAX; i++){
		for(j=0; j<JMAX; j++){
			p[i][j]=oldp[i][j];
		}
	}
	return(0);
}


int isoentro_rho(double rho[IMAX][JMAX], double rhon[IMAX][JMAX],
double p[IMAX][JMAX], double pn[IMAX][JMAX],
double phi[IMAX][JMAX],
double un[IMAX][JMAX], double vn[IMAX][JMAX]){
int i,j;
double cs2;
double pnl, pnr, pnd, pnu, pna;
double pl, pr, pd, pu, pa;
// --------------------------------------------------------------------------------
printf("(isoentro_rho) this function is NOT CONSIDERED to use....\n");
exit(-1);
// --------------------------------------------------------------------------------
for(i=0; i<IMAX; i++){
for(j=0; j<JMAX; j++){
// ----- shimizu2001 ---------------------------------------------------------------------
// cs2=sound_velocity(i, j, rho, p, phi);
// rhon[i][j]=rho[i][j]+1./cs2*(pn[i][j]-p[i][j]);
// ---------------------------------------------------------------------------------------
// ----- mizutani2002 --------------------------------------------------------------------
rhon[i][j]=rho[i][j]-rho[i][j]*((un[i+1][j]-un[i][j])/DX+(vn[i][j+1]-vn[i][j])/DY)*DT;
if(rhon[i][j]<=RHOair)rhon[i][j]=RHOair;
if(rhon[i][j]>=RHOwater)rhon[i][j]=RHOwater;
// ---------------------------------------------------------------------------------------
}
}
return(0);
}
int rho_from_phi(double rhon[IMAX][JMAX], double phin[IMAX][JMAX]){
int i,j;
double tmpphi[IMAX][JMAX];
for(i=0; i<IMAX; i++){
for(j=0; j<JMAX; j++){
rhon[i][j]=iphi(phin[i][j])*RHOwater + (1.-iphi(phin[i][j]))*RHOair;
if(rhon[i][j]<=0.){
printf("(rho_from_phi) Negative Density(=%8.4f) at (%d,%d), phin=%8.4f\n",
rhon[i][j], i, j, phin[i][j]);
exit(-1);
}
}
}
return(0);
}
int rho_from_phi_filtered(double rhon[IMAX][JMAX], double phin[IMAX][JMAX]){
int i,j, mgc=0, agc=0, wgc=0;
double tmpphi[IMAX][JMAX];
for(i=0; i<IMAX; i++){
for(j=0; j<JMAX; j++){
tmpphi[i][j]=0.2*(iphi(phin[i][j])+iphi(phin[i][j])+iphi(phin[i][j+1])+iphi(phin[i+1][j])+iphi(phin[i+1][j+1]));
}
}
for(i=0; i<IMAX; i++){
for(j=0; j<JMAX; j++){
if(tmpphi[i][j]<=0.2){
phin[i][j]=dphi(0.);
rhon[i][j]=RHOair;
agc++;
}
else if(tmpphi[i][j]>=0.8){
phin[i][j]=dphi(1.);
rhon[i][j]=RHOwater;
wgc++;
}
else{
rhon[i][j]=iphi(phin[i][j])*RHOwater + (1.-iphi(phin[i][j]))*RHOair;
phin[i][j]=dphi(0.5);
mgc++;
}
if(rhon[i][j]<=0.){
printf("(rho_from_phi) Negative Density(=%8.4f) at (%d,%d), phin=%8.4f\n",
rhon[i][j], i, j, phin[i][j]);
exit(-1);
}
}
}
printf("(rho_from_phi) Air=%d Water=%d Mid=%d\n",agc,wgc,mgc);
return(0);
}
int correct_phi(double phin[IMAX][JMAX]){
	/* original */
	int i,j,nom=0;
	double mixture[IMAX][JMAX], vlm=0., tvlm=0., difvlm, bi;
	for(i=0; i<IMAX; i++){
	for(j=0; j<JMAX; j++){
	if(iphi(phin[i][j])<0.){/* air cell */
	phin[i][j]=dphi(0.);
	}
	else if(iphi(phin[i][j])>1.){/* water cell */
	phin[i][j]=dphi(1.);
	}
	if(i>=2 && i<IMAX-1 && j>=2 && j<JMAX-1) vlm+=iphi(phin[i][j]); /* total volume */
	}
	}
	return(0);
}
int correct_phi_isoentro(double phin[IMAX][JMAX]){
int i,j;
int na=0, nf=0, nm=0;
for(i=0; i<IMAX; i++){
for(j=0; j<JMAX; j++){
if(iphi(phin[i][j])<=0.){
phin[i][j]=dphi(0.);
na++;
}
else if(iphi(phin[i][j])>=1.){
phin[i][j]=dphi(1.);
nf++;
}
else{
phin[i][j]=dphi(0.5);
nm++;
}
}
}
printf("Air:%3d(%8.4f%%) Water:%3d(%8.4f%%) Mix:%3d(%8.4f%%)\n",
na,(double)na/(double)(IMAX*JMAX)*100., nf,(double)nf/(double)(IMAX*JMAX)*100, nm,(double)nm/(double)(IMAX*JMAX)*100.);
return(0);
}
int correct_phi_adv(double phin[IMAX][JMAX]){
int i,j;
int na=0, nf=0, nm=0;
for(i=0; i<IMAX; i++){
for(j=0; j<JMAX; j++){
if(iphi(phin[i][j])<0.5){
phin[i][j]=dphi(0.);
na++;
}
else if(iphi(phin[i][j])>0.5){
phin[i][j]=dphi(1.);
nf++;
}
else{
phin[i][j]=dphi(0.5);
nm++;
}
}
}
printf("Air:%3d(%8.4f%%) Water:%3d(%8.4f%%) Mix:%3d(%8.4f%%)\n",
na,(double)na/(double)(IMAX*JMAX)*100., nf,(double)nf/(double)(IMAX*JMAX)*100, nm,(double)nm/(double)(IMAX*JMAX)*100.);
return(0);
}
int isoentro_phi(double phi[IMAX][JMAX], double phin[IMAX][JMAX], double rho[IMAX][JMAX],
double p[IMAX][JMAX], double pn[IMAX][JMAX],
double un[IMAX][JMAX], double vn[IMAX][JMAX]){
int i,j;
double cs2;
for(i=0; i<IMAX; i++){
for(j=0; j<JMAX; j++){
// ----- shimizu2001 ---------------------------------------------------------------------
// cs2=sound_velocity(i, j, rho, p, phi);
// phin[i][j]=phi[i][j]+phi[i][j]/(rho[i][j]*cs2)*(pn[i][j]-p[i][j]);
// ----- shimizu2001 ---------------------------------------------------------------------
// ----- mizutani2002 ---------------------------------------------------------------------
phin[i][j]=phi[i][j]-DT* phi[i][j]*((un[i+1][j]-un[i][j])/DX+ (vn[i][j+1]-vn[i][j])/DY);
// ----- mizutani2002 ---------------------------------------------------------------------
}
}
correct_phi(phin);
return(0);
}
int isoentro_update(double u[IMAX][JMAX], double un[IMAX][JMAX],
double v[IMAX][JMAX], double vn[IMAX][JMAX],
double rho[IMAX][JMAX], double rhon[IMAX][JMAX],
double phi[IMAX][JMAX], double phin[IMAX][JMAX]){
int i,j;
for(i=0; i<IMAX; i++){
for(j=0; j<JMAX; j++){
u[i][j]= un[i][j];
v[i][j]= vn[i][j];
rho[i][j]=rhon[i][j];
phi[i][j]=phin[i][j];
}
}
return(0);
}
int diffusion_update(double u[IMAX][JMAX], double un[IMAX][JMAX],
double v[IMAX][JMAX], double vn[IMAX][JMAX]){
int i,j;
for(i=0; i<IMAX; i++){
for(j=0; j<JMAX; j++){
u[i][j]=un[i][j];
v[i][j]=vn[i][j];
}
}
return(0);
}
int sc_update(double scn[IMAX][JMAX], double sc[IMAX][JMAX]){
int i,j;
for(i=0; i<IMAX; i++){
for(j=0; j<JMAX; j++){
sc[i][j]=scn[i][j];
}
}
return(0);
}
//---------- Boundary Condition ----------
int bc_scalar(double sc[IMAX][JMAX], double dvalue){
int i,j;
/***** LEFT SIDE *****/
/*... outlet ...*/
for(j=2; j<=OUTLET+1; j++){
/* left */
sc[0 ][j]=sc[2 ][j];
sc[1 ][j]=sc[2 ][j];
}
/*... wall ...*/
for(j=OUTLET+2; j<=JMAX-1; j++){
/* left */
sc[0 ][j]=sc[3 ][j];
sc[1 ][j]=sc[2 ][j];
}
/***** RIGHT SIDE *****/
/*... inlet ...*/
for(j=2; j<=INLET+1; j++){
sc[IMAX-1][j]=dvalue;
sc[IMAX-2][j]=dvalue;
}
/*... openboundary on inlet ...*/
for(j=INLET+2; j<=JMAX-1; j++){
sc[IMAX-1][j]=sc[IMAX-3][j];
sc[IMAX-2][j]=sc[IMAX-3][j];
}
/***** BOTTOM *****/
for(i=2; i<=IMAX-3; i++){
/* bottom */
sc[i][0]=sc[i][3];
sc[i][1]=sc[i][2];
}
/***** TOP (open boundary) *****/
for(i=2; i<=IMAX-3; i++){
/* top.. wall b.c. */
sc[i][JMAX-1]=sc[i][JMAX-3];
sc[i][JMAX-2]=sc[i][JMAX-3];
}
/* corner */
/* left & bottom
| | | | |
+---+---+---+---+--
3 | | | B | A |
+---+---+---+---+--
2 | | | D | C |
+---+---+---+---+--
1 | C | D | | |
+---+---+---+---+--
0 | A | B | | |