Athena: integrators/integrate_3d

Go to the documentation of this file.
00001 #include "../copyright.h"
00002 /*============================================================================*/
00003 /*! \file integrate_3d_ctu.c
00004  *  \brief Integrate MHD equations using 3D version of the directionally
00005  *   unsplit CTU integrator of Colella (1990). 
00006  *
00007  * PURPOSE: Integrate MHD equations using 3D version of the directionally
00008  *   unsplit CTU integrator of Colella (1990).  The variables updated are:
00009  *   -  U.[d,M1,M2,M3,E,B1c,B2c,B3c,s] -- where U is of type ConsS
00010  *   -  B1i, B2i, B3i  -- interface magnetic field
00011  *   Also adds gravitational source terms, self-gravity, optically thin cooling,
00012  *   shearing box source terms, and the H-correction of Sanders et al.
00013  *   - For adb hydro, requires (9*Cons1DS +  3*Real) = 48 3D arrays
00014  *   - For adb mhd, requires   (9*Cons1DS + 10*Real) = 73 3D arrays
00015  *   The H-correction of Sanders et al. adds another 3 arrays.  
00016  *
00017  * REFERENCES:
00018  * - P. Colella, "Multidimensional upwind methods for hyperbolic conservation
00019  *   laws", JCP, 87, 171 (1990)
00020  *
00021  * - T. Gardiner & J.M. Stone, "An unsplit Godunov method for ideal MHD via
00022  *   constrained transport in three dimensions", JCP, 227, 4123 (2008)
00023  *
00024  * - R. Sanders, E. Morano, & M.-C. Druguet, "Multidimensinal dissipation for
00025  *   upwind schemes: stability and applications to gas dynamics", JCP, 145, 511
00026  *   (1998)
00027  *
00028  * - J.M. Stone et al., "Athena: A new code for astrophysical MHD", ApJS,
00029  *   178, 137 (2008)
00030  *
00031  * CONTAINS PUBLIC FUNCTIONS: 
00032  * - integrate_3d_ctu()
00033  * - integrate_init_3d()
00034  * - integrate_destruct_3d() */
00035 /*============================================================================*/
00036 
00037 #include <math.h>
00038 #include <stdio.h>
00039 #include <stdlib.h>
00040 #include <string.h>
00041 #include "../defs.h"
00042 #include "../athena.h"
00043 #include "../globals.h"
00044 #include "prototypes.h"
00045 #include "../prototypes.h"
00046 #ifdef PARTICLES
00047 #include "../particles/particle.h"
00048 #endif
00049 
00050 #if defined(CTU_INTEGRATOR)
00051 #ifdef SPECIAL_RELATIVITY
00052 #error : The CTU integrator cannot be used for special relativity.
00053 #endif /* SPECIAL_RELATIVITY */
00054 
00055 /* The L/R states of conserved variables and fluxes at each cell face */
00056 static Cons1DS ***Ul_x1Face=NULL, ***Ur_x1Face=NULL;
00057 static Cons1DS ***Ul_x2Face=NULL, ***Ur_x2Face=NULL;
00058 static Cons1DS ***Ul_x3Face=NULL, ***Ur_x3Face=NULL;
00059 Cons1DS ***x1Flux=NULL, ***x2Flux=NULL, ***x3Flux=NULL;
00060 
00061 /* The interface magnetic fields and emfs */
00062 #ifdef MHD
00063 static Real ***B1_x1Face=NULL, ***B2_x2Face=NULL, ***B3_x3Face=NULL;
00064 Real ***emf1=NULL, ***emf2=NULL, ***emf3=NULL;
00065 static Real ***emf1_cc=NULL, ***emf2_cc=NULL, ***emf3_cc=NULL;
00066 #endif /* MHD */
00067 
00068 /* 1D scratch vectors used by lr_states and flux functions */
00069 static Real *Bxc=NULL, *Bxi=NULL;
00070 static Prim1DS *W=NULL, *Wl=NULL, *Wr=NULL;
00071 static Cons1DS *U1d=NULL;
00072 
00073 /* density and Pressure at t^{n+1/2} needed by MHD, cooling, and gravity */
00074 static Real ***dhalf = NULL, ***phalf=NULL;
00075 
00076 /* variables needed for H-correction of Sanders et al (1998) */
00077 extern Real etah;
00078 #ifdef H_CORRECTION
00079 static Real ***eta1=NULL, ***eta2=NULL, ***eta3=NULL;
00080 #endif
00081 
00082 /* variables needed to conserve net Bz in shearing box */
00083 #ifdef SHEARING_BOX
00084 static Real **remapEyiib=NULL, **remapEyoib=NULL;
00085 #endif
00086 
00087 /* VARIABLES NEED FOR CYLINDRICAL COORDINATES */
00088 #ifdef CYLINDRICAL
00089 static Real ***geom_src=NULL;
00090 #endif
00091 
00092 /*==============================================================================
00093  * PRIVATE FUNCTION PROTOTYPES: 
00094  *   integrate_emf1_corner() - the upwind CT method in GS05, for emf1
00095  *   integrate_emf2_corner() - the upwind CT method in GS05, for emf2
00096  *   integrate_emf3_corner() - the upwind CT method in GS05, for emf3
00097  *============================================================================*/
00098 
00099 #ifdef MHD
00100 static void integrate_emf1_corner(const GridS *pG);
00101 static void integrate_emf2_corner(const GridS *pG);
00102 static void integrate_emf3_corner(const GridS *pG);
00103 #endif /* MHD */
00104 
00105 /*=========================== PUBLIC FUNCTIONS ===============================*/
00106 /*----------------------------------------------------------------------------*/
00107 /*! \fn void integrate_3d_ctu(DomainS *pD)
00108  *  \brief 3D CTU integrator for MHD using 6-solve method */
00109 
00110 void integrate_3d_ctu(DomainS *pD)
00111 {
00112   GridS *pG=(pD->Grid);
00113   Real dtodx1=pG->dt/pG->dx1, dtodx2=pG->dt/pG->dx2, dtodx3=pG->dt/pG->dx3;
00114   Real hdt = 0.5*pG->dt, dx2=pG->dx2;
00115   Real q1 = 0.5*dtodx1, q2 = 0.5*dtodx2, q3 = 0.5*dtodx3;
00116   int i,il,iu, is = pG->is, ie = pG->ie;
00117   int j,jl,ju, js = pG->js, je = pG->je;
00118   int k,kl,ku, ks = pG->ks, ke = pG->ke;
00119   Real x1,x2,x3,phicl,phicr,phifc,phil,phir,phic,M1h,M2h,M3h,Bx=0.0;
00120 #ifndef BAROTROPIC
00121   Real coolfl,coolfr,coolf,Eh=0.0;
00122 #endif
00123 #ifdef MHD
00124   Real MHD_src_By,MHD_src_Bz,mdb1,mdb2,mdb3;
00125   Real db1,db2,db3,l1,l2,l3,B1,B2,B3,V1,V2,V3;
00126   Real B1ch,B2ch,B3ch;
00127 #endif
00128 #if defined(MHD) || defined(SELF_GRAVITY)
00129   Real dx1i=1.0/pG->dx1, dx2i=1.0/pG->dx2, dx3i=1.0/pG->dx3;
00130 #endif
00131 #ifdef H_CORRECTION
00132   Real cfr,cfl,lambdar,lambdal;
00133 #endif
00134 #if (NSCALARS > 0)
00135   int n;
00136 #endif
00137 #ifdef SELF_GRAVITY
00138   Real gxl,gxr,gyl,gyr,gzl,gzr,flx_m1l,flx_m1r,flx_m2l,flx_m2r,flx_m3l,flx_m3r;
00139 #endif
00140 #ifdef FEEDBACK
00141   Real dt1 = 1.0/pG->dt;
00142 #endif
00143 #ifdef SHEARING_BOX
00144   int my_iproc,my_jproc,my_kproc;
00145   Real M1n, dM2n; /* M1, dM2=(My+d*q*Omega_0*x) at time n */
00146   Real M1e, dM2e; /* M1, dM2 evolved by dt/2 */
00147   Real flx1_dM2, frx1_dM2, flx2_dM2, frx2_dM2, flx3_dM2, frx3_dM2;
00148   Real fact, qom, om_dt = Omega_0*pG->dt;
00149 #endif /* SHEARING_BOX */
00150 #ifdef STATIC_MESH_REFINEMENT
00151   int ncg,npg,dim;
00152   int ii,ics,ice,jj,jcs,jce,kk,kcs,kce,ips,ipe,jps,jpe,kps,kpe;
00153 #endif
00154 
00155   /* CYLINDRICAL VARIABLES */
00156 #ifdef CYLINDRICAL
00157 #ifndef ISOTHERMAL
00158   Real Pavgh;
00159 #endif
00160   Real g,gl,gr,rinv;
00161   Real geom_src_d, geom_src_Vx, geom_src_Vy, geom_src_P, geom_src_By, geom_src_Bz;
00162   const Real *r=pG->r, *ri=pG->ri;
00163 #ifdef FARGO
00164   Real Om, qshear, Mrn, Mpn, Mre, Mpe, Mrav, Mpav;
00165 #endif
00166 #endif /* CYLINDRICAL */
00167   Real lsf=1.0, rsf=1.0;
00168 
00169 /* With particles, one more ghost cell must be updated in predict step */
00170 #ifdef PARTICLES
00171   Real d1;
00172   il = is - 3;
00173   iu = ie + 3;
00174   jl = js - 3;
00175   ju = je + 3;
00176   kl = ks - 3;
00177   ku = ke + 3;
00178 #else
00179   il = is - 2;
00180   iu = ie + 2;
00181   jl = js - 2;
00182   ju = je + 2;
00183   kl = ks - 2;
00184   ku = ke + 2;
00185 #endif
00186 
00187 /* Set etah=0 so first calls to flux functions do not use H-correction */
00188   etah = 0.0;
00189 
00190 /* Compute predictor feedback from particle drag */
00191 #ifdef FEEDBACK
00192   feedback_predictor(pG);
00193 #endif
00194 
00195 /*=== STEP 1: Compute L/R x1-interface states and 1D x1-Fluxes ===============*/
00196 
00197 /*--- Step 1a ------------------------------------------------------------------
00198  * Load 1D vector of conserved variables;
00199  * U1d = (d, M1, M2, M3, E, B2c, B3c, s[n])
00200  */
00201 
00202   for (k=kl; k<=ku; k++) {
00203     for (j=jl; j<=ju; j++) {
00204       for (i=is-nghost; i<=ie+nghost; i++) {
00205         U1d[i].d  = pG->U[k][j][i].d;
00206         U1d[i].Mx = pG->U[k][j][i].M1;
00207         U1d[i].My = pG->U[k][j][i].M2;
00208         U1d[i].Mz = pG->U[k][j][i].M3;
00209 #ifndef BAROTROPIC
00210         U1d[i].E  = pG->U[k][j][i].E;
00211 #endif /* BAROTROPIC */
00212 #ifdef MHD
00213         U1d[i].By = pG->U[k][j][i].B2c;
00214         U1d[i].Bz = pG->U[k][j][i].B3c;
00215         Bxc[i] = pG->U[k][j][i].B1c;
00216         Bxi[i] = pG->B1i[k][j][i];
00217         B1_x1Face[k][j][i] = pG->B1i[k][j][i];
00218 #endif /* MHD */
00219 #if (NSCALARS > 0)
00220         for (n=0; n<NSCALARS; n++) U1d[i].s[n] = pG->U[k][j][i].s[n];
00221 #endif
00222       }
00223 
00224 /*--- Step 1b ------------------------------------------------------------------
00225  * Compute L and R states at X1-interfaces, add "MHD source terms" for 0.5*dt
00226  */
00227 
00228      for (i=is-nghost; i<=ie+nghost; i++) {
00229        W[i] = Cons1D_to_Prim1D(&U1d[i],&Bxc[i]);
00230 
00231         /* CALCULATE THE CELL-CENTERED GEOMETRIC SOURCE VECTOR NOW USING U^{n}
00232         * THIS WILL BE USED AT THE END OF THE INTEGRATION STEP AS A SOURCE TERM
00233         * FOR THE CELL-CENTERED CONSERVED VARIABLES (STEPS 6D,7D,8B) */
00234 #ifdef CYLINDRICAL
00235         geom_src[k][j][i]  = W[i].d*SQR(W[i].Vy);
00236 #ifdef MHD
00237         geom_src[k][j][i] += 0.5*(SQR(Bxc[i]) - SQR(W[i].By) + SQR(W[i].Bz));
00238 #endif
00239 #ifdef ISOTHERMAL
00240         geom_src[k][j][i] += Iso_csound2*W[i].d;
00241 #else
00242         geom_src[k][j][i] += W[i].P;
00243 #endif
00244         geom_src[k][j][i] /= x1vc(pG,i);
00245 #endif /* CYLINDRICAL */
00246      }
00247 
00248      lr_states(pG,W,Bxc,pG->dt,pG->dx1,il+1,iu-1,Wl,Wr,1);
00249 
00250 #ifdef MHD
00251       for (i=il+1; i<=iu; i++) {
00252 /* Source terms for left states in zone i-1 */
00253 #ifdef CYLINDRICAL
00254         rsf = ri[i]/r[i-1];  lsf = ri[i-1]/r[i-1];
00255         dx2i = 1.0/(r[i-1]*pG->dx2);
00256 #endif
00257         db1 = (rsf*pG->B1i[k  ][j  ][i  ] - lsf*pG->B1i[k][j][i-1])*dx1i;
00258         db2 = (    pG->B2i[k  ][j+1][i-1] -     pG->B2i[k][j][i-1])*dx2i;
00259         db3 = (    pG->B3i[k+1][j  ][i-1] -     pG->B3i[k][j][i-1])*dx3i;
00260 
00261         if(db1 >= 0.0){
00262           l3 = db1 < -db3 ? db1 : -db3;
00263           l3 = l3 > 0.0 ? l3 : 0.0;
00264 
00265           l2 = db1 < -db2 ? db1 : -db2;
00266           l2 = l2 > 0.0 ? l2 : 0.0;
00267         }
00268         else{
00269           l3 = db1 > -db3 ? db1 : -db3;
00270           l3 = l3 < 0.0 ? l3 : 0.0;
00271 
00272           l2 = db1 > -db2 ? db1 : -db2;
00273           l2 = l2 < 0.0 ? l2 : 0.0;
00274         }
00275 
00276         MHD_src_By = (pG->U[k][j][i-1].M2/pG->U[k][j][i-1].d)*l2;
00277         MHD_src_Bz = (pG->U[k][j][i-1].M3/pG->U[k][j][i-1].d)*l3;
00278 
00279         Wl[i].By += hdt*MHD_src_By;
00280         Wl[i].Bz += hdt*MHD_src_Bz;
00281 
00282 /* Source terms for right states in zone i */
00283 #ifdef CYLINDRICAL
00284         rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
00285         dx2i = 1.0/(r[i]*pG->dx2);
00286 #endif
00287         db1 = (rsf*pG->B1i[k  ][j  ][i+1] - lsf*pG->B1i[k][j][i])*dx1i;
00288         db2 = (    pG->B2i[k  ][j+1][i  ] -     pG->B2i[k][j][i])*dx2i;
00289         db3 = (    pG->B3i[k+1][j  ][i  ] -     pG->B3i[k][j][i])*dx3i;
00290 
00291         if(db1 >= 0.0){
00292           l3 = db1 < -db3 ? db1 : -db3;
00293           l3 = l3 > 0.0 ? l3 : 0.0;
00294 
00295           l2 = db1 < -db2 ? db1 : -db2;
00296           l2 = l2 > 0.0 ? l2 : 0.0;
00297         }
00298         else{
00299           l3 = db1 > -db3 ? db1 : -db3;
00300           l3 = l3 < 0.0 ? l3 : 0.0;
00301 
00302           l2 = db1 > -db2 ? db1 : -db2;
00303           l2 = l2 < 0.0 ? l2 : 0.0;
00304         }
00305 
00306         MHD_src_By = (pG->U[k][j][i].M2/pG->U[k][j][i].d)*l2;
00307         MHD_src_Bz = (pG->U[k][j][i].M3/pG->U[k][j][i].d)*l3;
00308 
00309         Wr[i].By += hdt*MHD_src_By;
00310         Wr[i].Bz += hdt*MHD_src_Bz;
00311       }
00312 #endif
00313 
00314 /*--- Step 1c ------------------------------------------------------------------
00315  * Add source terms from static gravitational potential for 0.5*dt to L/R states
00316  */
00317 
00318       if (StaticGravPot != NULL){
00319         for (i=il+1; i<=iu; i++) {
00320           cc_pos(pG,i,j,k,&x1,&x2,&x3);
00321 #ifdef CYLINDRICAL
00322           gl = (*x1GravAcc)(x1vc(pG,i-1),x2,x3);
00323           gr = (*x1GravAcc)(x1vc(pG,i),x2,x3);
00324 #ifdef FARGO
00325           /* Correct for force in the rotating frame */
00326                                         gl = gl - x1vc(pG,i-1)*SQR((*OrbitalProfile)(x1vc(pG,i-1)));
00327                                         gr = gr - x1vc(pG,i)*SQR((*OrbitalProfile)(x1vc(pG,i)));
00328 #endif
00329           /* APPLY GRAV. SOURCE TERMS TO VELOCITY USING ACCELERATION
00330           * FOR (dt/2) */
00331           Wl[i].Vx -= hdt*gl;
00332           Wr[i].Vx -= hdt*gr;
00333 #else
00334           phicr = (*StaticGravPot)( x1             ,x2,x3);
00335           phicl = (*StaticGravPot)((x1-    pG->dx1),x2,x3);
00336           phifc = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
00337 
00338           Wl[i].Vx -= dtodx1*(phifc - phicl);
00339           Wr[i].Vx -= dtodx1*(phicr - phifc);
00340 #endif
00341         }
00342       }
00343 
00344 /*--- Step 1c (cont) -----------------------------------------------------------
00345  * Add source terms for self-gravity for 0.5*dt to L/R states
00346  */
00347 
00348 #ifdef SELF_GRAVITY
00349       for (i=il+1; i<=iu; i++) {
00350         Wl[i].Vx -= q1*(pG->Phi[k][j][i] - pG->Phi[k][j][i-1]);
00351         Wr[i].Vx -= q1*(pG->Phi[k][j][i] - pG->Phi[k][j][i-1]);
00352       }
00353 #endif
00354 
00355 /*--- Step 1c (cont) -----------------------------------------------------------
00356  * Add source terms from optically-thin cooling for 0.5*dt to L/R states
00357  */
00358 
00359 #ifndef BAROTROPIC
00360       if (CoolingFunc != NULL){
00361         for (i=il+1; i<=iu; i++) {
00362           coolfl = (*CoolingFunc)(Wl[i].d,Wl[i].P,(0.5*pG->dt));
00363           coolfr = (*CoolingFunc)(Wr[i].d,Wr[i].P,(0.5*pG->dt));
00364 
00365           Wl[i].P -= 0.5*pG->dt*Gamma_1*coolfl;
00366           Wr[i].P -= 0.5*pG->dt*Gamma_1*coolfr;
00367         }
00368       }
00369 #endif /* BAROTROPIC */
00370 
00371 /*--- Step 1c (cont) -----------------------------------------------------------
00372  * Add source terms for shearing box (Coriolis forces) for 0.5*dt to L/R states
00373  * starting with tidal gravity terms added through the ShearingBoxPot
00374  *    Vx source term = (dt/2)*( 2 Omega_0 Vy)
00375  *    Vy source term = (dt/2)*(-2 Omega_0 Vx)
00376  *    Vy source term = (dt/2)*((q-2) Omega_0 Vx) (with FARGO)
00377  */
00378 
00379 #ifdef SHEARING_BOX
00380       if (ShearingBoxPot != NULL){
00381         for (i=il+1; i<=iu; i++) {
00382           cc_pos(pG,i,j,k,&x1,&x2,&x3);
00383           phicr = (*ShearingBoxPot)( x1             ,x2,x3);
00384           phicl = (*ShearingBoxPot)((x1-    pG->dx1),x2,x3);
00385           phifc = (*ShearingBoxPot)((x1-0.5*pG->dx1),x2,x3);
00386 
00387           Wl[i].Vx -= dtodx1*(phifc - phicl);
00388           Wr[i].Vx -= dtodx1*(phicr - phifc);
00389         }
00390       }
00391 
00392       for (i=il+1; i<=iu; i++) {
00393         Wl[i].Vx += pG->dt*Omega_0*W[i-1].Vy;
00394         Wr[i].Vx += pG->dt*Omega_0*W[i].Vy;
00395 #ifdef FARGO
00396         Wl[i].Vy += hdt*(qshear-2.)*Omega_0*W[i-1].Vx;
00397         Wr[i].Vy += hdt*(qshear-2.)*Omega_0*W[i].Vx;
00398 #else
00399         Wl[i].Vy -= pG->dt*Omega_0*W[i-1].Vx;
00400         Wr[i].Vy -= pG->dt*Omega_0*W[i].Vx;
00401 #endif
00402       }
00403 #endif /* SHEARING_BOX */
00404 
00405 #if defined(CYLINDRICAL) && defined(FARGO)
00406         for (i=il+1; i<=iu; i++) {
00407                 Om = (*OrbitalProfile)(x1vc(pG,i-1));
00408     qshear = (*ShearProfile)(x1vc(pG,i-1));
00409     Wl[i].Vx += (pG->dt)*Om*W[i-1].Vy;
00410     Wl[i].Vy += hdt*(qshear - 2.0)*Om*W[i-1].Vx;
00411 
00412                 Om = (*OrbitalProfile)(x1vc(pG,i));
00413     qshear = (*ShearProfile)(x1vc(pG,i));
00414     Wr[i].Vx += (pG->dt)*Om*W[i].Vy;
00415     Wr[i].Vy += hdt*(qshear - 2.0)*Om*W[i].Vx;
00416         }
00417 #endif /* Cylindrical + Fargo */
00418 
00419 /*--- Step 1c (cont) -----------------------------------------------------------
00420  * Add source terms for particle feedback for 0.5*dt to L/R states
00421  */
00422 
00423 #ifdef FEEDBACK
00424     for (i=il+1; i<=iu; i++) {
00425       d1 = 1.0/W[i-1].d;
00426       Wl[i].Vx -= pG->Coup[k][j][i-1].fb1*d1;
00427       Wl[i].Vy -= pG->Coup[k][j][i-1].fb2*d1;
00428       Wl[i].Vz -= pG->Coup[k][j][i-1].fb3*d1;
00429 
00430       d1 = 1.0/W[i].d;
00431       Wr[i].Vx -= pG->Coup[k][j][i].fb1*d1;
00432       Wr[i].Vy -= pG->Coup[k][j][i].fb2*d1;
00433       Wr[i].Vz -= pG->Coup[k][j][i].fb3*d1;
00434 
00435 #ifndef BAROTROPIC
00436       Wl[i].P += pG->Coup[k][j][i-1].Eloss*Gamma_1;
00437       Wr[i].P += pG->Coup[k][j][i].Eloss*Gamma_1;
00438 #endif
00439 
00440     }
00441 #endif /* FEEDBACK */
00442 
00443 /*--- Step 1c (cont) -----------------------------------------------------------
00444  * ADD THE GEOMETRIC SOURCE-TERMS NOW USING CELL-CENTERED PRIMITIVE 
00445  * VARIABLES AT TIME t^n
00446  */
00447 #ifdef CYLINDRICAL
00448       for (i=is-1; i<=ie+2; i++) {
00449 
00450         /* LEFT STATE GEOMETRIC SOURCE TERM (USES W[i-1]) */
00451         rinv = 1.0/x1vc(pG,i-1);
00452         geom_src_d  = -W[i-1].d*W[i-1].Vx*rinv;
00453         geom_src_Vx =  SQR(W[i-1].Vy);
00454         geom_src_Vy = -W[i-1].Vx*W[i-1].Vy;
00455 #ifdef MHD
00456         geom_src_Vx -= SQR(W[i-1].By)/W[i-1].d;
00457         geom_src_Vy += Bxc[i-1]*W[i-1].By/W[i-1].d;
00458         geom_src_By =  -W[i-1].Vy*Bxc[i-1]*rinv;
00459         geom_src_Bz =  -W[i-1].Vx*W[i-1].Bz*rinv;
00460 #endif /* MHD */
00461         geom_src_Vx *= rinv;
00462         geom_src_Vy *= rinv;
00463 #ifndef ISOTHERMAL
00464         geom_src_P  = -Gamma*W[i-1].P*W[i-1].Vx*rinv;
00465 #endif /* ISOTHERMAL */
00466 
00467         /* ADD SOURCE TERM TO LEFT STATE */
00468         Wl[i].d  += hdt*geom_src_d;
00469         Wl[i].Vx += hdt*geom_src_Vx;
00470         Wl[i].Vy += hdt*geom_src_Vy;
00471 #ifdef MHD
00472         Wl[i].By += hdt*geom_src_By;
00473         Wl[i].Bz += hdt*geom_src_Bz;
00474 #endif /* MHD */
00475 #ifndef ISOTHERMAL
00476         Wl[i].P  += hdt*geom_src_P;
00477 #endif /* ISOTHERMAL */
00478 
00479 
00480         /* RIGHT STATE GEOMETRIC SOURCE TERM (USES W[i]) */
00481         rinv = 1.0/x1vc(pG,i);
00482         geom_src_d  = -W[i].d*W[i].Vx*rinv;
00483         geom_src_Vx =  SQR(W[i].Vy);
00484         geom_src_Vy = -W[i].Vx*W[i].Vy;
00485 #ifdef MHD
00486         geom_src_Vx -= SQR(W[i].By)/W[i].d;
00487         geom_src_Vy += Bxc[i]*W[i].By/W[i].d;
00488         geom_src_By =  -W[i].Vy*Bxc[i]*rinv;
00489         geom_src_Bz =  -W[i].Vx*W[i].Bz*rinv;
00490 #endif /* MHD */
00491         geom_src_Vx *= rinv;
00492         geom_src_Vy *= rinv;
00493 #ifndef ISOTHERMAL
00494         geom_src_P  = -Gamma*W[i].P*W[i].Vx*rinv;
00495 #endif /* ISOTHERMAL */
00496 
00497         /* ADD SOURCE TERM TO RIGHT STATE */
00498         Wr[i].d  += hdt*geom_src_d;
00499         Wr[i].Vx += hdt*geom_src_Vx;
00500         Wr[i].Vy += hdt*geom_src_Vy;
00501 #ifdef MHD
00502         Wr[i].By += hdt*geom_src_By;
00503         Wr[i].Bz += hdt*geom_src_Bz;
00504 #endif /* MHD */
00505 #ifndef ISOTHERMAL
00506         Wr[i].P  += hdt*geom_src_P;
00507 #endif /* ISOTHERMAL */
00508       }
00509 #endif /* CYLINDRICAL */
00510 
00511 /*--- Step 1d ------------------------------------------------------------------
00512  * Compute 1D fluxes in x1-direction, storing into 3D array
00513  */
00514 
00515       for (i=il+1; i<=iu; i++) {
00516         Ul_x1Face[k][j][i] = Prim1D_to_Cons1D(&Wl[i],&Bxi[i]);
00517         Ur_x1Face[k][j][i] = Prim1D_to_Cons1D(&Wr[i],&Bxi[i]);
00518 
00519 #ifdef MHD
00520         Bx = B1_x1Face[k][j][i];
00521 #endif
00522         fluxes(Ul_x1Face[k][j][i],Ur_x1Face[k][j][i],Wl[i],Wr[i],Bx,
00523           &x1Flux[k][j][i]);
00524       }
00525     }
00526   }
00527 
00528 /*=== STEP 2: Compute L/R x2-interface states and 1D x2-Fluxes ===============*/
00529 
00530 /*--- Step 2a ------------------------------------------------------------------
00531  * Load 1D vector of conserved variables;
00532  * U1d = (d, M2, M3, M1, E, B3c, B1c, s[n])
00533  */
00534 
00535   for (k=kl; k<=ku; k++) {
00536     for (i=il; i<=iu; i++) {
00537 #ifdef CYLINDRICAL
00538       rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
00539       dx2 = r[i]*pG->dx2;
00540       dtodx2 = pG->dt/dx2;
00541 #endif
00542       for (j=js-nghost; j<=je+nghost; j++) {
00543         U1d[j].d  = pG->U[k][j][i].d;
00544         U1d[j].Mx = pG->U[k][j][i].M2;
00545         U1d[j].My = pG->U[k][j][i].M3;
00546         U1d[j].Mz = pG->U[k][j][i].M1;
00547 #ifndef BAROTROPIC
00548         U1d[j].E  = pG->U[k][j][i].E;
00549 #endif /* BAROTROPIC */
00550 #ifdef MHD
00551         U1d[j].By = pG->U[k][j][i].B3c;
00552         U1d[j].Bz = pG->U[k][j][i].B1c;
00553         Bxc[j] = pG->U[k][j][i].B2c;
00554         Bxi[j] = pG->B2i[k][j][i];
00555         B2_x2Face[k][j][i] = pG->B2i[k][j][i];
00556 #endif /* MHD */
00557 #if (NSCALARS > 0)
00558         for (n=0; n<NSCALARS; n++) U1d[j].s[n] = pG->U[k][j][i].s[n];
00559 #endif
00560       }
00561 
00562 /*--- Step 2b ------------------------------------------------------------------
00563  * Compute L and R states at X2-interfaces, add "MHD source terms" for 0.5*dt
00564  */
00565 
00566       for (j=js-nghost; j<=je+nghost; j++) {
00567         W[j] = Cons1D_to_Prim1D(&U1d[j],&Bxc[j]);
00568       }
00569 
00570       lr_states(pG,W,Bxc,pG->dt,dx2,jl+1,ju-1,Wl,Wr,2);
00571 
00572 #ifdef MHD
00573 #ifdef CYLINDRICAL
00574       dx2i = 1.0/dx2;
00575 #endif
00576       for (j=jl+1; j<=ju; j++) {
00577 /* Source terms for left states in zone j-1 */
00578         db1 = (rsf*pG->B1i[k  ][j-1][i+1] - lsf*pG->B1i[k][j-1][i])*dx1i;
00579         db2 = (    pG->B2i[k  ][j  ][i  ] -     pG->B2i[k][j-1][i])*dx2i;
00580         db3 = (    pG->B3i[k+1][j-1][i  ] -     pG->B3i[k][j-1][i])*dx3i;
00581 
00582         if(db2 >= 0.0){
00583           l1 = db2 < -db1 ? db2 : -db1;
00584           l1 = l1 > 0.0 ? l1 : 0.0;
00585 
00586           l3 = db2 < -db3 ? db2 : -db3;
00587           l3 = l3 > 0.0 ? l3 : 0.0;
00588         }
00589         else{
00590           l1 = db2 > -db1 ? db2 : -db1;
00591           l1 = l1 < 0.0 ? l1 : 0.0;
00592 
00593           l3 = db2 > -db3 ? db2 : -db3;
00594           l3 = l3 < 0.0 ? l3 : 0.0;
00595         }
00596 
00597         MHD_src_By = (pG->U[k][j-1][i].M3/pG->U[k][j-1][i].d)*l3;
00598         MHD_src_Bz = (pG->U[k][j-1][i].M1/pG->U[k][j-1][i].d)*l1;
00599 
00600         Wl[j].By += hdt*MHD_src_By;
00601         Wl[j].Bz += hdt*MHD_src_Bz;
00602 
00603 /* Source terms for right states in zone j */
00604         db1 = (rsf*pG->B1i[k  ][j  ][i+1] - lsf*pG->B1i[k][j][i])*dx1i;
00605         db2 = (    pG->B2i[k  ][j+1][i  ] -     pG->B2i[k][j][i])*dx2i;
00606         db3 = (    pG->B3i[k+1][j  ][i  ] -     pG->B3i[k][j][i])*dx3i;
00607 
00608         if(db2 >= 0.0){
00609           l1 = db2 < -db1 ? db2 : -db1;
00610           l1 = l1 > 0.0 ? l1 : 0.0;
00611 
00612           l3 = db2 < -db3 ? db2 : -db3;
00613           l3 = l3 > 0.0 ? l3 : 0.0;
00614         }
00615         else{
00616           l1 = db2 > -db1 ? db2 : -db1;
00617           l1 = l1 < 0.0 ? l1 : 0.0;
00618 
00619           l3 = db2 > -db3 ? db2 : -db3;
00620           l3 = l3 < 0.0 ? l3 : 0.0;
00621         }
00622 
00623         MHD_src_By = (pG->U[k][j][i].M3/pG->U[k][j][i].d)*l3;
00624         MHD_src_Bz = (pG->U[k][j][i].M1/pG->U[k][j][i].d)*l1;
00625 
00626         Wr[j].By += hdt*MHD_src_By;
00627         Wr[j].Bz += hdt*MHD_src_Bz;
00628       }
00629 #endif
00630 
00631 /*--- Step 2c ------------------------------------------------------------------
00632  * Add source terms from static gravitational potential for 0.5*dt to L/R states
00633  */
00634 
00635       if (StaticGravPot != NULL){
00636         for (j=jl+1; j<=ju; j++) {
00637           cc_pos(pG,i,j,k,&x1,&x2,&x3);
00638           phicr = (*StaticGravPot)(x1, x2             ,x3);
00639           phicl = (*StaticGravPot)(x1,(x2-    pG->dx2),x3);
00640           phifc = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
00641 
00642           Wl[j].Vx -= dtodx2*(phifc - phicl);
00643           Wr[j].Vx -= dtodx2*(phicr - phifc);
00644         }
00645       }
00646 
00647 /*--- Step 2c (cont) -----------------------------------------------------------
00648  * Add source terms for self-gravity for 0.5*dt to L/R states
00649  */
00650 
00651 #ifdef SELF_GRAVITY
00652       for (j=jl+1; j<=ju; j++) {
00653         Wl[j].Vx -= q2*(pG->Phi[k][j][i] - pG->Phi[k][j-1][i]);
00654         Wr[j].Vx -= q2*(pG->Phi[k][j][i] - pG->Phi[k][j-1][i]);
00655       }
00656 #endif
00657 
00658 /*--- Step 2c (cont) -----------------------------------------------------------
00659  * Add source terms from optically-thin cooling for 0.5*dt to L/R states
00660  */
00661 
00662 #ifndef BAROTROPIC
00663       if (CoolingFunc != NULL){
00664         for (j=jl+1; j<=ju; j++) {
00665           coolfl = (*CoolingFunc)(Wl[j].d,Wl[j].P,(0.5*pG->dt));
00666           coolfr = (*CoolingFunc)(Wr[j].d,Wr[j].P,(0.5*pG->dt));
00667 
00668           Wl[j].P -= 0.5*pG->dt*Gamma_1*coolfl;
00669           Wr[j].P -= 0.5*pG->dt*Gamma_1*coolfr;
00670         }
00671       }
00672 #endif /* BAROTROPIC */
00673 
00674 /*--- Step 2c (cont) -----------------------------------------------------------
00675  * Add source terms for particle feedback for 0.5*dt to L/R states
00676  */
00677 
00678 #ifdef FEEDBACK
00679    for (j=jl+1; j<=ju; j++) {
00680       d1 = 1.0/W[j-1].d;
00681       Wl[j].Vx -= pG->Coup[k][j-1][i].fb2*d1;
00682       Wl[j].Vy -= pG->Coup[k][j-1][i].fb3*d1;
00683       Wl[j].Vz -= pG->Coup[k][j-1][i].fb1*d1;
00684 
00685       d1 = 1.0/W[j].d;
00686       Wr[j].Vx -= pG->Coup[k][j][i].fb2*d1;
00687       Wr[j].Vy -= pG->Coup[k][j][i].fb3*d1;
00688       Wr[j].Vz -= pG->Coup[k][j][i].fb1*d1;
00689 
00690 #ifndef BAROTROPIC
00691       Wl[i].P += pG->Coup[k][j-1][i].Eloss*Gamma_1;
00692       Wr[i].P += pG->Coup[k][j][i].Eloss*Gamma_1;
00693 #endif
00694 
00695     }
00696 #endif /* FEEDBACK */
00697 
00698 
00699 /*--- Step 2d ------------------------------------------------------------------
00700  * Compute 1D fluxes in x2-direction, storing into 3D array
00701  */
00702 
00703       for (j=jl+1; j<=ju; j++) {
00704         Ul_x2Face[k][j][i] = Prim1D_to_Cons1D(&Wl[j],&Bxi[j]);
00705         Ur_x2Face[k][j][i] = Prim1D_to_Cons1D(&Wr[j],&Bxi[j]);
00706 
00707 #ifdef MHD
00708         Bx = B2_x2Face[k][j][i];
00709 #endif
00710         fluxes(Ul_x2Face[k][j][i],Ur_x2Face[k][j][i],Wl[j],Wr[j],Bx,
00711           &x2Flux[k][j][i]);
00712       }
00713     }
00714   }
00715 
00716 /*=== STEP 3: Compute L/R x3-interface states and 1D x3-Fluxes ===============*/
00717 
00718 /*--- Step 3a ------------------------------------------------------------------
00719  * Load 1D vector of conserved variables;
00720  * U1d = (d, M3, M1, M2, E, B1c, B2c, s[n])
00721  */
00722 
00723   for (j=jl; j<=ju; j++) {
00724     for (i=il; i<=iu; i++) {
00725       for (k=ks-nghost; k<=ke+nghost; k++) {
00726         U1d[k].d  = pG->U[k][j][i].d;
00727         U1d[k].Mx = pG->U[k][j][i].M3;
00728         U1d[k].My = pG->U[k][j][i].M1;
00729         U1d[k].Mz = pG->U[k][j][i].M2;
00730 #ifndef BAROTROPIC
00731         U1d[k].E  = pG->U[k][j][i].E;
00732 #endif /* BAROTROPIC */
00733 #ifdef MHD
00734         U1d[k].By = pG->U[k][j][i].B1c;
00735         U1d[k].Bz = pG->U[k][j][i].B2c;
00736         Bxc[k] = pG->U[k][j][i].B3c;
00737         Bxi[k] = pG->B3i[k][j][i];
00738         B3_x3Face[k][j][i] = pG->B3i[k][j][i];
00739 #endif /* MHD */
00740 #if (NSCALARS > 0)
00741         for (n=0; n<NSCALARS; n++) U1d[k].s[n] = pG->U[k][j][i].s[n];
00742 #endif
00743       }
00744 
00745 /*--- Step 3b ------------------------------------------------------------------
00746  * Compute L and R states at X3-interfaces, add "MHD source terms" for 0.5*dt
00747  */
00748 
00749       for (k=ks-nghost; k<=ke+nghost; k++) {
00750         W[k] = Cons1D_to_Prim1D(&U1d[k],&Bxc[k]);
00751       }
00752 
00753       lr_states(pG,W,Bxc,pG->dt,pG->dx3,kl+1,ku-1,Wl,Wr,3);
00754 
00755 #ifdef MHD
00756 #ifdef CYLINDRICAL
00757       rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
00758       dx2i = 1.0/(r[i]*pG->dx2);
00759 #endif
00760       for (k=kl+1; k<=ku; k++) {
00761 /* Source terms for left states in zone k-1 */
00762         db1 = (rsf*pG->B1i[k-1][j  ][i+1] - lsf*pG->B1i[k-1][j][i])*dx1i;
00763         db2 = (    pG->B2i[k-1][j+1][i  ] -     pG->B2i[k-1][j][i])*dx2i;
00764         db3 = (    pG->B3i[k  ][j  ][i  ] -     pG->B3i[k-1][j][i])*dx3i;
00765 
00766         if(db3 >= 0.0){
00767           l1 = db3 < -db1 ? db3 : -db1;
00768           l1 = l1 > 0.0 ? l1 : 0.0;
00769 
00770           l2 = db3 < -db2 ? db3 : -db2;
00771           l2 = l2 > 0.0 ? l2 : 0.0;
00772         }
00773         else{
00774           l1 = db3 > -db1 ? db3 : -db1;
00775           l1 = l1 < 0.0 ? l1 : 0.0;
00776 
00777           l2 = db3 > -db2 ? db3 : -db2;
00778           l2 = l2 < 0.0 ? l2 : 0.0;
00779         }
00780 
00781         MHD_src_By = (pG->U[k-1][j][i].M1/pG->U[k-1][j][i].d)*l1;
00782         MHD_src_Bz = (pG->U[k-1][j][i].M2/pG->U[k-1][j][i].d)*l2;
00783 
00784         Wl[k].By += hdt*MHD_src_By;
00785         Wl[k].Bz += hdt*MHD_src_Bz;
00786 
00787 /* Source terms for right states in zone k */
00788         db1 = (rsf*pG->B1i[k][j][i+1] - lsf*pG->B1i[k][j][i])*dx1i;
00789         db2 = (    pG->B2i[k][j+1][i] -     pG->B2i[k][j][i])*dx2i;
00790         db3 = (    pG->B3i[k+1][j][i] -     pG->B3i[k][j][i])*dx3i;
00791 
00792         if(db3 >= 0.0){
00793           l1 = db3 < -db1 ? db3 : -db1;
00794           l1 = l1 > 0.0 ? l1 : 0.0;
00795 
00796           l2 = db3 < -db2 ? db3 : -db2;
00797           l2 = l2 > 0.0 ? l2 : 0.0;
00798         }
00799         else{
00800           l1 = db3 > -db1 ? db3 : -db1;
00801           l1 = l1 < 0.0 ? l1 : 0.0;
00802 
00803           l2 = db3 > -db2 ? db3 : -db2;
00804           l2 = l2 < 0.0 ? l2 : 0.0;
00805         }
00806 
00807         MHD_src_By = (pG->U[k][j][i].M1/pG->U[k][j][i].d)*l1;
00808         MHD_src_Bz = (pG->U[k][j][i].M2/pG->U[k][j][i].d)*l2;
00809 
00810         Wr[k].By += hdt*MHD_src_By;
00811         Wr[k].Bz += hdt*MHD_src_Bz;
00812       }
00813 #endif
00814 
00815 /*--- Step 3c ------------------------------------------------------------------
00816  * Add source terms from static gravitational potential for 0.5*dt to L/R states
00817  */
00818 
00819       if (StaticGravPot != NULL){
00820         for (k=kl+1; k<=ku; k++) {
00821           cc_pos(pG,i,j,k,&x1,&x2,&x3);
00822           phicr = (*StaticGravPot)(x1,x2, x3             );
00823           phicl = (*StaticGravPot)(x1,x2,(x3-    pG->dx3));
00824           phifc = (*StaticGravPot)(x1,x2,(x3-0.5*pG->dx3));
00825 
00826           Wl[k].Vx -= dtodx3*(phifc - phicl);
00827           Wr[k].Vx -= dtodx3*(phicr - phifc);
00828         }
00829       }
00830 
00831 /*--- Step 3c (cont) -----------------------------------------------------------
00832  * Add source terms for self-gravity for 0.5*dt to L/R states
00833  */
00834 
00835 #ifdef SELF_GRAVITY
00836       for (k=kl+1; k<=ku; k++) {
00837         Wl[k].Vx -= q3*(pG->Phi[k][j][i] - pG->Phi[k-1][j][i]);
00838         Wr[k].Vx -= q3*(pG->Phi[k][j][i] - pG->Phi[k-1][j][i]);
00839       }
00840 #endif
00841 
00842 /*--- Step 3c (cont) -----------------------------------------------------------
00843  * Add source terms from optically-thin cooling for 0.5*dt to L/R states
00844  */
00845 
00846 #ifndef BAROTROPIC
00847       if (CoolingFunc != NULL){
00848         for (k=kl+1; k<=ku; k++) {
00849           coolfl = (*CoolingFunc)(Wl[k].d,Wl[k].P,(0.5*pG->dt));
00850           coolfr = (*CoolingFunc)(Wr[k].d,Wr[k].P,(0.5*pG->dt));
00851   
00852           Wl[k].P -= 0.5*pG->dt*Gamma_1*coolfl;
00853           Wr[k].P -= 0.5*pG->dt*Gamma_1*coolfr;
00854         }
00855       }
00856 #endif /* BAROTROPIC */
00857 
00858 /*--- Step 3c (cont) -----------------------------------------------------------
00859  * Add source terms for particle feedback for 0.5*dt to L/R states
00860  */
00861 
00862 #ifdef FEEDBACK
00863    for (k=kl+1; k<=ku; k++) {
00864       d1 = 1.0/W[k-1].d;
00865       Wl[k].Vx -= pG->Coup[k-1][j][i].fb3*d1;
00866       Wl[k].Vy -= pG->Coup[k-1][j][i].fb1*d1;
00867       Wl[k].Vz -= pG->Coup[k-1][j][i].fb2*d1;
00868 
00869       d1 = 1.0/W[k].d;
00870       Wr[k].Vx -= pG->Coup[k][j][i].fb3*d1;
00871       Wr[k].Vy -= pG->Coup[k][j][i].fb1*d1;
00872       Wr[k].Vz -= pG->Coup[k][j][i].fb2*d1;
00873 
00874 #ifndef BAROTROPIC
00875       Wl[i].P += pG->Coup[k-1][j][i].Eloss*Gamma_1;
00876       Wr[i].P += pG->Coup[k][j][i].Eloss*Gamma_1;
00877 #endif
00878     }
00879 #endif /* FEEDBACK */
00880 
00881 
00882 /*--- Step 3d ------------------------------------------------------------------
00883  * Compute 1D fluxes in x3-direction, storing into 3D array
00884  */
00885 
00886       for (k=kl+1; k<=ku; k++) {
00887         Ul_x3Face[k][j][i] = Prim1D_to_Cons1D(&Wl[k],&Bxi[k]);
00888         Ur_x3Face[k][j][i] = Prim1D_to_Cons1D(&Wr[k],&Bxi[k]);
00889 
00890 #ifdef MHD
00891         Bx = B3_x3Face[k][j][i];
00892 #endif
00893         fluxes(Ul_x3Face[k][j][i],Ur_x3Face[k][j][i],Wl[k],Wr[k],Bx,
00894           &x3Flux[k][j][i]);
00895       }
00896     }
00897   }
00898 
00899 /*=== STEP 4:  Update face-centered B for 0.5*dt =============================*/
00900 
00901 /*--- Step 4a ------------------------------------------------------------------
00902  * Calculate the cell centered value of emf1,2,3 at t^{n} and integrate
00903  * to corner.
00904  */
00905 
00906 #ifdef MHD
00907 /* emf1 */
00908   for (k=kl; k<=ku; k++) {
00909     for (j=jl; j<=ju; j++) {
00910       for (i=il; i<=iu; i++) {
00911         emf1_cc[k][j][i] = (pG->U[k][j][i].B2c*pG->U[k][j][i].M3 -
00912                             pG->U[k][j][i].B3c*pG->U[k][j][i].M2)
00913                               /pG->U[k][j][i].d;
00914         emf2_cc[k][j][i] = (pG->U[k][j][i].B3c*pG->U[k][j][i].M1 -
00915                             pG->U[k][j][i].B1c*pG->U[k][j][i].M3)
00916                               /pG->U[k][j][i].d;
00917         emf3_cc[k][j][i] = (pG->U[k][j][i].B1c*pG->U[k][j][i].M2 -
00918                             pG->U[k][j][i].B2c*pG->U[k][j][i].M1)
00919                               /pG->U[k][j][i].d;
00920       }
00921     }
00922   }
00923   integrate_emf1_corner(pG);
00924   integrate_emf2_corner(pG);
00925   integrate_emf3_corner(pG);
00926 
00927 /*--- Step 4b ------------------------------------------------------------------
00928  * Update the interface magnetic fields using CT for a half time step.
00929  */
00930 
00931   for (k=kl+1; k<=ku-1; k++) {
00932     for (j=jl+1; j<=ju-1; j++) {
00933       for (i=il+1; i<=iu-1; i++) {
00934 #ifdef CYLINDRICAL
00935         q2 = hdt/(ri[i]*pG->dx2);
00936 #endif
00937         B1_x1Face[k][j][i] += q3*(emf2[k+1][j  ][i  ] - emf2[k][j][i]) -
00938                               q2*(emf3[k  ][j+1][i  ] - emf3[k][j][i]);
00939         B2_x2Face[k][j][i] += q1*(emf3[k  ][j  ][i+1] - emf3[k][j][i]) -
00940                               q3*(emf1[k+1][j  ][i  ] - emf1[k][j][i]);
00941 #ifdef CYLINDRICAL
00942         rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
00943         q2 = hdt/(r[i]*pG->dx2);
00944 #endif
00945         B3_x3Face[k][j][i] += q2*(    emf1[k  ][j+1][i  ] -     emf1[k][j][i]) -
00946                               q1*(rsf*emf2[k  ][j  ][i+1] - lsf*emf2[k][j][i]);
00947       }
00948 #ifdef CYLINDRICAL
00949       q2 = hdt/(ri[iu]*pG->dx2);
00950 #endif
00951       B1_x1Face[k][j][iu] += q3*(emf2[k+1][j  ][iu]-emf2[k][j][iu]) -
00952                                q2*(emf3[k  ][j+1][iu]-emf3[k][j][iu]);
00953     }
00954     for (i=il+1; i<=iu-1; i++) {
00955       B2_x2Face[k][ju][i] += q1*(emf3[k  ][ju][i+1]-emf3[k][ju][i]) -
00956                                q3*(emf1[k+1][ju][i  ]-emf1[k][ju][i]);
00957     }
00958   }
00959   for (j=jl+1; j<=ju-1; j++) {
00960     for (i=il+1; i<=iu-1; i++) {
00961 #ifdef CYLINDRICAL
00962       rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
00963       q2 = hdt/(r[i]*pG->dx2);
00964 #endif
00965       B3_x3Face[ku][j][i] += q2*(    emf1[ku][j+1][i  ] -     emf1[ku][j][i]) -
00966                              q1*(rsf*emf2[ku][j  ][i+1] - lsf*emf2[ku][j][i]);
00967     }
00968   }
00969 #endif /* MHD */
00970 
00971 /*=== STEP 5: Correct x1-interface states with transverse flux gradients =====*/
00972 
00973 /*--- Step 5a ------------------------------------------------------------------
00974  * Correct x1-interface states using x2-fluxes computed in Step 2d.
00975  * Since the fluxes come from an x2-sweep, (x,y,z) on RHS -> (z,x,y) on LHS 
00976  */
00977 
00978   for (k=kl+1; k<=ku-1; k++) {
00979     for (j=jl+1; j<=ju-1; j++) {
00980       for (i=il+1; i<=iu; i++) {
00981 #ifdef CYLINDRICAL
00982         q2 = hdt/(r[i-1]*pG->dx2);
00983 #endif
00984         Ul_x1Face[k][j][i].d -=q2*(x2Flux[k][j+1][i-1].d -x2Flux[k][j][i-1].d );
00985         Ul_x1Face[k][j][i].Mx-=q2*(x2Flux[k][j+1][i-1].Mz-x2Flux[k][j][i-1].Mz);
00986         Ul_x1Face[k][j][i].My-=q2*(x2Flux[k][j+1][i-1].Mx-x2Flux[k][j][i-1].Mx);
00987         Ul_x1Face[k][j][i].Mz-=q2*(x2Flux[k][j+1][i-1].My-x2Flux[k][j][i-1].My);
00988 #ifndef BAROTROPIC
00989         Ul_x1Face[k][j][i].E -=q2*(x2Flux[k][j+1][i-1].E -x2Flux[k][j][i-1].E );
00990 #endif /* BAROTROPIC */
00991 #ifdef MHD
00992 /* Update B3 */
00993         Ul_x1Face[k][j][i].Bz+=q2*0.5*
00994           ((emf1[k  ][j+1][i-1] - emf1[k  ][j][i-1]) +
00995            (emf1[k+1][j+1][i-1] - emf1[k+1][j][i-1]));
00996 #endif
00997 
00998 #ifdef CYLINDRICAL
00999         q2 = hdt/(r[i]*pG->dx2);
01000 #endif
01001         Ur_x1Face[k][j][i].d -=q2*(x2Flux[k][j+1][i  ].d -x2Flux[k][j][i  ].d );
01002         Ur_x1Face[k][j][i].Mx-=q2*(x2Flux[k][j+1][i  ].Mz-x2Flux[k][j][i  ].Mz);
01003         Ur_x1Face[k][j][i].My-=q2*(x2Flux[k][j+1][i  ].Mx-x2Flux[k][j][i  ].Mx);
01004         Ur_x1Face[k][j][i].Mz-=q2*(x2Flux[k][j+1][i  ].My-x2Flux[k][j][i  ].My);
01005 #ifndef BAROTROPIC
01006         Ur_x1Face[k][j][i].E -=q2*(x2Flux[k][j+1][i  ].E -x2Flux[k][j][i  ].E );
01007 #endif /* BAROTROPIC */
01008 #ifdef MHD
01009 /* Update B3 */
01010         Ur_x1Face[k][j][i].Bz+=q2*0.5*
01011           ((emf1[k  ][j+1][i] - emf1[k  ][j][i]) +
01012            (emf1[k+1][j+1][i] - emf1[k+1][j][i]));
01013 #endif
01014 #if (NSCALARS > 0)
01015         for (n=0; n<NSCALARS; n++) {
01016           Ul_x1Face[k][j][i].s[n] -=
01017              q2*(x2Flux[k][j+1][i-1].s[n] - x2Flux[k][j][i-1].s[n]);
01018           Ur_x1Face[k][j][i].s[n] -=
01019              q2*(x2Flux[k][j+1][i  ].s[n] - x2Flux[k][j][i  ].s[n]);
01020         }
01021 #endif
01022 
01023 /*--- Step 5b ------------------------------------------------------------------
01024  * Correct x1-interface states using x3-fluxes computed in Step 3d.
01025  * Since the fluxes come from an x3-sweep, (x,y,z) on RHS -> (y,z,x) on LHS
01026  */
01027 
01028         Ul_x1Face[k][j][i].d -=q3*(x3Flux[k+1][j][i-1].d -x3Flux[k][j][i-1].d );
01029         Ul_x1Face[k][j][i].Mx-=q3*(x3Flux[k+1][j][i-1].My-x3Flux[k][j][i-1].My);
01030         Ul_x1Face[k][j][i].My-=q3*(x3Flux[k+1][j][i-1].Mz-x3Flux[k][j][i-1].Mz);
01031         Ul_x1Face[k][j][i].Mz-=q3*(x3Flux[k+1][j][i-1].Mx-x3Flux[k][j][i-1].Mx);
01032 #ifndef BAROTROPIC
01033         Ul_x1Face[k][j][i].E -=q3*(x3Flux[k+1][j][i-1].E -x3Flux[k][j][i-1].E );
01034 #endif /* BAROTROPIC */
01035 #ifdef MHD
01036 /* Update B2 */
01037         Ul_x1Face[k][j][i].By-=q3*0.5*
01038           ((emf1[k+1][j  ][i-1] - emf1[k][j  ][i-1]) +
01039            (emf1[k+1][j+1][i-1] - emf1[k][j+1][i-1]));
01040 #endif
01041 
01042         Ur_x1Face[k][j][i].d -=q3*(x3Flux[k+1][j][i  ].d -x3Flux[k][j][i  ].d );
01043         Ur_x1Face[k][j][i].Mx-=q3*(x3Flux[k+1][j][i  ].My-x3Flux[k][j][i  ].My);
01044         Ur_x1Face[k][j][i].My-=q3*(x3Flux[k+1][j][i  ].Mz-x3Flux[k][j][i  ].Mz);
01045         Ur_x1Face[k][j][i].Mz-=q3*(x3Flux[k+1][j][i  ].Mx-x3Flux[k][j][i  ].Mx);
01046 #ifndef BAROTROPIC
01047         Ur_x1Face[k][j][i].E -=q3*(x3Flux[k+1][j][i  ].E -x3Flux[k][j][i  ].E );
01048 #endif /* BAROTROPIC */
01049 #ifdef MHD
01050 /* Update B2 */
01051         Ur_x1Face[k][j][i].By-=q3*0.5*
01052           ((emf1[k+1][j  ][i] - emf1[k][j  ][i]) +
01053            (emf1[k+1][j+1][i] - emf1[k][j+1][i]));
01054 #endif
01055 #if (NSCALARS > 0)
01056         for (n=0; n<NSCALARS; n++) {
01057           Ul_x1Face[k][j][i].s[n] -=
01058              q3*(x3Flux[k+1][j][i-1].s[n] - x3Flux[k][j][i-1].s[n]);
01059           Ur_x1Face[k][j][i].s[n] -=
01060              q3*(x3Flux[k+1][j][i  ].s[n] - x3Flux[k][j][i  ].s[n]);
01061         }
01062 #endif
01063       }
01064     }
01065   }
01066 
01067 /*--- Step 5c ------------------------------------------------------------------
01068  * Add the "MHD source terms" from the x2- and x3-flux-gradients to the
01069  * conservative variables on the x1Face.  Limiting is used as in GS (2007)
01070  */
01071 
01072 #ifdef MHD
01073   for (k=kl+1; k<=ku-1; k++) {
01074     for (j=jl+1; j<=ju-1; j++) {
01075       for (i=il+1; i<=iu; i++) {
01076 #ifdef CYLINDRICAL
01077         rsf = ri[i]/r[i-1];  lsf = ri[i-1]/r[i-1];
01078         dx2i = 1.0/(r[i-1]*pG->dx2);
01079 #endif
01080         db1 = (rsf*pG->B1i[k  ][j  ][i  ] - lsf*pG->B1i[k][j][i-1])*dx1i;
01081         db2 = (    pG->B2i[k  ][j+1][i-1] -     pG->B2i[k][j][i-1])*dx2i;
01082         db3 = (    pG->B3i[k+1][j  ][i-1] -     pG->B3i[k][j][i-1])*dx3i;
01083         B1 = pG->U[k][j][i-1].B1c;
01084         B2 = pG->U[k][j][i-1].B2c;
01085         B3 = pG->U[k][j][i-1].B3c;
01086         V2 = pG->U[k][j][i-1].M2/pG->U[k][j][i-1].d;
01087         V3 = pG->U[k][j][i-1].M3/pG->U[k][j][i-1].d;
01088 
01089 /* Calculate mdb2 = min_mod(-db1,db2) */
01090         if(db1 > 0.0 && db2 < 0.0){
01091           mdb2 = db2 > -db1 ? db2 : -db1;
01092         }
01093         else if(db1 < 0.0 && db2 > 0.0){
01094           mdb2 = db2 < -db1 ? db2 : -db1;
01095         }
01096         else mdb2 = 0.0;
01097 
01098 /* Calculate mdb3 = min_mod(-db1,db3) */
01099         if(db1 > 0.0 && db3 < 0.0){
01100           mdb3 = db3 > -db1 ? db3 : -db1;
01101         }
01102         else if(db1 < 0.0 && db3 > 0.0){
01103           mdb3 = db3 < -db1 ? db3 : -db1;
01104         }
01105         else mdb3 = 0.0;
01106 
01107         Ul_x1Face[k][j][i].Mx += hdt*B1*db1;
01108         Ul_x1Face[k][j][i].My += hdt*B2*db1;
01109         Ul_x1Face[k][j][i].Mz += hdt*B3*db1;
01110         Ul_x1Face[k][j][i].By += hdt*V2*(-mdb3);
01111         Ul_x1Face[k][j][i].Bz += hdt*V3*(-mdb2);
01112 #ifndef BAROTROPIC
01113         Ul_x1Face[k][j][i].E  += hdt*(B2*V2*(-mdb3) + B3*V3*(-mdb2) );
01114 #endif /* BAROTROPIC */
01115 
01116 #ifdef CYLINDRICAL
01117         rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
01118         dx2i = 1.0/(r[i]*pG->dx2);
01119 #endif
01120         db1 = (rsf*pG->B1i[k  ][j  ][i+1] - lsf*pG->B1i[k][j][i])*dx1i;
01121         db2 = (    pG->B2i[k  ][j+1][i  ] -     pG->B2i[k][j][i])*dx2i;
01122         db3 = (    pG->B3i[k+1][j  ][i  ] -     pG->B3i[k][j][i])*dx3i;
01123         B1 = pG->U[k][j][i].B1c;
01124         B2 = pG->U[k][j][i].B2c;
01125         B3 = pG->U[k][j][i].B3c;
01126         V2 = pG->U[k][j][i].M2/pG->U[k][j][i].d;
01127         V3 = pG->U[k][j][i].M3/pG->U[k][j][i].d;
01128 
01129 /* Calculate mdb2 = min_mod(-db1,db2) */
01130         if(db1 > 0.0 && db2 < 0.0){
01131           mdb2 = db2 > -db1 ? db2 : -db1;
01132         }
01133         else if(db1 < 0.0 && db2 > 0.0){
01134           mdb2 = db2 < -db1 ? db2 : -db1;
01135         }
01136         else mdb2 = 0.0;
01137 
01138 /* Calculate mdb3 = min_mod(-db1,db3) */
01139         if(db1 > 0.0 && db3 < 0.0){
01140           mdb3 = db3 > -db1 ? db3 : -db1;
01141         }
01142         else if(db1 < 0.0 && db3 > 0.0){
01143           mdb3 = db3 < -db1 ? db3 : -db1;
01144         }
01145         else mdb3 = 0.0;
01146 
01147         Ur_x1Face[k][j][i].Mx += hdt*B1*db1;
01148         Ur_x1Face[k][j][i].My += hdt*B2*db1;
01149         Ur_x1Face[k][j][i].Mz += hdt*B3*db1;
01150         Ur_x1Face[k][j][i].By += hdt*V2*(-mdb3);
01151         Ur_x1Face[k][j][i].Bz += hdt*V3*(-mdb2);
01152 #ifndef BAROTROPIC
01153         Ur_x1Face[k][j][i].E  += hdt*(B2*V2*(-mdb3) + B3*V3*(-mdb2) );
01154 #endif /* BAROTROPIC */
01155       }
01156     }
01157   }
01158 #endif /* MHD */
01159 
01160 /*--- Step 5d ------------------------------------------------------------------
01161  * Add source terms for a static gravitational potential arising from x2-Flux
01162  * and x3-Flux gradients.  To improve conservation of total energy, average
01163  * the energy source term computed at cell faces.
01164  *    S_{M} = -(\rho) Grad(Phi);   S_{E} = -(\rho v) Grad{Phi}
01165  */
01166 
01167   if (StaticGravPot != NULL){
01168   for (k=kl+1; k<=ku-1; k++) {
01169     for (j=jl+1; j<=ju-1; j++) {
01170       for (i=il+1; i<=iu; i++) {
01171         cc_pos(pG,i,j,k,&x1,&x2,&x3);
01172         phic = (*StaticGravPot)(x1, x2             ,x3);
01173         phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
01174         phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
01175 
01176 /* correct right states; x2 and x3 gradients */
01177 #ifdef CYLINDRICAL
01178         q2 = hdt/(r[i]*pG->dx2);
01179 #endif
01180         Ur_x1Face[k][j][i].My -= q2*(phir-phil)*pG->U[k][j][i].d;
01181 #ifndef BAROTROPIC
01182         Ur_x1Face[k][j][i].E -= q2*(x2Flux[k][j  ][i  ].d*(phic - phil)
01183                                   + x2Flux[k][j+1][i  ].d*(phir - phic));
01184 #endif
01185 
01186         phir = (*StaticGravPot)(x1,x2,(x3+0.5*pG->dx3));
01187         phil = (*StaticGravPot)(x1,x2,(x3-0.5*pG->dx3));
01188         
01189         Ur_x1Face[k][j][i].Mz -= q3*(phir-phil)*pG->U[k][j][i].d;
01190 #ifndef BAROTROPIC
01191         Ur_x1Face[k][j][i].E -= q3*(x3Flux[k  ][j][i  ].d*(phic - phil)
01192                                   + x3Flux[k+1][j][i  ].d*(phir - phic));
01193 #endif
01194 
01195 /* correct left states; x2 and x3 gradients */
01196         phic = (*StaticGravPot)((x1-pG->dx1), x2             ,x3);
01197         phir = (*StaticGravPot)((x1-pG->dx1),(x2+0.5*pG->dx2),x3);
01198         phil = (*StaticGravPot)((x1-pG->dx1),(x2-0.5*pG->dx2),x3);
01199 
01200 #ifdef CYLINDRICAL
01201         q2 = hdt/(r[i-1]*pG->dx2);
01202 #endif
01203         Ul_x1Face[k][j][i].My -= q2*(phir-phil)*pG->U[k][j][i-1].d;
01204 #ifndef BAROTROPIC
01205         Ul_x1Face[k][j][i].E -= q2*(x2Flux[k][j  ][i-1].d*(phic - phil)
01206                                   + x2Flux[k][j+1][i-1].d*(phir - phic));
01207 #endif
01208 
01209         phir = (*StaticGravPot)((x1-pG->dx1),x2,(x3+0.5*pG->dx3));
01210         phil = (*StaticGravPot)((x1-pG->dx1),x2,(x3-0.5*pG->dx3));
01211         
01212         Ul_x1Face[k][j][i].Mz -= q3*(phir-phil)*pG->U[k][j][i-1].d;
01213 #ifndef BAROTROPIC
01214         Ul_x1Face[k][j][i].E -= q3*(x3Flux[k  ][j][i-1].d*(phic - phil)
01215                                   + x3Flux[k+1][j][i-1].d*(phir - phic));
01216 #endif
01217       }
01218     }
01219   }}
01220 
01221 /*--- Step 5d (cont) -----------------------------------------------------------
01222  * Add source terms for self gravity arising from x2-Flux and x3-Flux gradients
01223  *    S_{M} = -(\rho) Grad(Phi);   S_{E} = -(\rho v) Grad{Phi}
01224  */
01225 
01226 #ifdef SELF_GRAVITY
01227   for (k=kl+1; k<=ku-1; k++) {
01228     for (j=jl+1; j<=ju-1; j++) {
01229       for (i=il+1; i<=iu; i++) {
01230         phic = pG->Phi[k][j][i];
01231         phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j+1][i]);
01232         phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j-1][i]);
01233 
01234 /* correct right states; x2 and x3 gradients */
01235         Ur_x1Face[k][j][i].My -= q2*(phir-phil)*pG->U[k][j][i].d;
01236 #ifndef BAROTROPIC
01237         Ur_x1Face[k][j][i].E -= q2*(x2Flux[k][j  ][i  ].d*(phic - phil)
01238                                   + x2Flux[k][j+1][i  ].d*(phir - phic));
01239 #endif
01240 
01241         phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k+1][j][i]);
01242         phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k-1][j][i]);
01243 
01244         Ur_x1Face[k][j][i].Mz -= q3*(phir-phil)*pG->U[k][j][i].d;
01245 #ifndef BAROTROPIC
01246         Ur_x1Face[k][j][i].E -= q3*(x3Flux[k  ][j][i  ].d*(phic - phil)
01247                                   + x3Flux[k+1][j][i  ].d*(phir - phic));
01248 #endif
01249 
01250 /* correct left states; x2 and x3 gradients */
01251         phic = pG->Phi[k][j][i-1];
01252         phir = 0.5*(pG->Phi[k][j][i-1] + pG->Phi[k][j+1][i-1]);
01253         phil = 0.5*(pG->Phi[k][j][i-1] + pG->Phi[k][j-1][i-1]);
01254 
01255         Ul_x1Face[k][j][i].My -= q2*(phir-phil)*pG->U[k][j][i-1].d;
01256 #ifndef BAROTROPIC
01257         Ul_x1Face[k][j][i].E -= q2*(x2Flux[k][j  ][i-1].d*(phic - phil)
01258                                   + x2Flux[k][j+1][i-1].d*(phir - phic));
01259 #endif
01260 
01261         phir = 0.5*(pG->Phi[k][j][i-1] + pG->Phi[k+1][j][i-1]);
01262         phil = 0.5*(pG->Phi[k][j][i-1] + pG->Phi[k-1][j][i-1]);
01263 
01264         Ul_x1Face[k][j][i].Mz -= q3*(phir-phil)*pG->U[k][j][i-1].d;
01265 #ifndef BAROTROPIC
01266         Ul_x1Face[k][j][i].E -= q3*(x3Flux[k  ][j][i-1].d*(phic - phil)
01267                                   + x3Flux[k+1][j][i-1].d*(phir - phic));
01268 #endif
01269       }
01270     }
01271   }
01272 #endif /* SELF_GRAVITY */
01273 
01274 
01275 /*=== STEP 6: Correct x2-interface states with transverse flux gradients =====*/
01276 
01277 /*--- Step 6a ------------------------------------------------------------------
01278  * Correct x2-interface states using x1-fluxes computed in Step 1d.
01279  * Since the fluxes come from an x1-sweep, (x,y,z) on RHS -> (y,z,x) on LHS
01280  */
01281 
01282   for (k=kl+1; k<=ku-1; k++) {
01283     for (j=jl+1; j<=ju; j++) {
01284       for (i=il+1; i<=iu-1; i++) {
01285 #ifdef CYLINDRICAL
01286         rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
01287 #endif
01288         Ul_x2Face[k][j][i].d -=q1*(rsf*x1Flux[k][j-1][i+1].d -lsf*x1Flux[k][j-1][i].d );
01289         Ul_x2Face[k][j][i].Mx-=q1*(SQR(rsf)*x1Flux[k][j-1][i+1].My-SQR(lsf)*x1Flux[k][j-1][i].My);
01290         Ul_x2Face[k][j][i].My-=q1*(rsf*x1Flux[k][j-1][i+1].Mz-lsf*x1Flux[k][j-1][i].Mz);
01291         Ul_x2Face[k][j][i].Mz-=q1*(rsf*x1Flux[k][j-1][i+1].Mx-lsf*x1Flux[k][j-1][i].Mx);
01292 #ifndef BAROTROPIC
01293         Ul_x2Face[k][j][i].E -=q1*(rsf*x1Flux[k][j-1][i+1].E -lsf*x1Flux[k][j-1][i].E );
01294 #endif /* BAROTROPIC */
01295 #ifdef MHD
01296 /* Update B3 */
01297         Ul_x2Face[k][j][i].By-=q1*0.5*
01298           ((rsf*emf2[k  ][j-1][i+1] - lsf*emf2[k  ][j-1][i]) + 
01299            (rsf*emf2[k+1][j-1][i+1] - lsf*emf2[k+1][j-1][i]));
01300 #endif
01301 
01302         Ur_x2Face[k][j][i].d -=q1*(rsf*x1Flux[k][j  ][i+1].d -lsf*x1Flux[k][j  ][i].d );
01303         Ur_x2Face[k][j][i].Mx-=q1*(SQR(rsf)*x1Flux[k][j  ][i+1].My-SQR(lsf)*x1Flux[k][j  ][i].My);
01304         Ur_x2Face[k][j][i].My-=q1*(rsf*x1Flux[k][j  ][i+1].Mz-lsf*x1Flux[k][j  ][i].Mz);
01305         Ur_x2Face[k][j][i].Mz-=q1*(rsf*x1Flux[k][j  ][i+1].Mx-lsf*x1Flux[k][j  ][i].Mx);
01306 #ifndef BAROTROPIC
01307         Ur_x2Face[k][j][i].E -=q1*(rsf*x1Flux[k][j  ][i+1].E -lsf*x1Flux[k][j  ][i].E );
01308 #endif /* BAROTROPIC */
01309 #ifdef MHD
01310 /* Update B3 */
01311         Ur_x2Face[k][j][i].By-=q1*0.5*
01312           ((rsf*emf2[k  ][j][i+1] - lsf*emf2[k  ][j][i]) + 
01313            (rsf*emf2[k+1][j][i+1] - lsf*emf2[k+1][j][i]));
01314 #endif
01315 #if (NSCALARS > 0)
01316         for (n=0; n<NSCALARS; n++) {
01317           Ul_x2Face[k][j][i].s[n] -=
01318              q1*(rsf*x1Flux[k][j-1][i+1].s[n] - lsf*x1Flux[k][j-1][i].s[n]);
01319           Ur_x2Face[k][j][i].s[n] -=
01320              q1*(rsf*x1Flux[k][j  ][i+1].s[n] - lsf*x1Flux[k][j  ][i].s[n]);
01321         }
01322 #endif
01323 
01324 /*--- Step 6b ------------------------------------------------------------------
01325  * Correct x2-interface states using x3-fluxes computed in Step 3d.
01326  * Since the fluxes come from an x3-sweep, (x,y,z) on RHS -> (z,x,y) on LHS 
01327  */
01328 
01329         Ul_x2Face[k][j][i].d -=q3*(x3Flux[k+1][j-1][i].d -x3Flux[k][j-1][i].d );
01330         Ul_x2Face[k][j][i].Mx-=q3*(x3Flux[k+1][j-1][i].Mz-x3Flux[k][j-1][i].Mz);
01331         Ul_x2Face[k][j][i].My-=q3*(x3Flux[k+1][j-1][i].Mx-x3Flux[k][j-1][i].Mx);
01332         Ul_x2Face[k][j][i].Mz-=q3*(x3Flux[k+1][j-1][i].My-x3Flux[k][j-1][i].My);
01333 #ifndef BAROTROPIC
01334         Ul_x2Face[k][j][i].E -=q3*(x3Flux[k+1][j-1][i].E -x3Flux[k][j-1][i].E );
01335 #endif /* BAROTROPIC */
01336 #ifdef MHD
01337 /* Update B1 */
01338         Ul_x2Face[k][j][i].Bz+=q3*0.5*
01339           (lsf*(emf2[k+1][j-1][i  ] - emf2[k][j-1][i  ]) +
01340            rsf*(emf2[k+1][j-1][i+1] - emf2[k][j-1][i+1]));
01341 #endif
01342 
01343         Ur_x2Face[k][j][i].d -=q3*(x3Flux[k+1][j  ][i].d -x3Flux[k][j  ][i].d );
01344         Ur_x2Face[k][j][i].Mx-=q3*(x3Flux[k+1][j  ][i].Mz-x3Flux[k][j  ][i].Mz);
01345         Ur_x2Face[k][j][i].My-=q3*(x3Flux[k+1][j  ][i].Mx-x3Flux[k][j  ][i].Mx);
01346         Ur_x2Face[k][j][i].Mz-=q3*(x3Flux[k+1][j  ][i].My-x3Flux[k][j  ][i].My);
01347 #ifndef BAROTROPIC
01348         Ur_x2Face[k][j][i].E -=q3*(x3Flux[k+1][j  ][i].E -x3Flux[k][j  ][i].E );
01349 #endif /* BAROTROPIC */
01350 #ifdef MHD
01351 /* Update B1 */
01352         Ur_x2Face[k][j][i].Bz+=q3*0.5*
01353           (lsf*(emf2[k+1][j][i  ] - emf2[k][j][i  ]) +
01354            rsf*(emf2[k+1][j][i+1] - emf2[k][j][i+1]));
01355 #endif
01356 #if (NSCALARS > 0)
01357         for (n=0; n<NSCALARS; n++) {
01358           Ul_x2Face[k][j][i].s[n] -=
01359              q3*(x3Flux[k+1][j-1][i].s[n] - x3Flux[k][j-1][i].s[n]);
01360           Ur_x2Face[k][j][i].s[n] -=
01361              q3*(x3Flux[k+1][j  ][i].s[n] - x3Flux[k][j  ][i].s[n]);
01362         }
01363 #endif
01364       }
01365     }
01366   }
01367 
01368 /*--- Step 6c ------------------------------------------------------------------
01369  * Add the "MHD source terms" from the x1- and x3-flux-gradients to the
01370  * conservative variables on the x2Face.  Limiting is used as in GS (2007)
01371  */
01372 
01373 #ifdef MHD
01374   for (k=kl+1; k<=ku-1; k++) {
01375     for (j=jl+1; j<=ju; j++) {
01376       for (i=il+1; i<=iu-1; i++) {
01377 #ifdef CYLINDRICAL
01378         rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
01379 #endif
01380         db1 = (rsf*pG->B1i[k  ][j-1][i+1] - lsf*pG->B1i[k][j-1][i])*dx1i;
01381         db2 = (    pG->B2i[k  ][j  ][i  ] -     pG->B2i[k][j-1][i])*dx2i;
01382         db3 = (    pG->B3i[k+1][j-1][i  ] -     pG->B3i[k][j-1][i])*dx3i;
01383         B1 = pG->U[k][j-1][i].B1c;
01384         B2 = pG->U[k][j-1][i].B2c;
01385         B3 = pG->U[k][j-1][i].B3c;
01386         V1 = pG->U[k][j-1][i].M1/pG->U[k][j-1][i].d;
01387         V3 = pG->U[k][j-1][i].M3/pG->U[k][j-1][i].d;
01388 
01389 /* Calculate mdb1 = min_mod(-db2,db1) */
01390         if(db2 > 0.0 && db1 < 0.0){
01391           mdb1 = db1 > -db2 ? db1 : -db2;
01392         }
01393         else if(db2 < 0.0 && db1 > 0.0){
01394           mdb1 = db1 < -db2 ? db1 : -db2;
01395         }
01396         else mdb1 = 0.0;
01397 
01398 /* Calculate mdb3 = min_mod(-db2,db3) */
01399         if(db2 > 0.0 && db3 < 0.0){
01400           mdb3 = db3 > -db2 ? db3 : -db2;
01401         }
01402         else if(db2 < 0.0 && db3 > 0.0){
01403           mdb3 = db3 < -db2 ? db3 : -db2;
01404         }
01405         else mdb3 = 0.0;
01406 
01407         Ul_x2Face[k][j][i].Mz += hdt*B1*db2;
01408         Ul_x2Face[k][j][i].Mx += hdt*B2*db2;
01409         Ul_x2Face[k][j][i].My += hdt*B3*db2;
01410         Ul_x2Face[k][j][i].By += hdt*V3*(-mdb1);
01411         Ul_x2Face[k][j][i].Bz += hdt*V1*(-mdb3);
01412 #ifndef BAROTROPIC
01413         Ul_x2Face[k][j][i].E  += hdt*(B3*V3*(-mdb1) + B1*V1*(-mdb3) );
01414 #endif /* BAROTROPIC */
01415 
01416         db1 = (rsf*pG->B1i[k  ][j  ][i+1] - lsf*pG->B1i[k][j][i])*dx1i;
01417         db2 = (    pG->B2i[k  ][j+1][i  ] -     pG->B2i[k][j][i])*dx2i;
01418         db3 = (    pG->B3i[k+1][j  ][i  ] -     pG->B3i[k][j][i])*dx3i;
01419         B1 = pG->U[k][j][i].B1c;
01420         B2 = pG->U[k][j][i].B2c;
01421         B3 = pG->U[k][j][i].B3c;
01422         V1 = pG->U[k][j][i].M1/pG->U[k][j][i].d;
01423         V3 = pG->U[k][j][i].M3/pG->U[k][j][i].d;
01424 
01425 /* Calculate mdb1 = min_mod(-db2,db1) */
01426         if(db2 > 0.0 && db1 < 0.0){
01427           mdb1 = db1 > -db2 ? db1 : -db2;
01428         }
01429         else if(db2 < 0.0 && db1 > 0.0){
01430           mdb1 = db1 < -db2 ? db1 : -db2;
01431         }
01432         else mdb1 = 0.0;
01433 
01434 /* Calculate mdb3 = min_mod(-db2,db3) */
01435         if(db2 > 0.0 && db3 < 0.0){
01436           mdb3 = db3 > -db2 ? db3 : -db2;
01437         }
01438         else if(db2 < 0.0 && db3 > 0.0){
01439           mdb3 = db3 < -db2 ? db3 : -db2;
01440         }
01441         else mdb3 = 0.0;
01442 
01443         Ur_x2Face[k][j][i].Mz += hdt*B1*db2;
01444         Ur_x2Face[k][j][i].Mx += hdt*B2*db2;
01445         Ur_x2Face[k][j][i].My += hdt*B3*db2;
01446         Ur_x2Face[k][j][i].By += hdt*V3*(-mdb1);
01447         Ur_x2Face[k][j][i].Bz += hdt*V1*(-mdb3);
01448 #ifndef BAROTROPIC
01449         Ur_x2Face[k][j][i].E  += hdt*(B3*V3*(-mdb1) + B1*V1*(-mdb3) );
01450 #endif /* BAROTROPIC */
01451       }
01452     }
01453   }
01454 #endif /* MHD */
01455 
01456 /*--- Step 6d ------------------------------------------------------------------
01457  * Add source terms for a static gravitational potential arising from x1-Flux
01458  * and x3-Flux gradients. To improve conservation of total energy,
01459  * average the energy source term computed at cell faces.
01460  *    S_{M} = -(\rho) Grad(Phi);   S_{E} = -(\rho v) Grad{Phi}
01461  */
01462 
01463   if (StaticGravPot != NULL){
01464   for (k=kl+1; k<=ku-1; k++) {
01465     for (j=jl+1; j<=ju; j++) {
01466       for (i=il+1; i<=iu-1; i++) {
01467         cc_pos(pG,i,j,k,&x1,&x2,&x3);
01468         phic = (*StaticGravPot)((x1            ),x2,x3);
01469         phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
01470         phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
01471 
01472 /* correct right states; x1 and x3 gradients */
01473 #ifdef CYLINDRICAL
01474         g = (*x1GravAcc)(x1vc(pG,i),x2,x3);
01475 #ifdef FARGO
01476         g = g - x1vc(pG,i)*SQR((*OrbitalProfile)(x1vc(pG,i))); 
01477 #endif
01478         rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
01479         Ur_x2Face[k][j][i].Mz -= hdt*pG->U[k][j][i].d*g;
01480 #else
01481         Ur_x2Face[k][j][i].Mz -= q1*(phir-phil)*pG->U[k][j][i].d;
01482 #endif
01483 #ifndef BAROTROPIC
01484         Ur_x2Face[k][j][i].E -= q1*(lsf*x1Flux[k][j  ][i  ].d*(phic - phil)
01485                                   + rsf*x1Flux[k][j  ][i+1].d*(phir - phic));
01486 #endif
01487 
01488         phir = (*StaticGravPot)(x1,x2,(x3+0.5*pG->dx3));
01489         phil = (*StaticGravPot)(x1,x2,(x3-0.5*pG->dx3));
01490 
01491         Ur_x2Face[k][j][i].My -= q3*(phir-phil)*pG->U[k][j][i].d;
01492 #ifndef BAROTROPIC
01493         Ur_x2Face[k][j][i].E -= q3*(x3Flux[k  ][j  ][i].d*(phic - phil)
01494                                   + x3Flux[k+1][j  ][i].d*(phir - phic));
01495 #endif
01496 
01497 /* correct left states; x1 and x3 gradients */
01498         phic = (*StaticGravPot)((x1            ),(x2-pG->dx2),x3);
01499         phir = (*StaticGravPot)((x1+0.5*pG->dx1),(x2-pG->dx2),x3);
01500         phil = (*StaticGravPot)((x1-0.5*pG->dx1),(x2-pG->dx2),x3);
01501 
01502 #ifdef CYLINDRICAL
01503         g = (*x1GravAcc)(x1vc(pG,i),(x2-pG->dx2),x3);
01504 #ifdef FARGO
01505         g = g - x1vc(pG,i)*SQR((*OrbitalProfile)(x1vc(pG,i))); 
01506 #endif
01507 
01508         Ul_x2Face[k][j][i].Mz -= hdt*pG->U[k][j-1][i].d*g;
01509 #else
01510         Ul_x2Face[k][j][i].Mz -= q1*(phir-phil)*pG->U[k][j-1][i].d;
01511 #endif
01512 #ifndef BAROTROPIC
01513         Ul_x2Face[k][j][i].E -= q1*(lsf*x1Flux[k][j-1][i  ].d*(phic - phil)
01514                                   + rsf*x1Flux[k][j-1][i+1].d*(phir - phic));
01515 #endif
01516         phir = (*StaticGravPot)(x1,(x2-pG->dx2),(x3+0.5*pG->dx3));
01517         phil = (*StaticGravPot)(x1,(x2-pG->dx2),(x3-0.5*pG->dx3));
01518 
01519         Ul_x2Face[k][j][i].My -= q3*(phir-phil)*pG->U[k][j-1][i].d;
01520 #ifndef BAROTROPIC
01521         Ul_x2Face[k][j][i].E -= q3*(x3Flux[k  ][j-1][i].d*(phic - phil)
01522                                   + x3Flux[k+1][j-1][i].d*(phir - phic));
01523 #endif
01524       }
01525     }
01526   }}
01527 
01528 /*--- Step 6d (cont) -----------------------------------------------------------
01529  * Add source terms for self gravity arising from x1-Flux and x3-Flux gradients
01530  *    S_{M} = -(\rho) Grad(Phi);   S_{E} = -(\rho v) Grad{Phi}
01531  */
01532 
01533 #ifdef SELF_GRAVITY
01534   for (k=kl+1; k<=ku-1; k++) {
01535     for (j=jl+1; j<=ju; j++) {
01536       for (i=il+1; i<=iu-1; i++) {
01537         phic = pG->Phi[k][j][i];
01538         phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j][i+1]);
01539         phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j][i-1]);
01540 
01541 /* correct right states; x1 and x3 gradients */
01542         Ur_x2Face[k][j][i].Mz -= q1*(phir-phil)*pG->U[k][j][i].d;
01543 #ifndef BAROTROPIC
01544         Ur_x2Face[k][j][i].E -= q1*(x1Flux[k][j][i  ].d*(phic - phil)
01545                                   + x1Flux[k][j][i+1].d*(phir - phic));
01546 #endif
01547 
01548         phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k+1][j][i]);
01549         phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k-1][j][i]);
01550 
01551         Ur_x2Face[k][j][i].My -= q3*(phir-phil)*pG->U[k][j][i].d;
01552 #ifndef BAROTROPIC
01553         Ur_x2Face[k][j][i].E -= q3*(x3Flux[k  ][j][i].d*(phic - phil)
01554                                   + x3Flux[k+1][j][i].d*(phir - phic));
01555 #endif
01556 
01557 /* correct left states; x1 and x3 gradients */
01558         phic = pG->Phi[k][j-1][i];
01559         phir = 0.5*(pG->Phi[k][j-1][i] + pG->Phi[k][j-1][i+1]);
01560         phil = 0.5*(pG->Phi[k][j-1][i] + pG->Phi[k][j-1][i-1]);
01561 
01562         Ul_x2Face[k][j][i].Mz -= q1*(phir-phil)*pG->U[k][j-1][i].d;
01563 #ifndef BAROTROPIC
01564         Ul_x2Face[k][j][i].E -= q1*(x1Flux[k][j-1][i  ].d*(phic - phil)
01565                                   + x1Flux[k][j-1][i+1].d*(phir - phic));
01566 #endif
01567         phir = 0.5*(pG->Phi[k][j-1][i] + pG->Phi[k+1][j-1][i]);
01568         phil = 0.5*(pG->Phi[k][j-1][i] + pG->Phi[k-1][j-1][i]);
01569 
01570         Ul_x2Face[k][j][i].My -= q3*(phir-phil)*pG->U[k][j-1][i].d;
01571 #ifndef BAROTROPIC
01572         Ul_x2Face[k][j][i].E -= q3*(x3Flux[k  ][j-1][i].d*(phic - phil)
01573                                   + x3Flux[k+1][j-1][i].d*(phir - phic));
01574 #endif
01575       }
01576     }
01577   }
01578 #endif /* SELF_GRAVITY */
01579 
01580 /*--- Step 6d (cont) -----------------------------------------------------------
01581  * Add source terms for shearing box (Coriolis forces) for 0.5*dt arising from
01582  * x1-Flux gradient.  The tidal gravity terms are added via ShearingBoxPot
01583  *    Vx source term is (dt/2)( 2 Omega_0 Vy)
01584  *    Vy source term is (dt/2)(-2 Omega_0 Vx)
01585  *    Vy source term is (dt/2)((q-2) Omega_0 Vx) (with FARGO)
01586  */
01587 
01588 #ifdef SHEARING_BOX
01589   if (ShearingBoxPot != NULL){
01590   for (k=kl+1; k<=ku-1; k++) {
01591     for (j=jl+1; j<=ju; j++) {
01592       for (i=il+1; i<=iu-1; i++) {
01593         cc_pos(pG,i,j,k,&x1,&x2,&x3);
01594         phic = (*ShearingBoxPot)((x1            ),x2,x3);
01595         phir = (*ShearingBoxPot)((x1+0.5*pG->dx1),x2,x3);
01596         phil = (*ShearingBoxPot)((x1-0.5*pG->dx1),x2,x3);
01597 
01598 /* correct right states; x1 and x3 gradients */
01599         Ur_x2Face[k][j][i].Mz -= q1*(phir-phil)*pG->U[k][j][i].d;
01600 #ifndef BAROTROPIC
01601         Ur_x2Face[k][j][i].E -= q1*(x1Flux[k][j  ][i  ].d*(phic - phil)
01602                                   + x1Flux[k][j  ][i+1].d*(phir - phic));
01603 #endif
01604 
01605         phir = (*ShearingBoxPot)(x1,x2,(x3+0.5*pG->dx3));
01606         phil = (*ShearingBoxPot)(x1,x2,(x3-0.5*pG->dx3));
01607 
01608         Ur_x2Face[k][j][i].My -= q3*(phir-phil)*pG->U[k][j][i].d;
01609 #ifndef BAROTROPIC
01610         Ur_x2Face[k][j][i].E -= q3*(x3Flux[k  ][j  ][i].d*(phic - phil)
01611                                   + x3Flux[k+1][j  ][i].d*(phir - phic));
01612 #endif
01613 
01614 /* correct left states; x1 and x3 gradients */
01615         phic = (*ShearingBoxPot)((x1            ),(x2-pG->dx2),x3);
01616         phir = (*ShearingBoxPot)((x1+0.5*pG->dx1),(x2-pG->dx2),x3);
01617         phil = (*ShearingBoxPot)((x1-0.5*pG->dx1),(x2-pG->dx2),x3);
01618 
01619         Ul_x2Face[k][j][i].Mz -= q1*(phir-phil)*pG->U[k][j-1][i].d;
01620 #ifndef BAROTROPIC
01621         Ul_x2Face[k][j][i].E -= q1*(x1Flux[k][j-1][i  ].d*(phic - phil)
01622                                   + x1Flux[k][j-1][i+1].d*(phir - phic));
01623 #endif
01624         phir = (*ShearingBoxPot)(x1,(x2-pG->dx2),(x3+0.5*pG->dx3));
01625         phil = (*ShearingBoxPot)(x1,(x2-pG->dx2),(x3-0.5*pG->dx3));
01626 
01627         Ul_x2Face[k][j][i].My -= q3*(phir-phil)*pG->U[k][j-1][i].d;
01628 #ifndef BAROTROPIC
01629         Ul_x2Face[k][j][i].E -= q3*(x3Flux[k  ][j-1][i].d*(phic - phil)
01630                                   + x3Flux[k+1][j-1][i].d*(phir - phic));
01631 #endif
01632       }
01633     }
01634   }}
01635 
01636   for (k=kl+1; k<=ku-1; k++) {
01637     for (j=jl+1; j<=ju; j++) {
01638       for (i=il+1; i<=iu-1; i++) {
01639         Ur_x2Face[k][j][i].Mz += pG->dt*Omega_0*pG->U[k][j][i].M2;
01640         Ul_x2Face[k][j][i].Mz += pG->dt*Omega_0*pG->U[k][j-1][i].M2;
01641 #ifdef FARGO
01642         Ur_x2Face[k][j][i].Mx += hdt*(qshear-2.)*Omega_0*pG->U[k][j][i].M1;
01643         Ul_x2Face[k][j][i].Mx += hdt*(qshear-2.)*Omega_0*pG->U[k][j-1][i].M1;
01644 #else
01645         Ur_x2Face[k][j][i].Mx -= pG->dt*Omega_0*pG->U[k][j][i].M1;
01646         Ul_x2Face[k][j][i].Mx -= pG->dt*Omega_0*pG->U[k][j-1][i].M1;
01647 #endif
01648       }
01649     }
01650   }
01651 #endif /* SHEARING_BOX */
01652 
01653 #if defined(CYLINDRICAL) && defined(FARGO)
01654         for (k=kl+1; k<=ku-1; k++) {
01655                 for (j=jl+1; j<=ju; j++) {
01656                         for (i=il+1; i<=iu-1; i++) {
01657         Om = (*OrbitalProfile)(x1vc(pG,i));
01658                                 qshear = (*ShearProfile)(x1vc(pG,i));
01659                                 Ur_x2Face[k][j][i].Mz += pG->dt*Om*pG->U[k][j][i].M2;
01660                                 Ur_x2Face[k][j][i].Mx += hdt*(qshear-2.0)*Om*pG->U[k][j][i].M1;
01661 
01662                                 Ul_x2Face[k][j][i].Mz += pG->dt*Om*pG->U[k][j-1][i].M2;
01663                                 Ul_x2Face[k][j][i].Mx += hdt*(qshear-2.0)*Om*pG->U[k][j-1][i].M1;
01664                         }
01665                 }
01666         }
01667 #endif /* Cylindrical + Fargo */
01668 
01669 /*--- Step 6d (cont) -----------------------------------------------------------
01670  * ADD THE GEOMETRIC SOURCE-TERM IN THE X1-DIRECTION TO THE CORRECTED L/R 
01671  * STATES ON X2-FACES.  S_{M_R} = -(\rho V_\phi^2 - B_\phi^2)/R
01672  */
01673 #ifdef CYLINDRICAL
01674   for (k=kl+1; k<=ku-1; k++) {
01675     for (j=jl+1; j<=ju; j++) {
01676       for (i=il+1; i<=iu-1; i++) {
01677         Ur_x2Face[k][j][i].Mz += hdt*geom_src[k][j  ][i];
01678         Ul_x2Face[k][j][i].Mz += hdt*geom_src[k][j-1][i];
01679       }
01680     }
01681   }
01682 #endif /* CYLINDRICAL */
01683 
01684 /*=== STEP 7: Correct x3-interface states with transverse flux gradients =====*/
01685 
01686 /*--- Step 7a ------------------------------------------------------------------
01687  * Correct x3-interface states using x1-fluxes computed in Step 1d.
01688  * Since the fluxes come from an x1-sweep, (x,y,z) on RHS -> (z,x,y) on LHS 
01689  */
01690 
01691   for (k=kl+1; k<=ku; k++) {
01692     for (j=jl+1; j<=ju-1; j++) {
01693       for (i=il+1; i<=iu-1; i++) {
01694 #ifdef CYLINDRICAL
01695         rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
01696         q2 = hdt/(r[i]*pG->dx2);
01697 #endif
01698         Ul_x3Face[k][j][i].d -=q1*(rsf*x1Flux[k-1][j][i+1].d -lsf*x1Flux[k-1][j][i].d );
01699         Ul_x3Face[k][j][i].Mx-=q1*(rsf*x1Flux[k-1][j][i+1].Mz-lsf*x1Flux[k-1][j][i].Mz);
01700         Ul_x3Face[k][j][i].My-=q1*(rsf*x1Flux[k-1][j][i+1].Mx-lsf*x1Flux[k-1][j][i].Mx);
01701         Ul_x3Face[k][j][i].Mz-=q1*(SQR(rsf)*x1Flux[k-1][j][i+1].My-SQR(lsf)*x1Flux[k-1][j][i].My);
01702 #ifndef BAROTROPIC
01703         Ul_x3Face[k][j][i].E -=q1*(rsf*x1Flux[k-1][j][i+1].E -lsf*x1Flux[k-1][j][i].E );
01704 #endif /* BAROTROPIC */
01705 #ifdef MHD
01706 /* Update B2 */
01707         Ul_x3Face[k][j][i].Bz+=q1*0.5*
01708           ((emf3[k-1][j  ][i+1] - emf3[k-1][j  ][i]) +
01709            (emf3[k-1][j+1][i+1] - emf3[k-1][j+1][i]));
01710 #endif
01711 
01712         Ur_x3Face[k][j][i].d -=q1*(rsf*x1Flux[k  ][j][i+1].d -lsf*x1Flux[k  ][j][i].d );
01713         Ur_x3Face[k][j][i].Mx-=q1*(rsf*x1Flux[k  ][j][i+1].Mz-lsf*x1Flux[k  ][j][i].Mz);
01714         Ur_x3Face[k][j][i].My-=q1*(rsf*x1Flux[k  ][j][i+1].Mx-lsf*x1Flux[k  ][j][i].Mx);
01715         Ur_x3Face[k][j][i].Mz-=q1*(SQR(rsf)*x1Flux[k  ][j][i+1].My-SQR(lsf)*x1Flux[k  ][j][i].My);
01716 #ifndef BAROTROPIC
01717         Ur_x3Face[k][j][i].E -=q1*(rsf*x1Flux[k  ][j][i+1].E -lsf*x1Flux[k  ][j][i].E );
01718 #endif /* BAROTROPIC */
01719 #ifdef MHD
01720 /* Update B2 */
01721         Ur_x3Face[k][j][i].Bz+=q1*0.5*
01722           ((emf3[k][j  ][i+1] - emf3[k][j  ][i]) +
01723            (emf3[k][j+1][i+1] - emf3[k][j+1][i]));
01724 #endif
01725 #if (NSCALARS > 0)
01726         for (n=0; n<NSCALARS; n++) {
01727           Ul_x3Face[k][j][i].s[n] -=
01728              q1*(rsf*x1Flux[k-1][j][i+1].s[n] - lsf*x1Flux[k-1][j][i].s[n]);
01729           Ur_x3Face[k][j][i].s[n] -=
01730              q1*(rsf*x1Flux[k  ][j][i+1].s[n] - lsf*x1Flux[k  ][j][i].s[n]);
01731         }
01732 #endif
01733 
01734 /*--- Step 7b ------------------------------------------------------------------
01735  * Correct x3-interface states using x2-fluxes computed in Step 2d.
01736  * Since the fluxes come from an x2-sweep, (x,y,z) on RHS -> (y,z,x) on LHS 
01737  */
01738 
01739         Ul_x3Face[k][j][i].d -=q2*(x2Flux[k-1][j+1][i].d -x2Flux[k-1][j][i].d );
01740         Ul_x3Face[k][j][i].Mx-=q2*(x2Flux[k-1][j+1][i].My-x2Flux[k-1][j][i].My);
01741         Ul_x3Face[k][j][i].My-=q2*(x2Flux[k-1][j+1][i].Mz-x2Flux[k-1][j][i].Mz);
01742         Ul_x3Face[k][j][i].Mz-=q2*(x2Flux[k-1][j+1][i].Mx-x2Flux[k-1][j][i].Mx);
01743 #ifndef BAROTROPIC
01744         Ul_x3Face[k][j][i].E -=q2*(x2Flux[k-1][j+1][i].E -x2Flux[k-1][j][i].E );
01745 #endif /* BAROTROPIC */
01746 #ifdef MHD
01747 /* Update B1 */
01748         Ul_x3Face[k][j][i].By-=q2*0.5*
01749           ((emf3[k-1][j+1][i  ] - emf3[k-1][j][i  ]) +
01750            (emf3[k-1][j+1][i+1] - emf3[k-1][j][i+1]));
01751 #endif
01752 
01753         Ur_x3Face[k][j][i].d -=q2*(x2Flux[k  ][j+1][i].d -x2Flux[k  ][j][i].d );
01754         Ur_x3Face[k][j][i].Mx-=q2*(x2Flux[k  ][j+1][i].My-x2Flux[k  ][j][i].My);
01755         Ur_x3Face[k][j][i].My-=q2*(x2Flux[k  ][j+1][i].Mz-x2Flux[k  ][j][i].Mz);
01756         Ur_x3Face[k][j][i].Mz-=q2*(x2Flux[k  ][j+1][i].Mx-x2Flux[k  ][j][i].Mx);
01757 #ifndef BAROTROPIC
01758         Ur_x3Face[k][j][i].E -=q2*(x2Flux[k  ][j+1][i].E -x2Flux[k  ][j][i].E );
01759 #endif /* BAROTROPIC */
01760 #ifdef MHD
01761 /* Update B1 */
01762         Ur_x3Face[k][j][i].By-=q2*0.5*
01763           ((emf3[k][j+1][i  ] - emf3[k][j][i  ]) +
01764            (emf3[k][j+1][i+1] - emf3[k][j][i+1]));
01765 #endif
01766 #if (NSCALARS > 0)
01767         for (n=0; n<NSCALARS; n++) {
01768           Ul_x3Face[k][j][i].s[n] -=
01769              q2*(x2Flux[k-1][j+1][i].s[n] - x2Flux[k-1][j][i].s[n]);
01770           Ur_x3Face[k][j][i].s[n] -=
01771              q2*(x2Flux[k  ][j+1][i].s[n] - x2Flux[k  ][j][i].s[n]);
01772         }
01773 #endif
01774       }
01775     }
01776   }
01777 
01778 /*--- Step 7c ------------------------------------------------------------------
01779  * Add the "MHD source terms" from the x1- and x2-flux-gradients to the
01780  * conservative variables on the x3Face.  Limiting is used as in GS07.
01781  */
01782 
01783 #ifdef MHD
01784   for (k=kl+1; k<=ku; k++) {
01785     for (j=jl+1; j<=ju-1; j++) {
01786       for (i=il+1; i<=iu-1; i++) {
01787 #ifdef CYLINDRICAL
01788         rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
01789 #endif
01790         db1 = (rsf*pG->B1i[k-1][j  ][i+1] - lsf*pG->B1i[k-1][j][i])*dx1i;
01791         db2 = (    pG->B2i[k-1][j+1][i  ] -     pG->B2i[k-1][j][i])*dx2i;
01792         db3 = (    pG->B3i[k  ][j  ][i  ] -     pG->B3i[k-1][j][i])*dx3i;
01793         B1 = pG->U[k-1][j][i].B1c;
01794         B2 = pG->U[k-1][j][i].B2c;
01795         B3 = pG->U[k-1][j][i].B3c;
01796         V1 = pG->U[k-1][j][i].M1/pG->U[k-1][j][i].d;
01797         V2 = pG->U[k-1][j][i].M2/pG->U[k-1][j][i].d;
01798 
01799 /* Calculate mdb1 = min_mod(-db3,db1) */
01800         if(db3 > 0.0 && db1 < 0.0){
01801           mdb1 = db1 > -db3 ? db1 : -db3;
01802         }
01803         else if(db3 < 0.0 && db1 > 0.0){
01804           mdb1 = db1 < -db3 ? db1 : -db3;
01805         }
01806         else mdb1 = 0.0;
01807 
01808 /* Calculate mdb2 = min_mod(-db3,db2) */
01809         if(db3 > 0.0 && db2 < 0.0){
01810           mdb2 = db2 > -db3 ? db2 : -db3;
01811         }
01812         else if(db3 < 0.0 && db2 > 0.0){
01813           mdb2 = db2 < -db3 ? db2 : -db3;
01814         }
01815         else mdb2 = 0.0;
01816 
01817         Ul_x3Face[k][j][i].My += hdt*B1*db3;
01818         Ul_x3Face[k][j][i].Mz += hdt*B2*db3;
01819         Ul_x3Face[k][j][i].Mx += hdt*B3*db3;
01820         Ul_x3Face[k][j][i].By += hdt*V1*(-mdb2);
01821         Ul_x3Face[k][j][i].Bz += hdt*V2*(-mdb1);
01822 #ifndef BAROTROPIC
01823         Ul_x3Face[k][j][i].E  += hdt*(B1*V1*(-mdb2) + B2*V2*(-mdb1) );
01824 #endif /* BAROTROPIC */
01825 
01826         db1 = (rsf*pG->B1i[k  ][j  ][i+1] - lsf*pG->B1i[k][j][i])*dx1i;
01827         db2 = (    pG->B2i[k  ][j+1][i  ] -     pG->B2i[k][j][i])*dx2i;
01828         db3 = (    pG->B3i[k+1][j  ][i  ] -     pG->B3i[k][j][i])*dx3i;
01829         B1 = pG->U[k][j][i].B1c;
01830         B2 = pG->U[k][j][i].B2c;
01831         B3 = pG->U[k][j][i].B3c;
01832         V1 = pG->U[k][j][i].M1/pG->U[k][j][i].d;
01833         V2 = pG->U[k][j][i].M2/pG->U[k][j][i].d;
01834 
01835 /* Calculate mdb1 = min_mod(-db3,db1) */
01836         if(db3 > 0.0 && db1 < 0.0){
01837           mdb1 = db1 > -db3 ? db1 : -db3;
01838         }
01839         else if(db3 < 0.0 && db1 > 0.0){
01840           mdb1 = db1 < -db3 ? db1 : -db3;
01841         }
01842         else mdb1 = 0.0;
01843 
01844 /* Calculate mdb2 = min_mod(-db3,db2) */
01845         if(db3 > 0.0 && db2 < 0.0){
01846           mdb2 = db2 > -db3 ? db2 : -db3;
01847         }
01848         else if(db3 < 0.0 && db2 > 0.0){
01849           mdb2 = db2 < -db3 ? db2 : -db3;
01850         }
01851         else mdb2 = 0.0;
01852 
01853         Ur_x3Face[k][j][i].My += hdt*B1*db3;
01854         Ur_x3Face[k][j][i].Mz += hdt*B2*db3;
01855         Ur_x3Face[k][j][i].Mx += hdt*B3*db3;
01856         Ur_x3Face[k][j][i].By += hdt*V1*(-mdb2);
01857         Ur_x3Face[k][j][i].Bz += hdt*V2*(-mdb1);
01858 #ifndef BAROTROPIC
01859         Ur_x3Face[k][j][i].E  += hdt*(B1*V1*(-mdb2) + B2*V2*(-mdb1) );
01860 #endif /* BAROTROPIC */
01861       }
01862     }
01863   }
01864 #endif /* MHD */
01865 
01866 /*--- Step 7d ------------------------------------------------------------------
01867  * Add source terms for a static gravitational potential arising from x1-Flux
01868  * and x2-Flux gradients. To improve conservation of total energy,
01869  * average the energy source term computed at cell faces.
01870  *    S_{M} = -(\rho) Grad(Phi);   S_{E} = -(\rho v) Grad{Phi}
01871  */
01872 
01873   if (StaticGravPot != NULL){
01874   for (k=kl+1; k<=ku; k++) {
01875     for (j=jl+1; j<=ju-1; j++) {
01876       for (i=il+1; i<=iu-1; i++) {
01877         cc_pos(pG,i,j,k,&x1,&x2,&x3);
01878         phic = (*StaticGravPot)((x1            ),x2,x3);
01879         phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
01880         phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
01881 
01882 /* correct right states; x1 and x2 gradients */
01883 #ifdef CYLINDRICAL
01884         g = (*x1GravAcc)(x1vc(pG,i),x2,x3);
01885 #ifdef FARGO
01886         g = g - x1vc(pG,i)*SQR((*OrbitalProfile)(x1vc(pG,i)));
01887 #endif
01888         rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
01889         q2 = hdt/(r[i]*pG->dx2);
01890         Ur_x3Face[k][j][i].My -= hdt*pG->U[k][j][i].d*g;
01891 #else
01892         Ur_x3Face[k][j][i].My -= q1*(phir-phil)*pG->U[k][j][i].d;
01893 #endif
01894 #ifndef BAROTROPIC
01895         Ur_x3Face[k][j][i].E -= q1*(lsf*x1Flux[k  ][j][i  ].d*(phic - phil)
01896                                   + rsf*x1Flux[k  ][j][i+1].d*(phir - phic));
01897 #endif
01898 
01899         phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
01900         phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
01901 
01902         Ur_x3Face[k][j][i].Mz -= q2*(phir-phil)*pG->U[k][j][i].d;
01903 #ifndef BAROTROPIC
01904         Ur_x3Face[k][j][i].E -= q2*(x2Flux[k  ][j  ][i].d*(phic - phil)
01905                                   + x2Flux[k  ][j+1][i].d*(phir - phic));
01906 #endif
01907 
01908 /* correct left states; x1 and x2 gradients */
01909         phic = (*StaticGravPot)((x1            ),x2,(x3-pG->dx3));
01910         phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,(x3-pG->dx3));
01911         phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,(x3-pG->dx3));
01912 
01913 #ifdef CYLINDRICAL
01914         g = (*x1GravAcc)(x1vc(pG,i),x2,(x3-pG->dx3));
01915 #ifdef FARGO
01916         g = g - x1vc(pG,i)*SQR((*OrbitalProfile)(x1vc(pG,i)));
01917 #endif
01918         Ul_x3Face[k][j][i].My -= hdt*pG->U[k-1][j][i].d*g;
01919 #else
01920         Ul_x3Face[k][j][i].My -= q1*(phir-phil)*pG->U[k-1][j][i].d;
01921 #endif
01922 #ifndef BAROTROPIC
01923         Ul_x3Face[k][j][i].E -= q1*(lsf*x1Flux[k-1][j][i  ].d*(phic - phil)
01924                                   + rsf*x1Flux[k-1][j][i+1].d*(phir - phic));
01925 #endif
01926 
01927         phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),(x3-pG->dx3));
01928         phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),(x3-pG->dx3));
01929 
01930         Ul_x3Face[k][j][i].Mz -= q2*(phir-phil)*pG->U[k-1][j][i].d;
01931 #ifndef BAROTROPIC
01932         Ul_x3Face[k][j][i].E -= q2*(x2Flux[k-1][j  ][i].d*(phic - phil)
01933                                   + x2Flux[k-1][j+1][i].d*(phir - phic));
01934 #endif
01935       }
01936     }
01937   }}
01938 
01939 /*--- Step 7d (cont) -----------------------------------------------------------
01940  * Add source terms for self gravity arising from x1-Flux and x2-Flux gradients
01941  *    S_{M} = -(\rho) Grad(Phi);   S_{E} = -(\rho v) Grad{Phi}
01942  */
01943 
01944 #ifdef SELF_GRAVITY
01945   for (k=kl+1; k<=ku; k++) {
01946     for (j=jl+1; j<=ju-1; j++) {
01947       for (i=il+1; i<=iu-1; i++) {
01948         phic = pG->Phi[k][j][i];
01949         phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j][i+1]);
01950         phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j][i-1]);
01951 
01952 /* correct right states; x1 and x2 gradients */
01953         Ur_x3Face[k][j][i].My -= q1*(phir-phil)*pG->U[k][j][i].d;
01954 #ifndef BAROTROPIC
01955         Ur_x3Face[k][j][i].E -= q1*(x1Flux[k][j][i  ].d*(phic - phil)
01956                                   + x1Flux[k][j][i+1].d*(phir - phic));
01957 #endif
01958 
01959         phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j+1][i]);
01960         phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j-1][i]);
01961 
01962         Ur_x3Face[k][j][i].Mz -= q2*(phir-phil)*pG->U[k][j][i].d;
01963 #ifndef BAROTROPIC
01964         Ur_x3Face[k][j][i].E -= q2*(x2Flux[k][j  ][i].d*(phic - phil)
01965                                   + x2Flux[k][j+1][i].d*(phir - phic));
01966 #endif
01967 
01968 /* correct left states; x1 and x2 gradients */
01969         phic = pG->Phi[k-1][j][i];
01970         phir = 0.5*(pG->Phi[k-1][j][i] + pG->Phi[k-1][j][i+1]);
01971         phil = 0.5*(pG->Phi[k-1][j][i] + pG->Phi[k-1][j][i-1]);
01972 
01973         Ul_x3Face[k][j][i].My -= q1*(phir-phil)*pG->U[k-1][j][i].d;
01974 #ifndef BAROTROPIC
01975         Ul_x3Face[k][j][i].E -= q1*(x1Flux[k-1][j][i  ].d*(phic - phil)
01976                                   + x1Flux[k-1][j][i+1].d*(phir - phic));
01977 #endif
01978 
01979         phir = 0.5*(pG->Phi[k-1][j][i] + pG->Phi[k-1][j+1][i]);
01980         phil = 0.5*(pG->Phi[k-1][j][i] + pG->Phi[k-1][j-1][i]);
01981 
01982         Ul_x3Face[k][j][i].Mz -= q2*(phir-phil)*pG->U[k-1][j][i].d;
01983 #ifndef BAROTROPIC
01984         Ul_x3Face[k][j][i].E -= q2*(x2Flux[k-1][j  ][i].d*(phic - phil)
01985                                   + x2Flux[k-1][j+1][i].d*(phir - phic));
01986 #endif
01987       }
01988     }
01989   }
01990 #endif /* SELF_GRAVITY */
01991 
01992 /*--- Step 7d (cont) -----------------------------------------------------------
01993  * Add source terms for shearing box (Coriolis forces) for 0.5*dt arising from
01994  * x1-Flux gradient.  The tidal gravity terms are added via ShearingBoxPot
01995  *    Vx source term is (dt/2)( 2 Omega_0 Vy)
01996  *    Vy source term is (dt/2)(-2 Omega_0 Vx)
01997  *    Vy source term is (dt/2)((q-2) Omega_0 Vx) (with FARGO)
01998  */
01999 
02000 #ifdef SHEARING_BOX
02001   if (ShearingBoxPot != NULL){
02002   for (k=kl+1; k<=ku; k++) {
02003     for (j=jl+1; j<=ju-1; j++) {
02004       for (i=il+1; i<=iu-1; i++) {
02005         cc_pos(pG,i,j,k,&x1,&x2,&x3);
02006         phic = (*ShearingBoxPot)((x1            ),x2,x3);
02007         phir = (*ShearingBoxPot)((x1+0.5*pG->dx1),x2,x3);
02008         phil = (*ShearingBoxPot)((x1-0.5*pG->dx1),x2,x3);
02009 
02010 /* correct right states; x1 and x2 gradients */
02011         Ur_x3Face[k][j][i].My -= q1*(phir-phil)*pG->U[k][j][i].d;
02012 #ifndef BAROTROPIC
02013         Ur_x3Face[k][j][i].E -= q1*(x1Flux[k  ][j][i  ].d*(phic - phil)
02014                                   + x1Flux[k  ][j][i+1].d*(phir - phic));
02015 #endif
02016 
02017         phir = (*ShearingBoxPot)(x1,(x2+0.5*pG->dx2),x3);
02018         phil = (*ShearingBoxPot)(x1,(x2-0.5*pG->dx2),x3);
02019 
02020         Ur_x3Face[k][j][i].Mz -= q2*(phir-phil)*pG->U[k][j][i].d;
02021 #ifndef BAROTROPIC
02022         Ur_x3Face[k][j][i].E -= q2*(x2Flux[k  ][j  ][i].d*(phic - phil)
02023                                   + x2Flux[k  ][j+1][i].d*(phir - phic));
02024 #endif
02025 
02026 /* correct left states; x1 and x2 gradients */
02027         phic = (*ShearingBoxPot)((x1            ),x2,(x3-pG->dx3));
02028         phir = (*ShearingBoxPot)((x1+0.5*pG->dx1),x2,(x3-pG->dx3));
02029         phil = (*ShearingBoxPot)((x1-0.5*pG->dx1),x2,(x3-pG->dx3));
02030 
02031         Ul_x3Face[k][j][i].My -= q1*(phir-phil)*pG->U[k-1][j][i].d;
02032 #ifndef BAROTROPIC
02033         Ul_x3Face[k][j][i].E -= q1*(x1Flux[k-1][j][i  ].d*(phic - phil)
02034                                   + x1Flux[k-1][j][i+1].d*(phir - phic));
02035 #endif
02036 
02037         phir = (*ShearingBoxPot)(x1,(x2+0.5*pG->dx2),(x3-pG->dx3));
02038         phil = (*ShearingBoxPot)(x1,(x2-0.5*pG->dx2),(x3-pG->dx3));
02039 
02040         Ul_x3Face[k][j][i].Mz -= q2*(phir-phil)*pG->U[k-1][j][i].d;
02041 #ifndef BAROTROPIC
02042         Ul_x3Face[k][j][i].E -= q2*(x2Flux[k-1][j  ][i].d*(phic - phil)
02043                                   + x2Flux[k-1][j+1][i].d*(phir - phic));
02044 #endif
02045       }
02046     }
02047   }}
02048 
02049   for (k=kl+1; k<=ku; k++) {
02050     for (j=jl+1; j<=ju-1; j++) {
02051       for (i=il+1; i<=iu-1; i++) {
02052         Ur_x3Face[k][j][i].My += pG->dt*Omega_0*pG->U[k][j][i].M2;
02053         Ul_x3Face[k][j][i].My += pG->dt*Omega_0*pG->U[k-1][j][i].M2;
02054 #ifdef FARGO
02055         Ur_x3Face[k][j][i].Mz += hdt*(qshear-2.)*Omega_0*pG->U[k][j][i].M1;
02056         Ul_x3Face[k][j][i].Mz += hdt*(qshear-2.)*Omega_0*pG->U[k-1][j][i].M1;
02057 #else
02058         Ur_x3Face[k][j][i].Mz -= pG->dt*Omega_0*pG->U[k][j][i].M1;
02059         Ul_x3Face[k][j][i].Mz -= pG->dt*Omega_0*pG->U[k-1][j][i].M1;
02060 #endif
02061       }
02062     }
02063   }
02064 #endif /* SHEARING_BOX */
02065 
02066 #if defined(CYLINDRICAL) && defined(FARGO)
02067   for (k=kl+1; k<=ku; k++) {
02068     for (j=jl+1; j<=ju-1; j++) {
02069       for (i=il+1; i<=iu-1; i++) {
02070                                 Om = (*OrbitalProfile)(x1vc(pG,i));
02071                                 qshear = (*ShearProfile)(x1vc(pG,i));
02072 
02073                                 Ur_x3Face[k][j][i].My += pG->dt*Om*pG->U[k][j][i].M2;
02074                                 Ur_x3Face[k][j][i].Mz += hdt*(qshear-2.0)*Om*pG->U[k][j][i].M1;
02075 
02076                                 Ul_x3Face[k][j][i].My += pG->dt*Om*pG->U[k-1][j][i].M2;
02077                                 Ul_x3Face[k][j][i].Mz += hdt*(qshear-2.0)*Om*pG->U[k-1][j][i].M1;
02078                         }
02079                 }
02080         }
02081 #endif /* Cylindrical + Fargo */
02082 /*--- Step 7d (cont) -----------------------------------------------------------
02083  * ADD THE GEOMETRIC SOURCE-TERM IN THE X1-DIRECTION TO THE CORRECTED L/R 
02084  * STATES ON X3-FACES.  S_{M_R} = -(\rho V_\phi^2 - B_\phi^2)/R
02085  */
02086 #ifdef CYLINDRICAL
02087   for (k=kl+1; k<=ku; k++) {
02088     for (j=jl+1; j<=ju-1; j++) {
02089       for (i=il+1; i<=iu-1; i++) {
02090         Ul_x3Face[k][j][i].My += hdt*geom_src[k-1][j][i];
02091         Ur_x3Face[k][j][i].My += hdt*geom_src[k  ][j][i];
02092       }
02093     }
02094   }
02095 #endif /* CYLINDRICAL */
02096 
02097 /*=== STEP 8: Compute cell-centered values at n+1/2 ==========================*/
02098 
02099 /*--- Step 8a ------------------------------------------------------------------
02100  * Calculate d^{n+1/2} (needed with static potential, cooling, or MHD)
02101  */
02102 #ifndef MHD
02103 #ifndef PARTICLES
02104   if ((StaticGravPot != NULL) || (CoolingFunc != NULL))
02105 #endif
02106 #endif
02107   {
02108     for (k=kl+1; k<=ku-1; k++) {
02109       for (j=jl+1; j<=ju-1; j++) {
02110         for (i=il+1; i<=iu-1; i++) {
02111 #ifdef CYLINDRICAL
02112           rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
02113           q2 = hdt/(r[i]*pG->dx2);
02114 #endif
02115           dhalf[k][j][i] = pG->U[k][j][i].d 
02116             - q1*(rsf*x1Flux[k  ][j  ][i+1].d - lsf*x1Flux[k][j][i].d)
02117             - q2*(    x2Flux[k  ][j+1][i  ].d -     x2Flux[k][j][i].d)
02118             - q3*(    x3Flux[k+1][j  ][i  ].d -     x3Flux[k][j][i].d);
02119 #ifdef PARTICLES
02120           pG->Coup[k][j][i].grid_d = dhalf[k][j][i];
02121 #endif
02122         }
02123       }
02124     }
02125   }
02126 
02127 /*--- Step 8b ------------------------------------------------------------------
02128  * Calculate P^{n+1/2} (needed with cooling), and cell centered emf_cc^{n+1/2}
02129  */
02130 
02131 #ifndef MHD
02132 #ifndef PARTICLES
02133   if (CoolingFunc != NULL)
02134 #endif /* PARTICLES */
02135 #endif /* MHD */
02136   {
02137   for (k=kl+1; k<=ku-1; k++) {
02138     for (j=jl+1; j<=ju-1; j++) {
02139       for (i=il+1; i<=iu-1; i++) {
02140 #ifdef CYLINDRICAL
02141         rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
02142         q2 = hdt/(r[i]/pG->dx2);
02143 #endif
02144         M1h = pG->U[k][j][i].M1
02145            - q1*(rsf*x1Flux[k  ][j  ][i+1].Mx - lsf*x1Flux[k][j][i].Mx)
02146            - q2*(    x2Flux[k  ][j+1][i  ].Mz -     x2Flux[k][j][i].Mz)
02147            - q3*(    x3Flux[k+1][j  ][i  ].My -     x3Flux[k][j][i].My);
02148 
02149         M2h = pG->U[k][j][i].M2
02150            - q1*(SQR(rsf)*x1Flux[k  ][j  ][i+1].My - SQR(lsf)*x1Flux[k][j][i].My)
02151            - q2*(         x2Flux[k  ][j+1][i  ].Mx -          x2Flux[k][j][i].Mx)
02152            - q3*(         x3Flux[k+1][j  ][i  ].Mz -          x3Flux[k][j][i].Mz);
02153 
02154         M3h = pG->U[k][j][i].M3
02155            - q1*(lsf*x1Flux[k  ][j  ][i+1].Mz - rsf*x1Flux[k][j][i].Mz)
02156            - q2*(    x2Flux[k  ][j+1][i  ].My -     x2Flux[k][j][i].My)
02157            - q3*(    x3Flux[k+1][j  ][i  ].Mx -     x3Flux[k][j][i].Mx);
02158 
02159 #ifndef BAROTROPIC
02160         Eh = pG->U[k][j][i].E
02161            - q1*(rsf*x1Flux[k  ][j  ][i+1].E - lsf*x1Flux[k][j][i].E)
02162            - q2*(    x2Flux[k  ][j+1][i  ].E -     x2Flux[k][j][i].E)
02163            - q3*(    x3Flux[k+1][j  ][i  ].E -     x3Flux[k][j][i].E);
02164 #endif
02165 
02166 /* Add source terms for fixed gravitational potential */
02167         if (StaticGravPot != NULL){
02168           cc_pos(pG,i,j,k,&x1,&x2,&x3);
02169 #ifdef CYLINDRICAL
02170           g = (*x1GravAcc)(x1vc(pG,i),x2,x3);
02171 #ifdef FARGO
02172           g = g - x1vc(pG,i)*SQR((*OrbitalProfile)(x1vc(pG,i)));
02173 #endif
02174           M1h -= hdt*pG->U[k][j][i].d*g;
02175 #else
02176           phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
02177           phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
02178           M1h -= q1*(phir-phil)*pG->U[k][j][i].d;
02179 #endif
02180 
02181           phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
02182           phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
02183           M2h -= q2*(phir-phil)*pG->U[k][j][i].d;
02184 
02185           phir = (*StaticGravPot)(x1,x2,(x3+0.5*pG->dx3));
02186           phil = (*StaticGravPot)(x1,x2,(x3-0.5*pG->dx3));
02187           M3h -= q3*(phir-phil)*pG->U[k][j][i].d;
02188         }
02189 
02190 /* Add source terms due to self-gravity  */
02191 #ifdef SELF_GRAVITY
02192         phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j][i+1]);
02193         phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j][i-1]);
02194         M1h -= q1*(phir-phil)*pG->U[k][j][i].d;
02195 
02196         phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j+1][i]);
02197         phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j-1][i]);
02198         M2h -= q2*(phir-phil)*pG->U[k][j][i].d;
02199 
02200         phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k+1][j][i]);
02201         phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k-1][j][i]);
02202         M3h -= q3*(phir-phil)*pG->U[k][j][i].d;
02203 #endif /* SELF_GRAVITY */
02204 
02205 /* Add the tidal gravity and Coriolis terms for shearing box. */
02206 #ifdef SHEARING_BOX
02207         if (ShearingBoxPot != NULL){
02208           cc_pos(pG,i,j,k,&x1,&x2,&x3);
02209           phir = (*ShearingBoxPot)((x1+0.5*pG->dx1),x2,x3);
02210           phil = (*ShearingBoxPot)((x1-0.5*pG->dx1),x2,x3);
02211           M1h -= q1*(phir-phil)*pG->U[k][j][i].d;
02212 
02213           phir = (*ShearingBoxPot)(x1,(x2+0.5*pG->dx2),x3);
02214           phil = (*ShearingBoxPot)(x1,(x2-0.5*pG->dx2),x3);
02215           M2h -= q2*(phir-phil)*pG->U[k][j][i].d;
02216 
02217           phir = (*ShearingBoxPot)(x1,x2,(x3+0.5*pG->dx3));
02218           phil = (*ShearingBoxPot)(x1,x2,(x3-0.5*pG->dx3));
02219           M3h -= q3*(phir-phil)*pG->U[k][j][i].d;
02220         }
02221 
02222         M1h += pG->dt*Omega_0*pG->U[k][j][i].M2;
02223 #ifdef FARGO
02224         M2h += hdt*(qshear-2.)*Omega_0*pG->U[k][j][i].M1;
02225 #else
02226         M2h -= pG->dt*Omega_0*pG->U[k][j][i].M1;
02227 #endif
02228 #endif /* SHEARING_BOX */
02229 #if defined(CYLINDRICAL) && defined(FARGO)
02230         Om = (*OrbitalProfile)(x1vc(pG,i));
02231                                 qshear = (*ShearProfile)(x1vc(pG,i));
02232                           M1h += hdt*2.0*Om*pG->U[k][j][i].M2;
02233         M2h += hdt*Om*(qshear-2.0)*pG->U[k][j][i].M1;
02234 #endif /* Cylindrical + Fargo */
02235 
02236 /* Add the particle feedback terms */
02237 #ifdef FEEDBACK
02238       M1h -= pG->Coup[k][j][i].fb1;
02239       M2h -= pG->Coup[k][j][i].fb2;
02240       M3h -= pG->Coup[k][j][i].fb3;
02241 #endif /* FEEDBACK */
02242 
02243 /* Add the geometric source term */
02244 #ifdef CYLINDRICAL
02245       M1h += hdt*geom_src[k][j][i];
02246 #endif
02247 
02248 #ifndef BAROTROPIC
02249         phalf[k][j][i] = Eh - 0.5*(M1h*M1h + M2h*M2h + M3h*M3h)/dhalf[k][j][i];
02250 #endif
02251 
02252 #ifdef MHD
02253         B1ch = 0.5*(lsf*B1_x1Face[k][j][i] + rsf*B1_x1Face[k  ][j  ][i+1]);
02254         B2ch = 0.5*(    B2_x2Face[k][j][i] +     B2_x2Face[k  ][j+1][i  ]);
02255         B3ch = 0.5*(    B3_x3Face[k][j][i] +     B3_x3Face[k+1][j  ][i  ]);
02256         emf1_cc[k][j][i] = (B2ch*M3h - B3ch*M2h)/dhalf[k][j][i];
02257         emf2_cc[k][j][i] = (B3ch*M1h - B1ch*M3h)/dhalf[k][j][i];
02258         emf3_cc[k][j][i] = (B1ch*M2h - B2ch*M1h)/dhalf[k][j][i];
02259 #ifndef BAROTROPIC
02260         phalf[k][j][i] -= 0.5*(B1ch*B1ch + B2ch*B2ch + B3ch*B3ch);
02261 #endif
02262 #endif /* MHD */
02263 
02264 #ifndef BAROTROPIC
02265         phalf[k][j][i] *= Gamma_1;
02266 #endif
02267 
02268 #ifdef PARTICLES
02269       d1 = 1.0/dhalf[k][j][i];
02270       pG->Coup[k][j][i].grid_v1 = M1h*d1;
02271       pG->Coup[k][j][i].grid_v2 = M2h*d1;
02272       pG->Coup[k][j][i].grid_v3 = M3h*d1;
02273 #ifndef BAROTROPIC
02274       pG->Coup[k][j][i].grid_cs = sqrt(Gamma*phalf[k][j][i]*d1);
02275 #endif  /* BAROTROPIC */
02276 #endif /* PARTICLES */
02277 
02278       }
02279     }
02280   }
02281   }
02282 
02283 /*=== STEP 8.5: Integrate the particles, compute the feedback ================*/
02284 
02285 #ifdef PARTICLES
02286   Integrate_Particles(pG,pD);
02287 #ifdef FEEDBACK
02288   exchange_feedback(pG, pD);
02289 #endif
02290 #endif
02291 
02292 /*=== STEP 9: Compute 3D x1-Flux, x2-Flux, x3-Flux ===========================*/
02293 
02294 /*--- Step 9a ------------------------------------------------------------------
02295  * Compute maximum wavespeeds in multidimensions (eta in eq. 10 from Sanders et
02296  *  al. (1998)) for H-correction
02297  */
02298 
02299 #ifdef H_CORRECTION
02300   for (k=ks-1; k<=ke+1; k++) {
02301     for (j=js-1; j<=je+1; j++) {
02302       for (i=is-1; i<=ie+2; i++) {
02303 #ifdef MHD
02304         Bx = B1_x1Face[k][j][i];
02305 #endif
02306         cfr = cfast(&(Ur_x1Face[k][j][i]),&Bx);
02307         cfl = cfast(&(Ul_x1Face[k][j][i]),&Bx);
02308         lambdar = Ur_x1Face[k][j][i].Mx/Ur_x1Face[k][j][i].d + cfr;
02309         lambdal = Ul_x1Face[k][j][i].Mx/Ul_x1Face[k][j][i].d - cfl;
02310         eta1[k][j][i] = 0.5*fabs(lambdar - lambdal);
02311       }
02312     }
02313   }
02314 
02315   for (k=ks-1; k<=ke+1; k++) {
02316     for (j=js-1; j<=je+2; j++) {
02317       for (i=is-1; i<=ie+1; i++) {
02318 #ifdef MHD
02319         Bx = B2_x2Face[k][j][i];
02320 #endif
02321         cfr = cfast(&(Ur_x2Face[k][j][i]),&Bx);
02322         cfl = cfast(&(Ul_x2Face[k][j][i]),&Bx);
02323         lambdar = Ur_x2Face[k][j][i].Mx/Ur_x2Face[k][j][i].d + cfr;
02324         lambdal = Ul_x2Face[k][j][i].Mx/Ul_x2Face[k][j][i].d - cfl;
02325         eta2[k][j][i] = 0.5*fabs(lambdar - lambdal);
02326       }
02327     }
02328   }
02329 
02330   for (k=ks-1; k<=ke+2; k++) {
02331     for (j=js-1; j<=je+1; j++) {
02332       for (i=is-1; i<=ie+1; i++) {
02333 #ifdef MHD
02334         Bx = B3_x3Face[k][j][i];
02335 #endif
02336         cfr = cfast(&(Ur_x3Face[k][j][i]),&Bx);
02337         cfl = cfast(&(Ul_x3Face[k][j][i]),&Bx);
02338         lambdar = Ur_x3Face[k][j][i].Mx/Ur_x3Face[k][j][i].d + cfr;
02339         lambdal = Ul_x3Face[k][j][i].Mx/Ul_x3Face[k][j][i].d - cfl;
02340         eta3[k][j][i] = 0.5*fabs(lambdar - lambdal);
02341       }
02342     }
02343   }
02344 #endif /* H_CORRECTION */
02345 
02346 /*--- Step 9b ------------------------------------------------------------------
02347  * Compute 3D x1-fluxes from corrected L/R states.
02348  */
02349 
02350   for (k=ks-1; k<=ke+1; k++) {
02351     for (j=js-1; j<=je+1; j++) {
02352       for (i=is; i<=ie+1; i++) {
02353 #ifdef H_CORRECTION
02354         etah = MAX(eta2[k][j][i-1],eta2[k][j][i]);
02355         etah = MAX(etah,eta2[k][j+1][i-1]);
02356         etah = MAX(etah,eta2[k][j+1][i  ]);
02357 
02358         etah = MAX(etah,eta3[k  ][j][i-1]);
02359         etah = MAX(etah,eta3[k  ][j][i  ]);
02360         etah = MAX(etah,eta3[k+1][j][i-1]);
02361         etah = MAX(etah,eta3[k+1][j][i  ]);
02362 
02363         etah = MAX(etah,eta1[k  ][j][i  ]);
02364 #endif /* H_CORRECTION */
02365 #ifdef MHD
02366         Bx = B1_x1Face[k][j][i];
02367 #endif
02368         Wl[i] = Cons1D_to_Prim1D(&Ul_x1Face[k][j][i],&Bx);
02369         Wr[i] = Cons1D_to_Prim1D(&Ur_x1Face[k][j][i],&Bx);
02370 
02371         fluxes(Ul_x1Face[k][j][i],Ur_x1Face[k][j][i],Wl[i],Wr[i],Bx,
02372                &x1Flux[k][j][i]);
02373       }
02374     }
02375   }
02376 
02377 /*--- Step 9c ------------------------------------------------------------------
02378  * Compute 3D x2-fluxes from corrected L/R states.
02379  */
02380 
02381   for (k=ks-1; k<=ke+1; k++) {
02382     for (j=js; j<=je+1; j++) {
02383       for (i=is-1; i<=ie+1; i++) {
02384 #ifdef H_CORRECTION
02385         etah = MAX(eta1[k][j-1][i],eta1[k][j][i]);
02386         etah = MAX(etah,eta1[k][j-1][i+1]);
02387         etah = MAX(etah,eta1[k][j  ][i+1]);
02388 
02389         etah = MAX(etah,eta3[k  ][j-1][i]);
02390         etah = MAX(etah,eta3[k  ][j  ][i]);
02391         etah = MAX(etah,eta3[k+1][j-1][i]);
02392         etah = MAX(etah,eta3[k+1][j  ][i]);
02393 
02394         etah = MAX(etah,eta2[k  ][j  ][i]);
02395 #endif /* H_CORRECTION */
02396 #ifdef MHD
02397         Bx = B2_x2Face[k][j][i];
02398 #endif
02399         Wl[i] = Cons1D_to_Prim1D(&Ul_x2Face[k][j][i],&Bx);
02400         Wr[i] = Cons1D_to_Prim1D(&Ur_x2Face[k][j][i],&Bx);
02401 
02402         fluxes(Ul_x2Face[k][j][i],Ur_x2Face[k][j][i],Wl[i],Wr[i],Bx,
02403                &x2Flux[k][j][i]);
02404       }
02405     }
02406   }
02407 
02408 /*--- Step 9d ------------------------------------------------------------------
02409  * Compute 3D x3-fluxes from corrected L/R states.
02410  */
02411 
02412   for (k=ks; k<=ke+1; k++) {
02413     for (j=js-1; j<=je+1; j++) {
02414       for (i=is-1; i<=ie+1; i++) {
02415 #ifdef H_CORRECTION
02416         etah = MAX(eta1[k-1][j][i],eta1[k][j][i]);
02417         etah = MAX(etah,eta1[k-1][j][i+1]);
02418         etah = MAX(etah,eta1[k][j  ][i+1]);
02419 
02420         etah = MAX(etah,eta2[k-1][j  ][i]);
02421         etah = MAX(etah,eta2[k  ][j  ][i]);
02422         etah = MAX(etah,eta2[k-1][j+1][i]);
02423         etah = MAX(etah,eta2[k  ][j+1][i]);
02424 
02425         etah = MAX(etah,eta3[k  ][j  ][i]);
02426 #endif /* H_CORRECTION */
02427 #ifdef MHD
02428         Bx = B3_x3Face[k][j][i];
02429 #endif
02430         Wl[i] = Cons1D_to_Prim1D(&Ul_x3Face[k][j][i],&Bx);
02431         Wr[i] = Cons1D_to_Prim1D(&Ur_x3Face[k][j][i],&Bx);
02432 
02433         fluxes(Ul_x3Face[k][j][i],Ur_x3Face[k][j][i],Wl[i],Wr[i],Bx,
02434                &x3Flux[k][j][i]);
02435       }
02436     }
02437   }
02438 
02439 /*=== STEP 10: Update face-centered B for a full timestep ====================*/
02440 
02441 /*--- Step 10a -----------------------------------------------------------------
02442  * Integrate emf*^{n+1/2} to the grid cell corners
02443  */
02444 
02445 #ifdef MHD
02446   integrate_emf1_corner(pG);
02447   integrate_emf2_corner(pG);
02448   integrate_emf3_corner(pG);
02449 
02450 /* Remap Ey at is and ie+1 to conserve Bz in shearing box */
02451 #ifdef SHEARING_BOX
02452     get_myGridIndex(pD, myID_Comm_world, &my_iproc, &my_jproc, &my_kproc);
02453 
02454 /* compute remapped Ey from opposite side of grid */
02455 
02456     if (my_iproc == 0) {
02457       RemapEy_ix1(pD, emf2, remapEyiib);
02458     }
02459     if (my_iproc == (pD->NGrid[0]-1)) {
02460       RemapEy_ox1(pD, emf2, remapEyoib);
02461     }
02462 
02463 /* Now average Ey and remapped Ey */
02464 
02465     if (my_iproc == 0) {
02466       for(k=ks; k<=ke+1; k++) {
02467         for(j=js; j<=je; j++){
02468           emf2[k][j][is]  = 0.5*(emf2[k][j][is] + remapEyiib[k][j]);
02469         }
02470       }
02471     }
02472 
02473     if (my_iproc == (pD->NGrid[0]-1)) {
02474       for(k=ks; k<=ke+1; k++) {
02475         for(j=js; j<=je; j++){
02476           emf2[k][j][ie+1]  = 0.5*(emf2[k][j][ie+1] + remapEyoib[k][j]);
02477         }
02478       }
02479     }
02480 #endif /* SHEARING_BOX */
02481 
02482 /*--- Step 10b -----------------------------------------------------------------
02483  * Update the interface magnetic fields using CT for a full time step.
02484  */
02485 
02486   for (k=ks; k<=ke; k++) {
02487     for (j=js; j<=je; j++) {
02488       for (i=is; i<=ie; i++) {
02489 #ifdef CYLINDRICAL
02490         dtodx2 = pG->dt/(ri[i]*pG->dx2);
02491 #endif
02492         pG->B1i[k][j][i] += dtodx3*(emf2[k+1][j  ][i  ] - emf2[k][j][i]) -
02493                             dtodx2*(emf3[k  ][j+1][i  ] - emf3[k][j][i]);
02494         pG->B2i[k][j][i] += dtodx1*(emf3[k  ][j  ][i+1] - emf3[k][j][i]) -
02495                             dtodx3*(emf1[k+1][j  ][i  ] - emf1[k][j][i]);
02496 #ifdef CYLINDRICAL
02497         rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
02498         dtodx2 = pG->dt/(r[i]*pG->dx2);
02499 #endif
02500         pG->B3i[k][j][i] += dtodx2*(    emf1[k  ][j+1][i  ] -     emf1[k][j][i]) -
02501                             dtodx1*(rsf*emf2[k  ][j  ][i+1] - lsf*emf2[k][j][i]);
02502       }
02503 #ifdef CYLINDRICAL
02504       dtodx2 = pG->dt/(ri[ie+1]*pG->dx2);
02505 #endif
02506       pG->B1i[k][j][ie+1] +=
02507         dtodx3*(emf2[k+1][j  ][ie+1] - emf2[k][j][ie+1]) -
02508         dtodx2*(emf3[k  ][j+1][ie+1] - emf3[k][j][ie+1]);
02509     }
02510     for (i=is; i<=ie; i++) {
02511       pG->B2i[k][je+1][i] +=
02512         dtodx1*(emf3[k  ][je+1][i+1] - emf3[k][je+1][i]) -
02513         dtodx3*(emf1[k+1][je+1][i  ] - emf1[k][je+1][i]);
02514     }
02515   }
02516   for (j=js; j<=je; j++) {
02517     for (i=is; i<=ie; i++) {
02518 #ifdef CYLINDRICAL
02519       rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
02520       dtodx2 = pG->dt/(r[i]*pG->dx2);
02521 #endif
02522       pG->B3i[ke+1][j][i] += 
02523         dtodx2*(    emf1[ke+1][j+1][i  ] - emf1[ke+1][j][i]) -
02524         dtodx1*(rsf*emf2[ke+1][j  ][i+1] - lsf*emf2[ke+1][j][i]);
02525     }
02526   }
02527 #endif /* MHD */
02528 
02529 /*=== STEP 11: Add source terms for a full timestep using n+1/2 states =======*/
02530 
02531 /*--- Step 11a -----------------------------------------------------------------
02532  * ADD GEOMETRIC SOURCE TERMS.
02533  */
02534 #if defined(CYLINDRICAL) && !defined(FARGO)
02535   for (k=ks; k<=ke; k++) {
02536     for (j=js; j<=je; j++) {
02537       for (i=is; i<=ie; i++) {
02538         rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
02539         q2 = hdt/(r[i]*pG->dx2);
02540 
02541         /* CALCULATE d AT TIME n+1/2 */
02542         dhalf[k][j][i] = pG->U[k][j][i].d 
02543           - q1*(rsf*x1Flux[k  ][j  ][i+1].d - lsf*x1Flux[k][j][i].d)
02544           - q2*(    x2Flux[k  ][j+1][i  ].d -     x2Flux[k][j][i].d)
02545           - q3*(    x3Flux[k+1][j  ][i  ].d -     x3Flux[k][j][i].d);
02546 
02547         /* CALCULATE M2 AT TIME n+1/2 */
02548         M2h = pG->U[k][j][i].M2
02549           - q1*(SQR(rsf)*x1Flux[k  ][j  ][i+1].My - SQR(lsf)*x1Flux[k][j][i].My)
02550           - q2*(         x2Flux[k  ][j+1][i  ].Mx -          x2Flux[k][j][i].Mx)
02551           - q3*(         x3Flux[k+1][j  ][i  ].Mz -          x3Flux[k][j][i].Mz);
02552 
02553         /* ADD SOURCE TERM FOR FIXED GRAVITATIONAL POTENTIAL FOR 0.5*dt */
02554         if (StaticGravPot != NULL){
02555           phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
02556           phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
02557           M2h -= q2*(phir-phil)*pG->U[k][j][i].d;
02558         }
02559 
02560         /* COMPUTE GEOMETRIC SOURCE TERM AT TIME n+1/2 */
02561         geom_src[k][j][i] = SQR(M2h)/dhalf[k][j][i];
02562 #ifdef MHD
02563         B2ch = 0.5*(B2_x2Face[k][j][i] + B2_x2Face[k][j+1][i]);
02564         geom_src[k][j][i] -= SQR(B2ch);
02565 #endif
02566 #ifdef ISOTHERMAL
02567         geom_src[k][j][i] += Iso_csound2*dhalf[k][j][i];
02568 #ifdef MHD
02569         B1ch = 0.5*(lsf*B1_x1Face[k][j][i] + rsf*B1_x1Face[k  ][j  ][i+1]);
02570         B3ch = 0.5*(    B3_x3Face[k][j][i] +     B3_x3Face[k+1][j  ][i  ]);
02571         geom_src[k][j][i] += 0.5*(SQR(B1ch)+SQR(B2ch)+SQR(B3ch));
02572 #endif
02573 #else /* ISOTHERMAL */
02574         Pavgh = 0.5*(lsf*x1Flux[k][j][i  ].Pflux + rsf*x1Flux[k][j][i+1].Pflux);
02575         geom_src[k][j][i] += Pavgh;
02576 #endif
02577         geom_src[k][j][i] /= x1vc(pG,i);
02578 
02579         /* ADD TIME-CENTERED GEOMETRIC SOURCE TERM FOR FULL dt */
02580         pG->U[k][j][i].M1 += pG->dt*geom_src[k][j][i];
02581       }
02582     }
02583   }
02584 #endif /* CYLINDRICAL + !Fargo */
02585 
02586 // Add source terms using Heun's method for cylindrical fargo
02587 #if defined(CYLINDRICAL) && defined(FARGO)
02588   for (k=ks; k<=ke; k++) {
02589     for (j=js; j<=je; j++) {
02590       for (i=is; i<=ie; i++) {
02591                                 cc_pos(pG,i,j,k,&x1,&x2,&x3);
02592         rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
02593         q2 = hdt/(r[i]*pG->dx2);
02594 
02595         /* CALCULATE d AT TIME n+1/2 */
02596         dhalf[k][j][i] = pG->U[k][j][i].d 
02597           - q1*(rsf*x1Flux[k  ][j  ][i+1].d - lsf*x1Flux[k][j][i].d)
02598           - q2*(    x2Flux[k  ][j+1][i  ].d -     x2Flux[k][j][i].d)
02599           - q3*(    x3Flux[k+1][j  ][i  ].d -     x3Flux[k][j][i].d);
02600                                 /* Save current R/phi momenta */
02601                                 Mrn = pG->U[k][j][i].M1;
02602                                 Mpn = pG->U[k][j][i].M2;
02603 
02604                                 Om = (*OrbitalProfile)(x1vc(pG,i));
02605                                 qshear = (*ShearProfile)(x1vc(pG,i));
02606                                 g = (*x1GravAcc)(x1vc(pG,i),x2,x3) - x1vc(pG,i)*SQR((*OrbitalProfile)(x1vc(pG,i)));
02607                                 /* Use forward euler to approximate R/phi momenta at t^{n+1} */
02608                                 Mre = Mrn
02609                                         - (dtodx1/r[i])*(ri[i+1]*x1Flux[k][j][i+1].Mx - ri[i]*x1Flux[k][j][i].Mx)
02610                                         - (dtodx2/r[i])*(        x2Flux[k][j+1][i].Mz -       x2Flux[k][j][i].Mz)
02611                                         - (dtodx3)*(             x3Flux[k+1][j][i].My -       x3Flux[k][j][i].My);
02612                                 Mre += pG->dt*( 2.0*Om*Mpn + geom_src[k][j][i] - pG->U[k][j][i].d*g);
02613         
02614                                 Mpe = Mpn + pG->dt*Om*(qshear-2.0)*Mrn
02615                                         - (dtodx1/SQR(r[i]))*(SQR(ri[i+1])*x1Flux[k ][j ][i+1].My - SQR(ri[i])*x1Flux[k][j][i].My)
02616           - (dtodx2/r[i])*( x2Flux[k ][j+1][i ].Mx - x2Flux[k][j][i].Mx)
02617           - (dtodx3)*(      x3Flux[k+1][j ][i ].Mz - x3Flux[k][j][i].Mz);
02618         /* Average forward euler and current values to approximate at t^{n+1/2} */
02619                                 Mrav = 0.5*(Mrn+Mre);
02620                                 Mpav = 0.5*(Mpn+Mpe);
02621 
02622                                 /* Compute source terms at t^{n+1/2} */
02623                                 geom_src[k][j][i] = SQR(Mpav)/dhalf[k][j][i];
02624 #ifdef MHD
02625         B1ch = 0.5*(ri[i]*B1_x1Face[k][j][i] + ri[i+1]*B1_x1Face[k  ][j  ][i+1])/r[i];
02626         B2ch = 0.5*(B2_x2Face[k][j][i] + B2_x2Face[k][j+1][i]);
02627         B3ch = 0.5*(B3_x3Face[k][j][i] + B3_x3Face[k+1][j  ][i  ]);
02628         geom_src[k][j][i] += 0.5*( SQR(B1ch) - SQR(B2ch) + SQR(B3ch));
02629 #endif /* MHD */
02630 #ifdef ISOTHERMAL
02631                                 geom_src[k][j][i] += Iso_csound2*dhalf[k][j][i];
02632 #endif /* Isothermal */
02633                                 geom_src[k][j][i] /= x1vc(pG,i);
02634 
02635                                 /* Use average values to apply source terms for full time-step */
02636                                 pG->U[k][j][i].M1 += pG->dt*( 2.0*Om*Mpav + geom_src[k][j][i]);
02637         pG->U[k][j][i].M2 += pG->dt*( Om*(qshear-2.0)*Mrav);
02638 
02639                         }
02640                 }
02641         }
02642 
02643 #endif /* Cylindrical + Fargo */
02644 
02645 /*--- Step 11a -----------------------------------------------------------------
02646  * Add gravitational (or shearing box) source terms as a Static Potential.
02647  *   A Crank-Nicholson update is used for shearing box terms.
02648  *   The energy source terms computed at cell faces are averaged to improve
02649  * conservation of total energy.
02650  *    S_{M} = -(\rho)^{n+1/2} Grad(Phi);   S_{E} = -(\rho v)^{n+1/2} Grad{Phi}
02651  */
02652 
02653 #ifdef SHEARING_BOX
02654   fact = om_dt/(2. + (2.-qshear)*om_dt*om_dt);
02655   qom = qshear*Omega_0;
02656   for(k=ks; k<=ke; k++) {
02657     for(j=js; j<=je; j++) {
02658       for(i=is; i<=ie; i++) {
02659         cc_pos(pG,i,j,k,&x1,&x2,&x3);
02660 
02661 /* Store the current state */
02662         M1n  = pG->U[k][j][i].M1;
02663 #ifdef FARGO
02664         dM2n = pG->U[k][j][i].M2;
02665 #else
02666         dM2n = pG->U[k][j][i].M2 + qom*x1*pG->U[k][j][i].d;
02667 #endif
02668 
02669 /* Calculate the flux for the y-momentum fluctuation */
02670         frx1_dM2 = x1Flux[k][j][i+1].My;
02671         flx1_dM2 = x1Flux[k][j][i  ].My;
02672         frx2_dM2 = x2Flux[k][j+1][i].Mx;
02673         flx2_dM2 = x2Flux[k][j  ][i].Mx;
02674         frx3_dM2 = x3Flux[k+1][j][i].Mz;
02675         flx3_dM2 = x3Flux[k  ][j][i].Mz;
02676 #ifndef FARGO
02677         frx1_dM2 += qom*(x1+0.5*pG->dx1)*x1Flux[k][j][i+1].d;
02678         flx1_dM2 += qom*(x1-0.5*pG->dx1)*x1Flux[k][j][i  ].d;
02679         frx2_dM2 += qom*(x1            )*x2Flux[k][j+1][i].d;
02680         flx2_dM2 += qom*(x1            )*x2Flux[k][j  ][i].d;
02681         frx3_dM2 += qom*(x1            )*x3Flux[k+1][j][i].d;
02682         flx3_dM2 += qom*(x1            )*x3Flux[k  ][j][i].d;
02683 #endif
02684 
02685 /* Now evolve M1n and dM2n by dt/2 using Forward Euler */
02686         M1e = M1n - q1*(x1Flux[k][j][i+1].Mx - x1Flux[k][j][i].Mx)
02687                   - q2*(x2Flux[k][j+1][i].Mz - x2Flux[k][j][i].Mz)
02688                   - q3*(x3Flux[k+1][j][i].My - x3Flux[k][j][i].My);
02689 
02690         dM2e = dM2n - q1*(frx1_dM2 - flx1_dM2)
02691                     - q2*(frx2_dM2 - flx2_dM2) 
02692                     - q3*(frx3_dM2 - flx3_dM2);
02693 
02694 #ifdef FEEDBACK
02695       M1e -= 0.5*pG->Coup[k][j][i].fb1;
02696       dM2e -= 0.5*pG->Coup[k][j][i].fb2;
02697 #endif
02698 
02699 /* Update the 1- and 2-momentum for the Coriolis and tidal
02700  * potential momentum source terms using a Crank-Nicholson
02701  * discretization for the momentum fluctuation equation. */
02702 
02703         pG->U[k][j][i].M1 += (4.0*dM2e + 2.0*(qshear-2.)*om_dt*M1e)*fact;
02704         pG->U[k][j][i].M2 += 2.0*(qshear-2.)*(M1e + om_dt*dM2e)*fact;
02705 
02706 #ifndef FARGO
02707         pG->U[k][j][i].M2 -= 0.5*qshear*om_dt*
02708            (x1Flux[k][j][i].d + x1Flux[k][j][i+1].d);
02709 #endif
02710 
02711 /* Update the energy for a fixed potential.
02712  * This update is identical to non-SHEARING_BOX below  */
02713 
02714         phic = (*ShearingBoxPot)((x1            ),x2,x3);
02715         phir = (*ShearingBoxPot)((x1+0.5*pG->dx1),x2,x3);
02716         phil = (*ShearingBoxPot)((x1-0.5*pG->dx1),x2,x3);
02717 #ifndef BAROTROPIC
02718         pG->U[k][j][i].E -= dtodx1*(x1Flux[k][j][i  ].d*(phic - phil) +
02719                                     x1Flux[k][j][i+1].d*(phir - phic));
02720 #endif
02721 
02722         phir = (*ShearingBoxPot)(x1,(x2+0.5*pG->dx2),x3);
02723         phil = (*ShearingBoxPot)(x1,(x2-0.5*pG->dx2),x3);
02724 #ifndef BAROTROPIC
02725         pG->U[k][j][i].E -= dtodx2*(x2Flux[k][j  ][i].d*(phic - phil) +
02726                                     x2Flux[k][j+1][i].d*(phir - phic));
02727 #endif
02728 
02729         phir = (*ShearingBoxPot)(x1,x2,(x3+0.5*pG->dx3));
02730         phil = (*ShearingBoxPot)(x1,x2,(x3-0.5*pG->dx3));
02731 #ifndef BAROTROPIC
02732         pG->U[k][j][i].E -= dtodx3*(x3Flux[k  ][j][i].d*(phic - phil) +
02733                                     x3Flux[k+1][j][i].d*(phir - phic));
02734 #endif
02735       }
02736     }
02737   }
02738 
02739 #endif /* SHEARING_BOX */
02740 
02741   if (StaticGravPot != NULL){
02742     for (k=ks; k<=ke; k++) {
02743       for (j=js; j<=je; j++) {
02744         for (i=is; i<=ie; i++) {
02745           cc_pos(pG,i,j,k,&x1,&x2,&x3);
02746           phic = (*StaticGravPot)((x1            ),x2,x3);
02747           phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
02748           phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
02749 
02750 #ifdef CYLINDRICAL
02751           g = (*x1GravAcc)(x1vc(pG,i),x2,x3);
02752 #ifdef FARGO
02753                                         g = g - x1vc(pG,i)*SQR((*OrbitalProfile)(x1vc(pG,i)));
02754 #endif
02755           rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
02756           dtodx2 = pG->dt/(r[i]*pG->dx2);
02757           pG->U[k][j][i].M1 -= pG->dt*dhalf[k][j][i]*g;
02758 #else
02759           pG->U[k][j][i].M1 -= dtodx1*(phir-phil)*dhalf[k][j][i];
02760 #endif
02761 #ifndef BAROTROPIC
02762           pG->U[k][j][i].E -= dtodx1*(lsf*x1Flux[k][j][i  ].d*(phic - phil) +
02763                                       rsf*x1Flux[k][j][i+1].d*(phir - phic));
02764 #endif
02765           phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
02766           phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
02767           pG->U[k][j][i].M2 -= dtodx2*(phir-phil)*dhalf[k][j][i];
02768 #ifndef BAROTROPIC
02769           pG->U[k][j][i].E -= dtodx2*(x2Flux[k][j  ][i].d*(phic - phil) +
02770                                       x2Flux[k][j+1][i].d*(phir - phic));
02771 #endif
02772           phir = (*StaticGravPot)(x1,x2,(x3+0.5*pG->dx3));
02773           phil = (*StaticGravPot)(x1,x2,(x3-0.5*pG->dx3));
02774           pG->U[k][j][i].M3 -= dtodx3*(phir-phil)*dhalf[k][j][i];
02775 #ifndef BAROTROPIC
02776           pG->U[k][j][i].E -= dtodx3*(x3Flux[k  ][j][i].d*(phic - phil) +
02777                                       x3Flux[k+1][j][i].d*(phir - phic));
02778 #endif
02779         }
02780       }
02781     }
02782   }
02783 
02784 /*--- Step 11b -----------------------------------------------------------------
02785  * Add source terms for self-gravity.
02786  * A flux correction using Phi^{n+1} in the main loop is required to make
02787  * the source terms 2nd order: see selfg_flux_correction().
02788  */
02789 
02790 #ifdef SELF_GRAVITY
02791 /* Add fluxes and source terms due to (d/dx1) terms  */
02792 
02793   for (k=ks; k<=ke; k++){
02794     for (j=js; j<=je; j++){
02795       for (i=is; i<=ie; i++){
02796         phic = pG->Phi[k][j][i];
02797         phil = 0.5*(pG->Phi[k][j][i-1] + pG->Phi[k][j][i  ]);
02798         phir = 0.5*(pG->Phi[k][j][i  ] + pG->Phi[k][j][i+1]);
02799 
02800 /* gx, gy and gz centered at L and R x1-faces */
02801         gxl = (pG->Phi[k][j][i-1] - pG->Phi[k][j][i  ])*(dx1i);
02802         gxr = (pG->Phi[k][j][i  ] - pG->Phi[k][j][i+1])*(dx1i);
02803 
02804         gyl = 0.25*((pG->Phi[k][j-1][i-1] - pG->Phi[k][j+1][i-1]) +
02805                     (pG->Phi[k][j-1][i  ] - pG->Phi[k][j+1][i  ]) )*(dx2i);
02806         gyr = 0.25*((pG->Phi[k][j-1][i  ] - pG->Phi[k][j+1][i  ]) +
02807                     (pG->Phi[k][j-1][i+1] - pG->Phi[k][j+1][i+1]) )*(dx2i);
02808 
02809         gzl = 0.25*((pG->Phi[k-1][j][i-1] - pG->Phi[k+1][j][i-1]) +
02810                     (pG->Phi[k-1][j][i  ] - pG->Phi[k+1][j][i  ]) )*(dx3i);
02811         gzr = 0.25*((pG->Phi[k-1][j][i  ] - pG->Phi[k+1][j][i  ]) +
02812                     (pG->Phi[k-1][j][i+1] - pG->Phi[k+1][j][i+1]) )*(dx3i);
02813 
02814 /* momentum fluxes in x1.  2nd term is needed only if Jean's swindle used */
02815         flx_m1l = 0.5*(gxl*gxl-gyl*gyl-gzl*gzl)/four_pi_G + grav_mean_rho*phil;
02816         flx_m1r = 0.5*(gxr*gxr-gyr*gyr-gzr*gzr)/four_pi_G + grav_mean_rho*phir;
02817 
02818         flx_m2l = gxl*gyl/four_pi_G;
02819         flx_m2r = gxr*gyr/four_pi_G;
02820 
02821         flx_m3l = gxl*gzl/four_pi_G;
02822         flx_m3r = gxr*gzr/four_pi_G;
02823 
02824 /* Update momenta and energy with d/dx1 terms  */
02825         pG->U[k][j][i].M1 -= dtodx1*(flx_m1r - flx_m1l);
02826         pG->U[k][j][i].M2 -= dtodx1*(flx_m2r - flx_m2l);
02827         pG->U[k][j][i].M3 -= dtodx1*(flx_m3r - flx_m3l);
02828 #ifndef BAROTROPIC
02829         pG->U[k][j][i].E -= dtodx1*(x1Flux[k][j][i  ].d*(phic - phil) +
02830                                     x1Flux[k][j][i+1].d*(phir - phic));
02831 #endif /* BAROTROPIC */
02832       }
02833     }
02834   }
02835 
02836 /* Add fluxes and source terms due to (d/dx2) terms  */
02837 
02838   for (k=ks; k<=ke; k++){
02839     for (j=js; j<=je; j++){
02840       for (i=is; i<=ie; i++){
02841         phic = pG->Phi[k][j][i];
02842         phil = 0.5*(pG->Phi[k][j-1][i] + pG->Phi[k][j  ][i]);
02843         phir = 0.5*(pG->Phi[k][j  ][i] + pG->Phi[k][j+1][i]);
02844 
02845 /* gx, gy and gz centered at L and R x2-faces */
02846         gxl = 0.25*((pG->Phi[k][j-1][i-1] - pG->Phi[k][j-1][i+1]) +
02847                     (pG->Phi[k][j  ][i-1] - pG->Phi[k][j  ][i+1]) )*(dx1i);
02848         gxr = 0.25*((pG->Phi[k][j  ][i-1] - pG->Phi[k][j  ][i+1]) +
02849                     (pG->Phi[k][j+1][i-1] - pG->Phi[k][j+1][i+1]) )*(dx1i);
02850 
02851         gyl = (pG->Phi[k][j-1][i] - pG->Phi[k][j  ][i])*(dx2i);
02852         gyr = (pG->Phi[k][j  ][i] - pG->Phi[k][j+1][i])*(dx2i);
02853 
02854         gzl = 0.25*((pG->Phi[k-1][j-1][i] - pG->Phi[k+1][j-1][i]) +
02855                     (pG->Phi[k-1][j  ][i] - pG->Phi[k+1][j  ][i]) )*(dx3i);
02856         gzr = 0.25*((pG->Phi[k-1][j  ][i] - pG->Phi[k+1][j  ][i]) +
02857                     (pG->Phi[k-1][j+1][i] - pG->Phi[k+1][j+1][i]) )*(dx3i);
02858 
02859 /* momentum fluxes in x2.  2nd term is needed only if Jean's swindle used */
02860         flx_m1l = gyl*gxl/four_pi_G;
02861         flx_m1r = gyr*gxr/four_pi_G;
02862 
02863         flx_m2l = 0.5*(gyl*gyl-gxl*gxl-gzl*gzl)/four_pi_G + grav_mean_rho*phil;
02864         flx_m2r = 0.5*(gyr*gyr-gxr*gxr-gzr*gzr)/four_pi_G + grav_mean_rho*phir;
02865 
02866         flx_m3l = gyl*gzl/four_pi_G;
02867         flx_m3r = gyr*gzr/four_pi_G;
02868 
02869 /* Update momenta and energy with d/dx2 terms  */
02870         pG->U[k][j][i].M1 -= dtodx2*(flx_m1r - flx_m1l);
02871         pG->U[k][j][i].M2 -= dtodx2*(flx_m2r - flx_m2l);
02872         pG->U[k][j][i].M3 -= dtodx2*(flx_m3r - flx_m3l);
02873 #ifndef BAROTROPIC
02874         pG->U[k][j][i].E -= dtodx2*(x2Flux[k][j  ][i].d*(phic - phil) +
02875                                     x2Flux[k][j+1][i].d*(phir - phic));
02876 #endif /* BAROTROPIC */
02877       }
02878     }
02879   }
02880 
02881 /* Add fluxes and source terms due to (d/dx3) terms  */
02882 
02883   for (k=ks; k<=ke; k++){
02884     for (j=js; j<=je; j++){
02885       for (i=is; i<=ie; i++){
02886         phic = pG->Phi[k][j][i];
02887         phil = 0.5*(pG->Phi[k-1][j][i] + pG->Phi[k  ][j][i]);
02888         phir = 0.5*(pG->Phi[k  ][j][i] + pG->Phi[k+1][j][i]);
02889 
02890 /* gx, gy and gz centered at L and R x3-faces */
02891         gxl = 0.25*((pG->Phi[k-1][j][i-1] - pG->Phi[k-1][j][i+1]) +
02892                     (pG->Phi[k  ][j][i-1] - pG->Phi[k  ][j][i+1]) )*(dx1i);
02893         gxr = 0.25*((pG->Phi[k  ][j][i-1] - pG->Phi[k  ][j][i+1]) +
02894                     (pG->Phi[k+1][j][i-1] - pG->Phi[k+1][j][i+1]) )*(dx1i);
02895 
02896         gyl = 0.25*((pG->Phi[k-1][j-1][i] - pG->Phi[k-1][j+1][i]) +
02897                     (pG->Phi[k  ][j-1][i] - pG->Phi[k  ][j+1][i]) )*(dx2i);
02898         gyr = 0.25*((pG->Phi[k  ][j-1][i] - pG->Phi[k  ][j+1][i]) +
02899                     (pG->Phi[k+1][j-1][i] - pG->Phi[k+1][j+1][i]) )*(dx2i);
02900 
02901         gzl = (pG->Phi[k-1][j][i] - pG->Phi[k  ][j][i])*(dx3i);
02902         gzr = (pG->Phi[k  ][j][i] - pG->Phi[k+1][j][i])*(dx3i);
02903 
02904 /* momentum fluxes in x3.  2nd term is needed only if Jean's swindle used */
02905         flx_m1l = gzl*gxl/four_pi_G;
02906         flx_m1r = gzr*gxr/four_pi_G;
02907 
02908         flx_m2l = gzl*gyl/four_pi_G;
02909         flx_m2r = gzr*gyr/four_pi_G;
02910 
02911         flx_m3l = 0.5*(gzl*gzl-gxl*gxl-gyl*gyl)/four_pi_G + grav_mean_rho*phil;
02912         flx_m3r = 0.5*(gzr*gzr-gxr*gxr-gyr*gyr)/four_pi_G + grav_mean_rho*phir;
02913 
02914 /* Update momenta and energy with d/dx3 terms  */
02915         pG->U[k][j][i].M1 -= dtodx3*(flx_m1r - flx_m1l);
02916         pG->U[k][j][i].M2 -= dtodx3*(flx_m2r - flx_m2l);
02917         pG->U[k][j][i].M3 -= dtodx3*(flx_m3r - flx_m3l);
02918 #ifndef BAROTROPIC
02919         pG->U[k][j][i].E -= dtodx3*(x3Flux[k  ][j][i].d*(phic - phil) +
02920                                     x3Flux[k+1][j][i].d*(phir - phic));
02921 #endif /* BAROTROPIC */
02922       }
02923     }
02924   }
02925 
02926 /* Save mass fluxes in Grid structure for source term correction in main loop */
02927 
02928   for (k=ks; k<=ke+1; k++) {
02929     for (j=js; j<=je+1; j++) {
02930       for (i=is; i<=ie+1; i++) {
02931         pG->x1MassFlux[k][j][i] = x1Flux[k][j][i].d;
02932         pG->x2MassFlux[k][j][i] = x2Flux[k][j][i].d;
02933         pG->x3MassFlux[k][j][i] = x3Flux[k][j][i].d;
02934       }
02935     }
02936   }
02937 #endif /* SELF_GRAVITY */
02938 
02939 /*--- Step 11c -----------------------------------------------------------------
02940  * Add source terms for optically thin cooling
02941  */
02942 
02943 #ifndef BAROTROPIC
02944   if (CoolingFunc != NULL){
02945     for (k=ks; k<=ke; k++){
02946       for (j=js; j<=je; j++){
02947         for (i=is; i<=ie; i++){
02948           coolf = (*CoolingFunc)(dhalf[k][j][i],phalf[k][j][i],pG->dt);
02949           pG->U[k][j][i].E -= pG->dt*coolf;
02950         }
02951       }
02952     }
02953   }
02954 #endif /* BAROTROPIC */
02955 
02956 /*--- Step 11d -----------------------------------------------------------------
02957  * Add source terms for particle feedback
02958  */
02959 
02960 #ifdef FEEDBACK
02961   for (k=ks; k<=ke; k++)
02962     for (j=js; j<=je; j++)
02963       for (i=is; i<=ie; i++) {
02964       pG->U[k][j][i].M1 -= pG->Coup[k][j][i].fb1;
02965       pG->U[k][j][i].M2 -= pG->Coup[k][j][i].fb2;
02966       pG->U[k][j][i].M3 -= pG->Coup[k][j][i].fb3;
02967 #ifndef BAROTROPIC
02968       pG->U[k][j][i].E += pG->Coup[k][j][i].Eloss;
02969       pG->Coup[k][j][i].Eloss *= dt1; /* for history output purpose */
02970       pG->Eloss[k][j][i] *= dt1; /* for history output purpose */
02971 #endif
02972     }
02973 #endif
02974 
02975 /*=== STEP 12: Update cell-centered values for a full timestep ===============*/
02976 
02977 /*--- Step 12a -----------------------------------------------------------------
02978  * Update cell-centered variables in pG using 3D x1-Fluxes
02979  */
02980 
02981   for (k=ks; k<=ke; k++) {
02982     for (j=js; j<=je; j++) {
02983       for (i=is; i<=ie; i++) {
02984 #ifdef CYLINDRICAL
02985         rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
02986 #endif
02987         pG->U[k][j][i].d  -= dtodx1*(rsf*x1Flux[k][j][i+1].d -lsf*x1Flux[k][j][i].d );
02988         pG->U[k][j][i].M1 -= dtodx1*(rsf*x1Flux[k][j][i+1].Mx-lsf*x1Flux[k][j][i].Mx);
02989         pG->U[k][j][i].M2 -= dtodx1*(SQR(rsf)*x1Flux[k][j][i+1].My-SQR(lsf)*x1Flux[k][j][i].My);
02990         pG->U[k][j][i].M3 -= dtodx1*(rsf*x1Flux[k][j][i+1].Mz-lsf*x1Flux[k][j][i].Mz);
02991 #ifndef BAROTROPIC
02992         pG->U[k][j][i].E  -= dtodx1*(rsf*x1Flux[k][j][i+1].E -lsf*x1Flux[k][j][i].E );
02993 #endif /* BAROTROPIC */
02994 #if (NSCALARS > 0)
02995         for (n=0; n<NSCALARS; n++)
02996           pG->U[k][j][i].s[n] -= dtodx1*(rsf*x1Flux[k][j][i+1].s[n]
02997                                        - lsf*x1Flux[k][j][i  ].s[n]);
02998 #endif
02999       }
03000     }
03001   }
03002 
03003 /*--- Step 12b -----------------------------------------------------------------
03004  * Update cell-centered variables in pG using 3D x2-Fluxes
03005  */
03006 
03007   for (k=ks; k<=ke; k++) {
03008     for (j=js; j<=je; j++) {
03009       for (i=is; i<=ie; i++) {
03010 #ifdef CYLINDRICAL
03011         dtodx2 = pG->dt/(r[i]*pG->dx2);
03012 #endif
03013         pG->U[k][j][i].d  -= dtodx2*(x2Flux[k][j+1][i].d -x2Flux[k][j][i].d );
03014         pG->U[k][j][i].M1 -= dtodx2*(x2Flux[k][j+1][i].Mz-x2Flux[k][j][i].Mz);
03015         pG->U[k][j][i].M2 -= dtodx2*(x2Flux[k][j+1][i].Mx-x2Flux[k][j][i].Mx);
03016         pG->U[k][j][i].M3 -= dtodx2*(x2Flux[k][j+1][i].My-x2Flux[k][j][i].My);
03017 #ifndef BAROTROPIC
03018         pG->U[k][j][i].E -=dtodx2*(x2Flux[k][j+1][i].E -x2Flux[k][j][i].E );
03019 #endif /* BAROTROPIC */
03020 #if (NSCALARS > 0)
03021         for (n=0; n<NSCALARS; n++)
03022           pG->U[k][j][i].s[n] -= dtodx2*(x2Flux[k][j+1][i].s[n]
03023                                        - x2Flux[k][j  ][i].s[n]);
03024 #endif
03025       }
03026     }
03027   }
03028 
03029 /*--- Step 12c -----------------------------------------------------------------
03030  * Update cell-centered variables in pG using 3D x3-Fluxes
03031  */
03032 
03033   for (k=ks; k<=ke; k++) {
03034     for (j=js; j<=je; j++) {
03035       for (i=is; i<=ie; i++) {
03036         pG->U[k][j][i].d  -= dtodx3*(x3Flux[k+1][j][i].d -x3Flux[k][j][i].d );
03037         pG->U[k][j][i].M1 -= dtodx3*(x3Flux[k+1][j][i].My-x3Flux[k][j][i].My);
03038         pG->U[k][j][i].M2 -= dtodx3*(x3Flux[k+1][j][i].Mz-x3Flux[k][j][i].Mz);
03039         pG->U[k][j][i].M3 -= dtodx3*(x3Flux[k+1][j][i].Mx-x3Flux[k][j][i].Mx);
03040 #ifndef BAROTROPIC
03041         pG->U[k][j][i].E  -= dtodx3*(x3Flux[k+1][j][i].E -x3Flux[k][j][i].E );
03042 #endif /* BAROTROPIC */
03043 #if (NSCALARS > 0)
03044         for (n=0; n<NSCALARS; n++)
03045           pG->U[k][j][i].s[n] -= dtodx3*(x3Flux[k+1][j][i].s[n]
03046                                        - x3Flux[k  ][j][i].s[n]);
03047 #endif
03048       }
03049     }
03050   }
03051 
03052 /*--- Step 12d -----------------------------------------------------------------
03053  * Set cell centered magnetic fields to average of updated face centered fields.
03054  */
03055 
03056 #ifdef MHD
03057   for (k=ks; k<=ke; k++) {
03058     for (j=js; j<=je; j++) {
03059       for (i=is; i<=ie; i++) {
03060 #ifdef CYLINDRICAL
03061         rsf = ri[i+1]/r[i];  lsf = ri[i]/r[i];
03062 #endif
03063         pG->U[k][j][i].B1c = 0.5*(lsf*pG->B1i[k][j][i] + rsf*pG->B1i[k][j][i+1]);
03064         pG->U[k][j][i].B2c = 0.5*(    pG->B2i[k][j][i] +     pG->B2i[k][j+1][i]);
03065         pG->U[k][j][i].B3c = 0.5*(    pG->B3i[k][j][i] +     pG->B3i[k+1][j][i]);
03066       }
03067     }
03068   }
03069 #endif /* MHD */
03070 
03071 #ifdef STATIC_MESH_REFINEMENT
03072 /*--- Step 12e -----------------------------------------------------------------
03073  * With SMR, store fluxes at boundaries of child and parent grids.  */
03074 /* Loop over all child grids ------------------*/
03075 
03076   for (ncg=0; ncg<pG->NCGrid; ncg++) {
03077 
03078 /* x1-boundaries of child Grids (interior to THIS Grid) */
03079 
03080     for (dim=0; dim<2; dim++){
03081       if (pG->CGrid[ncg].myFlx[dim] != NULL) {
03082 
03083         if (dim==0) i = pG->CGrid[ncg].ijks[0];
03084         if (dim==1) i = pG->CGrid[ncg].ijke[0] + 1;
03085         jcs = pG->CGrid[ncg].ijks[1];
03086         jce = pG->CGrid[ncg].ijke[1];
03087         kcs = pG->CGrid[ncg].ijks[2];
03088         kce = pG->CGrid[ncg].ijke[2];
03089 
03090         for (k=kcs, kk=0; k<=kce; k++, kk++){
03091           for (j=jcs, jj=0; j<=jce; j++, jj++){
03092             pG->CGrid[ncg].myFlx[dim][kk][jj].d  = x1Flux[k][j][i].d; 
03093             pG->CGrid[ncg].myFlx[dim][kk][jj].M1 = x1Flux[k][j][i].Mx; 
03094             pG->CGrid[ncg].myFlx[dim][kk][jj].M2 = x1Flux[k][j][i].My;
03095             pG->CGrid[ncg].myFlx[dim][kk][jj].M3 = x1Flux[k][j][i].Mz; 
03096 #ifndef BAROTROPIC
03097             pG->CGrid[ncg].myFlx[dim][kk][jj].E  = x1Flux[k][j][i].E; 
03098 #endif /* BAROTROPIC */
03099 #ifdef MHD
03100             pG->CGrid[ncg].myFlx[dim][kk][jj].B1c = 0.0;
03101             pG->CGrid[ncg].myFlx[dim][kk][jj].B2c = x1Flux[k][j][i].By; 
03102             pG->CGrid[ncg].myFlx[dim][kk][jj].B3c = x1Flux[k][j][i].Bz; 
03103 #endif /* MHD */
03104 #if (NSCALARS > 0)
03105             for (n=0; n<NSCALARS; n++)
03106               pG->CGrid[ncg].myFlx[dim][kk][jj].s[n]  = x1Flux[k][j][i].s[n]; 
03107 #endif
03108           }
03109         }
03110 #ifdef MHD
03111         for (k=kcs, kk=0; k<=kce+1; k++, kk++){
03112           for (j=jcs, jj=0; j<=jce; j++, jj++){
03113             pG->CGrid[ncg].myEMF2[dim][kk][jj] = emf2[k][j][i];
03114           }
03115         }
03116         for (k=kcs, kk=0; k<=kce; k++, kk++){
03117           for (j=jcs, jj=0; j<=jce+1; j++, jj++){
03118             pG->CGrid[ncg].myEMF3[dim][kk][jj] = emf3[k][j][i];
03119           }
03120         }
03121 #endif /* MHD */
03122       }
03123     }
03124 
03125 /* x2-boundaries of child Grids (interior to THIS Grid) */
03126 
03127     for (dim=2; dim<4; dim++){
03128       if (pG->CGrid[ncg].myFlx[dim] != NULL) {
03129 
03130         ics = pG->CGrid[ncg].ijks[0];
03131         ice = pG->CGrid[ncg].ijke[0];
03132         if (dim==2) j = pG->CGrid[ncg].ijks[1];
03133         if (dim==3) j = pG->CGrid[ncg].ijke[1] + 1;
03134         kcs = pG->CGrid[ncg].ijks[2];
03135         kce = pG->CGrid[ncg].ijke[2];
03136 
03137         for (k=kcs, kk=0; k<=kce; k++, kk++){
03138           for (i=ics, ii=0; i<=ice; i++, ii++){
03139             pG->CGrid[ncg].myFlx[dim][kk][ii].d  = x2Flux[k][j][i].d; 
03140             pG->CGrid[ncg].myFlx[dim][kk][ii].M1 = x2Flux[k][j][i].Mz; 
03141             pG->CGrid[ncg].myFlx[dim][kk][ii].M2 = x2Flux[k][j][i].Mx;
03142             pG->CGrid[ncg].myFlx[dim][kk][ii].M3 = x2Flux[k][j][i].My; 
03143 #ifndef BAROTROPIC
03144             pG->CGrid[ncg].myFlx[dim][kk][ii].E  = x2Flux[k][j][i].E; 
03145 #endif /* BAROTROPIC */
03146 #ifdef MHD
03147             pG->CGrid[ncg].myFlx[dim][kk][ii].B1c = x2Flux[k][j][i].Bz; 
03148             pG->CGrid[ncg].myFlx[dim][kk][ii].B2c = 0.0;
03149             pG->CGrid[ncg].myFlx[dim][kk][ii].B3c = x2Flux[k][j][i].By; 
03150 #endif /* MHD */
03151 #if (NSCALARS > 0)
03152             for (n=0; n<NSCALARS; n++)
03153               pG->CGrid[ncg].myFlx[dim][kk][ii].s[n]  = x2Flux[k][j][i].s[n]; 
03154 #endif
03155           }
03156         }
03157 #ifdef MHD
03158         for (k=kcs, kk=0; k<=kce+1; k++, kk++){
03159           for (i=ics, ii=0; i<=ice; i++, ii++){
03160             pG->CGrid[ncg].myEMF1[dim][kk][ii] = emf1[k][j][i];
03161           }
03162         }
03163         for (k=kcs, kk=0; k<=kce; k++, kk++){
03164           for (i=ics, ii=0; i<=ice+1; i++, ii++){
03165             pG->CGrid[ncg].myEMF3[dim][kk][ii] = emf3[k][j][i];
03166           }
03167         }
03168 #endif /* MHD */
03169       }
03170     }
03171 
03172 /* x3-boundaries of child Grids (interior to THIS Grid) */
03173 
03174     for (dim=4; dim<6; dim++){
03175       if (pG->CGrid[ncg].myFlx[dim] != NULL) {
03176 
03177         ics = pG->CGrid[ncg].ijks[0];
03178         ice = pG->CGrid[ncg].ijke[0];
03179         jcs = pG->CGrid[ncg].ijks[1];
03180         jce = pG->CGrid[ncg].ijke[1];
03181         if (dim==4) k = pG->CGrid[ncg].ijks[2];
03182         if (dim==5) k = pG->CGrid[ncg].ijke[2] + 1;
03183 
03184         for (j=jcs, jj=0; j<=jce; j++, jj++){
03185           for (i=ics, ii=0; i<=ice; i++, ii++){
03186             pG->CGrid[ncg].myFlx[dim][jj][ii].d  = x3Flux[k][j][i].d; 
03187             pG->CGrid[ncg].myFlx[dim][jj][ii].M1 = x3Flux[k][j][i].My; 
03188             pG->CGrid[ncg].myFlx[dim][jj][ii].M2 = x3Flux[k][j][i].Mz;
03189             pG->CGrid[ncg].myFlx[dim][jj][ii].M3 = x3Flux[k][j][i].Mx; 
03190 #ifndef BAROTROPIC
03191             pG->CGrid[ncg].myFlx[dim][jj][ii].E  = x3Flux[k][j][i].E; 
03192 #endif /* BAROTROPIC */
03193 #ifdef MHD
03194             pG->CGrid[ncg].myFlx[dim][jj][ii].B1c = x3Flux[k][j][i].By; 
03195             pG->CGrid[ncg].myFlx[dim][jj][ii].B2c = x3Flux[k][j][i].Bz; 
03196             pG->CGrid[ncg].myFlx[dim][jj][ii].B3c = 0.0;
03197 #endif /* MHD */
03198 #if (NSCALARS > 0)
03199             for (n=0; n<NSCALARS; n++)
03200               pG->CGrid[ncg].myFlx[dim][jj][ii].s[n]  = x3Flux[k][j][i].s[n]; 
03201 #endif
03202           }
03203         }
03204 #ifdef MHD
03205         for (j=jcs, jj=0; j<=jce+1; j++, jj++){
03206           for (i=ics, ii=0; i<=ice; i++, ii++){
03207             pG->CGrid[ncg].myEMF1[dim][jj][ii] = emf1[k][j][i];
03208           }
03209         }
03210         for (j=jcs, jj=0; j<=jce; j++, jj++){
03211           for (i=ics, ii=0; i<=ice+1; i++, ii++){
03212             pG->CGrid[ncg].myEMF2[dim][jj][ii] = emf2[k][j][i];
03213           }
03214         }
03215 #endif /* MHD */
03216       }
03217     }
03218   } /* end loop over child Grids */
03219 
03220 /* Loop over all parent grids ------------------*/
03221 
03222   for (npg=0; npg<pG->NPGrid; npg++) {
03223 
03224 /* x1-boundaries of parent Grids (at boundaries of THIS Grid)  */
03225 
03226     for (dim=0; dim<2; dim++){
03227       if (pG->PGrid[npg].myFlx[dim] != NULL) {
03228 
03229         if (dim==0) i = pG->PGrid[npg].ijks[0];
03230         if (dim==1) i = pG->PGrid[npg].ijke[0] + 1;
03231         jps = pG->PGrid[npg].ijks[1];
03232         jpe = pG->PGrid[npg].ijke[1];
03233         kps = pG->PGrid[npg].ijks[2];
03234         kpe = pG->PGrid[npg].ijke[2];
03235 
03236         for (k=kps, kk=0; k<=kpe; k++, kk++){
03237           for (j=jps, jj=0; j<=jpe; j++, jj++){
03238             pG->PGrid[npg].myFlx[dim][kk][jj].d  = x1Flux[k][j][i].d; 
03239             pG->PGrid[npg].myFlx[dim][kk][jj].M1 = x1Flux[k][j][i].Mx; 
03240             pG->PGrid[npg].myFlx[dim][kk][jj].M2 = x1Flux[k][j][i].My;
03241             pG->PGrid[npg].myFlx[dim][kk][jj].M3 = x1Flux[k][j][i].Mz; 
03242 #ifndef BAROTROPIC
03243             pG->PGrid[npg].myFlx[dim][kk][jj].E  = x1Flux[k][j][i].E; 
03244 #endif /* BAROTROPIC */
03245 #ifdef MHD
03246             pG->PGrid[npg].myFlx[dim][kk][jj].B1c = 0.0;
03247             pG->PGrid[npg].myFlx[dim][kk][jj].B2c = x1Flux[k][j][i].By; 
03248             pG->PGrid[npg].myFlx[dim][kk][jj].B3c = x1Flux[k][j][i].Bz; 
03249 #endif /* MHD */
03250 #if (NSCALARS > 0)
03251             for (n=0; n<NSCALARS; n++)
03252               pG->PGrid[npg].myFlx[dim][kk][jj].s[n]  = x1Flux[k][j][i].s[n]; 
03253 #endif
03254           }
03255         }
03256 #ifdef MHD
03257         for (k=kps, kk=0; k<=kpe+1; k++, kk++){
03258           for (j=jps, jj=0; j<=jpe; j++, jj++){
03259             pG->PGrid[npg].myEMF2[dim][kk][jj] = emf2[k][j][i];
03260           }
03261         }
03262         for (k=kps, kk=0; k<=kpe; k++, kk++){
03263           for (j=jps, jj=0; j<=jpe+1; j++, jj++){
03264             pG->PGrid[npg].myEMF3[dim][kk][jj] = emf3[k][j][i];
03265           }
03266         }
03267 #endif /* MHD */
03268       }
03269     }
03270 
03271 /* x2-boundaries of parent Grids (at boundaries of THIS Grid)  */
03272 
03273     for (dim=2; dim<4; dim++){
03274       if (pG->PGrid[npg].myFlx[dim] != NULL) {
03275 
03276         ips = pG->PGrid[npg].ijks[0];
03277         ipe = pG->PGrid[npg].ijke[0];
03278         if (dim==2) j = pG->PGrid[npg].ijks[1];
03279         if (dim==3) j = pG->PGrid[npg].ijke[1] + 1;
03280         kps = pG->PGrid[npg].ijks[2];
03281         kpe = pG->PGrid[npg].ijke[2];
03282 
03283         for (k=kps, kk=0; k<=kpe; k++, kk++){
03284           for (i=ips, ii=0; i<=ipe; i++, ii++){
03285             pG->PGrid[npg].myFlx[dim][kk][ii].d  = x2Flux[k][j][i].d; 
03286             pG->PGrid[npg].myFlx[dim][kk][ii].M1 = x2Flux[k][j][i].Mz; 
03287             pG->PGrid[npg].myFlx[dim][kk][ii].M2 = x2Flux[k][j][i].Mx;
03288             pG->PGrid[npg].myFlx[dim][kk][ii].M3 = x2Flux[k][j][i].My; 
03289 #ifndef BAROTROPIC
03290             pG->PGrid[npg].myFlx[dim][kk][ii].E  = x2Flux[k][j][i].E; 
03291 #endif /* BAROTROPIC */
03292 #ifdef MHD
03293             pG->PGrid[npg].myFlx[dim][kk][ii].B1c = x2Flux[k][j][i].Bz; 
03294             pG->PGrid[npg].myFlx[dim][kk][ii].B2c = 0.0;
03295             pG->PGrid[npg].myFlx[dim][kk][ii].B3c = x2Flux[k][j][i].By; 
03296 #endif /* MHD */
03297 #if (NSCALARS > 0)
03298             for (n=0; n<NSCALARS; n++)
03299               pG->PGrid[npg].myFlx[dim][kk][ii].s[n]  = x2Flux[k][j][i].s[n]; 
03300 #endif
03301           }
03302         }
03303 #ifdef MHD
03304         for (k=kps, kk=0; k<=kpe+1; k++, kk++){
03305           for (i=ips, ii=0; i<=ipe; i++, ii++){
03306             pG->PGrid[npg].myEMF1[dim][kk][ii] = emf1[k][j][i];
03307           }
03308         }
03309         for (k=kps, kk=0; k<=kpe; k++, kk++){
03310           for (i=ips, ii=0; i<=ipe+1; i++, ii++){
03311             pG->PGrid[npg].myEMF3[dim][kk][ii] = emf3[k][j][i];
03312           }
03313         }
03314 #endif /* MHD */
03315       }
03316     }
03317 
03318 /* x3-boundaries of parent Grids (at boundaries of THIS Grid)  */
03319 
03320     for (dim=4; dim<6; dim++){
03321       if (pG->PGrid[npg].myFlx[dim] != NULL) {
03322 
03323         ips = pG->PGrid[npg].ijks[0];
03324         ipe = pG->PGrid[npg].ijke[0];
03325         jps = pG->PGrid[npg].ijks[1];
03326         jpe = pG->PGrid[npg].ijke[1];
03327         if (dim==4) k = pG->PGrid[npg].ijks[2];
03328         if (dim==5) k = pG->PGrid[npg].ijke[2] + 1;
03329 
03330         for (j=jps, jj=0; j<=jpe; j++, jj++){
03331           for (i=ips, ii=0; i<=ipe; i++, ii++){
03332             pG->PGrid[npg].myFlx[dim][jj][ii].d  = x3Flux[k][j][i].d; 
03333             pG->PGrid[npg].myFlx[dim][jj][ii].M1 = x3Flux[k][j][i].My; 
03334             pG->PGrid[npg].myFlx[dim][jj][ii].M2 = x3Flux[k][j][i].Mz;
03335             pG->PGrid[npg].myFlx[dim][jj][ii].M3 = x3Flux[k][j][i].Mx; 
03336 #ifndef BAROTROPIC
03337             pG->PGrid[npg].myFlx[dim][jj][ii].E  = x3Flux[k][j][i].E; 
03338 #endif /* BAROTROPIC */
03339 #ifdef MHD
03340             pG->PGrid[npg].myFlx[dim][jj][ii].B1c = x3Flux[k][j][i].By; 
03341             pG->PGrid[npg].myFlx[dim][jj][ii].B2c = x3Flux[k][j][i].Bz; 
03342             pG->PGrid[npg].myFlx[dim][jj][ii].B3c = 0.0;
03343 #endif /* MHD */
03344 #if (NSCALARS > 0)
03345             for (n=0; n<NSCALARS; n++)
03346               pG->PGrid[npg].myFlx[dim][jj][ii].s[n]  = x3Flux[k][j][i].s[n]; 
03347 #endif
03348           }
03349         }
03350 #ifdef MHD
03351         for (j=jps, jj=0; j<=jpe+1; j++, jj++){
03352           for (i=ips, ii=0; i<=ipe; i++, ii++){
03353             pG->PGrid[npg].myEMF1[dim][jj][ii] = emf1[k][j][i];
03354           }
03355         }
03356         for (j=jps, jj=0; j<=jpe; j++, jj++){
03357           for (i=ips, ii=0; i<=ipe+1; i++, ii++){
03358             pG->PGrid[npg].myEMF2[dim][jj][ii] = emf2[k][j][i];
03359           }
03360         }
03361 #endif /* MHD */
03362       }
03363     }
03364   }
03365 
03366 #endif /* STATIC_MESH_REFINEMENT */
03367   return;
03368 }
03369 
03370 /*----------------------------------------------------------------------------*/
03371 /*! \fn void integrate_init_3d(MeshS *pM)
03372  *  \brief Allocate temporary integration arrays 
03373 */
03374 void integrate_init_3d(MeshS *pM)
03375 {
03376   int nmax,size1=0,size2=0,size3=0,nl,nd;
03377 
03378 /* Cycle over all Grids on this processor to find maximum Nx1, Nx2, Nx3 */
03379   for (nl=0; nl<(pM->NLevels); nl++){
03380     for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
03381       if (pM->Domain[nl][nd].Grid != NULL) {
03382         if (pM->Domain[nl][nd].Grid->Nx[0] > size1){
03383           size1 = pM->Domain[nl][nd].Grid->Nx[0];
03384         }
03385         if (pM->Domain[nl][nd].Grid->Nx[1] > size2){
03386           size2 = pM->Domain[nl][nd].Grid->Nx[1];
03387         }
03388         if (pM->Domain[nl][nd].Grid->Nx[2] > size3){
03389           size3 = pM->Domain[nl][nd].Grid->Nx[2];
03390         }
03391       }
03392     }
03393   }
03394 
03395   size1 = size1 + 2*nghost;
03396   size2 = size2 + 2*nghost;
03397   size3 = size3 + 2*nghost;
03398   nmax = MAX((MAX(size1,size2)),size3);
03399 
03400 #ifdef MHD
03401   if ((emf1 = (Real***)calloc_3d_array(size3,size2,size1,sizeof(Real)))==NULL)
03402     goto on_error;
03403   if ((emf2 = (Real***)calloc_3d_array(size3,size2,size1,sizeof(Real)))==NULL)
03404     goto on_error;
03405   if ((emf3 = (Real***)calloc_3d_array(size3,size2,size1,sizeof(Real)))==NULL)
03406     goto on_error;
03407 
03408   if ((emf1_cc=(Real***)calloc_3d_array(size3,size2,size1,sizeof(Real)))==NULL)
03409     goto on_error;
03410   if ((emf2_cc=(Real***)calloc_3d_array(size3,size2,size1,sizeof(Real)))==NULL)
03411     goto on_error;
03412   if ((emf3_cc=(Real***)calloc_3d_array(size3,size2,size1,sizeof(Real)))==NULL)
03413     goto on_error;
03414 #endif /* MHD */
03415 #ifdef H_CORRECTION
03416   if ((eta1 = (Real***)calloc_3d_array(size3,size2,size1,sizeof(Real)))==NULL)
03417     goto on_error;
03418   if ((eta2 = (Real***)calloc_3d_array(size3,size2,size1,sizeof(Real)))==NULL)
03419     goto on_error;
03420   if ((eta3 = (Real***)calloc_3d_array(size3,size2,size1,sizeof(Real)))==NULL)
03421     goto on_error;
03422 #endif /* H_CORRECTION */
03423 
03424   if ((Bxc = (Real*)malloc(nmax*sizeof(Real))) == NULL) goto on_error;
03425   if ((Bxi = (Real*)malloc(nmax*sizeof(Real))) == NULL) goto on_error;
03426 
03427 #ifdef MHD
03428   if ((B1_x1Face = (Real***)calloc_3d_array(size3,size2,size1, sizeof(Real)))
03429     == NULL) goto on_error;
03430   if ((B2_x2Face = (Real***)calloc_3d_array(size3,size2,size1, sizeof(Real)))
03431     == NULL) goto on_error;
03432   if ((B3_x3Face = (Real***)calloc_3d_array(size3,size2,size1, sizeof(Real)))
03433     == NULL) goto on_error;
03434 #endif /* MHD */
03435 
03436   if ((U1d=(Cons1DS*)malloc(nmax*sizeof(Cons1DS))) == NULL) goto on_error;
03437   if ((W  =(Prim1DS*)malloc(nmax*sizeof(Prim1DS))) == NULL) goto on_error;
03438   if ((Wl =(Prim1DS*)malloc(nmax*sizeof(Prim1DS))) == NULL) goto on_error;
03439   if ((Wr =(Prim1DS*)malloc(nmax*sizeof(Prim1DS))) == NULL) goto on_error;
03440 
03441   if ((Ul_x1Face=(Cons1DS***)calloc_3d_array(size3,size2,size1,sizeof(Cons1DS)))
03442     == NULL) goto on_error;
03443   if ((Ur_x1Face=(Cons1DS***)calloc_3d_array(size3,size2,size1,sizeof(Cons1DS)))
03444     == NULL) goto on_error;
03445   if ((Ul_x2Face=(Cons1DS***)calloc_3d_array(size3,size2,size1,sizeof(Cons1DS)))
03446     == NULL) goto on_error;
03447   if ((Ur_x2Face=(Cons1DS***)calloc_3d_array(size3,size2,size1,sizeof(Cons1DS)))
03448     == NULL) goto on_error;
03449   if ((Ul_x3Face=(Cons1DS***)calloc_3d_array(size3,size2,size1,sizeof(Cons1DS)))
03450     == NULL) goto on_error;
03451   if ((Ur_x3Face=(Cons1DS***)calloc_3d_array(size3,size2,size1,sizeof(Cons1DS)))
03452     == NULL) goto on_error;
03453   if ((x1Flux   =(Cons1DS***)calloc_3d_array(size3,size2,size1,sizeof(Cons1DS)))
03454     == NULL) goto on_error;
03455   if ((x2Flux   =(Cons1DS***)calloc_3d_array(size3,size2,size1,sizeof(Cons1DS)))
03456     == NULL) goto on_error;
03457   if ((x3Flux   =(Cons1DS***)calloc_3d_array(size3,size2,size1,sizeof(Cons1DS)))
03458     == NULL) goto on_error;
03459 
03460 #ifdef CYLINDRICAL
03461 #ifndef MHD
03462 #ifndef PARTICLES
03463   if((StaticGravPot != NULL) || (CoolingFunc != NULL))
03464 #endif
03465 #endif
03466 #endif
03467   {
03468   if ((dhalf = (Real***)calloc_3d_array(size3,size2,size1,sizeof(Real)))==NULL)
03469     goto on_error;
03470   if ((phalf = (Real***)calloc_3d_array(size3,size2,size1,sizeof(Real)))==NULL)
03471     goto on_error;
03472   }
03473 
03474 #ifdef SHEARING_BOX
03475   if ((remapEyiib = (Real**)calloc_2d_array(size3,size2, sizeof(Real))) == NULL)
03476     goto on_error;
03477   if ((remapEyoib = (Real**)calloc_2d_array(size3,size2, sizeof(Real))) == NULL)
03478     goto on_error;
03479 #endif
03480 
03481   /* DATA STRUCTURES FOR CYLINDRICAL COORDINATES */
03482 #ifdef CYLINDRICAL
03483   if ((geom_src = (Real***)calloc_3d_array(size3, size2, size1, sizeof(Real))) == NULL)
03484     goto on_error;
03485 #endif
03486 
03487   return;
03488 
03489   on_error:
03490     integrate_destruct();
03491     ath_error("[integrate_init]: malloc returned a NULL pointer\n");
03492 }
03493 
03494 /*----------------------------------------------------------------------------*/
03495 /*! \fn void integrate_destruct_3d(void)
03496  *  \brief Free temporary integration arrays 
03497  */
03498 void integrate_destruct_3d(void)
03499 {
03500 
03501 #ifdef MHD
03502   if (emf1    != NULL) free_3d_array(emf1);
03503   if (emf2    != NULL) free_3d_array(emf2);
03504   if (emf3    != NULL) free_3d_array(emf3);
03505   if (emf1_cc != NULL) free_3d_array(emf1_cc);
03506   if (emf2_cc != NULL) free_3d_array(emf2_cc);
03507   if (emf3_cc != NULL) free_3d_array(emf3_cc);
03508 #endif /* MHD */
03509 #ifdef H_CORRECTION
03510   if (eta1 != NULL) free_3d_array(eta1);
03511   if (eta2 != NULL) free_3d_array(eta2);
03512   if (eta3 != NULL) free_3d_array(eta3);
03513 #endif /* H_CORRECTION */
03514 
03515   if (Bxc != NULL) free(Bxc);
03516   if (Bxi != NULL) free(Bxi);
03517 #ifdef MHD
03518   if (B1_x1Face != NULL) free_3d_array(B1_x1Face);
03519   if (B2_x2Face != NULL) free_3d_array(B2_x2Face);
03520   if (B3_x3Face != NULL) free_3d_array(B3_x3Face);
03521 #endif /* MHD */
03522 
03523   if (U1d      != NULL) free(U1d);
03524   if (W        != NULL) free(W);
03525   if (Wl       != NULL) free(Wl);
03526   if (Wr       != NULL) free(Wr);
03527 
03528   if (Ul_x1Face != NULL) free_3d_array(Ul_x1Face);
03529   if (Ur_x1Face != NULL) free_3d_array(Ur_x1Face);
03530   if (Ul_x2Face != NULL) free_3d_array(Ul_x2Face);
03531   if (Ur_x2Face != NULL) free_3d_array(Ur_x2Face);
03532   if (Ul_x3Face != NULL) free_3d_array(Ul_x3Face);
03533   if (Ur_x3Face != NULL) free_3d_array(Ur_x3Face);
03534   if (x1Flux    != NULL) free_3d_array(x1Flux);
03535   if (x2Flux    != NULL) free_3d_array(x2Flux);
03536   if (x3Flux    != NULL) free_3d_array(x3Flux);
03537   if (dhalf     != NULL) free_3d_array(dhalf);
03538   if (phalf     != NULL) free_3d_array(phalf);
03539 #ifdef SHEARING_BOX
03540   if (remapEyiib != NULL) free_2d_array(remapEyiib);
03541   if (remapEyoib != NULL) free_2d_array(remapEyoib);
03542 #endif
03543 
03544   /* DATA STRUCTURES FOR CYLINDRICAL COORDINATES */
03545 #ifdef CYLINDRICAL
03546   if (geom_src  != NULL) free_3d_array(geom_src);
03547 #endif
03548 
03549   return;
03550 }
03551 
03552 /*=========================== PRIVATE FUNCTIONS ==============================*/
03553 
03554 /*----------------------------------------------------------------------------*/
03555 /*! \fn static void integrate_emf1_corner(const GridS *pG)
03556  *  \brief Integrates face centered B-fluxes to compute corner EMFs.  
03557  *
03558  *  Note: 
03559  * - x1Flux.By = VxBy - BxVy = v1*b2-b1*v2 = -EMFZ
03560  * - x1Flux.Bz = VxBz - BxVz = v1*b3-b1*v3 = EMFY
03561  * - x2Flux.By = VxBy - BxVy = v2*b3-b2*v3 = -EMFX
03562  * - x2Flux.Bz = VxBz - BxVz = v2*b1-b2*v1 = EMFZ
03563  * - x3Flux.By = VxBy - BxVy = v3*b1-b3*v1 = -EMFY
03564  * - x3Flux.Bz = VxBz - BxVz = v3*b2-b3*v2 = EMFX 
03565  */
03566 #ifdef MHD
03567 static void integrate_emf1_corner(const GridS *pG)
03568 {
03569   int i, is = pG->is, ie = pG->ie;
03570   int j, js = pG->js, je = pG->je;
03571   int k, ks = pG->ks, ke = pG->ke;
03572   Real de1_l2, de1_r2, de1_l3, de1_r3;
03573 
03574   for (k=ks-1; k<=ke+2; k++) {
03575     for (j=js-1; j<=je+2; j++) {
03576       for (i=is-2; i<=ie+2; i++) {
03577 /* NOTE: The x2-Flux of By is -E1. */
03578 /*       The x3-Flux of Bz is +E1. */
03579         if (x2Flux[k-1][j][i].d > 0.0)
03580           de1_l3 = x3Flux[k][j-1][i].Bz - emf1_cc[k-1][j-1][i];
03581         else if (x2Flux[k-1][j][i].d < 0.0)
03582           de1_l3 = x3Flux[k][j][i].Bz - emf1_cc[k-1][j][i];
03583         else {
03584           de1_l3 = 0.5*(x3Flux[k][j-1][i].Bz - emf1_cc[k-1][j-1][i] +
03585                         x3Flux[k][j  ][i].Bz - emf1_cc[k-1][j  ][i] );
03586         }
03587 
03588         if (x2Flux[k][j][i].d > 0.0)
03589           de1_r3 = x3Flux[k][j-1][i].Bz - emf1_cc[k][j-1][i];
03590         else if (x2Flux[k][j][i].d < 0.0)
03591           de1_r3 = x3Flux[k][j][i].Bz - emf1_cc[k][j][i];
03592         else {
03593           de1_r3 = 0.5*(x3Flux[k][j-1][i].Bz - emf1_cc[k][j-1][i] +
03594                         x3Flux[k][j  ][i].Bz - emf1_cc[k][j  ][i] );
03595         }
03596 
03597         if (x3Flux[k][j-1][i].d > 0.0)
03598           de1_l2 = -x2Flux[k-1][j][i].By - emf1_cc[k-1][j-1][i];
03599         else if (x3Flux[k][j-1][i].d < 0.0)
03600           de1_l2 = -x2Flux[k][j][i].By - emf1_cc[k][j-1][i];
03601         else {
03602           de1_l2 = 0.5*(-x2Flux[k-1][j][i].By - emf1_cc[k-1][j-1][i]
03603                         -x2Flux[k  ][j][i].By - emf1_cc[k  ][j-1][i] );
03604         }
03605 
03606         if (x3Flux[k][j][i].d > 0.0)
03607           de1_r2 = -x2Flux[k-1][j][i].By - emf1_cc[k-1][j][i];
03608         else if (x3Flux[k][j][i].d < 0.0)
03609           de1_r2 = -x2Flux[k][j][i].By - emf1_cc[k][j][i];
03610         else {
03611           de1_r2 = 0.5*(-x2Flux[k-1][j][i].By - emf1_cc[k-1][j][i]
03612                         -x2Flux[k  ][j][i].By - emf1_cc[k  ][j][i] );
03613         }
03614 
03615         emf1[k][j][i] = 0.25*(  x3Flux[k][j][i].Bz + x3Flux[k][j-1][i].Bz
03616                               - x2Flux[k][j][i].By - x2Flux[k-1][j][i].By 
03617                               + de1_l2 + de1_r2 + de1_l3 + de1_r3);
03618       }
03619     }
03620   }
03621 
03622   return;
03623 }
03624 
03625 /*! \fn static void integrate_emf2_corner(const GridS *pG)
03626  *  \brief Integrates face centered B-fluxes to compute corner EMFs.  
03627  *
03628  *  Note: 
03629  * - x1Flux.By = VxBy - BxVy = v1*b2-b1*v2 = -EMFZ
03630  * - x1Flux.Bz = VxBz - BxVz = v1*b3-b1*v3 = EMFY
03631  * - x2Flux.By = VxBy - BxVy = v2*b3-b2*v3 = -EMFX
03632  * - x2Flux.Bz = VxBz - BxVz = v2*b1-b2*v1 = EMFZ
03633  * - x3Flux.By = VxBy - BxVy = v3*b1-b3*v1 = -EMFY
03634  * - x3Flux.Bz = VxBz - BxVz = v3*b2-b3*v2 = EMFX 
03635  */
03636 static void integrate_emf2_corner(const GridS *pG)
03637 {
03638   int i, is = pG->is, ie = pG->ie;
03639   int j, js = pG->js, je = pG->je;
03640   int k, ks = pG->ks, ke = pG->ke;
03641   Real de2_l1, de2_r1, de2_l3, de2_r3;
03642 
03643   for (k=ks-1; k<=ke+2; k++) {
03644     for (j=js-2; j<=je+2; j++) {
03645       for (i=is-1; i<=ie+2; i++) {
03646 /* NOTE: The x1-Flux of Bz is +E2. */
03647 /*       The x3-Flux of By is -E2. */
03648         if (x1Flux[k-1][j][i].d > 0.0)
03649           de2_l3 = -x3Flux[k][j][i-1].By - emf2_cc[k-1][j][i-1];
03650         else if (x1Flux[k-1][j][i].d < 0.0)
03651           de2_l3 = -x3Flux[k][j][i].By - emf2_cc[k-1][j][i];
03652         else {
03653           de2_l3 = 0.5*(-x3Flux[k][j][i-1].By - emf2_cc[k-1][j][i-1] 
03654                         -x3Flux[k][j][i  ].By - emf2_cc[k-1][j][i  ] );
03655         }
03656 
03657         if (x1Flux[k][j][i].d > 0.0)
03658           de2_r3 = -x3Flux[k][j][i-1].By - emf2_cc[k][j][i-1];
03659         else if (x1Flux[k][j][i].d < 0.0)
03660           de2_r3 = -x3Flux[k][j][i].By - emf2_cc[k][j][i];
03661         else {
03662           de2_r3 = 0.5*(-x3Flux[k][j][i-1].By - emf2_cc[k][j][i-1] 
03663                         -x3Flux[k][j][i  ].By - emf2_cc[k][j][i  ] );
03664         }
03665 
03666         if (x3Flux[k][j][i-1].d > 0.0)
03667           de2_l1 = x1Flux[k-1][j][i].Bz - emf2_cc[k-1][j][i-1];
03668         else if (x3Flux[k][j][i-1].d < 0.0)
03669           de2_l1 = x1Flux[k][j][i].Bz - emf2_cc[k][j][i-1];
03670         else {
03671           de2_l1 = 0.5*(x1Flux[k-1][j][i].Bz - emf2_cc[k-1][j][i-1] +
03672                         x1Flux[k  ][j][i].Bz - emf2_cc[k  ][j][i-1] );
03673         }
03674 
03675         if (x3Flux[k][j][i].d > 0.0)
03676           de2_r1 = x1Flux[k-1][j][i].Bz - emf2_cc[k-1][j][i];
03677         else if (x3Flux[k][j][i].d < 0.0)
03678           de2_r1 = x1Flux[k][j][i].Bz - emf2_cc[k][j][i];
03679         else {
03680           de2_r1 = 0.5*(x1Flux[k-1][j][i].Bz - emf2_cc[k-1][j][i] +
03681                         x1Flux[k  ][j][i].Bz - emf2_cc[k  ][j][i] );
03682         }
03683 
03684         emf2[k][j][i] = 0.25*(  x1Flux[k][j][i].Bz + x1Flux[k-1][j][i  ].Bz
03685                               - x3Flux[k][j][i].By - x3Flux[k  ][j][i-1].By
03686                               + de2_l1 + de2_r1 + de2_l3 + de2_r3);
03687       }
03688     }
03689   }
03690 
03691   return;
03692 }
03693 
03694 /*! \fn static void integrate_emf3_corner(const GridS *pG)
03695  *  \brief Integrates face centered B-fluxes to compute corner EMFs.  
03696  *
03697  *  Note: 
03698  * - x1Flux.By = VxBy - BxVy = v1*b2-b1*v2 = -EMFZ
03699  * - x1Flux.Bz = VxBz - BxVz = v1*b3-b1*v3 = EMFY
03700  * - x2Flux.By = VxBy - BxVy = v2*b3-b2*v3 = -EMFX
03701  * - x2Flux.Bz = VxBz - BxVz = v2*b1-b2*v1 = EMFZ
03702  * - x3Flux.By = VxBy - BxVy = v3*b1-b3*v1 = -EMFY
03703  * - x3Flux.Bz = VxBz - BxVz = v3*b2-b3*v2 = EMFX 
03704  */
03705 static void integrate_emf3_corner(const GridS *pG)
03706 {
03707   int i, is = pG->is, ie = pG->ie;
03708   int j, js = pG->js, je = pG->je;
03709   int k, ks = pG->ks, ke = pG->ke;
03710   Real de3_l1, de3_r1, de3_l2, de3_r2;
03711   Real rsf=1.0,lsf=1.0;
03712 
03713   for (k=ks-2; k<=ke+2; k++) {
03714     for (j=js-1; j<=je+2; j++) {
03715       for (i=is-1; i<=ie+2; i++) {
03716 /* NOTE: The x1-Flux of By is -E3. */
03717 /*       The x2-Flux of Bx is +E3. */
03718 #ifdef CYLINDRICAL
03719         rsf = pG->ri[i]/pG->r[i];  lsf = pG->ri[i]/pG->r[i-1];
03720 #endif
03721         if (x1Flux[k][j-1][i].d > 0.0)
03722           de3_l2 = (x2Flux[k][j][i-1].Bz - emf3_cc[k][j-1][i-1])*lsf;
03723         else if (x1Flux[k][j-1][i].d < 0.0)
03724           de3_l2 = (x2Flux[k][j][i].Bz - emf3_cc[k][j-1][i])*rsf;
03725         else {
03726           de3_l2 = 0.5*((x2Flux[k][j][i-1].Bz - emf3_cc[k][j-1][i-1])*lsf + 
03727                         (x2Flux[k][j][i  ].Bz - emf3_cc[k][j-1][i  ])*rsf );
03728         }
03729 
03730         if (x1Flux[k][j][i].d > 0.0)
03731           de3_r2 = (x2Flux[k][j][i-1].Bz - emf3_cc[k][j][i-1])*lsf;
03732         else if (x1Flux[k][j][i].d < 0.0)
03733           de3_r2 = (x2Flux[k][j][i].Bz - emf3_cc[k][j][i])*rsf;
03734         else {
03735           de3_r2 = 0.5*((x2Flux[k][j][i-1].Bz - emf3_cc[k][j][i-1])*lsf + 
03736                         (x2Flux[k][j][i  ].Bz - emf3_cc[k][j][i  ])*rsf );
03737         }
03738 
03739         if (x2Flux[k][j][i-1].d > 0.0)
03740           de3_l1 = -x1Flux[k][j-1][i].By - emf3_cc[k][j-1][i-1];
03741         else if (x2Flux[k][j][i-1].d < 0.0)
03742           de3_l1 = -x1Flux[k][j][i].By - emf3_cc[k][j][i-1];
03743         else {
03744           de3_l1 = 0.5*(-x1Flux[k][j-1][i].By - emf3_cc[k][j-1][i-1]
03745                         -x1Flux[k][j  ][i].By - emf3_cc[k][j  ][i-1] );
03746         }
03747 
03748         if (x2Flux[k][j][i].d > 0.0)
03749           de3_r1 = -x1Flux[k][j-1][i].By - emf3_cc[k][j-1][i];
03750         else if (x2Flux[k][j][i].d < 0.0)
03751           de3_r1 = -x1Flux[k][j][i].By - emf3_cc[k][j][i];
03752         else {
03753           de3_r1 = 0.5*(-x1Flux[k][j-1][i].By - emf3_cc[k][j-1][i]
03754                         -x1Flux[k][j  ][i].By - emf3_cc[k][j  ][i] );
03755         }
03756 
03757         emf3[k][j][i] = 0.25*(  x2Flux[k][j  ][i-1].Bz + x2Flux[k][j][i].Bz
03758                               - x1Flux[k][j-1][i  ].By - x1Flux[k][j][i].By
03759                               + de3_l1 + de3_r1 + de3_l2 + de3_r2);
03760       }
03761     }
03762   }
03763 
03764   return;
03765 }
03766 #endif /* MHD */
03767 
03768 #endif /* CTU_INTEGRATOR */
integrators/integrate_3d_ctu.c
