00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 #include <math.h>
00031 #include <stdio.h>
00032 #include <stdlib.h>
00033 #include "defs.h"
00034 #include "athena.h"
00035 #include "globals.h"
00036 #include "prototypes.h"
00037
00038
00039 static ConsS ***RootSoln=NULL;
00040 static int wave_flag;
00041
00042
00043
00044
00045 void problem(DomainS *pDomain)
00046 {
00047 GridS *pGrid=(pDomain->Grid);
00048 ConsS ***Soln;
00049 int i=0,j=0,k=0;
00050 int is,ie,js,je,ks,ke,n,m,nx1,nx2,nx3,Nx1,Nx2;
00051 Real amp,vflow,angle;
00052 Real d0,p0,u0,v0,w0,h0;
00053 Real x1,x2,x3,r,ev[NWAVE],rem[NWAVE][NWAVE],lem[NWAVE][NWAVE];
00054 Real x1size,x2size,sin_a,cos_a,lambda;
00055 #ifdef MHD
00056 Real bx0,by0,bz0,xfact,yfact;
00057 Real **az;
00058 #endif
00059 is = pGrid->is; ie = pGrid->ie;
00060 js = pGrid->js; je = pGrid->je;
00061 ks = pGrid->ks; ke = pGrid->ke;
00062 nx1 = (ie-is)+1 + 2*nghost;
00063 nx2 = (je-js)+1 + 2*nghost;
00064 nx3 = (ke-ks)+1 + 2*nghost;
00065
00066
00067
00068 Nx1 = pDomain->Nx[0];
00069 Nx2 = pDomain->Nx[1];
00070 if (Nx1 == 1 || Nx2 == 1) {
00071 ath_error("[problem]: this test only works with Nx1 & Nx2 > 1\n");
00072 }
00073
00074
00075
00076 #ifdef MHD
00077 if ((az = (Real**)calloc_2d_array(nx2, nx1, sizeof(Real))) == NULL)
00078 ath_error("[problem]: Error allocating memory for \"az\"\n");
00079 #endif
00080
00081 if ((Soln = (ConsS***)calloc_3d_array(nx3,nx2,nx1,sizeof(ConsS)))==NULL)
00082 ath_error("[problem]: Error allocating memory\n");
00083 if (pDomain->Level == 0){
00084 if ((RootSoln = (ConsS***)calloc_3d_array(nx3,nx2,nx1,sizeof(ConsS)))
00085 == NULL) ath_error("[problem]: Error alloc memory for RootSoln\n");
00086 }
00087
00088
00089
00090 wave_flag = par_geti("problem","wave_flag");
00091 amp = par_getd("problem","amp");
00092 vflow = par_getd("problem","vflow");
00093
00094
00095
00096 x1size = pDomain->RootMaxX[0] - pDomain->RootMinX[0];
00097 x2size = pDomain->RootMaxX[1] - pDomain->RootMinX[1];
00098 angle = atan(x1size/x2size);
00099 sin_a = sin(angle);
00100 cos_a = cos(angle);
00101
00102
00103 if (cos_a >= sin_a) {
00104 lambda = x1size*cos_a;
00105 } else {
00106 lambda = x2size*sin_a;
00107 }
00108
00109
00110
00111
00112 d0 = 1.0;
00113 #ifndef ISOTHERMAL
00114 p0 = 1.0/Gamma;
00115 u0 = vflow*sqrt(Gamma*p0/d0);
00116 #else
00117 u0 = vflow*Iso_csound;
00118 #endif
00119 v0 = 0.0;
00120 w0 = 0.0;
00121 #ifdef MHD
00122 bx0 = 1.0;
00123 by0 = sqrt(2.0);
00124 bz0 = 0.5;
00125 xfact = 0.0;
00126 yfact = 1.0;
00127 #endif
00128
00129 for (n=0; n<NWAVE; n++) {
00130 for (m=0; m<NWAVE; m++) {
00131 rem[n][m] = 0.0;
00132 lem[n][m] = 0.0;
00133 }
00134 }
00135
00136 #ifdef HYDRO
00137 #ifdef ISOTHERMAL
00138 esys_roe_iso_hyd(u0,v0,w0, ev,rem,lem);
00139 #else
00140 h0 = ((p0/Gamma_1 + 0.5*d0*(u0*u0+v0*v0+w0*w0)) + p0)/d0;
00141 esys_roe_adb_hyd(u0,v0,w0,h0,ev,rem,lem);
00142 ath_pout(0,"Ux - Cs = %e, %e\n",ev[0],rem[0][wave_flag]);
00143 ath_pout(0,"Ux = %e, %e\n",ev[1],rem[1][wave_flag]);
00144 ath_pout(0,"Ux + Cs = %e, %e\n",ev[4],rem[4][wave_flag]);
00145 #endif
00146 #endif
00147
00148 #ifdef MHD
00149 #if defined(ISOTHERMAL)
00150 esys_roe_iso_mhd(d0,u0,v0,w0,bx0,by0,bz0,xfact,yfact,ev,rem,lem);
00151 ath_pout(0,"Ux - Cf = %e, %e\n",ev[0],rem[0][wave_flag]);
00152 ath_pout(0,"Ux - Ca = %e, %e\n",ev[1],rem[1][wave_flag]);
00153 ath_pout(0,"Ux - Cs = %e, %e\n",ev[2],rem[2][wave_flag]);
00154 ath_pout(0,"Ux + Cs = %e, %e\n",ev[3],rem[3][wave_flag]);
00155 ath_pout(0,"Ux + Ca = %e, %e\n",ev[4],rem[4][wave_flag]);
00156 ath_pout(0,"Ux + Cf = %e, %e\n",ev[5],rem[5][wave_flag]);
00157 #else
00158 h0 = ((p0/Gamma_1+0.5*(bx0*bx0+by0*by0+bz0*bz0)+0.5*d0*(u0*u0+v0*v0+w0*w0))
00159 + (p0+0.5*(bx0*bx0+by0*by0+bz0*bz0)))/d0;
00160 esys_roe_adb_mhd(d0,u0,v0,w0,h0,bx0,by0,bz0,xfact,yfact,ev,rem,lem);
00161 ath_pout(0,"Ux - Cf = %e, %e\n",ev[0],rem[0][wave_flag]);
00162 ath_pout(0,"Ux - Ca = %e, %e\n",ev[1],rem[1][wave_flag]);
00163 ath_pout(0,"Ux - Cs = %e, %e\n",ev[2],rem[2][wave_flag]);
00164 ath_pout(0,"Ux = %e, %e\n",ev[3],rem[3][wave_flag]);
00165 ath_pout(0,"Ux + Cs = %e, %e\n",ev[4],rem[4][wave_flag]);
00166 ath_pout(0,"Ux + Ca = %e, %e\n",ev[5],rem[5][wave_flag]);
00167 ath_pout(0,"Ux + Cf = %e, %e\n",ev[6],rem[6][wave_flag]);
00168 #endif
00169 #endif
00170
00171
00172
00173
00174 #ifdef MHD
00175 for (j=js; j<=je+1; j++) {
00176 for (i=is; i<=ie+1; i++) {
00177 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00178 x1 = x1 - 0.5*pGrid->dx1;
00179 x2 = x2 - 0.5*pGrid->dx2;
00180 r = (x1*cos_a + x2*sin_a)/lambda;
00181
00182 az[j][i] = -bx0*(x1*sin_a - x2*cos_a)
00183 - by0*(x1*cos_a + x2*sin_a)
00184 + amp*lambda*cos(2.0*PI*r)*rem[NWAVE-2][wave_flag]/(2.0*PI);
00185 }
00186 }
00187 #endif
00188
00189
00190
00191
00192
00193 for (k=ks; k<=ke; k++) {
00194 for (j=js; j<=je; j++) {
00195 for (i=is; i<=ie; i++) {
00196 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00197 r = (x1*cos_a + x2*sin_a)/lambda;
00198
00199 Soln[k][j][i].d = d0 + amp*sin(2.0*PI*r)*rem[0][wave_flag];
00200
00201 #ifndef ISOTHERMAL
00202 #ifdef HYDRO
00203 Soln[k][j][i].E = p0/Gamma_1 + 0.5*d0*u0*u0
00204 + amp*sin(2.0*PI*r)*rem[4][wave_flag];
00205 #else
00206 Soln[k][j][i].E = p0/Gamma_1 + 0.5*d0*u0*u0 + 0.5*(bx0*bx0+by0*by0+bz0*bz0)
00207 + amp*sin(2.0*PI*r)*rem[4][wave_flag];
00208 #endif
00209 #endif
00210
00211 Soln[k][j][i].M1 = d0*vflow*cos_a
00212 + amp*sin(2.0*PI*r)*rem[1][wave_flag]*cos_a
00213 - amp*sin(2.0*PI*r)*rem[2][wave_flag]*sin_a;
00214 Soln[k][j][i].M2 = d0*vflow*sin_a
00215 + amp*sin(2.0*PI*r)*rem[1][wave_flag]*sin_a
00216 + amp*sin(2.0*PI*r)*rem[2][wave_flag]*cos_a;
00217 Soln[k][j][i].M3 = amp*sin(2.0*PI*r)*rem[3][wave_flag];
00218
00219 #ifdef MHD
00220 Soln[k][j][i].B3c = bz0 + amp*sin(2.0*PI*r)*rem[NWAVE-1][wave_flag];
00221 #endif
00222
00223 #if (NSCALARS > 0)
00224 for (n=0; n<NSCALARS; n++)
00225 Soln[k][j][i].s[n] = amp*(1.0 + sin(2.0*PI*r));
00226 #endif
00227
00228 }}}
00229
00230 #ifdef MHD
00231
00232
00233 for (k=ks; k<=ke; k++) {
00234 for (j=js; j<=je; j++) {
00235 for (i=is; i<=ie+1; i++) {
00236 pGrid->B1i[k][j][i] = (az[j+1][i] - az[j][i])/pGrid->dx2;
00237 }}}
00238 for (k=ks; k<=ke; k++) {
00239 for (j=js; j<=je+1; j++) {
00240 for (i=is; i<=ie; i++) {
00241 pGrid->B2i[k][j][i] =-(az[j][i+1] - az[j][i])/pGrid->dx1;
00242 }}}
00243 for (k=ks; k<=ke; k++) {
00244 for (j=js; j<=je; j++) {
00245 for (i=is; i<=ie; i++) {
00246 pGrid->B3i[k][j][i] = Soln[k][j][i].B3c;
00247 }}}
00248 if (pGrid->Nx[2] > 1) {
00249 for (j=js; j<=je; j++) {
00250 for (i=is; i<=ie; i++) {
00251 pGrid->B3i[ke+1][j][i] = Soln[ke][j][i].B3c;
00252 }}
00253 }
00254
00255
00256
00257 for (k=ks; k<=ke; k++) {
00258 for (j=js; j<=je; j++) {
00259 for (i=is; i<=ie; i++) {
00260 Soln[k][j][i].B1c = 0.5*(pGrid->B1i[k][j][i] + pGrid->B1i[k][j][i+1]);
00261 Soln[k][j][i].B2c = 0.5*(pGrid->B2i[k][j][i] + pGrid->B2i[k][j+1][i]);
00262 }}}
00263 free_2d_array(az);
00264 #endif
00265
00266
00267
00268 for (k=ks; k<=ke; k++) {
00269 for (j=js; j<=je; j++) {
00270 for (i=is; i<=ie; i++) {
00271 pGrid->U[k][j][i].d = Soln[k][j][i].d;
00272 #ifndef ISOTHERMAL
00273 pGrid->U[k][j][i].E = Soln[k][j][i].E;
00274 #endif
00275 pGrid->U[k][j][i].M1 = Soln[k][j][i].M1;
00276 pGrid->U[k][j][i].M2 = Soln[k][j][i].M2;
00277 pGrid->U[k][j][i].M3 = Soln[k][j][i].M3;
00278 #ifdef MHD
00279 pGrid->U[k][j][i].B1c = Soln[k][j][i].B1c;
00280 pGrid->U[k][j][i].B2c = Soln[k][j][i].B2c;
00281 pGrid->U[k][j][i].B3c = Soln[k][j][i].B3c;
00282 #endif
00283 #if (NSCALARS > 0)
00284 for (n=0; n<NSCALARS; n++)
00285 pGrid->U[k][j][i].s[n] = Soln[k][j][i].s[n];
00286 #endif
00287 }}}
00288
00289
00290
00291 #ifdef SELF_GRAVITY
00292 four_pi_G = par_getd("problem","four_pi_G");
00293 grav_mean_rho = d0;
00294 #endif
00295
00296
00297
00298 #ifdef RESISTIVITY
00299 eta_Ohm = par_getd("problem","eta_O");
00300 Q_AD = par_getd_def("problem","Q_AD",0.0);
00301 d_ind = par_getd_def("problem","d_ind",0.0);
00302 #endif
00303 #ifdef NAVIER_STOKES
00304 nu_V = par_getd("problem","nu");
00305 #endif
00306
00307
00308
00309 if (pDomain->Level == 0) {
00310 for (k=ks; k<=ke; k++) {
00311 for (j=js; j<=je; j++) {
00312 for (i=is; i<=ie; i++) {
00313 RootSoln[k][j][i].d = Soln[k][j][i].d ;
00314 RootSoln[k][j][i].M1 = Soln[k][j][i].M1;
00315 RootSoln[k][j][i].M2 = Soln[k][j][i].M2;
00316 RootSoln[k][j][i].M3 = Soln[k][j][i].M3;
00317 #ifndef ISOTHERMAL
00318 RootSoln[k][j][i].E = Soln[k][j][i].E ;
00319 #endif
00320 #ifdef MHD
00321 RootSoln[k][j][i].B1c = Soln[k][j][i].B1c;
00322 RootSoln[k][j][i].B2c = Soln[k][j][i].B2c;
00323 RootSoln[k][j][i].B3c = Soln[k][j][i].B3c;
00324 #endif
00325 #if (NSCALARS > 0)
00326 for (n=0; n<NSCALARS; n++)
00327 RootSoln[k][j][i].s[n] = Soln[k][j][i].s[n];
00328 #endif
00329 }}}
00330 }
00331
00332 free_3d_array(Soln);
00333
00334 return;
00335 }
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349 void problem_write_restart(MeshS *pM, FILE *fp)
00350 {
00351 return;
00352 }
00353
00354 void problem_read_restart(MeshS *pM, FILE *fp)
00355 {
00356 return;
00357 }
00358
00359 #if (NSCALARS > 0)
00360
00361
00362
00363 static Real color(const GridS *pG, const int i, const int j, const int k)
00364 {
00365 return pG->U[k][j][i].s[0]/pG->U[k][j][i].d;
00366 }
00367 #endif
00368
00369 ConsFun_t get_usr_expr(const char *expr)
00370 {
00371 #if (NSCALARS > 0)
00372 if(strcmp(expr,"color")==0) return color;
00373 #endif
00374 return NULL;
00375 }
00376
00377 VOutFun_t get_usr_out_fun(const char *name){
00378 return NULL;
00379 }
00380
00381 #ifdef RESISTIVITY
00382 void get_eta_user(GridS *pG, int i, int j, int k,
00383 Real *eta_O, Real *eta_H, Real *eta_A)
00384 {
00385 *eta_O = 0.0;
00386 *eta_H = 0.0;
00387 *eta_A = 0.0;
00388
00389 return;
00390 }
00391 #endif
00392
00393 void Userwork_in_loop(MeshS *pM)
00394 {
00395 }
00396
00397
00398
00399
00400
00401
00402
00403 void Userwork_after_loop(MeshS *pM)
00404 {
00405 GridS *pGrid;
00406 int i=0,j=0,k=0;
00407 #if (NSCALARS > 0)
00408 int n;
00409 #endif
00410 int is,ie,js,je,ks,ke;
00411 Real rms_error=0.0;
00412 ConsS error,total_error;
00413 FILE *fp;
00414 char *fname;
00415 int Nx1, Nx2, Nx3, count;
00416 #if defined MPI_PARALLEL
00417 double err[8+NSCALARS], tot_err[8+NSCALARS];
00418 int ierr,myID;
00419 #endif
00420
00421 total_error.d = 0.0;
00422 total_error.M1 = 0.0;
00423 total_error.M2 = 0.0;
00424 total_error.M3 = 0.0;
00425 #ifdef MHD
00426 total_error.B1c = 0.0;
00427 total_error.B2c = 0.0;
00428 total_error.B3c = 0.0;
00429 #endif
00430 #ifndef ISOTHERMAL
00431 total_error.E = 0.0;
00432 #endif
00433 #if (NSCALARS > 0)
00434 for (n=0; n<NSCALARS; n++) total_error.s[n] = 0.0;
00435 #endif
00436
00437
00438
00439 pGrid=pM->Domain[0][0].Grid;
00440 if (pGrid == NULL) return;
00441
00442
00443
00444 is = pGrid->is; ie = pGrid->ie;
00445 js = pGrid->js; je = pGrid->je;
00446 ks = pGrid->ks; ke = pGrid->ke;
00447 for (k=ks; k<=ke; k++) {
00448 for (j=js; j<=je; j++) {
00449 error.d = 0.0;
00450 error.M1 = 0.0;
00451 error.M2 = 0.0;
00452 error.M3 = 0.0;
00453 #ifdef MHD
00454 error.B1c = 0.0;
00455 error.B2c = 0.0;
00456 error.B3c = 0.0;
00457 #endif
00458 #ifndef ISOTHERMAL
00459 error.E = 0.0;
00460 #endif
00461 #if (NSCALARS > 0)
00462 for (n=0; n<NSCALARS; n++) error.s[n] = 0.0;
00463 #endif
00464
00465 for (i=is; i<=ie; i++) {
00466 error.d += fabs(pGrid->U[k][j][i].d - RootSoln[k][j][i].d );
00467 error.M1 += fabs(pGrid->U[k][j][i].M1 - RootSoln[k][j][i].M1);
00468 error.M2 += fabs(pGrid->U[k][j][i].M2 - RootSoln[k][j][i].M2);
00469 error.M3 += fabs(pGrid->U[k][j][i].M3 - RootSoln[k][j][i].M3);
00470 #ifdef MHD
00471 error.B1c += fabs(pGrid->U[k][j][i].B1c - RootSoln[k][j][i].B1c);
00472 error.B2c += fabs(pGrid->U[k][j][i].B2c - RootSoln[k][j][i].B2c);
00473 error.B3c += fabs(pGrid->U[k][j][i].B3c - RootSoln[k][j][i].B3c);
00474 #endif
00475 #ifndef ISOTHERMAL
00476 error.E += fabs(pGrid->U[k][j][i].E - RootSoln[k][j][i].E );
00477 #endif
00478 #if (NSCALARS > 0)
00479 for (n=0; n<NSCALARS; n++)
00480 error.s[n] += fabs(pGrid->U[k][j][i].s[n] - RootSoln[k][j][i].s[n]);;
00481 #endif
00482 }
00483
00484 total_error.d += error.d;
00485 total_error.M1 += error.M1;
00486 total_error.M2 += error.M2;
00487 total_error.M3 += error.M3;
00488 #ifdef MHD
00489 total_error.B1c += error.B1c;
00490 total_error.B2c += error.B2c;
00491 total_error.B3c += error.B3c;
00492 #endif
00493 #ifndef ISOTHERMAL
00494 total_error.E += error.E;
00495 #endif
00496 #if (NSCALARS > 0)
00497 for (n=0; n<NSCALARS; n++) total_error.s[n] += error.s[n];
00498 #endif
00499 }}
00500
00501 #ifdef MPI_PARALLEL
00502 Nx1 = pM->Domain[0][0].Nx[0];
00503 Nx2 = pM->Domain[0][0].Nx[1];
00504 Nx3 = pM->Domain[0][0].Nx[2];
00505 #else
00506 Nx1 = ie - is + 1;
00507 Nx2 = je - js + 1;
00508 Nx3 = ke - ks + 1;
00509 #endif
00510 count = Nx1*Nx2*Nx3;
00511
00512 #ifdef MPI_PARALLEL
00513
00514
00515
00516 err[0] = total_error.d;
00517 err[1] = total_error.M1;
00518 err[2] = total_error.M2;
00519 err[3] = total_error.M3;
00520 #ifdef MHD
00521 err[4] = total_error.B1c;
00522 err[5] = total_error.B2c;
00523 err[6] = total_error.B3c;
00524 #endif
00525 #ifndef ISOTHERMAL
00526 err[7] = total_error.E;
00527 #endif
00528 #if (NSCALARS > 0)
00529 for (n=0; n<NSCALARS; n++) err[8+n] = total_error.s[n];
00530 #endif
00531
00532
00533 ierr = MPI_Reduce(err,tot_err,(8+NSCALARS),MPI_DOUBLE,MPI_SUM,0,
00534 pM->Domain[0][0].Comm_Domain);
00535
00536
00537
00538 ierr = MPI_Comm_rank(pM->Domain[0][0].Comm_Domain, &myID);
00539 if(myID == 0){
00540 total_error.d = tot_err[0];
00541 total_error.M1 = tot_err[1];
00542 total_error.M2 = tot_err[2];
00543 total_error.M3 = tot_err[3];
00544 #ifdef MHD
00545 total_error.B1c = tot_err[4];
00546 total_error.B2c = tot_err[5];
00547 total_error.B3c = tot_err[6];
00548 #endif
00549 #ifndef ISOTHERMAL
00550 total_error.E = tot_err[7];
00551 #endif
00552 #if (NSCALARS > 0)
00553 for (n=0; n<NSCALARS; n++) total_error.s[n] = err[8+n];
00554 #endif
00555
00556 }
00557 else return;
00558
00559 #endif
00560
00561
00562
00563 rms_error = SQR(total_error.d) + SQR(total_error.M1) + SQR(total_error.M2)
00564 + SQR(total_error.M3);
00565 #ifdef MHD
00566 rms_error += SQR(total_error.B1c) + SQR(total_error.B2c)
00567 + SQR(total_error.B3c);
00568 #endif
00569 #ifndef ISOTHERMAL
00570 rms_error += SQR(total_error.E);
00571 #endif
00572 rms_error = sqrt(rms_error)/(double)count;
00573
00574
00575
00576
00577 #ifdef MPI_PARALLEL
00578 fname = ath_fname("../","LinWave-errors",NULL,NULL,1,wave_flag,NULL,"dat");
00579 #else
00580 fname = ath_fname(NULL,"LinWave-errors",NULL,NULL,1,wave_flag,NULL,"dat");
00581 #endif
00582
00583
00584 if((fp=fopen(fname,"r")) != NULL){
00585 if((fp = freopen(fname,"a",fp)) == NULL){
00586 ath_error("[Userwork_after_loop]: Unable to reopen file.\n");
00587 free(fname);
00588 return;
00589 }
00590 }
00591
00592 else{
00593 if((fp = fopen(fname,"w")) == NULL){
00594 ath_error("[Userwork_after_loop]: Unable to open file.\n");
00595 free(fname);
00596 return;
00597 }
00598
00599 fprintf(fp,"# Nx1 Nx2 Nx3 RMS-Error d M1 M2 M3");
00600 #ifndef ISOTHERMAL
00601 fprintf(fp," E");
00602 #endif
00603 #ifdef MHD
00604 fprintf(fp," B1c B2c B3c");
00605 #endif
00606 #if (NSCALARS > 0)
00607 for (n=0; n<NSCALARS; n++) {
00608 fprintf(fp," S[ %d ]",n);
00609 }
00610 #endif
00611 fprintf(fp,"\n#\n");
00612 }
00613
00614 fprintf(fp,"%d %d %d %e",Nx1,Nx2,Nx3,rms_error);
00615
00616 fprintf(fp," %e %e %e %e",
00617 (total_error.d/(double) count),
00618 (total_error.M1/(double)count),
00619 (total_error.M2/(double)count),
00620 (total_error.M3/(double)count));
00621
00622 #ifndef ISOTHERMAL
00623 fprintf(fp," %e",(total_error.E/(double)count));
00624 #endif
00625
00626 #ifdef MHD
00627 fprintf(fp," %e %e %e",
00628 (total_error.B1c/(double)count),
00629 (total_error.B2c/(double)count),
00630 (total_error.B3c/(double)count));
00631 #endif
00632 #if (NSCALARS > 0)
00633 for (n=0; n<NSCALARS; n++) {
00634 fprintf(fp," %e",total_error.s[n]/(double)count);
00635 }
00636 #endif
00637
00638 fprintf(fp,"\n");
00639
00640 fclose(fp);
00641 free(fname);
00642
00643 return;
00644 }