source: Athena/tags/version-3.1/src/integrate_3d_ctu.c @ 616

Revision 616, 80.1 KB checked in by jstone, 6 years ago (diff)

Fixed bug in V3.1; self-gravity not included in computation of cell-centered
emfs at n+1/2 in CTU+CT integrators.

Line 
1#include "copyright.h"
2/*==============================================================================
3 * FILE: integrate_3d_ctu.c
4 *
5 * PURPOSE: Updates the input Grid structure pointed to by *pG by one
6 *   timestep using directionally unsplit CTU method of Colella (1990).  The
7 *   variables updated are:
8 *      U.[d,M1,M2,M3,E,B1c,B2c,B3c,s] -- where U is of type Gas
9 *      B1i, B2i, B3i  -- interface magnetic field
10 *   Also adds gravitational source terms, self-gravity, and the H-correction
11 *   of Sanders et al.
12 *     For adb hydro, requires (9*Cons1D +  3*Real) = 48 3D arrays
13 *     For adb mhd, requires   (9*Cons1D + 10*Real) = 73 3D arrays
14 *   The H-correction of Sanders et al. adds another 3 arrays. 
15 *
16 * REFERENCES:
17 *   P. Colella, "Multidimensional upwind methods for hyperbolic conservation
18 *   laws", JCP, 87, 171 (1990)
19 *
20 *   T. Gardiner & J.M. Stone, "An unsplit Godunov method for ideal MHD via
21 *   constrained transport", JCP, 205, 509 (2005)
22 *
23 *   T. Gardiner & J.M. Stone, "An unsplit Godunov method for ideal MHD via
24 *   constrained transport in Three Dimensions", JCP, ***, *** (2007)
25 *
26 *   R. Sanders, E. Morano, & M.-C. Druguet, "Multidimensinal dissipation for
27 *   upwind schemes: stability and applications to gas dynamics", JCP, 145, 511
28 *   (1998)
29 *
30 * CONTAINS PUBLIC FUNCTIONS:
31 *   integrate_3d_ctu()
32 *   integrate_init_3d()
33 *   integrate_destruct_3d()
34 *============================================================================*/
35
36#include <math.h>
37#include <stdio.h>
38#include <stdlib.h>
39#include <string.h>
40#include "defs.h"
41#include "athena.h"
42#include "globals.h"
43#include "prototypes.h"
44
45/* The L/R states of conserved variables and fluxes at each cell face */
46static Cons1D ***Ul_x1Face=NULL, ***Ur_x1Face=NULL;
47static Cons1D ***Ul_x2Face=NULL, ***Ur_x2Face=NULL;
48static Cons1D ***Ul_x3Face=NULL, ***Ur_x3Face=NULL;
49static Cons1D ***x1Flux=NULL, ***x2Flux=NULL, ***x3Flux=NULL;
50
51/* The interface magnetic fields and emfs */
52static Real ***B1_x1Face=NULL, ***B2_x2Face=NULL, ***B3_x3Face=NULL;
53#ifdef MHD
54static Real ***emf1=NULL, ***emf2=NULL, ***emf3=NULL;
55static Real ***emf1_cc=NULL, ***emf2_cc=NULL, ***emf3_cc=NULL;
56#endif /* MHD */
57
58/* 1D scratch vectors used by lr_states and flux functions */
59static Real *Bxc=NULL, *Bxi=NULL;
60static Prim1D *W=NULL, *Wl=NULL, *Wr=NULL;
61static Cons1D *U1d=NULL;
62
63/* density at t^{n+1/2} needed by both MHD and to make gravity 2nd order */
64static Real ***dhalf = NULL;
65
66/* variables needed for H-correction of Sanders et al (1998) */
67extern Real etah;
68#ifdef H_CORRECTION
69static Real ***eta1=NULL, ***eta2=NULL, ***eta3=NULL;
70#endif
71
72/*==============================================================================
73 * PRIVATE FUNCTION PROTOTYPES:
74 *   integrate_emf1_corner() - the upwind CT method in GS05, for emf1
75 *   integrate_emf2_corner() - the upwind CT method in GS05, for emf2
76 *   integrate_emf3_corner() - the upwind CT method in GS05, for emf3
77 *============================================================================*/
78
79#ifdef MHD
80static void integrate_emf1_corner(const Grid *pG);
81static void integrate_emf2_corner(const Grid *pG);
82static void integrate_emf3_corner(const Grid *pG);
83#endif /* MHD */
84
85/*=========================== PUBLIC FUNCTIONS ===============================*/
86/*----------------------------------------------------------------------------*/
87/* integrate_3d: 3D CTU integrator for MHD using 6-solve method */
88
89void integrate_3d_ctu(Grid *pG)
90{
91  Real dtodx1=pG->dt/pG->dx1, dtodx2=pG->dt/pG->dx2, dtodx3=pG->dt/pG->dx3;
92  Real dx1i=1.0/pG->dx1, dx2i=1.0/pG->dx2, dx3i=1.0/pG->dx3;
93  Real q1 = 0.5*dtodx1, q2 = 0.5*dtodx2, q3 = 0.5*dtodx3;
94  int i, is = pG->is, ie = pG->ie;
95  int j, js = pG->js, je = pG->je;
96  int k, ks = pG->ks, ke = pG->ke;
97#ifdef MHD
98  Real MHD_src_By,MHD_src_Bz,mdb1,mdb2,mdb3;
99  Real db1,db2,db3,l1,l2,l3,B1,B2,B3,V1,V2,V3;
100  Real d, M1, M2, M3, B1c, B2c, B3c;
101  Real hdt = 0.5*pG->dt;
102#endif
103#ifdef H_CORRECTION
104  Real cfr,cfl,lambdar,lambdal;
105#endif
106#if (NSCALARS > 0)
107  int n;
108#endif
109  Real x1,x2,x3,phicl,phicr,phifc,phil,phir,phic;
110#ifdef SELF_GRAVITY
111  Real gxl,gxr,gyl,gyr,gzl,gzr,flx_m1l,flx_m1r,flx_m2l,flx_m2r,flx_m3l,flx_m3r;
112#endif
113
114/*--- Step 1a ------------------------------------------------------------------
115 * Load 1D vector of conserved variables;
116 * U1d = (d, M1, M2, M3, E, B2c, B3c, s[n])
117 */
118
119  for (k=ks-2; k<=ke+2; k++) {
120    for (j=js-2; j<=je+2; j++) {
121      for (i=is-nghost; i<=ie+nghost; i++) {
122        U1d[i].= pG->U[k][j][i].d;
123        U1d[i].Mx = pG->U[k][j][i].M1;
124        U1d[i].My = pG->U[k][j][i].M2;
125        U1d[i].Mz = pG->U[k][j][i].M3;
126#ifndef ISOTHERMAL
127        U1d[i].= pG->U[k][j][i].E;
128#endif /* ISOTHERMAL */
129#ifdef MHD
130        U1d[i].By = pG->U[k][j][i].B2c;
131        U1d[i].Bz = pG->U[k][j][i].B3c;
132        Bxc[i] = pG->U[k][j][i].B1c;
133        Bxi[i] = pG->B1i[k][j][i];
134        B1_x1Face[k][j][i] = pG->B1i[k][j][i];
135#endif /* MHD */
136#if (NSCALARS > 0)
137        for (n=0; n<NSCALARS; n++) U1d[i].s[n] = pG->U[k][j][i].s[n];
138#endif
139      }
140
141/*--- Step 1b ------------------------------------------------------------------
142 * Compute L and R states at X1-interfaces.
143 */
144
145     for (i=is-nghost; i<=ie+nghost; i++) {
146       Cons1D_to_Prim1D(&U1d[i],&W[i],&Bxc[i]);
147     }
148
149     lr_states(W,Bxc,pG->dt,dtodx1,is-1,ie+1,Wl,Wr);
150
151/* Add "MHD source terms" for 0.5*dt */
152
153#ifdef MHD
154      for (i=is-1; i<=ie+2; i++) {
155
156/* Source terms for left states in zone i-1 */
157        db1 = (pG->B1i[][][] - pG->B1i[k][j][i-1])*dx1i;
158        db2 = (pG->B2i[][j+1][i-1] - pG->B2i[k][j][i-1])*dx2i;
159        db3 = (pG->B3i[k+1][][i-1] - pG->B3i[k][j][i-1])*dx3i;
160
161        if(db1 >= 0.0){
162          l3 = db1 < -db3 ? db1 : -db3;
163          l3 = l3 > 0.0 ? l3 : 0.0;
164
165          l2 = db1 < -db2 ? db1 : -db2;
166          l2 = l2 > 0.0 ? l2 : 0.0;
167        }
168        else{
169          l3 = db1 > -db3 ? db1 : -db3;
170          l3 = l3 < 0.0 ? l3 : 0.0;
171
172          l2 = db1 > -db2 ? db1 : -db2;
173          l2 = l2 < 0.0 ? l2 : 0.0;
174        }
175
176        MHD_src_By = (pG->U[k][j][i-1].M2/pG->U[k][j][i-1].d)*l2;
177        MHD_src_Bz = (pG->U[k][j][i-1].M3/pG->U[k][j][i-1].d)*l3;
178
179        Wl[i].By += hdt*MHD_src_By;
180        Wl[i].Bz += hdt*MHD_src_Bz;
181
182/* Source terms for right states in zone i */
183        db1 = (pG->B1i[][][i+1] - pG->B1i[k][j][i])*dx1i;
184        db2 = (pG->B2i[][j+1][] - pG->B2i[k][j][i])*dx2i;
185        db3 = (pG->B3i[k+1][][] - pG->B3i[k][j][i])*dx3i;
186
187        if(db1 >= 0.0){
188          l3 = db1 < -db3 ? db1 : -db3;
189          l3 = l3 > 0.0 ? l3 : 0.0;
190
191          l2 = db1 < -db2 ? db1 : -db2;
192          l2 = l2 > 0.0 ? l2 : 0.0;
193        }
194        else{
195          l3 = db1 > -db3 ? db1 : -db3;
196          l3 = l3 < 0.0 ? l3 : 0.0;
197
198          l2 = db1 > -db2 ? db1 : -db2;
199          l2 = l2 < 0.0 ? l2 : 0.0;
200        }
201
202        MHD_src_By = (pG->U[k][j][i].M2/pG->U[k][j][i].d)*l2;
203        MHD_src_Bz = (pG->U[k][j][i].M3/pG->U[k][j][i].d)*l3;
204
205        Wr[i].By += hdt*MHD_src_By;
206        Wr[i].Bz += hdt*MHD_src_Bz;
207      }
208#endif
209
210/*--- Step 1c ------------------------------------------------------------------
211 * Add gravitational source terms from static potential for 0.5*dt to L/R states
212 */
213
214      if (StaticGravPot != NULL){
215        for (i=is-1; i<=ie+2; i++) {
216          cc_pos(pG,i,j,k,&x1,&x2,&x3);
217          phicr = (*StaticGravPot)( x1                ,x2,x3);
218          phicl = (*StaticGravPot)((x1-    pG->dx1),x2,x3);
219          phifc = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
220
221/* Apply gravitational source terms to velocity using gradient of potential
222 * for (dt/2).   S_{V} = -Grad(Phi) */
223          Wl[i].Vx -= dtodx1*(phifc - phicl);
224          Wr[i].Vx -= dtodx1*(phicr - phifc);
225        }
226      }
227
228/*--- Step 1d ------------------------------------------------------------------
229 * Add gravitational source terms for self-gravity for dt/2 to L/R states
230 */
231
232#ifdef SELF_GRAVITY
233      for (i=is-1; i<=ie+2; i++) {
234        Wl[i].Vx -= q1*(pG->Phi[k][j][i] - pG->Phi[k][j][i-1]);
235        Wr[i].Vx -= q1*(pG->Phi[k][j][i] - pG->Phi[k][j][i-1]);
236      }
237#endif
238
239
240/*--- Step 1e ------------------------------------------------------------------
241 * Compute 1D fluxes in x1-direction, storing into 3D array
242 */
243
244      for (i=is-1; i<=ie+2; i++) {
245        Prim1D_to_Cons1D(&Ul_x1Face[k][j][i],&Wl[i],&Bxi[i]);
246        Prim1D_to_Cons1D(&Ur_x1Face[k][j][i],&Wr[i],&Bxi[i]);
247
248        GET_FLUXES(Ul_x1Face[k][j][i],Ur_x1Face[k][j][i],Wl[i],Wr[i],
249                   B1_x1Face[k][j][i],&x1Flux[k][j][i]);
250      }
251    }
252  }
253
254/*--- Step 2a ------------------------------------------------------------------
255 * Load 1D vector of conserved variables;
256 * U1d = (d, M2, M3, M1, E, B3c, B1c, s[n])
257 */
258
259  for (k=ks-2; k<=ke+2; k++) {
260    for (i=is-2; i<=ie+2; i++) {
261      for (j=js-nghost; j<=je+nghost; j++) {
262        U1d[j].= pG->U[k][j][i].d;
263        U1d[j].Mx = pG->U[k][j][i].M2;
264        U1d[j].My = pG->U[k][j][i].M3;
265        U1d[j].Mz = pG->U[k][j][i].M1;
266#ifndef ISOTHERMAL
267        U1d[j].= pG->U[k][j][i].E;
268#endif /* ISOTHERMAL */
269#ifdef MHD
270        U1d[j].By = pG->U[k][j][i].B3c;
271        U1d[j].Bz = pG->U[k][j][i].B1c;
272        Bxc[j] = pG->U[k][j][i].B2c;
273        Bxi[j] = pG->B2i[k][j][i];
274        B2_x2Face[k][j][i] = pG->B2i[k][j][i];
275#endif /* MHD */
276#if (NSCALARS > 0)
277        for (n=0; n<NSCALARS; n++) U1d[j].s[n] = pG->U[k][j][i].s[n];
278#endif
279      }
280
281/*--- Step 2b ------------------------------------------------------------------
282 * Compute L and R states at X2-interfaces.
283 */
284
285      for (j=js-nghost; j<=je+nghost; j++) {
286        Cons1D_to_Prim1D(&U1d[j],&W[j],&Bxc[j]);
287      }
288
289      lr_states(W,Bxc,pG->dt,dtodx2,js-1,je+1,Wl,Wr);
290
291/* Add "MHD source terms" for 0.5*dt */
292
293#ifdef MHD
294      for (j=js-1; j<=je+2; j++) {
295/* Source terms for left states in zone j-1 */
296        db1 = (pG->B1i[][j-1][i+1] - pG->B1i[k][j-1][i])*dx1i;
297        db2 = (pG->B2i[][][] - pG->B2i[k][j-1][i])*dx2i;
298        db3 = (pG->B3i[k+1][j-1][] - pG->B3i[k][j-1][i])*dx3i;
299
300        if(db2 >= 0.0){
301          l1 = db2 < -db1 ? db2 : -db1;
302          l1 = l1 > 0.0 ? l1 : 0.0;
303
304          l3 = db2 < -db3 ? db2 : -db3;
305          l3 = l3 > 0.0 ? l3 : 0.0;
306        }
307        else{
308          l1 = db2 > -db1 ? db2 : -db1;
309          l1 = l1 < 0.0 ? l1 : 0.0;
310
311          l3 = db2 > -db3 ? db2 : -db3;
312          l3 = l3 < 0.0 ? l3 : 0.0;
313        }
314
315        MHD_src_By = (pG->U[k][j-1][i].M3/pG->U[k][j-1][i].d)*l3;
316        MHD_src_Bz = (pG->U[k][j-1][i].M1/pG->U[k][j-1][i].d)*l1;
317
318        Wl[j].By += hdt*MHD_src_By;
319        Wl[j].Bz += hdt*MHD_src_Bz;
320
321/* Source terms for right states in zone j */
322        db1 = (pG->B1i[][][i+1] - pG->B1i[k][j][i])*dx1i;
323        db2 = (pG->B2i[][j+1][] - pG->B2i[k][j][i])*dx2i;
324        db3 = (pG->B3i[k+1][][] - pG->B3i[k][j][i])*dx3i;
325
326        if(db2 >= 0.0){
327          l1 = db2 < -db1 ? db2 : -db1;
328          l1 = l1 > 0.0 ? l1 : 0.0;
329
330          l3 = db2 < -db3 ? db2 : -db3;
331          l3 = l3 > 0.0 ? l3 : 0.0;
332        }
333        else{
334          l1 = db2 > -db1 ? db2 : -db1;
335          l1 = l1 < 0.0 ? l1 : 0.0;
336
337          l3 = db2 > -db3 ? db2 : -db3;
338          l3 = l3 < 0.0 ? l3 : 0.0;
339        }
340
341        MHD_src_By = (pG->U[k][j][i].M3/pG->U[k][j][i].d)*l3;
342        MHD_src_Bz = (pG->U[k][j][i].M1/pG->U[k][j][i].d)*l1;
343
344        Wr[j].By += hdt*MHD_src_By;
345        Wr[j].Bz += hdt*MHD_src_Bz;
346      }
347#endif
348
349/*--- Step 2c ------------------------------------------------------------------
350 * Add gravitational source terms from static potential for 0.5*dt to L/R states
351 */
352
353      if (StaticGravPot != NULL){
354        for (j=js-1; j<=je+2; j++) {
355          cc_pos(pG,i,j,k,&x1,&x2,&x3);
356          phicr = (*StaticGravPot)(x1, x2                ,x3);
357          phicl = (*StaticGravPot)(x1,(x2-    pG->dx2),x3);
358          phifc = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
359
360/* Apply gravitational source terms to velocity using gradient of potential.
361 * for (dt/2).   S_{V} = -Grad(Phi) */
362          Wl[j].Vx -= dtodx2*(phifc - phicl);
363          Wr[j].Vx -= dtodx2*(phicr - phifc);
364        }
365      }
366
367/*--- Step 2d ------------------------------------------------------------------
368 * Add gravitational source terms for self-gravity for dt/2 to L/R states
369 */
370
371#ifdef SELF_GRAVITY
372      for (j=js-1; j<=je+2; j++) {
373        Wl[j].Vx -= q2*(pG->Phi[k][j][i] - pG->Phi[k][j-1][i]);
374        Wr[j].Vx -= q2*(pG->Phi[k][j][i] - pG->Phi[k][j-1][i]);
375      }
376#endif
377
378
379/*--- Step 2e ------------------------------------------------------------------
380 * Compute 1D fluxes in x2-direction, storing into 3D array
381 */
382
383      for (j=js-1; j<=je+2; j++) {
384        Prim1D_to_Cons1D(&Ul_x2Face[k][j][i],&Wl[j],&Bxi[j]);
385        Prim1D_to_Cons1D(&Ur_x2Face[k][j][i],&Wr[j],&Bxi[j]);
386
387        GET_FLUXES(Ul_x2Face[k][j][i],Ur_x2Face[k][j][i],Wl[j],Wr[j],
388                   B2_x2Face[k][j][i],&x2Flux[k][j][i]);
389      }
390    }
391  }
392
393
394/*--- Step 3a ------------------------------------------------------------------
395 * Load 1D vector of conserved variables;
396 * U1d = (d, M3, M1, M2, E, B1c, B2c, s[n])
397 */
398
399  for (j=js-2; j<=je+2; j++) {
400    for (i=is-2; i<=ie+2; i++) {
401      for (k=ks-nghost; k<=ke+nghost; k++) {
402        U1d[k].= pG->U[k][j][i].d;
403        U1d[k].Mx = pG->U[k][j][i].M3;
404        U1d[k].My = pG->U[k][j][i].M1;
405        U1d[k].Mz = pG->U[k][j][i].M2;
406#ifndef ISOTHERMAL
407        U1d[k].= pG->U[k][j][i].E;
408#endif /* ISOTHERMAL */
409#ifdef MHD
410        U1d[k].By = pG->U[k][j][i].B1c;
411        U1d[k].Bz = pG->U[k][j][i].B2c;
412        Bxc[k] = pG->U[k][j][i].B3c;
413        Bxi[k] = pG->B3i[k][j][i];
414        B3_x3Face[k][j][i] = pG->B3i[k][j][i];
415#endif /* MHD */
416#if (NSCALARS > 0)
417        for (n=0; n<NSCALARS; n++) U1d[k].s[n] = pG->U[k][j][i].s[n];
418#endif
419      }
420
421/*--- Step 3b ------------------------------------------------------------------
422 * Compute L and R states at X3-interfaces.
423 */
424
425      for (k=ks-nghost; k<=ke+nghost; k++) {
426        Cons1D_to_Prim1D(&U1d[k],&W[k],&Bxc[k]);
427      }
428
429      lr_states(W,Bxc,pG->dt,dtodx3,ks-1,ke+1,Wl,Wr);
430
431/* Add "MHD source terms" for 0.5*dt */
432
433#ifdef MHD
434      for (k=ks-1; k<=ke+2; k++) {
435/* Source terms for left states in zone k-1 */
436        db1 = (pG->B1i[k-1][][i+1] - pG->B1i[k-1][j][i])*dx1i;
437        db2 = (pG->B2i[k-1][j+1][] - pG->B2i[k-1][j][i])*dx2i;
438        db3 = (pG->B3i[][][] - pG->B3i[k-1][j][i])*dx3i;
439
440        if(db3 >= 0.0){
441          l1 = db3 < -db1 ? db3 : -db1;
442          l1 = l1 > 0.0 ? l1 : 0.0;
443
444          l2 = db3 < -db2 ? db3 : -db2;
445          l2 = l2 > 0.0 ? l2 : 0.0;
446        }
447        else{
448          l1 = db3 > -db1 ? db3 : -db1;
449          l1 = l1 < 0.0 ? l1 : 0.0;
450
451          l2 = db3 > -db2 ? db3 : -db2;
452          l2 = l2 < 0.0 ? l2 : 0.0;
453        }
454
455        MHD_src_By = (pG->U[k-1][j][i].M1/pG->U[k-1][j][i].d)*l1;
456        MHD_src_Bz = (pG->U[k-1][j][i].M2/pG->U[k-1][j][i].d)*l2;
457
458        Wl[k].By += hdt*MHD_src_By;
459        Wl[k].Bz += hdt*MHD_src_Bz;
460
461/* Source terms for right states in zone k */
462        db1 = (pG->B1i[k][j][i+1] - pG->B1i[k][j][i])*dx1i;
463        db2 = (pG->B2i[k][j+1][i] - pG->B2i[k][j][i])*dx2i;
464        db3 = (pG->B3i[k+1][j][i] - pG->B3i[k][j][i])*dx3i;
465
466        if(db3 >= 0.0){
467          l1 = db3 < -db1 ? db3 : -db1;
468          l1 = l1 > 0.0 ? l1 : 0.0;
469
470          l2 = db3 < -db2 ? db3 : -db2;
471          l2 = l2 > 0.0 ? l2 : 0.0;
472        }
473        else{
474          l1 = db3 > -db1 ? db3 : -db1;
475          l1 = l1 < 0.0 ? l1 : 0.0;
476
477          l2 = db3 > -db2 ? db3 : -db2;
478          l2 = l2 < 0.0 ? l2 : 0.0;
479        }
480
481        MHD_src_By = (pG->U[k][j][i].M1/pG->U[k][j][i].d)*l1;
482        MHD_src_Bz = (pG->U[k][j][i].M2/pG->U[k][j][i].d)*l2;
483
484        Wr[k].By += hdt*MHD_src_By;
485        Wr[k].Bz += hdt*MHD_src_Bz;
486      }
487#endif
488
489/*--- Step 3c ------------------------------------------------------------------
490 * Add gravitational source terms from static potential for 0.5*dt to L/R states
491 */
492
493      if (StaticGravPot != NULL){
494        for (k=ks-1; k<=ke+2; k++) {
495          cc_pos(pG,i,j,k,&x1,&x2,&x3);
496          phicr = (*StaticGravPot)(x1,x2, x3                );
497          phicl = (*StaticGravPot)(x1,x2,(x3-    pG->dx3));
498          phifc = (*StaticGravPot)(x1,x2,(x3-0.5*pG->dx3));
499
500/* Apply gravitational source terms to velocity using gradient of potential.
501 * for (dt/2).   S_{V} = -Grad(Phi) */
502          Wl[k].Vx -= dtodx3*(phifc - phicl);
503          Wr[k].Vx -= dtodx3*(phicr - phifc);
504        }
505      }
506
507/*--- Step 3d ------------------------------------------------------------------
508 * Add gravitational source terms for self-gravity for dt/2 to L/R states
509 */
510
511#ifdef SELF_GRAVITY
512      for (k=ks-1; k<=ke+2; k++) {
513        Wl[k].Vx -= q3*(pG->Phi[k][j][i] - pG->Phi[k-1][j][i]);
514        Wr[k].Vx -= q3*(pG->Phi[k][j][i] - pG->Phi[k-1][j][i]);
515      }
516#endif
517
518
519/*--- Step 3e ------------------------------------------------------------------
520 * Compute 1D fluxes in x3-direction, storing into 3D array
521 */
522
523      for (k=ks-1; k<=ke+2; k++) {
524        Prim1D_to_Cons1D(&Ul_x3Face[k][j][i],&Wl[k],&Bxi[k]);
525        Prim1D_to_Cons1D(&Ur_x3Face[k][j][i],&Wr[k],&Bxi[k]);
526
527        GET_FLUXES(Ul_x3Face[k][j][i],Ur_x3Face[k][j][i],Wl[k],Wr[k],
528                   B3_x3Face[k][j][i],&x3Flux[k][j][i]);
529      }
530    }
531  }
532
533/*--- Step 4 -------------------------------------------------------------------
534 * Calculate the cell centered value of emf1,2,3 at t^{n} and integrate
535 * to corner.
536 */
537
538#ifdef MHD
539/* emf1 */
540  for (k=ks-2; k<=ke+2; k++) {
541    for (j=js-2; j<=je+2; j++) {
542      for (i=is-2; i<=ie+2; i++) {
543        emf1_cc[k][j][i] = (pG->U[k][j][i].B2c*pG->U[k][j][i].M3 -
544                            pG->U[k][j][i].B3c*pG->U[k][j][i].M2)
545                              /pG->U[k][j][i].d;
546        emf2_cc[k][j][i] = (pG->U[k][j][i].B3c*pG->U[k][j][i].M1 -
547                            pG->U[k][j][i].B1c*pG->U[k][j][i].M3)
548                              /pG->U[k][j][i].d;
549        emf3_cc[k][j][i] = (pG->U[k][j][i].B1c*pG->U[k][j][i].M2 -
550                            pG->U[k][j][i].B2c*pG->U[k][j][i].M1)
551                              /pG->U[k][j][i].d;
552      }
553    }
554  }
555  integrate_emf1_corner(pG);
556  integrate_emf2_corner(pG);
557  integrate_emf3_corner(pG);
558#endif /* MHD */
559
560
561/*--- Step 5 ------------------------------------------------------------------
562 * Update the interface magnetic fields using CT for a half time step.
563 */
564
565#ifdef MHD
566  for (k=ks-1; k<=ke+1; k++) {
567    for (j=js-1; j<=je+1; j++) {
568      for (i=is-1; i<=ie+1; i++) {
569        B1_x1Face[k][j][i] += q3*(emf2[k+1][][] - emf2[k][j][i]) -
570                              q2*(emf3[][j+1][] - emf3[k][j][i]);
571        B2_x2Face[k][j][i] += q1*(emf3[][][i+1] - emf3[k][j][i]) -
572                              q3*(emf1[k+1][][] - emf1[k][j][i]);
573        B3_x3Face[k][j][i] += q2*(emf1[][j+1][] - emf1[k][j][i]) -
574                              q1*(emf2[][][i+1] - emf2[k][j][i]);
575      }
576      B1_x1Face[k][j][ie+2] += q3*(emf2[k+1][][ie+2]-emf2[k][j][ie+2]) -
577                               q2*(emf3[][j+1][ie+2]-emf3[k][j][ie+2]);
578    }
579    for (i=is-1; i<=ie+1; i++) {
580      B2_x2Face[k][je+2][i] += q1*(emf3[][je+2][i+1]-emf3[k][je+2][i]) -
581                               q3*(emf1[k+1][je+2][]-emf1[k][je+2][i]);
582    }
583  }
584  for (j=js-1; j<=je+1; j++) {
585    for (i=is-1; i<=ie+1; i++) {
586      B3_x3Face[ke+2][j][i] += q2*(emf1[ke+2][j+1][]-emf1[ke+2][j][i]) -
587                               q1*(emf2[ke+2][][i+1]-emf2[ke+2][j][i]);
588    }
589  }
590#endif
591
592/*--- Step 6a ------------------------------------------------------------------
593 * Correct the L/R states at x1-interfaces using x2-fluxes computed in Step 2e.
594 * Since the fluxes come from an x2-sweep, (x,y,z) on RHS -> (z,x,y) on LHS
595 */
596
597  for (k=ks-1; k<=ke+1; k++) {
598    for (j=js-1; j<=je+1; j++) {
599      for (i=is-1; i<=ie+2; i++) {
600        Ul_x1Face[k][j][i].d -=q2*(x2Flux[k][j+1][i-1].d -x2Flux[k][j][i-1].d );
601        Ul_x1Face[k][j][i].Mx-=q2*(x2Flux[k][j+1][i-1].Mz-x2Flux[k][j][i-1].Mz);
602        Ul_x1Face[k][j][i].My-=q2*(x2Flux[k][j+1][i-1].Mx-x2Flux[k][j][i-1].Mx);
603        Ul_x1Face[k][j][i].Mz-=q2*(x2Flux[k][j+1][i-1].My-x2Flux[k][j][i-1].My);
604#ifndef ISOTHERMAL
605        Ul_x1Face[k][j][i].E -=q2*(x2Flux[k][j+1][i-1].E -x2Flux[k][j][i-1].E );
606#endif /* ISOTHERMAL */
607#ifdef MHD
608/* Update B3 */
609        Ul_x1Face[k][j][i].Bz+=q2*0.5*
610          ((emf1[][j+1][i-1] - emf1[][j][i-1]) +
611           (emf1[k+1][j+1][i-1] - emf1[k+1][j][i-1]));
612#endif
613
614        Ur_x1Face[k][j][i].d -=q2*(x2Flux[k][j+1][].d -x2Flux[k][j][].d );
615        Ur_x1Face[k][j][i].Mx-=q2*(x2Flux[k][j+1][].Mz-x2Flux[k][j][].Mz);
616        Ur_x1Face[k][j][i].My-=q2*(x2Flux[k][j+1][].Mx-x2Flux[k][j][].Mx);
617        Ur_x1Face[k][j][i].Mz-=q2*(x2Flux[k][j+1][].My-x2Flux[k][j][].My);
618#ifndef ISOTHERMAL
619        Ur_x1Face[k][j][i].E -=q2*(x2Flux[k][j+1][].E -x2Flux[k][j][].E );
620#endif /* ISOTHERMAL */
621#ifdef MHD
622/* Update B3 */
623        Ur_x1Face[k][j][i].Bz+=q2*0.5*
624          ((emf1[][j+1][i] - emf1[][j][i]) +
625           (emf1[k+1][j+1][i] - emf1[k+1][j][i]));
626#endif
627#if (NSCALARS > 0)
628        for (n=0; n<NSCALARS; n++) {
629          Ul_x1Face[k][j][i].s[n] -=
630             q2*(x2Flux[k][j+1][i-1].s[n] - x2Flux[k][j][i-1].s[n]);
631          Ur_x1Face[k][j][i].s[n] -=
632             q2*(x2Flux[k][j+1][].s[n] - x2Flux[k][j][].s[n]);
633        }
634#endif
635
636
637
638/*--- Step 6b ------------------------------------------------------------------
639 * Correct the L/R states at x1-interfaces using x3-fluxes computed in Step 3e.
640 * Since the fluxes come from an x3-sweep, (x,y,z) on RHS -> (y,z,x) on LHS */
641
642        Ul_x1Face[k][j][i].d -=q3*(x3Flux[k+1][j][i-1].d -x3Flux[k][j][i-1].d );
643        Ul_x1Face[k][j][i].Mx-=q3*(x3Flux[k+1][j][i-1].My-x3Flux[k][j][i-1].My);
644        Ul_x1Face[k][j][i].My-=q3*(x3Flux[k+1][j][i-1].Mz-x3Flux[k][j][i-1].Mz);
645        Ul_x1Face[k][j][i].Mz-=q3*(x3Flux[k+1][j][i-1].Mx-x3Flux[k][j][i-1].Mx);
646#ifndef ISOTHERMAL
647        Ul_x1Face[k][j][i].E -=q3*(x3Flux[k+1][j][i-1].E -x3Flux[k][j][i-1].E );
648#endif /* ISOTHERMAL */
649#ifdef MHD
650/* Update B2 */
651        Ul_x1Face[k][j][i].By-=q3*0.5*
652          ((emf1[k+1][][i-1] - emf1[k][][i-1]) +
653           (emf1[k+1][j+1][i-1] - emf1[k][j+1][i-1]));
654#endif
655
656        Ur_x1Face[k][j][i].d -=q3*(x3Flux[k+1][j][].d -x3Flux[k][j][].d );
657        Ur_x1Face[k][j][i].Mx-=q3*(x3Flux[k+1][j][].My-x3Flux[k][j][].My);
658        Ur_x1Face[k][j][i].My-=q3*(x3Flux[k+1][j][].Mz-x3Flux[k][j][].Mz);
659        Ur_x1Face[k][j][i].Mz-=q3*(x3Flux[k+1][j][].Mx-x3Flux[k][j][].Mx);
660#ifndef ISOTHERMAL
661        Ur_x1Face[k][j][i].E -=q3*(x3Flux[k+1][j][].E -x3Flux[k][j][].E );
662#endif /* ISOTHERMAL */
663#ifdef MHD
664/* Update B2 */
665        Ur_x1Face[k][j][i].By-=q3*0.5*
666          ((emf1[k+1][][i] - emf1[k][][i]) +
667           (emf1[k+1][j+1][i] - emf1[k][j+1][i]));
668#endif
669#if (NSCALARS > 0)
670        for (n=0; n<NSCALARS; n++) {
671          Ul_x1Face[k][j][i].s[n] -=
672             q3*(x3Flux[k+1][j][i-1].s[n] - x3Flux[k][j][i-1].s[n]);
673          Ur_x1Face[k][j][i].s[n] -=
674             q3*(x3Flux[k+1][j][].s[n] - x3Flux[k][j][].s[n]);
675        }
676#endif
677      }
678    }
679  }
680
681/*--- Step 6c ------------------------------------------------------------------
682 * Add the "MHD source terms" from the x2- and x3-flux-gradients to the
683 * conservative variables on the x1Face.  Limiting is used as in GS (2007)
684 */
685#ifdef MHD
686  for (k=ks-1; k<=ke+1; k++) {
687    for (j=js-1; j<=je+1; j++) {
688      for (i=is-1; i<=ie+2; i++) {
689        db1 = (pG->B1i[][][] - pG->B1i[k][j][i-1])*dx1i;
690        db2 = (pG->B2i[][j+1][i-1] - pG->B2i[k][j][i-1])*dx2i;
691        db3 = (pG->B3i[k+1][][i-1] - pG->B3i[k][j][i-1])*dx3i;
692        B1 = pG->U[k][j][i-1].B1c;
693        B2 = pG->U[k][j][i-1].B2c;
694        B3 = pG->U[k][j][i-1].B3c;
695        V2 = pG->U[k][j][i-1].M2/pG->U[k][j][i-1].d;
696        V3 = pG->U[k][j][i-1].M3/pG->U[k][j][i-1].d;
697
698/* Calculate mdb2 = min_mod(-db1,db2) */
699        if(db1 > 0.0 && db2 < 0.0){
700          mdb2 = db2 > -db1 ? db2 : -db1;
701        }
702        else if(db1 < 0.0 && db2 > 0.0){
703          mdb2 = db2 < -db1 ? db2 : -db1;
704        }
705        else mdb2 = 0.0;
706
707/* Calculate mdb3 = min_mod(-db1,db3) */
708        if(db1 > 0.0 && db3 < 0.0){
709          mdb3 = db3 > -db1 ? db3 : -db1;
710        }
711        else if(db1 < 0.0 && db3 > 0.0){
712          mdb3 = db3 < -db1 ? db3 : -db1;
713        }
714        else mdb3 = 0.0;
715
716        Ul_x1Face[k][j][i].Mx += hdt*B1*db1;
717        Ul_x1Face[k][j][i].My += hdt*B2*db1;
718        Ul_x1Face[k][j][i].Mz += hdt*B3*db1;
719        Ul_x1Face[k][j][i].By += hdt*V2*(-mdb3);
720        Ul_x1Face[k][j][i].Bz += hdt*V3*(-mdb2);
721#ifndef ISOTHERMAL
722        Ul_x1Face[k][j][i].+= hdt*(B2*V2*(-mdb3) + B3*V3*(-mdb2) );
723#endif /* ISOTHERMAL */
724
725        db1 = (pG->B1i[][][i+1] - pG->B1i[k][j][i])*dx1i;
726        db2 = (pG->B2i[][j+1][] - pG->B2i[k][j][i])*dx2i;
727        db3 = (pG->B3i[k+1][][] - pG->B3i[k][j][i])*dx3i;
728        B1 = pG->U[k][j][i].B1c;
729        B2 = pG->U[k][j][i].B2c;
730        B3 = pG->U[k][j][i].B3c;
731        V2 = pG->U[k][j][i].M2/pG->U[k][j][i].d;
732        V3 = pG->U[k][j][i].M3/pG->U[k][j][i].d;
733
734/* Calculate mdb2 = min_mod(-db1,db2) */
735        if(db1 > 0.0 && db2 < 0.0){
736          mdb2 = db2 > -db1 ? db2 : -db1;
737        }
738        else if(db1 < 0.0 && db2 > 0.0){
739          mdb2 = db2 < -db1 ? db2 : -db1;
740        }
741        else mdb2 = 0.0;
742
743/* Calculate mdb3 = min_mod(-db1,db3) */
744        if(db1 > 0.0 && db3 < 0.0){
745          mdb3 = db3 > -db1 ? db3 : -db1;
746        }
747        else if(db1 < 0.0 && db3 > 0.0){
748          mdb3 = db3 < -db1 ? db3 : -db1;
749        }
750        else mdb3 = 0.0;
751
752        Ur_x1Face[k][j][i].Mx += hdt*B1*db1;
753        Ur_x1Face[k][j][i].My += hdt*B2*db1;
754        Ur_x1Face[k][j][i].Mz += hdt*B3*db1;
755        Ur_x1Face[k][j][i].By += hdt*V2*(-mdb3);
756        Ur_x1Face[k][j][i].Bz += hdt*V3*(-mdb2);
757#ifndef ISOTHERMAL
758        Ur_x1Face[k][j][i].+= hdt*(B2*V2*(-mdb3) + B3*V3*(-mdb2) );
759#endif /* ISOTHERMAL */
760      }
761    }
762  }
763#endif /* MHD */
764
765/*--- Step 6d ------------------------------------------------------------------
766 * Add gravitational source terms in x2-direction and x3-direction to
767 * corrected L/R states on x1-faces.  To improve conservation of total energy,
768 * we average the energy source term computed at cell faces.
769 *    S_{M} = -(\rho) Grad(Phi);   S_{E} = -(\rho v) Grad{Phi}
770 */
771
772  if (StaticGravPot != NULL){
773  for (k=ks-1; k<=ke+1; k++) {
774    for (j=js-1; j<=je+1; j++) {
775      for (i=is-1; i<=ie+2; i++) {
776        cc_pos(pG,i,j,k,&x1,&x2,&x3);
777        phic = (*StaticGravPot)(x1, x2                ,x3);
778        phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
779        phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
780
781/* correct right states */
782        Ur_x1Face[k][j][i].My -= q2*(phir-phil)*pG->U[k][j][i].d;
783#ifndef ISOTHERMAL
784        Ur_x1Face[k][j][i].E -= q2*(x2Flux[k][][].d*(phic - phil)
785                                  + x2Flux[k][j+1][].d*(phir - phic));
786#endif
787
788        phir = (*StaticGravPot)(x1,x2,(x3+0.5*pG->dx3));
789        phil = (*StaticGravPot)(x1,x2,(x3-0.5*pG->dx3));
790       
791        Ur_x1Face[k][j][i].Mz -= q3*(phir-phil)*pG->U[k][j][i].d;
792#ifndef ISOTHERMAL
793        Ur_x1Face[k][j][i].E -= q3*(x3Flux[][j][].d*(phic - phil)
794                                  + x3Flux[k+1][j][].d*(phir - phic));
795#endif
796
797/* correct left states */
798        phic = (*StaticGravPot)((x1-pG->dx1), x2                ,x3);
799        phir = (*StaticGravPot)((x1-pG->dx1),(x2+0.5*pG->dx2),x3);
800        phil = (*StaticGravPot)((x1-pG->dx1),(x2-0.5*pG->dx2),x3);
801
802        Ul_x1Face[k][j][i].My -= q2*(phir-phil)*pG->U[k][j][i-1].d;
803#ifndef ISOTHERMAL
804        Ul_x1Face[k][j][i].E -= q2*(x2Flux[k][][i-1].d*(phic - phil)
805                                  + x2Flux[k][j+1][i-1].d*(phir - phic));
806#endif
807
808        phir = (*StaticGravPot)((x1-pG->dx1),x2,(x3+0.5*pG->dx3));
809        phil = (*StaticGravPot)((x1-pG->dx1),x2,(x3-0.5*pG->dx3));
810       
811        Ul_x1Face[k][j][i].Mz -= q3*(phir-phil)*pG->U[k][j][i-1].d;
812#ifndef ISOTHERMAL
813        Ul_x1Face[k][j][i].E -= q3*(x3Flux[][j][i-1].d*(phic - phil)
814                                  + x3Flux[k+1][j][i-1].d*(phir - phic));
815#endif
816      }
817    }
818  }}
819
820/*--- Step 6e ------------------------------------------------------------------
821 * Add source terms for self gravity to L/R states.
822 *    S_{M} = -(\rho) Grad(Phi);   S_{E} = -(\rho v) Grad{Phi}
823 */
824
825#ifdef SELF_GRAVITY
826  for (k=ks-1; k<=ke+1; k++) {
827    for (j=js-1; j<=je+1; j++) {
828      for (i=is-1; i<=ie+2; i++) {
829        phic = pG->Phi[k][j][i];
830        phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j+1][i]);
831        phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j-1][i]);
832
833/* correct right states */
834        Ur_x1Face[k][j][i].My -= q2*(phir-phil)*pG->U[k][j][i].d;
835#ifndef ISOTHERMAL
836        Ur_x1Face[k][j][i].E -= q2*(x2Flux[k][][].d*(phic - phil)
837                                  + x2Flux[k][j+1][].d*(phir - phic));
838#endif
839
840        phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k+1][j][i]);
841        phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k-1][j][i]);
842
843        Ur_x1Face[k][j][i].Mz -= q3*(phir-phil)*pG->U[k][j][i].d;
844#ifndef ISOTHERMAL
845        Ur_x1Face[k][j][i].E -= q3*(x3Flux[][j][].d*(phic - phil)
846                                  + x3Flux[k+1][j][].d*(phir - phic));
847#endif
848
849/* correct left states */
850        phic = pG->Phi[k][j][i-1];
851        phir = 0.5*(pG->Phi[k][j][i-1] + pG->Phi[k][j+1][i-1]);
852        phil = 0.5*(pG->Phi[k][j][i-1] + pG->Phi[k][j-1][i-1]);
853
854        Ul_x1Face[k][j][i].My -= q2*(phir-phil)*pG->U[k][j][i-1].d;
855#ifndef ISOTHERMAL
856        Ul_x1Face[k][j][i].E -= q2*(x2Flux[k][][i-1].d*(phic - phil)
857                                  + x2Flux[k][j+1][i-1].d*(phir - phic));
858#endif
859
860        phir = 0.5*(pG->Phi[k][j][i-1] + pG->Phi[k+1][j][i-1]);
861        phil = 0.5*(pG->Phi[k][j][i-1] + pG->Phi[k-1][j][i-1]);
862
863        Ul_x1Face[k][j][i].Mz -= q3*(phir-phil)*pG->U[k][j][i-1].d;
864#ifndef ISOTHERMAL
865        Ul_x1Face[k][j][i].E -= q3*(x3Flux[][j][i-1].d*(phic - phil)
866                                  + x3Flux[k+1][j][i-1].d*(phir - phic));
867#endif
868      }
869    }
870  }
871#endif /* SELF_GRAVITY */
872
873/*--- Step 7a ------------------------------------------------------------------
874 * Correct the L/R states at x2-interfaces using x1-fluxes computed in Step 1e.
875 * Since the fluxes come from an x1-sweep, (x,y,z) on RHS -> (y,z,x) on LHS
876 */
877
878  for (k=ks-1; k<=ke+1; k++) {
879    for (j=js-1; j<=je+2; j++) {
880      for (i=is-1; i<=ie+1; i++) {
881        Ul_x2Face[k][j][i].d -=q1*(x1Flux[k][j-1][i+1].d -x1Flux[k][j-1][i].d );
882        Ul_x2Face[k][j][i].Mx-=q1*(x1Flux[k][j-1][i+1].My-x1Flux[k][j-1][i].My);
883        Ul_x2Face[k][j][i].My-=q1*(x1Flux[k][j-1][i+1].Mz-x1Flux[k][j-1][i].Mz);
884        Ul_x2Face[k][j][i].Mz-=q1*(x1Flux[k][j-1][i+1].Mx-x1Flux[k][j-1][i].Mx);
885#ifndef ISOTHERMAL
886        Ul_x2Face[k][j][i].E -=q1*(x1Flux[k][j-1][i+1].E -x1Flux[k][j-1][i].E );
887#endif /* ISOTHERMAL */
888#ifdef MHD
889/* Update B3 */
890        Ul_x2Face[k][j][i].By-=q1*0.5*
891          ((emf2[][j-1][i+1] - emf2[][j-1][i]) + 
892           (emf2[k+1][j-1][i+1] - emf2[k+1][j-1][i]));
893#endif
894
895        Ur_x2Face[k][j][i].d -=q1*(x1Flux[k][][i+1].d -x1Flux[k][][i].d );
896        Ur_x2Face[k][j][i].Mx-=q1*(x1Flux[k][][i+1].My-x1Flux[k][][i].My);
897        Ur_x2Face[k][j][i].My-=q1*(x1Flux[k][][i+1].Mz-x1Flux[k][][i].Mz);
898        Ur_x2Face[k][j][i].Mz-=q1*(x1Flux[k][][i+1].Mx-x1Flux[k][][i].Mx);
899#ifndef ISOTHERMAL
900        Ur_x2Face[k][j][i].E -=q1*(x1Flux[k][][i+1].E -x1Flux[k][][i].E );
901#endif /* ISOTHERMAL */
902#ifdef MHD
903/* Update B3 */
904        Ur_x2Face[k][j][i].By-=q1*0.5*
905          ((emf2[][j][i+1] - emf2[][j][i]) + 
906           (emf2[k+1][j][i+1] - emf2[k+1][j][i]));
907#endif
908#if (NSCALARS > 0)
909        for (n=0; n<NSCALARS; n++) {
910          Ul_x2Face[k][j][i].s[n] -=
911             q1*(x1Flux[k][j-1][i+1].s[n] - x1Flux[k][j-1][i].s[n]);
912          Ur_x2Face[k][j][i].s[n] -=
913             q1*(x1Flux[k][][i+1].s[n] - x1Flux[k][][i].s[n]);
914        }
915#endif
916
917/*--- Step 7b ------------------------------------------------------------------
918 * Correct the L/R states at x2-interfaces using x3-fluxes computed in Step 3e.
919 * Since the fluxes come from an x3-sweep, (x,y,z) on RHS -> (z,x,y) on LHS
920 */
921
922        Ul_x2Face[k][j][i].d -=q3*(x3Flux[k+1][j-1][i].d -x3Flux[k][j-1][i].d );
923        Ul_x2Face[k][j][i].Mx-=q3*(x3Flux[k+1][j-1][i].Mz-x3Flux[k][j-1][i].Mz);
924        Ul_x2Face[k][j][i].My-=q3*(x3Flux[k+1][j-1][i].Mx-x3Flux[k][j-1][i].Mx);
925        Ul_x2Face[k][j][i].Mz-=q3*(x3Flux[k+1][j-1][i].My-x3Flux[k][j-1][i].My);
926#ifndef ISOTHERMAL
927        Ul_x2Face[k][j][i].E -=q3*(x3Flux[k+1][j-1][i].E -x3Flux[k][j-1][i].E );
928#endif /* ISOTHERMAL */
929#ifdef MHD
930/* Update B1 */
931        Ul_x2Face[k][j][i].Bz+=q3*0.5*
932          ((emf2[k+1][j-1][] - emf2[k][j-1][]) +
933           (emf2[k+1][j-1][i+1] - emf2[k][j-1][i+1]));
934#endif
935
936        Ur_x2Face[k][j][i].d -=q3*(x3Flux[k+1][][i].d -x3Flux[k][][i].d );
937        Ur_x2Face[k][j][i].Mx-=q3*(x3Flux[k+1][][i].Mz-x3Flux[k][][i].Mz);
938        Ur_x2Face[k][j][i].My-=q3*(x3Flux[k+1][][i].Mx-x3Flux[k][][i].Mx);
939        Ur_x2Face[k][j][i].Mz-=q3*(x3Flux[k+1][][i].My-x3Flux[k][][i].My);
940#ifndef ISOTHERMAL
941        Ur_x2Face[k][j][i].E -=q3*(x3Flux[k+1][][i].E -x3Flux[k][][i].E );
942#endif /* ISOTHERMAL */
943#ifdef MHD
944/* Update B1 */
945        Ur_x2Face[k][j][i].Bz+=q3*0.5*
946          ((emf2[k+1][j][] - emf2[k][j][]) +
947           (emf2[k+1][j][i+1] - emf2[k][j][i+1]));
948#endif
949#if (NSCALARS > 0)
950        for (n=0; n<NSCALARS; n++) {
951          Ul_x2Face[k][j][i].s[n] -=
952             q3*(x3Flux[k+1][j-1][i].s[n] - x3Flux[k][j-1][i].s[n]);
953          Ur_x2Face[k][j][i].s[n] -=
954             q3*(x3Flux[k+1][][i].s[n] - x3Flux[k][][i].s[n]);
955        }
956#endif
957      }
958    }
959  }
960
961/*--- Step 7c ------------------------------------------------------------------
962 * Add the "MHD source terms" from the x1- and x3-flux-gradients to the
963 * conservative variables on the x2Face.  Limiting is used as in GS (2007)
964 */
965#ifdef MHD
966  for (k=ks-1; k<=ke+1; k++) {
967    for (j=js-1; j<=je+2; j++) {
968      for (i=is-1; i<=ie+1; i++) {
969        db1 = (pG->B1i[][j-1][i+1] - pG->B1i[k][j-1][i])*dx1i;
970        db2 = (pG->B2i[][][] - pG->B2i[k][j-1][i])*dx2i;
971        db3 = (pG->B3i[k+1][j-1][] - pG->B3i[k][j-1][i])*dx3i;
972        B1 = pG->U[k][j-1][i].B1c;
973        B2 = pG->U[k][j-1][i].B2c;
974        B3 = pG->U[k][j-1][i].B3c;
975        V1 = pG->U[k][j-1][i].M1/pG->U[k][j-1][i].d;
976        V3 = pG->U[k][j-1][i].M3/pG->U[k][j-1][i].d;
977
978/* Calculate mdb1 = min_mod(-db2,db1) */
979        if(db2 > 0.0 && db1 < 0.0){
980          mdb1 = db1 > -db2 ? db1 : -db2;
981        }
982        else if(db2 < 0.0 && db1 > 0.0){
983          mdb1 = db1 < -db2 ? db1 : -db2;
984        }
985        else mdb1 = 0.0;
986
987/* Calculate mdb3 = min_mod(-db2,db3) */
988        if(db2 > 0.0 && db3 < 0.0){
989          mdb3 = db3 > -db2 ? db3 : -db2;
990        }
991        else if(db2 < 0.0 && db3 > 0.0){
992          mdb3 = db3 < -db2 ? db3 : -db2;
993        }
994        else mdb3 = 0.0;
995
996        Ul_x2Face[k][j][i].Mz += hdt*B1*db2;
997        Ul_x2Face[k][j][i].Mx += hdt*B2*db2;
998        Ul_x2Face[k][j][i].My += hdt*B3*db2;
999        Ul_x2Face[k][j][i].By += hdt*V3*(-mdb1);
1000        Ul_x2Face[k][j][i].Bz += hdt*V1*(-mdb3);
1001#ifndef ISOTHERMAL
1002        Ul_x2Face[k][j][i].+= hdt*(B3*V3*(-mdb1) + B1*V1*(-mdb3) );
1003#endif /* ISOTHERMAL */
1004
1005        db1 = (pG->B1i[][][i+1] - pG->B1i[k][j][i])*dx1i;
1006        db2 = (pG->B2i[][j+1][] - pG->B2i[k][j][i])*dx2i;
1007        db3 = (pG->B3i[k+1][][] - pG->B3i[k][j][i])*dx3i;
1008        B1 = pG->U[k][j][i].B1c;
1009        B2 = pG->U[k][j][i].B2c;
1010        B3 = pG->U[k][j][i].B3c;
1011        V1 = pG->U[k][j][i].M1/pG->U[k][j][i].d;
1012        V3 = pG->U[k][j][i].M3/pG->U[k][j][i].d;
1013
1014/* Calculate mdb1 = min_mod(-db2,db1) */
1015        if(db2 > 0.0 && db1 < 0.0){
1016          mdb1 = db1 > -db2 ? db1 : -db2;
1017        }
1018        else if(db2 < 0.0 && db1 > 0.0){
1019          mdb1 = db1 < -db2 ? db1 : -db2;
1020        }
1021        else mdb1 = 0.0;
1022
1023/* Calculate mdb3 = min_mod(-db2,db3) */
1024        if(db2 > 0.0 && db3 < 0.0){
1025          mdb3 = db3 > -db2 ? db3 : -db2;
1026        }
1027        else if(db2 < 0.0 && db3 > 0.0){
1028          mdb3 = db3 < -db2 ? db3 : -db2;
1029        }
1030        else mdb3 = 0.0;
1031
1032        Ur_x2Face[k][j][i].Mz += hdt*B1*db2;
1033        Ur_x2Face[k][j][i].Mx += hdt*B2*db2;
1034        Ur_x2Face[k][j][i].My += hdt*B3*db2;
1035        Ur_x2Face[k][j][i].By += hdt*V3*(-mdb1);
1036        Ur_x2Face[k][j][i].Bz += hdt*V1*(-mdb3);
1037#ifndef ISOTHERMAL
1038        Ur_x2Face[k][j][i].+= hdt*(B3*V3*(-mdb1) + B1*V1*(-mdb3) );
1039#endif /* ISOTHERMAL */
1040      }
1041    }
1042  }
1043#endif /* MHD */
1044
1045/*--- Step 7d ------------------------------------------------------------------
1046 * Add gravitational source terms in x1-direction and x3-direction to
1047 * corrected L/R states on x2-faces. To improve conservation of total energy,
1048 * we average the energy source term computed at cell faces.
1049 *    S_{M} = -(\rho) Grad(Phi);   S_{E} = -(\rho v) Grad{Phi}
1050 */
1051
1052  if (StaticGravPot != NULL){
1053  for (k=ks-1; k<=ke+1; k++) {
1054    for (j=js-1; j<=je+2; j++) {
1055      for (i=is-1; i<=ie+1; i++) {
1056        cc_pos(pG,i,j,k,&x1,&x2,&x3);
1057        phic = (*StaticGravPot)((x1               ),x2,x3);
1058        phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
1059        phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
1060
1061/* correct right states */
1062        Ur_x2Face[k][j][i].Mz -= q1*(phir-phil)*pG->U[k][j][i].d;
1063#ifndef ISOTHERMAL
1064        Ur_x2Face[k][j][i].E -= q1*(x1Flux[k][][].d*(phic - phil)
1065                                  + x1Flux[k][][i+1].d*(phir - phic));
1066#endif
1067
1068        phir = (*StaticGravPot)(x1,x2,(x3+0.5*pG->dx3));
1069        phil = (*StaticGravPot)(x1,x2,(x3-0.5*pG->dx3));
1070
1071        Ur_x2Face[k][j][i].My -= q3*(phir-phil)*pG->U[k][j][i].d;
1072#ifndef ISOTHERMAL
1073        Ur_x2Face[k][j][i].E -= q3*(x3Flux[][][i].d*(phic - phil)
1074                                  + x3Flux[k+1][][i].d*(phir - phic));
1075#endif
1076
1077/* correct left states */
1078        phic = (*StaticGravPot)((x1               ),(x2-pG->dx2),x3);
1079        phir = (*StaticGravPot)((x1+0.5*pG->dx1),(x2-pG->dx2),x3);
1080        phil = (*StaticGravPot)((x1-0.5*pG->dx1),(x2-pG->dx2),x3);
1081
1082        Ul_x2Face[k][j][i].Mz -= q1*(phir-phil)*pG->U[k][j-1][i].d;
1083#ifndef ISOTHERMAL
1084        Ul_x2Face[k][j][i].E -= q1*(x1Flux[k][j-1][].d*(phic - phil)
1085                                  + x1Flux[k][j-1][i+1].d*(phir - phic));
1086#endif
1087        phir = (*StaticGravPot)(x1,(x2-pG->dx2),(x3+0.5*pG->dx3));
1088        phil = (*StaticGravPot)(x1,(x2-pG->dx2),(x3-0.5*pG->dx3));
1089
1090        Ul_x2Face[k][j][i].My -= q3*(phir-phil)*pG->U[k][j-1][i].d;
1091#ifndef ISOTHERMAL
1092        Ul_x2Face[k][j][i].E -= q3*(x3Flux[][j-1][i].d*(phic - phil)
1093                                  + x3Flux[k+1][j-1][i].d*(phir - phic));
1094#endif
1095      }
1096    }
1097  }}
1098
1099/*--- Step 7e ------------------------------------------------------------------
1100 * Add source terms for self gravity to L/R states.
1101 *    S_{M} = -(\rho) Grad(Phi);   S_{E} = -(\rho v) Grad{Phi}
1102 */
1103
1104#ifdef SELF_GRAVITY
1105  for (k=ks-1; k<=ke+1; k++) {
1106    for (j=js-1; j<=je+2; j++) {
1107      for (i=is-1; i<=ie+1; i++) {
1108        phic = pG->Phi[k][j][i];
1109        phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j][i+1]);
1110        phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j][i-1]);
1111
1112/* correct right states */
1113        Ur_x2Face[k][j][i].Mz -= q1*(phir-phil)*pG->U[k][j][i].d;
1114#ifndef ISOTHERMAL
1115        Ur_x2Face[k][j][i].E -= q1*(x1Flux[k][j][].d*(phic - phil)
1116                                  + x1Flux[k][j][i+1].d*(phir - phic));
1117#endif
1118
1119        phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k+1][j][i]);
1120        phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k-1][j][i]);
1121
1122        Ur_x2Face[k][j][i].My -= q3*(phir-phil)*pG->U[k][j][i].d;
1123#ifndef ISOTHERMAL
1124        Ur_x2Face[k][j][i].E -= q3*(x3Flux[][j][i].d*(phic - phil)
1125                                  + x3Flux[k+1][j][i].d*(phir - phic));
1126#endif
1127
1128/* correct left states */
1129        phic = pG->Phi[k][j-1][i];
1130        phir = 0.5*(pG->Phi[k][j-1][i] + pG->Phi[k][j-1][i+1]);
1131        phil = 0.5*(pG->Phi[k][j-1][i] + pG->Phi[k][j-1][i-1]);
1132
1133        Ul_x2Face[k][j][i].Mz -= q1*(phir-phil)*pG->U[k][j-1][i].d;
1134#ifndef ISOTHERMAL
1135        Ul_x2Face[k][j][i].E -= q1*(x1Flux[k][j-1][].d*(phic - phil)
1136                                  + x1Flux[k][j-1][i+1].d*(phir - phic));
1137#endif
1138        phir = 0.5*(pG->Phi[k][j-1][i] + pG->Phi[k+1][j-1][i]);
1139        phil = 0.5*(pG->Phi[k][j-1][i] + pG->Phi[k-1][j-1][i]);
1140
1141        Ul_x2Face[k][j][i].My -= q3*(phir-phil)*pG->U[k][j-1][i].d;
1142#ifndef ISOTHERMAL
1143        Ul_x2Face[k][j][i].E -= q3*(x3Flux[][j-1][i].d*(phic - phil)
1144                                  + x3Flux[k+1][j-1][i].d*(phir - phic));
1145#endif
1146      }
1147    }
1148  }
1149#endif /* SELF_GRAVITY */
1150
1151
1152/*--- Step 8a ------------------------------------------------------------------
1153 * Correct the L/R states at x3-interfaces using x1-fluxes computed in Step 1e.
1154 * Since the fluxes come from an x1-sweep, (x,y,z) on RHS -> (z,x,y) on LHS
1155 */
1156
1157  for (k=ks-1; k<=ke+2; k++) {
1158    for (j=js-1; j<=je+1; j++) {
1159      for (i=is-1; i<=ie+1; i++) {
1160        Ul_x3Face[k][j][i].d -=q1*(x1Flux[k-1][j][i+1].d -x1Flux[k-1][j][i].d );
1161        Ul_x3Face[k][j][i].Mx-=q1*(x1Flux[k-1][j][i+1].Mz-x1Flux[k-1][j][i].Mz);
1162        Ul_x3Face[k][j][i].My-=q1*(x1Flux[k-1][j][i+1].Mx-x1Flux[k-1][j][i].Mx);
1163        Ul_x3Face[k][j][i].Mz-=q1*(x1Flux[k-1][j][i+1].My-x1Flux[k-1][j][i].My);
1164#ifndef ISOTHERMAL
1165        Ul_x3Face[k][j][i].E -=q1*(x1Flux[k-1][j][i+1].E -x1Flux[k-1][j][i].E );
1166#endif /* ISOTHERMAL */
1167#ifdef MHD
1168/* Update B2 */
1169        Ul_x3Face[k][j][i].Bz+=q1*0.5*
1170          ((emf3[k-1][][i+1] - emf3[k-1][][i]) +
1171           (emf3[k-1][j+1][i+1] - emf3[k-1][j+1][i]));
1172#endif
1173
1174        Ur_x3Face[k][j][i].d -=q1*(x1Flux[][j][i+1].d -x1Flux[][j][i].d );
1175        Ur_x3Face[k][j][i].Mx-=q1*(x1Flux[][j][i+1].Mz-x1Flux[][j][i].Mz);
1176        Ur_x3Face[k][j][i].My-=q1*(x1Flux[][j][i+1].Mx-x1Flux[][j][i].Mx);
1177        Ur_x3Face[k][j][i].Mz-=q1*(x1Flux[][j][i+1].My-x1Flux[][j][i].My);
1178#ifndef ISOTHERMAL
1179        Ur_x3Face[k][j][i].E -=q1*(x1Flux[][j][i+1].E -x1Flux[][j][i].E );
1180#endif /* ISOTHERMAL */
1181#ifdef MHD
1182/* Update B2 */
1183        Ur_x3Face[k][j][i].Bz+=q1*0.5*
1184          ((emf3[k][][i+1] - emf3[k][][i]) +
1185           (emf3[k][j+1][i+1] - emf3[k][j+1][i]));
1186#endif
1187#if (NSCALARS > 0)
1188        for (n=0; n<NSCALARS; n++) {
1189          Ul_x3Face[k][j][i].s[n] -=
1190             q1*(x1Flux[k-1][j][i+1].s[n] - x1Flux[k-1][j][i].s[n]);
1191          Ur_x3Face[k][j][i].s[n] -=
1192             q1*(x1Flux[][j][i+1].s[n] - x1Flux[][j][i].s[n]);
1193        }
1194#endif
1195
1196/*--- Step 8b ------------------------------------------------------------------
1197 * Correct the L/R states at x3-interfaces using x2-fluxes computed in Step 2e.
1198 * Since the fluxes come from an x2-sweep, (x,y,z) on RHS -> (y,z,x) on LHS
1199 */
1200
1201        Ul_x3Face[k][j][i].d -=q2*(x2Flux[k-1][j+1][i].d -x2Flux[k-1][j][i].d );
1202        Ul_x3Face[k][j][i].Mx-=q2*(x2Flux[k-1][j+1][i].My-x2Flux[k-1][j][i].My);
1203        Ul_x3Face[k][j][i].My-=q2*(x2Flux[k-1][j+1][i].Mz-x2Flux[k-1][j][i].Mz);
1204        Ul_x3Face[k][j][i].Mz-=q2*(x2Flux[k-1][j+1][i].Mx-x2Flux[k-1][j][i].Mx);
1205#ifndef ISOTHERMAL
1206        Ul_x3Face[k][j][i].E -=q2*(x2Flux[k-1][j+1][i].E -x2Flux[k-1][j][i].E );
1207#endif /* ISOTHERMAL */
1208#ifdef MHD
1209/* Update B1 */
1210        Ul_x3Face[k][j][i].By-=q2*0.5*
1211          ((emf3[k-1][j+1][] - emf3[k-1][j][]) +
1212           (emf3[k-1][j+1][i+1] - emf3[k-1][j][i+1]));
1213#endif
1214
1215        Ur_x3Face[k][j][i].d -=q2*(x2Flux[][j+1][i].d -x2Flux[][j][i].d );
1216        Ur_x3Face[k][j][i].Mx-=q2*(x2Flux[][j+1][i].My-x2Flux[][j][i].My);
1217        Ur_x3Face[k][j][i].My-=q2*(x2Flux[][j+1][i].Mz-x2Flux[][j][i].Mz);
1218        Ur_x3Face[k][j][i].Mz-=q2*(x2Flux[][j+1][i].Mx-x2Flux[][j][i].Mx);
1219#ifndef ISOTHERMAL
1220        Ur_x3Face[k][j][i].E -=q2*(x2Flux[][j+1][i].E -x2Flux[][j][i].E );
1221#endif /* ISOTHERMAL */
1222#ifdef MHD
1223/* Update B1 */
1224        Ur_x3Face[k][j][i].By-=q2*0.5*
1225          ((emf3[k][j+1][] - emf3[k][j][]) +
1226           (emf3[k][j+1][i+1] - emf3[k][j][i+1]));
1227#endif
1228#if (NSCALARS > 0)
1229        for (n=0; n<NSCALARS; n++) {
1230          Ul_x3Face[k][j][i].s[n] -=
1231             q2*(x2Flux[k-1][j+1][i].s[n] - x2Flux[k-1][j][i].s[n]);
1232          Ur_x3Face[k][j][i].s[n] -=
1233             q2*(x2Flux[][j+1][i].s[n] - x2Flux[][j][i].s[n]);
1234        }
1235#endif
1236      }
1237    }
1238  }
1239
1240/*--- Step 8c ------------------------------------------------------------------
1241 * Add the "MHD source terms" from the x1- and x2-flux-gradients to the
1242 * conservative variables on the x3Face.  Limiting is used as in GS07.
1243 */
1244#ifdef MHD
1245  for (k=ks-1; k<=ke+2; k++) {
1246    for (j=js-1; j<=je+1; j++) {
1247      for (i=is-1; i<=ie+1; i++) {
1248        db1 = (pG->B1i[k-1][][i+1] - pG->B1i[k-1][j][i])*dx1i;
1249        db2 = (pG->B2i[k-1][j+1][] - pG->B2i[k-1][j][i])*dx2i;
1250        db3 = (pG->B3i[][][] - pG->B3i[k-1][j][i])*dx3i;
1251        B1 = pG->U[k-1][j][i].B1c;
1252        B2 = pG->U[k-1][j][i].B2c;
1253        B3 = pG->U[k-1][j][i].B3c;
1254        V1 = pG->U[k-1][j][i].M1/pG->U[k-1][j][i].d;
1255        V2 = pG->U[k-1][j][i].M2/pG->U[k-1][j][i].d;
1256
1257/* Calculate mdb1 = min_mod(-db3,db1) */
1258        if(db3 > 0.0 && db1 < 0.0){
1259          mdb1 = db1 > -db3 ? db1 : -db3;
1260        }
1261        else if(db3 < 0.0 && db1 > 0.0){
1262          mdb1 = db1 < -db3 ? db1 : -db3;
1263        }
1264        else mdb1 = 0.0;
1265
1266/* Calculate mdb2 = min_mod(-db3,db2) */
1267        if(db3 > 0.0 && db2 < 0.0){
1268          mdb2 = db2 > -db3 ? db2 : -db3;
1269        }
1270        else if(db3 < 0.0 && db2 > 0.0){
1271          mdb2 = db2 < -db3 ? db2 : -db3;
1272        }
1273        else mdb2 = 0.0;
1274
1275        Ul_x3Face[k][j][i].My += hdt*B1*db3;
1276        Ul_x3Face[k][j][i].Mz += hdt*B2*db3;
1277        Ul_x3Face[k][j][i].Mx += hdt*B3*db3;
1278        Ul_x3Face[k][j][i].By += hdt*V1*(-mdb2);
1279        Ul_x3Face[k][j][i].Bz += hdt*V2*(-mdb1);
1280#ifndef ISOTHERMAL
1281        Ul_x3Face[k][j][i].+= hdt*(B1*V1*(-mdb2) + B2*V2*(-mdb1) );
1282#endif /* ISOTHERMAL */
1283
1284        db1 = (pG->B1i[][][i+1] - pG->B1i[k][j][i])*dx1i;
1285        db2 = (pG->B2i[][j+1][] - pG->B2i[k][j][i])*dx2i;
1286        db3 = (pG->B3i[k+1][][] - pG->B3i[k][j][i])*dx3i;
1287        B1 = pG->U[k][j][i].B1c;
1288        B2 = pG->U[k][j][i].B2c;
1289        B3 = pG->U[k][j][i].B3c;
1290        V1 = pG->U[k][j][i].M1/pG->U[k][j][i].d;
1291        V2 = pG->U[k][j][i].M2/pG->U[k][j][i].d;
1292
1293/* Calculate mdb1 = min_mod(-db3,db1) */
1294        if(db3 > 0.0 && db1 < 0.0){
1295          mdb1 = db1 > -db3 ? db1 : -db3;
1296        }
1297        else if(db3 < 0.0 && db1 > 0.0){
1298          mdb1 = db1 < -db3 ? db1 : -db3;
1299        }
1300        else mdb1 = 0.0;
1301
1302/* Calculate mdb2 = min_mod(-db3,db2) */
1303        if(db3 > 0.0 && db2 < 0.0){
1304          mdb2 = db2 > -db3 ? db2 : -db3;
1305        }
1306        else if(db3 < 0.0 && db2 > 0.0){
1307          mdb2 = db2 < -db3 ? db2 : -db3;
1308        }
1309        else mdb2 = 0.0;
1310
1311        Ur_x3Face[k][j][i].My += hdt*B1*db3;
1312        Ur_x3Face[k][j][i].Mz += hdt*B2*db3;
1313        Ur_x3Face[k][j][i].Mx += hdt*B3*db3;
1314        Ur_x3Face[k][j][i].By += hdt*V1*(-mdb2);
1315        Ur_x3Face[k][j][i].Bz += hdt*V2*(-mdb1);
1316#ifndef ISOTHERMAL
1317        Ur_x3Face[k][j][i].+= hdt*(B1*V1*(-mdb2) + B2*V2*(-mdb1) );
1318#endif /* ISOTHERMAL */
1319      }
1320    }
1321  }
1322#endif /* MHD */
1323
1324/*--- Step 8d ------------------------------------------------------------------
1325 * Add gravitational source terms in x1-direction and x2-direction to
1326 * corrected L/R states on x3-faces.  To improve conservation of total energy,
1327 * we average the energy source term computed at cell faces.
1328 *    S_{M} = -(\rho) Grad(Phi);   S_{E} = -(\rho v) Grad{Phi}
1329 */
1330
1331  if (StaticGravPot != NULL){
1332  for (k=ks-1; k<=ke+2; k++) {
1333    for (j=js-1; j<=je+1; j++) {
1334      for (i=is-1; i<=ie+1; i++) {
1335        cc_pos(pG,i,j,k,&x1,&x2,&x3);
1336        phic = (*StaticGravPot)((x1               ),x2,x3);
1337        phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
1338        phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
1339
1340/* correct right states */
1341        Ur_x3Face[k][j][i].My -= q1*(phir-phil)*pG->U[k][j][i].d;
1342#ifndef ISOTHERMAL
1343        Ur_x3Face[k][j][i].E -= q1*(x1Flux[][j][].d*(phic - phil)
1344                                  + x1Flux[][j][i+1].d*(phir - phic));
1345#endif
1346
1347        phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
1348        phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
1349
1350        Ur_x3Face[k][j][i].Mz -= q2*(phir-phil)*pG->U[k][j][i].d;
1351#ifndef ISOTHERMAL
1352        Ur_x3Face[k][j][i].E -= q2*(x2Flux[][][i].d*(phic - phil)
1353                                  + x2Flux[][j+1][i].d*(phir - phic));
1354#endif
1355
1356/* correct left states */
1357        phic = (*StaticGravPot)((x1               ),x2,(x3-pG->dx3));
1358        phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,(x3-pG->dx3));
1359        phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,(x3-pG->dx3));
1360
1361        Ul_x3Face[k][j][i].My -= q1*(phir-phil)*pG->U[k-1][j][i].d;
1362#ifndef ISOTHERMAL
1363        Ul_x3Face[k][j][i].E -= q1*(x1Flux[k-1][j][].d*(phic - phil)
1364                                  + x1Flux[k-1][j][i+1].d*(phir - phic));
1365#endif
1366
1367        phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),(x3-pG->dx3));
1368        phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),(x3-pG->dx3));
1369
1370        Ul_x3Face[k][j][i].Mz -= q2*(phir-phil)*pG->U[k-1][j][i].d;
1371#ifndef ISOTHERMAL
1372        Ul_x3Face[k][j][i].E -= q2*(x2Flux[k-1][][i].d*(phic - phil)
1373                                  + x2Flux[k-1][j+1][i].d*(phir - phic));
1374#endif
1375      }
1376    }
1377  }}
1378
1379/*--- Step 8e ------------------------------------------------------------------
1380 * Add source terms for self gravity to L/R states.
1381 *    S_{M} = -(\rho) Grad(Phi);   S_{E} = -(\rho v) Grad{Phi}
1382 */
1383
1384#ifdef SELF_GRAVITY
1385  for (k=ks-1; k<=ke+2; k++) {
1386    for (j=js-1; j<=je+1; j++) {
1387      for (i=is-1; i<=ie+1; i++) {
1388        phic = pG->Phi[k][j][i];
1389        phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j][i+1]);
1390        phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j][i-1]);
1391
1392/* correct right states */
1393        Ur_x3Face[k][j][i].My -= q1*(phir-phil)*pG->U[k][j][i].d;
1394#ifndef ISOTHERMAL
1395        Ur_x3Face[k][j][i].E -= q1*(x1Flux[k][j][].d*(phic - phil)
1396                                  + x1Flux[k][j][i+1].d*(phir - phic));
1397#endif
1398
1399        phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j+1][i]);
1400        phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j-1][i]);
1401
1402        Ur_x3Face[k][j][i].Mz -= q2*(phir-phil)*pG->U[k][j][i].d;
1403#ifndef ISOTHERMAL
1404        Ur_x3Face[k][j][i].E -= q2*(x2Flux[k][][i].d*(phic - phil)
1405                                  + x2Flux[k][j+1][i].d*(phir - phic));
1406#endif
1407
1408/* correct left states */
1409        phic = pG->Phi[k-1][j][i];
1410        phir = 0.5*(pG->Phi[k-1][j][i] + pG->Phi[k-1][j][i+1]);
1411        phil = 0.5*(pG->Phi[k-1][j][i] + pG->Phi[k-1][j][i-1]);
1412
1413        Ul_x3Face[k][j][i].My -= q1*(phir-phil)*pG->U[k-1][j][i].d;
1414#ifndef ISOTHERMAL
1415        Ul_x3Face[k][j][i].E -= q1*(x1Flux[k-1][j][].d*(phic - phil)
1416                                  + x1Flux[k-1][j][i+1].d*(phir - phic));
1417#endif
1418
1419        phir = 0.5*(pG->Phi[k-1][j][i] + pG->Phi[k-1][j+1][i]);
1420        phil = 0.5*(pG->Phi[k-1][j][i] + pG->Phi[k-1][j-1][i]);
1421
1422        Ul_x3Face[k][j][i].Mz -= q2*(phir-phil)*pG->U[k-1][j][i].d;
1423#ifndef ISOTHERMAL
1424        Ul_x3Face[k][j][i].E -= q2*(x2Flux[k-1][][i].d*(phic - phil)
1425                                  + x2Flux[k-1][j+1][i].d*(phir - phic));
1426#endif
1427      }
1428    }
1429  }
1430#endif /* SELF_GRAVITY */
1431
1432/*--- Step 9 -------------------------------------------------------------------
1433 * Calculate the cell centered value of emf1,2,3 at the half-time-step
1434 */
1435
1436  if (dhalf != NULL){
1437    for (k=ks-1; k<=ke+1; k++) {
1438      for (j=js-1; j<=je+1; j++) {
1439        for (i=is-1; i<=ie+1; i++) {
1440          dhalf[k][j][i] = pG->U[k][j][i].d
1441            - q1*(x1Flux[][][i+1].d - x1Flux[k][j][i].d)
1442            - q2*(x2Flux[][j+1][].d - x2Flux[k][j][i].d)
1443            - q3*(x3Flux[k+1][][].d - x3Flux[k][j][i].d);
1444        }
1445      }
1446    }
1447  }
1448
1449#ifdef MHD
1450  for (k=ks-1; k<=ke+1; k++) {
1451    for (j=js-1; j<=je+1; j++) {
1452      for (i=is-1; i<=ie+1; i++) {
1453        cc_pos(pG,i,j,k,&x1,&x2,&x3);
1454
1455        d  = dhalf[k][j][i];
1456
1457        M1 = pG->U[k][j][i].M1
1458           - q1*(x1Flux[][][i+1].Mx - x1Flux[k][j][i].Mx)
1459           - q2*(x2Flux[][j+1][].Mz - x2Flux[k][j][i].Mz)
1460           - q3*(x3Flux[k+1][][].My - x3Flux[k][j][i].My);
1461
1462        M2 = pG->U[k][j][i].M2
1463           - q1*(x1Flux[][][i+1].My - x1Flux[k][j][i].My)
1464           - q2*(x2Flux[][j+1][].Mx - x2Flux[k][j][i].Mx)
1465           - q3*(x3Flux[k+1][][].Mz - x3Flux[k][j][i].Mz);
1466
1467        M3 = pG->U[k][j][i].M3
1468           - q1*(x1Flux[][][i+1].Mz - x1Flux[k][j][i].Mz)
1469           - q2*(x2Flux[][j+1][].My - x2Flux[k][j][i].My)
1470           - q3*(x3Flux[k+1][][].Mx - x3Flux[k][j][i].Mx);
1471
1472/* Add source terms for fixed gravitational potential */
1473        if (StaticGravPot != NULL){
1474          phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
1475          phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
1476          M1 -= q1*(phir-phil)*pG->U[k][j][i].d;
1477
1478          phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
1479          phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
1480          M2 -= q2*(phir-phil)*pG->U[k][j][i].d;
1481
1482          phir = (*StaticGravPot)(x1,x2,(x3+0.5*pG->dx3));
1483          phil = (*StaticGravPot)(x1,x2,(x3-0.5*pG->dx3));
1484          M3 -= q3*(phir-phil)*pG->U[k][j][i].d;
1485        }
1486
1487/* Add source terms due to self-gravity  */
1488#ifdef SELF_GRAVITY
1489        phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j][i+1]);
1490        phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j][i-1]);
1491        M1 -= q1*(phir-phil)*pG->U[k][j][i].d;
1492
1493        phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j+1][i]);
1494        phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k][j-1][i]);
1495        M2 -= q2*(phir-phil)*pG->U[k][j][i].d;
1496
1497        phir = 0.5*(pG->Phi[k][j][i] + pG->Phi[k+1][j][i]);
1498        phil = 0.5*(pG->Phi[k][j][i] + pG->Phi[k-1][j][i]);
1499        M3 -= q3*(phir-phil)*pG->U[k][j][i].d;
1500#endif /* SELF_GRAVITY */
1501
1502        B1c = 0.5*(B1_x1Face[k][j][i] + B1_x1Face[][][i+1]);
1503        B2c = 0.5*(B2_x2Face[k][j][i] + B2_x2Face[][j+1][]);
1504        B3c = 0.5*(B3_x3Face[k][j][i] + B3_x3Face[k+1][][]);
1505
1506        emf1_cc[k][j][i] = (B2c*M3 - B3c*M2)/d;
1507
1508        emf2_cc[k][j][i] = (B3c*M1 - B1c*M3)/d;
1509
1510        emf3_cc[k][j][i] = (B1c*M2 - B2c*M1)/d;
1511      }
1512    }
1513  }
1514#endif
1515
1516/*--- Step 10a -----------------------------------------------------------------
1517 * Compute maximum wavespeeds in multidimensions (eta in eq. 10 from Sanders et
1518 *  al. (1998)) for H-correction
1519 */
1520
1521#ifdef H_CORRECTION
1522  for (k=ks-1; k<=ke+1; k++) {
1523    for (j=js-1; j<=je+1; j++) {
1524      for (i=is-1; i<=ie+2; i++) {
1525        cfr = cfast(&(Ur_x1Face[k][j][i]), &(B1_x1Face[k][j][i]));
1526        cfl = cfast(&(Ul_x1Face[k][j][i]), &(B1_x1Face[k][j][i]));
1527        lambdar = Ur_x1Face[k][j][i].Mx/Ur_x1Face[k][j][i].d + cfr;
1528        lambdal = Ul_x1Face[k][j][i].Mx/Ul_x1Face[k][j][i].d - cfl;
1529        eta1[k][j][i] = 0.5*fabs(lambdar - lambdal);
1530      }
1531    }
1532  }
1533
1534  for (k=ks-1; k<=ke+1; k++) {
1535    for (j=js-1; j<=je+2; j++) {
1536      for (i=is-1; i<=ie+1; i++) {
1537        cfr = cfast(&(Ur_x2Face[k][j][i]), &(B2_x2Face[k][j][i]));
1538        cfl = cfast(&(Ul_x2Face[k][j][i]), &(B2_x2Face[k][j][i]));
1539        lambdar = Ur_x2Face[k][j][i].Mx/Ur_x2Face[k][j][i].d + cfr;
1540        lambdal = Ul_x2Face[k][j][i].Mx/Ul_x2Face[k][j][i].d - cfl;
1541        eta2[k][j][i] = 0.5*fabs(lambdar - lambdal);
1542      }
1543    }
1544  }
1545
1546  for (k=ks-1; k<=ke+2; k++) {
1547    for (j=js-1; j<=je+1; j++) {
1548      for (i=is-1; i<=ie+1; i++) {
1549        cfr = cfast(&(Ur_x3Face[k][j][i]), &(B3_x3Face[k][j][i]));
1550        cfl = cfast(&(Ul_x3Face[k][j][i]), &(B3_x3Face[k][j][i]));
1551        lambdar = Ur_x3Face[k][j][i].Mx/Ur_x3Face[k][j][i].d + cfr;
1552        lambdal = Ul_x3Face[k][j][i].Mx/Ul_x3Face[k][j][i].d - cfl;
1553        eta3[k][j][i] = 0.5*fabs(lambdar - lambdal);
1554      }
1555    }
1556  }
1557#endif /* H_CORRECTION */
1558
1559/*--- Step 10b -----------------------------------------------------------------
1560 * Compute x1-fluxes from corrected L/R states.
1561 */
1562
1563  for (k=ks-1; k<=ke+1; k++) {
1564    for (j=js-1; j<=je+1; j++) {
1565      for (i=is; i<=ie+1; i++) {
1566#ifdef H_CORRECTION
1567        etah = MAX(eta2[k][j][i-1],eta2[k][j][i]);
1568        etah = MAX(etah,eta2[k][j+1][i-1]);
1569        etah = MAX(etah,eta2[k][j+1][]);
1570
1571        etah = MAX(etah,eta3[][j][i-1]);
1572        etah = MAX(etah,eta3[][j][]);
1573        etah = MAX(etah,eta3[k+1][j][i-1]);
1574        etah = MAX(etah,eta3[k+1][j][]);
1575
1576        etah = MAX(etah,eta1[][j][]);
1577#endif /* H_CORRECTION */
1578        Cons1D_to_Prim1D(&Ul_x1Face[k][j][i],&Wl[i],&B1_x1Face[k][j][i]);
1579        Cons1D_to_Prim1D(&Ur_x1Face[k][j][i],&Wr[i],&B1_x1Face[k][j][i]);
1580
1581        GET_FLUXES(Ul_x1Face[k][j][i],Ur_x1Face[k][j][i],Wl[i],Wr[i],
1582                   B1_x1Face[k][j][i],&x1Flux[k][j][i]);
1583      }
1584    }
1585  }
1586
1587/*--- Step 10c -----------------------------------------------------------------
1588 * Compute x2-fluxes from corrected L/R states.
1589 */
1590
1591  for (k=ks-1; k<=ke+1; k++) {
1592    for (j=js; j<=je+1; j++) {
1593      for (i=is-1; i<=ie+1; i++) {
1594#ifdef H_CORRECTION
1595        etah = MAX(eta1[k][j-1][i],eta1[k][j][i]);
1596        etah = MAX(etah,eta1[k][j-1][i+1]);
1597        etah = MAX(etah,eta1[k][][i+1]);
1598
1599        etah = MAX(etah,eta3[][j-1][i]);
1600        etah = MAX(etah,eta3[][][i]);
1601        etah = MAX(etah,eta3[k+1][j-1][i]);
1602        etah = MAX(etah,eta3[k+1][][i]);
1603
1604        etah = MAX(etah,eta2[][][i]);
1605#endif /* H_CORRECTION */
1606        Cons1D_to_Prim1D(&Ul_x2Face[k][j][i],&Wl[i],&B2_x2Face[k][j][i]);
1607        Cons1D_to_Prim1D(&Ur_x2Face[k][j][i],&Wr[i],&B2_x2Face[k][j][i]);
1608
1609        GET_FLUXES(Ul_x2Face[k][j][i],Ur_x2Face[k][j][i],Wl[i],Wr[i],
1610                   B2_x2Face[k][j][i],&x2Flux[k][j][i]);
1611      }
1612    }
1613  }
1614
1615/*--- Step 10d -----------------------------------------------------------------
1616 * Compute x3-fluxes from corrected L/R states.
1617 */
1618
1619  for (k=ks; k<=ke+1; k++) {
1620    for (j=js-1; j<=je+1; j++) {
1621      for (i=is-1; i<=ie+1; i++) {
1622#ifdef H_CORRECTION
1623        etah = MAX(eta1[k-1][j][i],eta1[k][j][i]);
1624        etah = MAX(etah,eta1[k-1][j][i+1]);
1625        etah = MAX(etah,eta1[k][][i+1]);
1626
1627        etah = MAX(etah,eta2[k-1][][i]);
1628        etah = MAX(etah,eta2[][][i]);
1629        etah = MAX(etah,eta2[k-1][j+1][i]);
1630        etah = MAX(etah,eta2[][j+1][i]);
1631
1632        etah = MAX(etah,eta3[][][i]);
1633#endif /* H_CORRECTION */
1634        Cons1D_to_Prim1D(&Ul_x3Face[k][j][i],&Wl[i],&B3_x3Face[k][j][i]);
1635        Cons1D_to_Prim1D(&Ur_x3Face[k][j][i],&Wr[i],&B3_x3Face[k][j][i]);
1636
1637        GET_FLUXES(Ul_x3Face[k][j][i],Ur_x3Face[k][j][i],Wl[i],Wr[i],
1638                   B3_x3Face[k][j][i],&x3Flux[k][j][i]);
1639      }
1640    }
1641  }
1642
1643/*--- Step 11 ------------------------------------------------------------------
1644 * Integrate emf*^{n+1/2} to the grid cell corners and then update the
1645 * interface magnetic fields using CT for a full time step.
1646 */
1647
1648#ifdef MHD
1649  integrate_emf1_corner(pG);
1650  integrate_emf2_corner(pG);
1651  integrate_emf3_corner(pG);
1652
1653  for (k=ks; k<=ke; k++) {
1654    for (j=js; j<=je; j++) {
1655      for (i=is; i<=ie; i++) {
1656        pG->B1i[k][j][i] += dtodx3*(emf2[k+1][][] - emf2[k][j][i]) -
1657                            dtodx2*(emf3[][j+1][] - emf3[k][j][i]);
1658        pG->B2i[k][j][i] += dtodx1*(emf3[][][i+1] - emf3[k][j][i]) -
1659                            dtodx3*(emf1[k+1][][] - emf1[k][j][i]);
1660        pG->B3i[k][j][i] += dtodx2*(emf1[][j+1][] - emf1[k][j][i]) -
1661                            dtodx1*(emf2[][][i+1] - emf2[k][j][i]);
1662      }
1663      pG->B1i[k][j][ie+1] +=
1664        dtodx3*(emf2[k+1][][ie+1] - emf2[k][j][ie+1]) -
1665        dtodx2*(emf3[][j+1][ie+1] - emf3[k][j][ie+1]);
1666    }
1667    for (i=is; i<=ie; i++) {
1668      pG->B2i[k][je+1][i] +=
1669        dtodx1*(emf3[][je+1][i+1] - emf3[k][je+1][i]) -
1670        dtodx3*(emf1[k+1][je+1][] - emf1[k][je+1][i]);
1671    }
1672  }
1673  for (j=js; j<=je; j++) {
1674    for (i=is; i<=ie; i++) {
1675      pG->B3i[ke+1][j][i] += 
1676        dtodx2*(emf1[ke+1][j+1][] - emf1[ke+1][j][i]) -
1677        dtodx1*(emf2[ke+1][][i+1] - emf2[ke+1][j][i]);
1678    }
1679  }
1680#endif
1681
1682/*--- Step 12a -----------------------------------------------------------------
1683 * Add the gravitational source terms due to a Static Potential, and the
1684 * shearing box source terms (written in terms of a Static Potential).
1685 *    S_{M} = -(\rho)^{n+1/2} Grad(Phi);   S_{E} = -(\rho v)^{n+1/2} Grad{Phi}
1686 */
1687
1688  if (StaticGravPot != NULL){
1689    for (k=ks; k<=ke; k++) {
1690      for (j=js; j<=je; j++) {
1691        for (i=is; i<=ie; i++) {
1692          cc_pos(pG,i,j,k,&x1,&x2,&x3);
1693          phic = (*StaticGravPot)(x1,x2,x3);
1694          phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
1695          phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
1696
1697          pG->U[k][j][i].M1 -= dtodx1*(phir-phil)*dhalf[k][j][i];
1698#ifndef ISOTHERMAL
1699          pG->U[k][j][i].E -= dtodx1*(x1Flux[k][j][].d*(phic - phil) +
1700                                      x1Flux[k][j][i+1].d*(phir - phic));
1701#endif
1702
1703          phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
1704          phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
1705
1706          pG->U[k][j][i].M2 -= dtodx2*(phir-phil)*dhalf[k][j][i];
1707#ifndef ISOTHERMAL
1708          pG->U[k][j][i].E -= dtodx2*(x2Flux[k][][i].d*(phic - phil) +
1709                                      x2Flux[k][j+1][i].d*(phir - phic));
1710#endif
1711
1712          phir = (*StaticGravPot)(x1,x2,(x3+0.5*pG->dx3));
1713          phil = (*StaticGravPot)(x1,x2,(x3-0.5*pG->dx3));
1714
1715          pG->U[k][j][i].M3 -= dtodx3*(phir-phil)*dhalf[k][j][i];
1716#ifndef ISOTHERMAL
1717          pG->U[k][j][i].E -= dtodx3*(x3Flux[][j][i].d*(phic - phil) +
1718                                      x3Flux[k+1][j][i].d*(phir - phic));
1719#endif
1720        }
1721      }
1722    }
1723  }
1724
1725/*--- Step 12b -----------------------------------------------------------------
1726 * Add gravitational source terms for self-gravity.
1727 * A flux correction using Phi^{n+1} in the main loop is required to make
1728 * the source terms 2nd order: see selfg_flux_correction().
1729 */
1730#ifdef SELF_GRAVITY
1731/* Add fluxes and source terms due to (d/dx1) terms  */
1732
1733  for (k=ks; k<=ke; k++){
1734    for (j=js; j<=je; j++){
1735      for (i=is; i<=ie; i++){
1736        phic = pG->Phi[k][j][i];
1737        phil = 0.5*(pG->Phi[k][j][i-1] + pG->Phi[k][j][]);
1738        phir = 0.5*(pG->Phi[k][j][] + pG->Phi[k][j][i+1]);
1739
1740/* gx, gy and gz centered at L and R x1-faces */
1741        gxl = (pG->Phi[k][j][i-1] - pG->Phi[k][j][])*(dx1i);
1742        gxr = (pG->Phi[k][j][] - pG->Phi[k][j][i+1])*(dx1i);
1743
1744        gyl = 0.25*((pG->Phi[k][j-1][i-1] - pG->Phi[k][j+1][i-1]) +
1745                    (pG->Phi[k][j-1][] - pG->Phi[k][j+1][]) )*(dx2i);
1746        gyr = 0.25*((pG->Phi[k][j-1][] - pG->Phi[k][j+1][]) +
1747                    (pG->Phi[k][j-1][i+1] - pG->Phi[k][j+1][i+1]) )*(dx2i);
1748
1749        gzl = 0.25*((pG->Phi[k-1][j][i-1] - pG->Phi[k+1][j][i-1]) +
1750                    (pG->Phi[k-1][j][] - pG->Phi[k+1][j][]) )*(dx3i);
1751        gzr = 0.25*((pG->Phi[k-1][j][] - pG->Phi[k+1][j][]) +
1752                    (pG->Phi[k-1][j][i+1] - pG->Phi[k+1][j][i+1]) )*(dx3i);
1753
1754/* momentum fluxes in x1.  2nd term is needed only if Jean's swindle used */
1755        flx_m1l = 0.5*(gxl*gxl-gyl*gyl-gzl*gzl)/four_pi_G + grav_mean_rho*phil;
1756        flx_m1r = 0.5*(gxr*gxr-gyr*gyr-gzr*gzr)/four_pi_G + grav_mean_rho*phir;
1757
1758        flx_m2l = gxl*gyl/four_pi_G;
1759        flx_m2r = gxr*gyr/four_pi_G;
1760
1761        flx_m3l = gxl*gzl/four_pi_G;
1762        flx_m3r = gxr*gzr/four_pi_G;
1763
1764/* Update momenta and energy with d/dx1 terms  */
1765        pG->U[k][j][i].M1 -= dtodx1*(flx_m1r - flx_m1l);
1766        pG->U[k][j][i].M2 -= dtodx1*(flx_m2r - flx_m2l);
1767        pG->U[k][j][i].M3 -= dtodx1*(flx_m3r - flx_m3l);
1768#ifdef ADIABATIC
1769        pG->U[k][j][i].E -= dtodx1*(x1Flux[k][j][].d*(phic - phil) +
1770                                    x1Flux[k][j][i+1].d*(phir - phic));
1771#endif /* ADIABATIC */
1772      }
1773    }
1774  }
1775
1776/* Add fluxes and source terms due to (d/dx2) terms  */
1777
1778  for (k=ks; k<=ke; k++){
1779    for (j=js; j<=je; j++){
1780      for (i=is; i<=ie; i++){
1781        phic = pG->Phi[k][j][i];
1782        phil = 0.5*(pG->Phi[k][j-1][i] + pG->Phi[k][][i]);
1783        phir = 0.5*(pG->Phi[k][][i] + pG->Phi[k][j+1][i]);
1784
1785/* gx, gy and gz centered at L and R x2-faces */
1786        gxl = 0.25*((pG->Phi[k][j-1][i-1] - pG->Phi[k][j-1][i+1]) +
1787                    (pG->Phi[k][][i-1] - pG->Phi[k][][i+1]) )*(dx1i);
1788        gxr = 0.25*((pG->Phi[k][][i-1] - pG->Phi[k][][i+1]) +
1789                    (pG->Phi[k][j+1][i-1] - pG->Phi[k][j+1][i+1]) )*(dx1i);
1790
1791        gyl = (pG->Phi[k][j-1][i] - pG->Phi[k][][i])*(dx2i);
1792        gyr = (pG->Phi[k][][i] - pG->Phi[k][j+1][i])*(dx2i);
1793
1794        gzl = 0.25*((pG->Phi[k-1][j-1][i] - pG->Phi[k+1][j-1][i]) +
1795                    (pG->Phi[k-1][][i] - pG->Phi[k+1][][i]) )*(dx3i);
1796        gzr = 0.25*((pG->Phi[k-1][][i] - pG->Phi[k+1][][i]) +
1797                    (pG->Phi[k-1][j+1][i] - pG->Phi[k+1][j+1][i]) )*(dx3i);
1798
1799/* momentum fluxes in x2.  2nd term is needed only if Jean's swindle used */
1800        flx_m1l = gyl*gxl/four_pi_G;
1801        flx_m1r = gyr*gxr/four_pi_G;
1802
1803        flx_m2l = 0.5*(gyl*gyl-gxl*gxl-gzl*gzl)/four_pi_G + grav_mean_rho*phil;
1804        flx_m2r = 0.5*(gyr*gyr-gxr*gxr-gzr*gzr)/four_pi_G + grav_mean_rho*phir;
1805
1806        flx_m3l = gyl*gzl/four_pi_G;
1807        flx_m3r = gyr*gzr/four_pi_G;
1808
1809/* Update momenta and energy with d/dx2 terms  */
1810        pG->U[k][j][i].M1 -= dtodx2*(flx_m1r - flx_m1l);
1811        pG->U[k][j][i].M2 -= dtodx2*(flx_m2r - flx_m2l);
1812        pG->U[k][j][i].M3 -= dtodx2*(flx_m3r - flx_m3l);
1813#ifdef ADIABATIC
1814        pG->U[k][j][i].E -= dtodx2*(x2Flux[k][][i].d*(phic - phil) +
1815                                    x2Flux[k][j+1][i].d*(phir - phic));
1816#endif /* ADIABATIC */
1817      }
1818    }
1819  }
1820
1821/* Add fluxes and source terms due to (d/dx3) terms  */
1822
1823  for (k=ks; k<=ke; k++){
1824    for (j=js; j<=je; j++){
1825      for (i=is; i<=ie; i++){
1826        phic = pG->Phi[k][j][i];
1827        phil = 0.5*(pG->Phi[k-1][j][i] + pG->Phi[][j][i]);
1828        phir = 0.5*(pG->Phi[][j][i] + pG->Phi[k+1][j][i]);
1829
1830/* gx, gy and gz centered at L and R x3-faces */
1831        gxl = 0.25*((pG->Phi[k-1][j][i-1] - pG->Phi[k-1][j][i+1]) +
1832                    (pG->Phi[][j][i-1] - pG->Phi[][j][i+1]) )*(dx1i);
1833        gxr = 0.25*((pG->Phi[][j][i-1] - pG->Phi[][j][i+1]) +
1834                    (pG->Phi[k+1][j][i-1] - pG->Phi[k+1][j][i+1]) )*(dx1i);
1835
1836        gyl = 0.25*((pG->Phi[k-1][j-1][i] - pG->Phi[k-1][j+1][i]) +
1837                    (pG->Phi[][j-1][i] - pG->Phi[][j+1][i]) )*(dx2i);
1838        gyr = 0.25*((pG->Phi[][j-1][i] - pG->Phi[][j+1][i]) +
1839                    (pG->Phi[k+1][j-1][i] - pG->Phi[k+1][j+1][i]) )*(dx2i);
1840
1841        gzl = (pG->Phi[k-1][j][i] - pG->Phi[][j][i])*(dx3i);
1842        gzr = (pG->Phi[][j][i] - pG->Phi[k+1][j][i])*(dx3i);
1843
1844/* momentum fluxes in x3.  2nd term is needed only if Jean's swindle used */
1845        flx_m1l = gzl*gxl/four_pi_G;
1846        flx_m1r = gzr*gxr/four_pi_G;
1847
1848        flx_m2l = gzl*gyl/four_pi_G;
1849        flx_m2r = gzr*gyr/four_pi_G;
1850
1851        flx_m3l = 0.5*(gzl*gzl-gxl*gxl-gyl*gyl)/four_pi_G + grav_mean_rho*phil;
1852        flx_m3r = 0.5*(gzr*gzr-gxr*gxr-gyr*gyr)/four_pi_G + grav_mean_rho*phir;
1853
1854/* Update momenta and energy with d/dx3 terms  */
1855        pG->U[k][j][i].M1 -= dtodx3*(flx_m1r - flx_m1l);
1856        pG->U[k][j][i].M2 -= dtodx3*(flx_m2r - flx_m2l);
1857        pG->U[k][j][i].M3 -= dtodx3*(flx_m3r - flx_m3l);
1858#ifdef ADIABATIC
1859        pG->U[k][j][i].E -= dtodx3*(x3Flux[][j][i].d*(phic - phil) +
1860                                    x3Flux[k+1][j][i].d*(phir - phic));
1861#endif /* ADIABATIC */
1862      }
1863    }
1864  }
1865
1866/* Save mass fluxes in Grid structure for source term correction in main loop */
1867
1868  for (k=ks; k<=ke+1; k++) {
1869    for (j=js; j<=je+1; j++) {
1870      for (i=is; i<=ie+1; i++) {
1871        pG->x1MassFlux[k][j][i] = x1Flux[k][j][i].d;
1872        pG->x2MassFlux[k][j][i] = x2Flux[k][j][i].d;
1873        pG->x3MassFlux[k][j][i] = x3Flux[k][j][i].d;
1874      }
1875    }
1876  }
1877#endif /* SELF_GRAVITY */
1878
1879/*--- Step 13a -----------------------------------------------------------------
1880 * Update cell-centered variables in pG using x1-fluxes
1881 */
1882
1883  for (k=ks; k<=ke; k++) {
1884    for (j=js; j<=je; j++) {
1885      for (i=is; i<=ie; i++) {
1886        pG->U[k][j][i].-= dtodx1*(x1Flux[k][j][i+1].d -x1Flux[k][j][i].d );
1887        pG->U[k][j][i].M1 -= dtodx1*(x1Flux[k][j][i+1].Mx-x1Flux[k][j][i].Mx);
1888        pG->U[k][j][i].M2 -= dtodx1*(x1Flux[k][j][i+1].My-x1Flux[k][j][i].My);
1889        pG->U[k][j][i].M3 -= dtodx1*(x1Flux[k][j][i+1].Mz-x1Flux[k][j][i].Mz);
1890#ifndef ISOTHERMAL
1891        pG->U[k][j][i].-= dtodx1*(x1Flux[k][j][i+1].E -x1Flux[k][j][i].E );
1892#endif /* ISOTHERMAL */
1893#ifdef MHD
1894        pG->U[k][j][i].B2c -= dtodx1*(x1Flux[k][j][i+1].By-x1Flux[k][j][i].By);
1895        pG->U[k][j][i].B3c -= dtodx1*(x1Flux[k][j][i+1].Bz-x1Flux[k][j][i].Bz);
1896#endif /* MHD */
1897#if (NSCALARS > 0)
1898        for (n=0; n<NSCALARS; n++)
1899          pG->U[k][j][i].s[n] -= dtodx1*(x1Flux[k][j][i+1].s[n]
1900                                       - x1Flux[k][j][].s[n]);
1901#endif
1902      }
1903    }
1904  }
1905
1906/*--- Step 13b -----------------------------------------------------------------
1907 * Update cell-centered variables in pG using x2-fluxes
1908 */
1909
1910  for (k=ks; k<=ke; k++) {
1911    for (j=js; j<=je; j++) {
1912      for (i=is; i<=ie; i++) {
1913        pG->U[k][j][i].-= dtodx2*(x2Flux[k][j+1][i].d -x2Flux[k][j][i].d );
1914        pG->U[k][j][i].M1 -= dtodx2*(x2Flux[k][j+1][i].Mz-x2Flux[k][j][i].Mz);
1915        pG->U[k][j][i].M2 -= dtodx2*(x2Flux[k][j+1][i].Mx-x2Flux[k][j][i].Mx);
1916        pG->U[k][j][i].M3 -= dtodx2*(x2Flux[k][j+1][i].My-x2Flux[k][j][i].My);
1917#ifndef ISOTHERMAL
1918        pG->U[k][j][i].E -=dtodx2*(x2Flux[k][j+1][i].E -x2Flux[k][j][i].E );
1919#endif /* ISOTHERMAL */
1920#ifdef MHD
1921        pG->U[k][j][i].B3c -= dtodx2*(x2Flux[k][j+1][i].By-x2Flux[k][j][i].By);
1922        pG->U[k][j][i].B1c -= dtodx2*(x2Flux[k][j+1][i].Bz-x2Flux[k][j][i].Bz);
1923#endif /* MHD */
1924#if (NSCALARS > 0)
1925        for (n=0; n<NSCALARS; n++)
1926          pG->U[k][j][i].s[n] -= dtodx2*(x2Flux[k][j+1][i].s[n]
1927                                       - x2Flux[k][][i].s[n]);
1928#endif
1929      }
1930    }
1931  }
1932
1933/*--- Step 13c -----------------------------------------------------------------
1934 * Update cell-centered variables in pG using x3-fluxes
1935 */
1936
1937  for (k=ks; k<=ke; k++) {
1938    for (j=js; j<=je; j++) {
1939      for (i=is; i<=ie; i++) {
1940        pG->U[k][j][i].-= dtodx3*(x3Flux[k+1][j][i].d -x3Flux[k][j][i].d );
1941        pG->U[k][j][i].M1 -= dtodx3*(x3Flux[k+1][j][i].My-x3Flux[k][j][i].My);
1942        pG->U[k][j][i].M2 -= dtodx3*(x3Flux[k+1][j][i].Mz-x3Flux[k][j][i].Mz);
1943        pG->U[k][j][i].M3 -= dtodx3*(x3Flux[k+1][j][i].Mx-x3Flux[k][j][i].Mx);
1944#ifndef ISOTHERMAL
1945        pG->U[k][j][i].-= dtodx3*(x3Flux[k+1][j][i].E -x3Flux[k][j][i].E );
1946#endif /* ISOTHERMAL */
1947#ifdef MHD
1948        pG->U[k][j][i].B1c -= dtodx3*(x3Flux[k+1][j][i].By-x3Flux[k][j][i].By);
1949        pG->U[k][j][i].B2c -= dtodx3*(x3Flux[k+1][j][i].Bz-x3Flux[k][j][i].Bz);
1950#endif /* MHD */
1951#if (NSCALARS > 0)
1952        for (n=0; n<NSCALARS; n++)
1953          pG->U[k][j][i].s[n] -= dtodx3*(x3Flux[k+1][j][i].s[n]
1954                                       - x3Flux[][j][i].s[n]);
1955#endif
1956      }
1957    }
1958  }
1959
1960/*--- Step 14 ------------------------------------------------------------------
1961 * LAST STEP!
1962 * Set cell centered magnetic fields to average of updated face centered fields.
1963 */
1964
1965#ifdef MHD
1966  for (k=ks; k<=ke; k++) {
1967    for (j=js; j<=je; j++) {
1968      for (i=is; i<=ie; i++) {
1969        pG->U[k][j][i].B1c = 0.5*(pG->B1i[k][j][i]+pG->B1i[k][j][i+1]);
1970        pG->U[k][j][i].B2c = 0.5*(pG->B2i[k][j][i]+pG->B2i[k][j+1][i]);
1971        pG->U[k][j][i].B3c = 0.5*(pG->B3i[k][j][i]+pG->B3i[k+1][j][i]);
1972      }
1973    }
1974  }
1975#endif /* MHD */
1976
1977  return;
1978}
1979
1980/*----------------------------------------------------------------------------*/
1981/* integrate_destruct_3d:  Free temporary integration arrays
1982 */
1983
1984void integrate_destruct_3d(void)
1985{
1986
1987#ifdef MHD
1988  if (emf1    != NULL) free_3d_array(emf1);
1989  if (emf2    != NULL) free_3d_array(emf2);
1990  if (emf3    != NULL) free_3d_array(emf3);
1991  if (emf1_cc != NULL) free_3d_array(emf1_cc);
1992  if (emf2_cc != NULL) free_3d_array(emf2_cc);
1993  if (emf3_cc != NULL) free_3d_array(emf3_cc);
1994#endif /* MHD */
1995#ifdef H_CORRECTION
1996  if (eta1 != NULL) free_3d_array(eta1);
1997  if (eta2 != NULL) free_3d_array(eta2);
1998  if (eta3 != NULL) free_3d_array(eta3);
1999#endif /* H_CORRECTION */
2000
2001  if (Bxc != NULL) free(Bxc);
2002  if (Bxi != NULL) free(Bxi);
2003  if (B1_x1Face != NULL) free_3d_array(B1_x1Face);
2004  if (B2_x2Face != NULL) free_3d_array(B2_x2Face);
2005  if (B3_x3Face != NULL) free_3d_array(B3_x3Face);
2006
2007  if (U1d      != NULL) free(U1d);
2008  if (W        != NULL) free(W);
2009  if (Wl       != NULL) free(Wl);
2010  if (Wr       != NULL) free(Wr);
2011
2012  if (Ul_x1Face != NULL) free_3d_array(Ul_x1Face);
2013  if (Ur_x1Face != NULL) free_3d_array(Ur_x1Face);
2014  if (Ul_x2Face != NULL) free_3d_array(Ul_x2Face);
2015  if (Ur_x2Face != NULL) free_3d_array(Ur_x2Face);
2016  if (Ul_x3Face != NULL) free_3d_array(Ul_x3Face);
2017  if (Ur_x3Face != NULL) free_3d_array(Ur_x3Face);
2018  if (x1Flux    != NULL) free_3d_array(x1Flux);
2019  if (x2Flux    != NULL) free_3d_array(x2Flux);
2020  if (x3Flux    != NULL) free_3d_array(x3Flux);
2021  if (dhalf     != NULL) free_3d_array(dhalf);
2022
2023  return;
2024}
2025
2026/*----------------------------------------------------------------------------*/
2027/* integrate_init_3d: Allocate temporary integration arrays
2028*/
2029
2030void integrate_init_3d(int nx1, int nx2, int nx3)
2031{
2032  int nmax;
2033  int Nx1 = nx1 + 2*nghost;
2034  int Nx2 = nx2 + 2*nghost;
2035  int Nx3 = nx3 + 2*nghost;
2036  nmax = MAX(MAX(Nx1,Nx2),Nx3);
2037
2038#ifdef MHD
2039  if ((emf1 = (Real***)calloc_3d_array(Nx3, Nx2, Nx1, sizeof(Real))) == NULL)
2040    goto on_error;
2041  if ((emf2 = (Real***)calloc_3d_array(Nx3, Nx2, Nx1, sizeof(Real))) == NULL)
2042    goto on_error;
2043  if ((emf3 = (Real***)calloc_3d_array(Nx3, Nx2, Nx1, sizeof(Real))) == NULL)
2044    goto on_error;
2045
2046  if ((emf1_cc = (Real***)calloc_3d_array(Nx3, Nx2, Nx1, sizeof(Real))) == NULL)
2047    goto on_error;
2048  if ((emf2_cc = (Real***)calloc_3d_array(Nx3, Nx2, Nx1, sizeof(Real))) == NULL)
2049    goto on_error;
2050  if ((emf3_cc = (Real***)calloc_3d_array(Nx3, Nx2, Nx1, sizeof(Real))) == NULL)
2051    goto on_error;
2052#endif /* MHD */
2053#ifdef H_CORRECTION
2054  if ((eta1 = (Real***)calloc_3d_array(Nx3, Nx2, Nx1, sizeof(Real))) == NULL)
2055    goto on_error;
2056  if ((eta2 = (Real***)calloc_3d_array(Nx3, Nx2, Nx1, sizeof(Real))) == NULL)
2057    goto on_error;
2058  if ((eta3 = (Real***)calloc_3d_array(Nx3, Nx2, Nx1, sizeof(Real))) == NULL)
2059    goto on_error;
2060#endif /* H_CORRECTION */
2061
2062  if ((Bxc = (Real*)malloc(nmax*sizeof(Real))) == NULL) goto on_error;
2063  if ((Bxi = (Real*)malloc(nmax*sizeof(Real))) == NULL) goto on_error;
2064
2065  if ((B1_x1Face = (Real***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Real))) == NULL)
2066    goto on_error;
2067  if ((B2_x2Face = (Real***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Real))) == NULL)
2068    goto on_error;
2069  if ((B3_x3Face = (Real***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Real))) == NULL)
2070    goto on_error;
2071
2072  if ((U1d =      (Cons1D*)malloc(nmax*sizeof(Cons1D))) == NULL) goto on_error;
2073  if ((W   =      (Prim1D*)malloc(nmax*sizeof(Prim1D))) == NULL) goto on_error;
2074  if ((Wl  =      (Prim1D*)malloc(nmax*sizeof(Prim1D))) == NULL) goto on_error;
2075  if ((Wr  =      (Prim1D*)malloc(nmax*sizeof(Prim1D))) == NULL) goto on_error;
2076
2077  if ((Ul_x1Face = (Cons1D***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Cons1D)))
2078    == NULL) goto on_error;
2079  if ((Ur_x1Face = (Cons1D***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Cons1D)))
2080    == NULL) goto on_error;
2081  if ((Ul_x2Face = (Cons1D***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Cons1D)))
2082    == NULL) goto on_error;
2083  if ((Ur_x2Face = (Cons1D***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Cons1D)))
2084    == NULL) goto on_error;
2085  if ((Ul_x3Face = (Cons1D***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Cons1D))) 
2086    == NULL) goto on_error;
2087  if ((Ur_x3Face = (Cons1D***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Cons1D))) 
2088    == NULL) goto on_error;
2089  if ((x1Flux    = (Cons1D***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Cons1D))) 
2090    == NULL) goto on_error;
2091  if ((x2Flux    = (Cons1D***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Cons1D))) 
2092    == NULL) goto on_error;
2093  if ((x3Flux    = (Cons1D***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Cons1D))) 
2094    == NULL) goto on_error;
2095
2096
2097#if defined MHD
2098  if ((dhalf = (Real***)calloc_3d_array(Nx3, Nx2, Nx1, sizeof(Real))) == NULL)
2099    goto on_error;
2100#else
2101  if(StaticGravPot != NULL){
2102    if ((dhalf = (Real***)calloc_3d_array(Nx3, Nx2, Nx1, sizeof(Real))) == NULL)
2103      goto on_error;
2104  }
2105#endif
2106
2107  return;
2108
2109  on_error:
2110  integrate_destruct();
2111  ath_error("[integrate_init]: malloc returned a NULL pointer\n");
2112}
2113
2114
2115/*=========================== PRIVATE FUNCTIONS ==============================*/
2116
2117/*----------------------------------------------------------------------------*/
2118/* integrate_emf1_corner
2119 * integrate_emf2_corner
2120 * integrate_emf3_corner
2121 *   Integrates face centered B-fluxes to compute corner EMFs.  Note:
2122 *   x1Flux.By = VxBy - BxVy = v1*b2-b1*v2 = -EMFZ
2123 *   x1Flux.Bz = VxBz - BxVz = v1*b3-b1*v3 = EMFY
2124 *   x2Flux.By = VxBy - BxVy = v2*b3-b2*v3 = -EMFX
2125 *   x2Flux.Bz = VxBz - BxVz = v2*b1-b2*v1 = EMFZ
2126 *   x3Flux.By = VxBy - BxVy = v3*b1-b3*v1 = -EMFY
2127 *   x3Flux.Bz = VxBz - BxVz = v3*b2-b3*v2 = EMFX
2128 */
2129
2130#ifdef MHD
2131static void integrate_emf1_corner(const Grid *pG)
2132{
2133  int i, is = pG->is, ie = pG->ie;
2134  int j, js = pG->js, je = pG->je;
2135  int k, ks = pG->ks, ke = pG->ke;
2136  Real de1_l2, de1_r2, de1_l3, de1_r3;
2137
2138  for (k=ks-1; k<=ke+2; k++) {
2139    for (j=js-1; j<=je+2; j++) {
2140      for (i=is-2; i<=ie+2; i++) {
2141/* NOTE: The x2-Flux of By is -E1. */
2142/*       The x3-Flux of Bz is +E1. */
2143        if (x2Flux[k-1][j][i].d > 0.0)
2144          de1_l3 = x3Flux[k][j-1][i].Bz - emf1_cc[k-1][j-1][i];
2145        else if (x2Flux[k-1][j][i].d < 0.0)
2146          de1_l3 = x3Flux[k][j][i].Bz - emf1_cc[k-1][j][i];
2147        else {
2148          de1_l3 = 0.5*(x3Flux[k][j-1][i].Bz - emf1_cc[k-1][j-1][i] +
2149                        x3Flux[k][][i].Bz - emf1_cc[k-1][][i] );
2150        }
2151
2152        if (x2Flux[k][j][i].d > 0.0)
2153          de1_r3 = x3Flux[k][j-1][i].Bz - emf1_cc[k][j-1][i];
2154        else if (x2Flux[k][j][i].d < 0.0)
2155          de1_r3 = x3Flux[k][j][i].Bz - emf1_cc[k][j][i];
2156        else {
2157          de1_r3 = 0.5*(x3Flux[k][j-1][i].Bz - emf1_cc[k][j-1][i] +
2158                        x3Flux[k][][i].Bz - emf1_cc[k][][i] );
2159        }
2160
2161        if (x3Flux[k][j-1][i].d > 0.0)
2162          de1_l2 = -x2Flux[k-1][j][i].By - emf1_cc[k-1][j-1][i];
2163        else if (x3Flux[k][j-1][i].d < 0.0)
2164          de1_l2 = -x2Flux[k][j][i].By - emf1_cc[k][j-1][i];
2165        else {
2166          de1_l2 = 0.5*(-x2Flux[k-1][j][i].By - emf1_cc[k-1][j-1][i]
2167                        -x2Flux[][j][i].By - emf1_cc[][j-1][i] );
2168        }
2169
2170        if (x3Flux[k][j][i].d > 0.0)
2171          de1_r2 = -x2Flux[k-1][j][i].By - emf1_cc[k-1][j][i];
2172        else if (x3Flux[k][j][i].d < 0.0)
2173          de1_r2 = -x2Flux[k][j][i].By - emf1_cc[k][j][i];
2174        else {
2175          de1_r2 = 0.5*(-x2Flux[k-1][j][i].By - emf1_cc[k-1][j][i]
2176                        -x2Flux[][j][i].By - emf1_cc[][j][i] );
2177        }
2178
2179        emf1[k][j][i] = 0.25*(  x3Flux[k][j][i].Bz + x3Flux[k][j-1][i].Bz
2180                              - x2Flux[k][j][i].By - x2Flux[k-1][j][i].By
2181                              + de1_l2 + de1_r2 + de1_l3 + de1_r3);
2182      }
2183    }
2184  }
2185
2186  return;
2187}
2188
2189
2190static void integrate_emf2_corner(const Grid *pG)
2191{
2192  int i, is = pG->is, ie = pG->ie;
2193  int j, js = pG->js, je = pG->je;
2194  int k, ks = pG->ks, ke = pG->ke;
2195  Real de2_l1, de2_r1, de2_l3, de2_r3;
2196
2197  for (k=ks-1; k<=ke+2; k++) {
2198    for (j=js-2; j<=je+2; j++) {
2199      for (i=is-1; i<=ie+2; i++) {
2200/* NOTE: The x1-Flux of Bz is +E2. */
2201/*       The x3-Flux of By is -E2. */
2202        if (x1Flux[k-1][j][i].d > 0.0)
2203          de2_l3 = -x3Flux[k][j][i-1].By - emf2_cc[k-1][j][i-1];
2204        else if (x1Flux[k-1][j][i].d < 0.0)
2205          de2_l3 = -x3Flux[k][j][i].By - emf2_cc[k-1][j][i];
2206        else {
2207          de2_l3 = 0.5*(-x3Flux[k][j][i-1].By - emf2_cc[k-1][j][i-1] 
2208                        -x3Flux[k][j][].By - emf2_cc[k-1][j][] );
2209        }
2210
2211        if (x1Flux[k][j][i].d > 0.0)
2212          de2_r3 = -x3Flux[k][j][i-1].By - emf2_cc[k][j][i-1];
2213        else if (x1Flux[k][j][i].d < 0.0)
2214          de2_r3 = -x3Flux[k][j][i].By - emf2_cc[k][j][i];
2215        else {
2216          de2_r3 = 0.5*(-x3Flux[k][j][i-1].By - emf2_cc[k][j][i-1] 
2217                        -x3Flux[k][j][].By - emf2_cc[k][j][] );
2218        }
2219
2220        if (x3Flux[k][j][i-1].d > 0.0)
2221          de2_l1 = x1Flux[k-1][j][i].Bz - emf2_cc[k-1][j][i-1];
2222        else if (x3Flux[k][j][i-1].d < 0.0)
2223          de2_l1 = x1Flux[k][j][i].Bz - emf2_cc[k][j][i-1];
2224        else {
2225          de2_l1 = 0.5*(x1Flux[k-1][j][i].Bz - emf2_cc[k-1][j][i-1] +
2226                        x1Flux[][j][i].Bz - emf2_cc[][j][i-1] );
2227        }
2228
2229        if (x3Flux[k][j][i].d > 0.0)
2230          de2_r1 = x1Flux[k-1][j][i].Bz - emf2_cc[k-1][j][i];
2231        else if (x3Flux[k][j][i].d < 0.0)
2232          de2_r1 = x1Flux[k][j][i].Bz - emf2_cc[k][j][i];
2233        else {
2234          de2_r1 = 0.5*(x1Flux[k-1][j][i].Bz - emf2_cc[k-1][j][i] +
2235                        x1Flux[][j][i].Bz - emf2_cc[k-1][j][i] );
2236        }
2237
2238        emf2[k][j][i] = 0.25*(  x1Flux[k][j][i].Bz + x1Flux[k-1][j][].Bz
2239                              - x3Flux[k][j][i].By - x3Flux[][j][i-1].By
2240                              + de2_l1 + de2_r1 + de2_l3 + de2_r3);
2241      }
2242    }
2243  }
2244
2245  return;
2246}
2247
2248static void integrate_emf3_corner(const Grid *pG)
2249{
2250  int i, is = pG->is, ie = pG->ie;
2251  int j, js = pG->js, je = pG->je;
2252  int k, ks = pG->ks, ke = pG->ke;
2253  Real de3_l1, de3_r1, de3_l2, de3_r2;
2254
2255  for (k=ks-2; k<=ke+2; k++) {
2256    for (j=js-1; j<=je+2; j++) {
2257      for (i=is-1; i<=ie+2; i++) {
2258/* NOTE: The x1-Flux of By is -E3. */
2259/*       The x2-Flux of Bx is +E3. */
2260        if (x1Flux[k][j-1][i].d > 0.0)
2261          de3_l2 = x2Flux[k][j][i-1].Bz - emf3_cc[k][j-1][i-1];
2262        else if (x1Flux[k][j-1][i].d < 0.0)
2263          de3_l2 = x2Flux[k][j][i].Bz - emf3_cc[k][j-1][i];
2264        else {
2265          de3_l2 = 0.5*(x2Flux[k][j][i-1].Bz - emf3_cc[k][j-1][i-1] + 
2266                        x2Flux[k][j][].Bz - emf3_cc[k][j-1][] );
2267        }
2268
2269        if (x1Flux[k][j][i].d > 0.0)
2270          de3_r2 = x2Flux[k][j][i-1].Bz - emf3_cc[k][j][i-1];
2271        else if (x1Flux[k][j][i].d < 0.0)
2272          de3_r2 = x2Flux[k][j][i].Bz - emf3_cc[k][j][i];
2273        else {
2274          de3_r2 = 0.5*(x2Flux[k][j][i-1].Bz - emf3_cc[k][j][i-1] + 
2275                        x2Flux[k][j][].Bz - emf3_cc[k][j][] );
2276        }
2277
2278        if (x2Flux[k][j][i-1].d > 0.0)
2279          de3_l1 = -x1Flux[k][j-1][i].By - emf3_cc[k][j-1][i-1];
2280        else if (x2Flux[k][j][i-1].d < 0.0)
2281          de3_l1 = -x1Flux[k][j][i].By - emf3_cc[k][j][i-1];
2282        else {
2283          de3_l1 = 0.5*(-x1Flux[k][j-1][i].By - emf3_cc[k][j-1][i-1]
2284                        -x1Flux[k][][i].By - emf3_cc[k][][i-1] );
2285        }
2286
2287        if (x2Flux[k][j][i].d > 0.0)
2288          de3_r1 = -x1Flux[k][j-1][i].By - emf3_cc[k][j-1][i];
2289        else if (x2Flux[k][j][i].d < 0.0)
2290          de3_r1 = -x1Flux[k][j][i].By - emf3_cc[k][j][i];
2291        else {
2292          de3_r1 = 0.5*(-x1Flux[k][j-1][i].By - emf3_cc[k][j-1][i]
2293                        -x1Flux[k][][i].By - emf3_cc[k][][i] );
2294        }
2295
2296        emf3[k][j][i] = 0.25*(  x2Flux[k][][i-1].Bz + x2Flux[k][j][i].Bz
2297                              - x1Flux[k][j-1][].By - x1Flux[k][j][i].By
2298                              + de3_l1 + de3_r1 + de3_l2 + de3_r2);
2299      }
2300    }
2301  }
2302
2303  return;
2304}
2305#endif /* MHD */
Note: See TracBrowser for help on using the repository browser.