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

Revision 607, 79.5 KB checked in by jstone, 6 years ago (diff)

Removed shearing_box_evolution, since do not want to distribute.

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        if (StaticGravPot != NULL){
1462          phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
1463          phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
1464          M1 -= q1*(phir-phil)*pG->U[k][j][i].d;
1465        }
1466
1467        M2 = pG->U[k][j][i].M2
1468           - q1*(x1Flux[][][i+1].My - x1Flux[k][j][i].My)
1469           - q2*(x2Flux[][j+1][].Mx - x2Flux[k][j][i].Mx)
1470           - q3*(x3Flux[k+1][][].Mz - x3Flux[k][j][i].Mz);
1471        if (StaticGravPot != NULL){
1472          phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
1473          phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
1474          M2 -= q2*(phir-phil)*pG->U[k][j][i].d;
1475        }
1476
1477        M3 = pG->U[k][j][i].M3
1478           - q1*(x1Flux[][][i+1].Mz - x1Flux[k][j][i].Mz)
1479           - q2*(x2Flux[][j+1][].My - x2Flux[k][j][i].My)
1480           - q3*(x3Flux[k+1][][].Mx - x3Flux[k][j][i].Mx);
1481        if (StaticGravPot != NULL){
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        B1c = 0.5*(B1_x1Face[k][j][i] + B1_x1Face[][][i+1]);
1488        B2c = 0.5*(B2_x2Face[k][j][i] + B2_x2Face[][j+1][]);
1489        B3c = 0.5*(B3_x3Face[k][j][i] + B3_x3Face[k+1][][]);
1490
1491        emf1_cc[k][j][i] = (B2c*M3 - B3c*M2)/d;
1492
1493        emf2_cc[k][j][i] = (B3c*M1 - B1c*M3)/d;
1494
1495        emf3_cc[k][j][i] = (B1c*M2 - B2c*M1)/d;
1496      }
1497    }
1498  }
1499#endif
1500
1501/*--- Step 10a -----------------------------------------------------------------
1502 * Compute maximum wavespeeds in multidimensions (eta in eq. 10 from Sanders et
1503 *  al. (1998)) for H-correction
1504 */
1505
1506#ifdef H_CORRECTION
1507  for (k=ks-1; k<=ke+1; k++) {
1508    for (j=js-1; j<=je+1; j++) {
1509      for (i=is-1; i<=ie+2; i++) {
1510        cfr = cfast(&(Ur_x1Face[k][j][i]), &(B1_x1Face[k][j][i]));
1511        cfl = cfast(&(Ul_x1Face[k][j][i]), &(B1_x1Face[k][j][i]));
1512        lambdar = Ur_x1Face[k][j][i].Mx/Ur_x1Face[k][j][i].d + cfr;
1513        lambdal = Ul_x1Face[k][j][i].Mx/Ul_x1Face[k][j][i].d - cfl;
1514        eta1[k][j][i] = 0.5*fabs(lambdar - lambdal);
1515      }
1516    }
1517  }
1518
1519  for (k=ks-1; k<=ke+1; k++) {
1520    for (j=js-1; j<=je+2; j++) {
1521      for (i=is-1; i<=ie+1; i++) {
1522        cfr = cfast(&(Ur_x2Face[k][j][i]), &(B2_x2Face[k][j][i]));
1523        cfl = cfast(&(Ul_x2Face[k][j][i]), &(B2_x2Face[k][j][i]));
1524        lambdar = Ur_x2Face[k][j][i].Mx/Ur_x2Face[k][j][i].d + cfr;
1525        lambdal = Ul_x2Face[k][j][i].Mx/Ul_x2Face[k][j][i].d - cfl;
1526        eta2[k][j][i] = 0.5*fabs(lambdar - lambdal);
1527      }
1528    }
1529  }
1530
1531  for (k=ks-1; k<=ke+2; k++) {
1532    for (j=js-1; j<=je+1; j++) {
1533      for (i=is-1; i<=ie+1; i++) {
1534        cfr = cfast(&(Ur_x3Face[k][j][i]), &(B3_x3Face[k][j][i]));
1535        cfl = cfast(&(Ul_x3Face[k][j][i]), &(B3_x3Face[k][j][i]));
1536        lambdar = Ur_x3Face[k][j][i].Mx/Ur_x3Face[k][j][i].d + cfr;
1537        lambdal = Ul_x3Face[k][j][i].Mx/Ul_x3Face[k][j][i].d - cfl;
1538        eta3[k][j][i] = 0.5*fabs(lambdar - lambdal);
1539      }
1540    }
1541  }
1542#endif /* H_CORRECTION */
1543
1544/*--- Step 10b -----------------------------------------------------------------
1545 * Compute x1-fluxes from corrected L/R states.
1546 */
1547
1548  for (k=ks-1; k<=ke+1; k++) {
1549    for (j=js-1; j<=je+1; j++) {
1550      for (i=is; i<=ie+1; i++) {
1551#ifdef H_CORRECTION
1552        etah = MAX(eta2[k][j][i-1],eta2[k][j][i]);
1553        etah = MAX(etah,eta2[k][j+1][i-1]);
1554        etah = MAX(etah,eta2[k][j+1][]);
1555
1556        etah = MAX(etah,eta3[][j][i-1]);
1557        etah = MAX(etah,eta3[][j][]);
1558        etah = MAX(etah,eta3[k+1][j][i-1]);
1559        etah = MAX(etah,eta3[k+1][j][]);
1560
1561        etah = MAX(etah,eta1[][j][]);
1562#endif /* H_CORRECTION */
1563        Cons1D_to_Prim1D(&Ul_x1Face[k][j][i],&Wl[i],&B1_x1Face[k][j][i]);
1564        Cons1D_to_Prim1D(&Ur_x1Face[k][j][i],&Wr[i],&B1_x1Face[k][j][i]);
1565
1566        GET_FLUXES(Ul_x1Face[k][j][i],Ur_x1Face[k][j][i],Wl[i],Wr[i],
1567                   B1_x1Face[k][j][i],&x1Flux[k][j][i]);
1568      }
1569    }
1570  }
1571
1572/*--- Step 10c -----------------------------------------------------------------
1573 * Compute x2-fluxes from corrected L/R states.
1574 */
1575
1576  for (k=ks-1; k<=ke+1; k++) {
1577    for (j=js; j<=je+1; j++) {
1578      for (i=is-1; i<=ie+1; i++) {
1579#ifdef H_CORRECTION
1580        etah = MAX(eta1[k][j-1][i],eta1[k][j][i]);
1581        etah = MAX(etah,eta1[k][j-1][i+1]);
1582        etah = MAX(etah,eta1[k][][i+1]);
1583
1584        etah = MAX(etah,eta3[][j-1][i]);
1585        etah = MAX(etah,eta3[][][i]);
1586        etah = MAX(etah,eta3[k+1][j-1][i]);
1587        etah = MAX(etah,eta3[k+1][][i]);
1588
1589        etah = MAX(etah,eta2[][][i]);
1590#endif /* H_CORRECTION */
1591        Cons1D_to_Prim1D(&Ul_x2Face[k][j][i],&Wl[i],&B2_x2Face[k][j][i]);
1592        Cons1D_to_Prim1D(&Ur_x2Face[k][j][i],&Wr[i],&B2_x2Face[k][j][i]);
1593
1594        GET_FLUXES(Ul_x2Face[k][j][i],Ur_x2Face[k][j][i],Wl[i],Wr[i],
1595                   B2_x2Face[k][j][i],&x2Flux[k][j][i]);
1596      }
1597    }
1598  }
1599
1600/*--- Step 10d -----------------------------------------------------------------
1601 * Compute x3-fluxes from corrected L/R states.
1602 */
1603
1604  for (k=ks; k<=ke+1; k++) {
1605    for (j=js-1; j<=je+1; j++) {
1606      for (i=is-1; i<=ie+1; i++) {
1607#ifdef H_CORRECTION
1608        etah = MAX(eta1[k-1][j][i],eta1[k][j][i]);
1609        etah = MAX(etah,eta1[k-1][j][i+1]);
1610        etah = MAX(etah,eta1[k][][i+1]);
1611
1612        etah = MAX(etah,eta2[k-1][][i]);
1613        etah = MAX(etah,eta2[][][i]);
1614        etah = MAX(etah,eta2[k-1][j+1][i]);
1615        etah = MAX(etah,eta2[][j+1][i]);
1616
1617        etah = MAX(etah,eta3[][][i]);
1618#endif /* H_CORRECTION */
1619        Cons1D_to_Prim1D(&Ul_x3Face[k][j][i],&Wl[i],&B3_x3Face[k][j][i]);
1620        Cons1D_to_Prim1D(&Ur_x3Face[k][j][i],&Wr[i],&B3_x3Face[k][j][i]);
1621
1622        GET_FLUXES(Ul_x3Face[k][j][i],Ur_x3Face[k][j][i],Wl[i],Wr[i],
1623                   B3_x3Face[k][j][i],&x3Flux[k][j][i]);
1624      }
1625    }
1626  }
1627
1628/*--- Step 11 ------------------------------------------------------------------
1629 * Integrate emf*^{n+1/2} to the grid cell corners and then update the
1630 * interface magnetic fields using CT for a full time step.
1631 */
1632
1633#ifdef MHD
1634  integrate_emf1_corner(pG);
1635  integrate_emf2_corner(pG);
1636  integrate_emf3_corner(pG);
1637
1638  for (k=ks; k<=ke; k++) {
1639    for (j=js; j<=je; j++) {
1640      for (i=is; i<=ie; i++) {
1641        pG->B1i[k][j][i] += dtodx3*(emf2[k+1][][] - emf2[k][j][i]) -
1642                            dtodx2*(emf3[][j+1][] - emf3[k][j][i]);
1643        pG->B2i[k][j][i] += dtodx1*(emf3[][][i+1] - emf3[k][j][i]) -
1644                            dtodx3*(emf1[k+1][][] - emf1[k][j][i]);
1645        pG->B3i[k][j][i] += dtodx2*(emf1[][j+1][] - emf1[k][j][i]) -
1646                            dtodx1*(emf2[][][i+1] - emf2[k][j][i]);
1647      }
1648      pG->B1i[k][j][ie+1] +=
1649        dtodx3*(emf2[k+1][][ie+1] - emf2[k][j][ie+1]) -
1650        dtodx2*(emf3[][j+1][ie+1] - emf3[k][j][ie+1]);
1651    }
1652    for (i=is; i<=ie; i++) {
1653      pG->B2i[k][je+1][i] +=
1654        dtodx1*(emf3[][je+1][i+1] - emf3[k][je+1][i]) -
1655        dtodx3*(emf1[k+1][je+1][] - emf1[k][je+1][i]);
1656    }
1657  }
1658  for (j=js; j<=je; j++) {
1659    for (i=is; i<=ie; i++) {
1660      pG->B3i[ke+1][j][i] += 
1661        dtodx2*(emf1[ke+1][j+1][] - emf1[ke+1][j][i]) -
1662        dtodx1*(emf2[ke+1][][i+1] - emf2[ke+1][j][i]);
1663    }
1664  }
1665#endif
1666
1667/*--- Step 12a -----------------------------------------------------------------
1668 * Add the gravitational source terms due to a Static Potential, and the
1669 * shearing box source terms (written in terms of a Static Potential).
1670 *    S_{M} = -(\rho)^{n+1/2} Grad(Phi);   S_{E} = -(\rho v)^{n+1/2} Grad{Phi}
1671 */
1672
1673  if (StaticGravPot != NULL){
1674    for (k=ks; k<=ke; k++) {
1675      for (j=js; j<=je; j++) {
1676        for (i=is; i<=ie; i++) {
1677          cc_pos(pG,i,j,k,&x1,&x2,&x3);
1678          phic = (*StaticGravPot)(x1,x2,x3);
1679          phir = (*StaticGravPot)((x1+0.5*pG->dx1),x2,x3);
1680          phil = (*StaticGravPot)((x1-0.5*pG->dx1),x2,x3);
1681
1682          pG->U[k][j][i].M1 -= dtodx1*(phir-phil)*dhalf[k][j][i];
1683#ifndef ISOTHERMAL
1684          pG->U[k][j][i].E -= dtodx1*(x1Flux[k][j][].d*(phic - phil) +
1685                                      x1Flux[k][j][i+1].d*(phir - phic));
1686#endif
1687
1688          phir = (*StaticGravPot)(x1,(x2+0.5*pG->dx2),x3);
1689          phil = (*StaticGravPot)(x1,(x2-0.5*pG->dx2),x3);
1690
1691          pG->U[k][j][i].M2 -= dtodx2*(phir-phil)*dhalf[k][j][i];
1692#ifndef ISOTHERMAL
1693          pG->U[k][j][i].E -= dtodx2*(x2Flux[k][][i].d*(phic - phil) +
1694                                      x2Flux[k][j+1][i].d*(phir - phic));
1695#endif
1696
1697          phir = (*StaticGravPot)(x1,x2,(x3+0.5*pG->dx3));
1698          phil = (*StaticGravPot)(x1,x2,(x3-0.5*pG->dx3));
1699
1700          pG->U[k][j][i].M3 -= dtodx3*(phir-phil)*dhalf[k][j][i];
1701#ifndef ISOTHERMAL
1702          pG->U[k][j][i].E -= dtodx3*(x3Flux[][j][i].d*(phic - phil) +
1703                                      x3Flux[k+1][j][i].d*(phir - phic));
1704#endif
1705        }
1706      }
1707    }
1708  }
1709
1710/*--- Step 12b -----------------------------------------------------------------
1711 * Add gravitational source terms for self-gravity.
1712 * A flux correction using Phi^{n+1} in the main loop is required to make
1713 * the source terms 2nd order: see selfg_flux_correction().
1714 */
1715#ifdef SELF_GRAVITY
1716/* Add fluxes and source terms due to (d/dx1) terms  */
1717
1718  for (k=ks; k<=ke; k++){
1719    for (j=js; j<=je; j++){
1720      for (i=is; i<=ie; i++){
1721        phic = pG->Phi[k][j][i];
1722        phil = 0.5*(pG->Phi[k][j][i-1] + pG->Phi[k][j][]);
1723        phir = 0.5*(pG->Phi[k][j][] + pG->Phi[k][j][i+1]);
1724
1725/* gx, gy and gz centered at L and R x1-faces */
1726        gxl = (pG->Phi[k][j][i-1] - pG->Phi[k][j][])*(dx1i);
1727        gxr = (pG->Phi[k][j][] - pG->Phi[k][j][i+1])*(dx1i);
1728
1729        gyl = 0.25*((pG->Phi[k][j-1][i-1] - pG->Phi[k][j+1][i-1]) +
1730                    (pG->Phi[k][j-1][] - pG->Phi[k][j+1][]) )*(dx2i);
1731        gyr = 0.25*((pG->Phi[k][j-1][] - pG->Phi[k][j+1][]) +
1732                    (pG->Phi[k][j-1][i+1] - pG->Phi[k][j+1][i+1]) )*(dx2i);
1733
1734        gzl = 0.25*((pG->Phi[k-1][j][i-1] - pG->Phi[k+1][j][i-1]) +
1735                    (pG->Phi[k-1][j][] - pG->Phi[k+1][j][]) )*(dx3i);
1736        gzr = 0.25*((pG->Phi[k-1][j][] - pG->Phi[k+1][j][]) +
1737                    (pG->Phi[k-1][j][i+1] - pG->Phi[k+1][j][i+1]) )*(dx3i);
1738
1739/* momentum fluxes in x1.  2nd term is needed only if Jean's swindle used */
1740        flx_m1l = 0.5*(gxl*gxl-gyl*gyl-gzl*gzl)/four_pi_G + grav_mean_rho*phil;
1741        flx_m1r = 0.5*(gxr*gxr-gyr*gyr-gzr*gzr)/four_pi_G + grav_mean_rho*phir;
1742
1743        flx_m2l = gxl*gyl/four_pi_G;
1744        flx_m2r = gxr*gyr/four_pi_G;
1745
1746        flx_m3l = gxl*gzl/four_pi_G;
1747        flx_m3r = gxr*gzr/four_pi_G;
1748
1749/* Update momenta and energy with d/dx1 terms  */
1750        pG->U[k][j][i].M1 -= dtodx1*(flx_m1r - flx_m1l);
1751        pG->U[k][j][i].M2 -= dtodx1*(flx_m2r - flx_m2l);
1752        pG->U[k][j][i].M3 -= dtodx1*(flx_m3r - flx_m3l);
1753#ifdef ADIABATIC
1754        pG->U[k][j][i].E -= dtodx1*(x1Flux[k][j][].d*(phic - phil) +
1755                                    x1Flux[k][j][i+1].d*(phir - phic));
1756#endif /* ADIABATIC */
1757      }
1758    }
1759  }
1760
1761/* Add fluxes and source terms due to (d/dx2) terms  */
1762
1763  for (k=ks; k<=ke; k++){
1764    for (j=js; j<=je; j++){
1765      for (i=is; i<=ie; i++){
1766        phic = pG->Phi[k][j][i];
1767        phil = 0.5*(pG->Phi[k][j-1][i] + pG->Phi[k][][i]);
1768        phir = 0.5*(pG->Phi[k][][i] + pG->Phi[k][j+1][i]);
1769
1770/* gx, gy and gz centered at L and R x2-faces */
1771        gxl = 0.25*((pG->Phi[k][j-1][i-1] - pG->Phi[k][j-1][i+1]) +
1772                    (pG->Phi[k][][i-1] - pG->Phi[k][][i+1]) )*(dx1i);
1773        gxr = 0.25*((pG->Phi[k][][i-1] - pG->Phi[k][][i+1]) +
1774                    (pG->Phi[k][j+1][i-1] - pG->Phi[k][j+1][i+1]) )*(dx1i);
1775
1776        gyl = (pG->Phi[k][j-1][i] - pG->Phi[k][][i])*(dx2i);
1777        gyr = (pG->Phi[k][][i] - pG->Phi[k][j+1][i])*(dx2i);
1778
1779        gzl = 0.25*((pG->Phi[k-1][j-1][i] - pG->Phi[k+1][j-1][i]) +
1780                    (pG->Phi[k-1][][i] - pG->Phi[k+1][][i]) )*(dx3i);
1781        gzr = 0.25*((pG->Phi[k-1][][i] - pG->Phi[k+1][][i]) +
1782                    (pG->Phi[k-1][j+1][i] - pG->Phi[k+1][j+1][i]) )*(dx3i);
1783
1784/* momentum fluxes in x2.  2nd term is needed only if Jean's swindle used */
1785        flx_m1l = gyl*gxl/four_pi_G;
1786        flx_m1r = gyr*gxr/four_pi_G;
1787
1788        flx_m2l = 0.5*(gyl*gyl-gxl*gxl-gzl*gzl)/four_pi_G + grav_mean_rho*phil;
1789        flx_m2r = 0.5*(gyr*gyr-gxr*gxr-gzr*gzr)/four_pi_G + grav_mean_rho*phir;
1790
1791        flx_m3l = gyl*gzl/four_pi_G;
1792        flx_m3r = gyr*gzr/four_pi_G;
1793
1794/* Update momenta and energy with d/dx2 terms  */
1795        pG->U[k][j][i].M1 -= dtodx2*(flx_m1r - flx_m1l);
1796        pG->U[k][j][i].M2 -= dtodx2*(flx_m2r - flx_m2l);
1797        pG->U[k][j][i].M3 -= dtodx2*(flx_m3r - flx_m3l);
1798#ifdef ADIABATIC
1799        pG->U[k][j][i].E -= dtodx2*(x2Flux[k][][i].d*(phic - phil) +
1800                                    x2Flux[k][j+1][i].d*(phir - phic));
1801#endif /* ADIABATIC */
1802      }
1803    }
1804  }
1805
1806/* Add fluxes and source terms due to (d/dx3) terms  */
1807
1808  for (k=ks; k<=ke; k++){
1809    for (j=js; j<=je; j++){
1810      for (i=is; i<=ie; i++){
1811        phic = pG->Phi[k][j][i];
1812        phil = 0.5*(pG->Phi[k-1][j][i] + pG->Phi[][j][i]);
1813        phir = 0.5*(pG->Phi[][j][i] + pG->Phi[k+1][j][i]);
1814
1815/* gx, gy and gz centered at L and R x3-faces */
1816        gxl = 0.25*((pG->Phi[k-1][j][i-1] - pG->Phi[k-1][j][i+1]) +
1817                    (pG->Phi[][j][i-1] - pG->Phi[][j][i+1]) )*(dx1i);
1818        gxr = 0.25*((pG->Phi[][j][i-1] - pG->Phi[][j][i+1]) +
1819                    (pG->Phi[k+1][j][i-1] - pG->Phi[k+1][j][i+1]) )*(dx1i);
1820
1821        gyl = 0.25*((pG->Phi[k-1][j-1][i] - pG->Phi[k-1][j+1][i]) +
1822                    (pG->Phi[][j-1][i] - pG->Phi[][j+1][i]) )*(dx2i);
1823        gyr = 0.25*((pG->Phi[][j-1][i] - pG->Phi[][j+1][i]) +
1824                    (pG->Phi[k+1][j-1][i] - pG->Phi[k+1][j+1][i]) )*(dx2i);
1825
1826        gzl = (pG->Phi[k-1][j][i] - pG->Phi[][j][i])*(dx3i);
1827        gzr = (pG->Phi[][j][i] - pG->Phi[k+1][j][i])*(dx3i);
1828
1829/* momentum fluxes in x3.  2nd term is needed only if Jean's swindle used */
1830        flx_m1l = gzl*gxl/four_pi_G;
1831        flx_m1r = gzr*gxr/four_pi_G;
1832
1833        flx_m2l = gzl*gyl/four_pi_G;
1834        flx_m2r = gzr*gyr/four_pi_G;
1835
1836        flx_m3l = 0.5*(gzl*gzl-gxl*gxl-gyl*gyl)/four_pi_G + grav_mean_rho*phil;
1837        flx_m3r = 0.5*(gzr*gzr-gxr*gxr-gyr*gyr)/four_pi_G + grav_mean_rho*phir;
1838
1839/* Update momenta and energy with d/dx3 terms  */
1840        pG->U[k][j][i].M1 -= dtodx3*(flx_m1r - flx_m1l);
1841        pG->U[k][j][i].M2 -= dtodx3*(flx_m2r - flx_m2l);
1842        pG->U[k][j][i].M3 -= dtodx3*(flx_m3r - flx_m3l);
1843#ifdef ADIABATIC
1844        pG->U[k][j][i].E -= dtodx3*(x3Flux[][j][i].d*(phic - phil) +
1845                                    x3Flux[k+1][j][i].d*(phir - phic));
1846#endif /* ADIABATIC */
1847      }
1848    }
1849  }
1850
1851/* Save mass fluxes in Grid structure for source term correction in main loop */
1852
1853  for (k=ks; k<=ke+1; k++) {
1854    for (j=js; j<=je+1; j++) {
1855      for (i=is; i<=ie+1; i++) {
1856        pG->x1MassFlux[k][j][i] = x1Flux[k][j][i].d;
1857        pG->x2MassFlux[k][j][i] = x2Flux[k][j][i].d;
1858        pG->x3MassFlux[k][j][i] = x3Flux[k][j][i].d;
1859      }
1860    }
1861  }
1862#endif /* SELF_GRAVITY */
1863
1864/*--- Step 13a -----------------------------------------------------------------
1865 * Update cell-centered variables in pG using x1-fluxes
1866 */
1867
1868  for (k=ks; k<=ke; k++) {
1869    for (j=js; j<=je; j++) {
1870      for (i=is; i<=ie; i++) {
1871        pG->U[k][j][i].-= dtodx1*(x1Flux[k][j][i+1].d -x1Flux[k][j][i].d );
1872        pG->U[k][j][i].M1 -= dtodx1*(x1Flux[k][j][i+1].Mx-x1Flux[k][j][i].Mx);
1873        pG->U[k][j][i].M2 -= dtodx1*(x1Flux[k][j][i+1].My-x1Flux[k][j][i].My);
1874        pG->U[k][j][i].M3 -= dtodx1*(x1Flux[k][j][i+1].Mz-x1Flux[k][j][i].Mz);
1875#ifndef ISOTHERMAL
1876        pG->U[k][j][i].-= dtodx1*(x1Flux[k][j][i+1].E -x1Flux[k][j][i].E );
1877#endif /* ISOTHERMAL */
1878#ifdef MHD
1879        pG->U[k][j][i].B2c -= dtodx1*(x1Flux[k][j][i+1].By-x1Flux[k][j][i].By);
1880        pG->U[k][j][i].B3c -= dtodx1*(x1Flux[k][j][i+1].Bz-x1Flux[k][j][i].Bz);
1881#endif /* MHD */
1882#if (NSCALARS > 0)
1883        for (n=0; n<NSCALARS; n++)
1884          pG->U[k][j][i].s[n] -= dtodx1*(x1Flux[k][j][i+1].s[n]
1885                                       - x1Flux[k][j][].s[n]);
1886#endif
1887      }
1888    }
1889  }
1890
1891/*--- Step 13b -----------------------------------------------------------------
1892 * Update cell-centered variables in pG using x2-fluxes
1893 */
1894
1895  for (k=ks; k<=ke; k++) {
1896    for (j=js; j<=je; j++) {
1897      for (i=is; i<=ie; i++) {
1898        pG->U[k][j][i].-= dtodx2*(x2Flux[k][j+1][i].d -x2Flux[k][j][i].d );
1899        pG->U[k][j][i].M1 -= dtodx2*(x2Flux[k][j+1][i].Mz-x2Flux[k][j][i].Mz);
1900        pG->U[k][j][i].M2 -= dtodx2*(x2Flux[k][j+1][i].Mx-x2Flux[k][j][i].Mx);
1901        pG->U[k][j][i].M3 -= dtodx2*(x2Flux[k][j+1][i].My-x2Flux[k][j][i].My);
1902#ifndef ISOTHERMAL
1903        pG->U[k][j][i].E -=dtodx2*(x2Flux[k][j+1][i].E -x2Flux[k][j][i].E );
1904#endif /* ISOTHERMAL */
1905#ifdef MHD
1906        pG->U[k][j][i].B3c -= dtodx2*(x2Flux[k][j+1][i].By-x2Flux[k][j][i].By);
1907        pG->U[k][j][i].B1c -= dtodx2*(x2Flux[k][j+1][i].Bz-x2Flux[k][j][i].Bz);
1908#endif /* MHD */
1909#if (NSCALARS > 0)
1910        for (n=0; n<NSCALARS; n++)
1911          pG->U[k][j][i].s[n] -= dtodx2*(x2Flux[k][j+1][i].s[n]
1912                                       - x2Flux[k][][i].s[n]);
1913#endif
1914      }
1915    }
1916  }
1917
1918/*--- Step 13c -----------------------------------------------------------------
1919 * Update cell-centered variables in pG using x3-fluxes
1920 */
1921
1922  for (k=ks; k<=ke; k++) {
1923    for (j=js; j<=je; j++) {
1924      for (i=is; i<=ie; i++) {
1925        pG->U[k][j][i].-= dtodx3*(x3Flux[k+1][j][i].d -x3Flux[k][j][i].d );
1926        pG->U[k][j][i].M1 -= dtodx3*(x3Flux[k+1][j][i].My-x3Flux[k][j][i].My);
1927        pG->U[k][j][i].M2 -= dtodx3*(x3Flux[k+1][j][i].Mz-x3Flux[k][j][i].Mz);
1928        pG->U[k][j][i].M3 -= dtodx3*(x3Flux[k+1][j][i].Mx-x3Flux[k][j][i].Mx);
1929#ifndef ISOTHERMAL
1930        pG->U[k][j][i].-= dtodx3*(x3Flux[k+1][j][i].E -x3Flux[k][j][i].E );
1931#endif /* ISOTHERMAL */
1932#ifdef MHD
1933        pG->U[k][j][i].B1c -= dtodx3*(x3Flux[k+1][j][i].By-x3Flux[k][j][i].By);
1934        pG->U[k][j][i].B2c -= dtodx3*(x3Flux[k+1][j][i].Bz-x3Flux[k][j][i].Bz);
1935#endif /* MHD */
1936#if (NSCALARS > 0)
1937        for (n=0; n<NSCALARS; n++)
1938          pG->U[k][j][i].s[n] -= dtodx3*(x3Flux[k+1][j][i].s[n]
1939                                       - x3Flux[][j][i].s[n]);
1940#endif
1941      }
1942    }
1943  }
1944
1945/*--- Step 14 ------------------------------------------------------------------
1946 * LAST STEP!
1947 * Set cell centered magnetic fields to average of updated face centered fields.
1948 */
1949
1950#ifdef MHD
1951  for (k=ks; k<=ke; k++) {
1952    for (j=js; j<=je; j++) {
1953      for (i=is; i<=ie; i++) {
1954        pG->U[k][j][i].B1c = 0.5*(pG->B1i[k][j][i]+pG->B1i[k][j][i+1]);
1955        pG->U[k][j][i].B2c = 0.5*(pG->B2i[k][j][i]+pG->B2i[k][j+1][i]);
1956        pG->U[k][j][i].B3c = 0.5*(pG->B3i[k][j][i]+pG->B3i[k+1][j][i]);
1957      }
1958    }
1959  }
1960#endif /* MHD */
1961
1962  return;
1963}
1964
1965/*----------------------------------------------------------------------------*/
1966/* integrate_destruct_3d:  Free temporary integration arrays
1967 */
1968
1969void integrate_destruct_3d(void)
1970{
1971
1972#ifdef MHD
1973  if (emf1    != NULL) free_3d_array(emf1);
1974  if (emf2    != NULL) free_3d_array(emf2);
1975  if (emf3    != NULL) free_3d_array(emf3);
1976  if (emf1_cc != NULL) free_3d_array(emf1_cc);
1977  if (emf2_cc != NULL) free_3d_array(emf2_cc);
1978  if (emf3_cc != NULL) free_3d_array(emf3_cc);
1979#endif /* MHD */
1980#ifdef H_CORRECTION
1981  if (eta1 != NULL) free_3d_array(eta1);
1982  if (eta2 != NULL) free_3d_array(eta2);
1983  if (eta3 != NULL) free_3d_array(eta3);
1984#endif /* H_CORRECTION */
1985
1986  if (Bxc != NULL) free(Bxc);
1987  if (Bxi != NULL) free(Bxi);
1988  if (B1_x1Face != NULL) free_3d_array(B1_x1Face);
1989  if (B2_x2Face != NULL) free_3d_array(B2_x2Face);
1990  if (B3_x3Face != NULL) free_3d_array(B3_x3Face);
1991
1992  if (U1d      != NULL) free(U1d);
1993  if (W        != NULL) free(W);
1994  if (Wl       != NULL) free(Wl);
1995  if (Wr       != NULL) free(Wr);
1996
1997  if (Ul_x1Face != NULL) free_3d_array(Ul_x1Face);
1998  if (Ur_x1Face != NULL) free_3d_array(Ur_x1Face);
1999  if (Ul_x2Face != NULL) free_3d_array(Ul_x2Face);
2000  if (Ur_x2Face != NULL) free_3d_array(Ur_x2Face);
2001  if (Ul_x3Face != NULL) free_3d_array(Ul_x3Face);
2002  if (Ur_x3Face != NULL) free_3d_array(Ur_x3Face);
2003  if (x1Flux    != NULL) free_3d_array(x1Flux);
2004  if (x2Flux    != NULL) free_3d_array(x2Flux);
2005  if (x3Flux    != NULL) free_3d_array(x3Flux);
2006  if (dhalf     != NULL) free_3d_array(dhalf);
2007
2008  return;
2009}
2010
2011/*----------------------------------------------------------------------------*/
2012/* integrate_init_3d: Allocate temporary integration arrays
2013*/
2014
2015void integrate_init_3d(int nx1, int nx2, int nx3)
2016{
2017  int nmax;
2018  int Nx1 = nx1 + 2*nghost;
2019  int Nx2 = nx2 + 2*nghost;
2020  int Nx3 = nx3 + 2*nghost;
2021  nmax = MAX(MAX(Nx1,Nx2),Nx3);
2022
2023#ifdef MHD
2024  if ((emf1 = (Real***)calloc_3d_array(Nx3, Nx2, Nx1, sizeof(Real))) == NULL)
2025    goto on_error;
2026  if ((emf2 = (Real***)calloc_3d_array(Nx3, Nx2, Nx1, sizeof(Real))) == NULL)
2027    goto on_error;
2028  if ((emf3 = (Real***)calloc_3d_array(Nx3, Nx2, Nx1, sizeof(Real))) == NULL)
2029    goto on_error;
2030
2031  if ((emf1_cc = (Real***)calloc_3d_array(Nx3, Nx2, Nx1, sizeof(Real))) == NULL)
2032    goto on_error;
2033  if ((emf2_cc = (Real***)calloc_3d_array(Nx3, Nx2, Nx1, sizeof(Real))) == NULL)
2034    goto on_error;
2035  if ((emf3_cc = (Real***)calloc_3d_array(Nx3, Nx2, Nx1, sizeof(Real))) == NULL)
2036    goto on_error;
2037#endif /* MHD */
2038#ifdef H_CORRECTION
2039  if ((eta1 = (Real***)calloc_3d_array(Nx3, Nx2, Nx1, sizeof(Real))) == NULL)
2040    goto on_error;
2041  if ((eta2 = (Real***)calloc_3d_array(Nx3, Nx2, Nx1, sizeof(Real))) == NULL)
2042    goto on_error;
2043  if ((eta3 = (Real***)calloc_3d_array(Nx3, Nx2, Nx1, sizeof(Real))) == NULL)
2044    goto on_error;
2045#endif /* H_CORRECTION */
2046
2047  if ((Bxc = (Real*)malloc(nmax*sizeof(Real))) == NULL) goto on_error;
2048  if ((Bxi = (Real*)malloc(nmax*sizeof(Real))) == NULL) goto on_error;
2049
2050  if ((B1_x1Face = (Real***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Real))) == NULL)
2051    goto on_error;
2052  if ((B2_x2Face = (Real***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Real))) == NULL)
2053    goto on_error;
2054  if ((B3_x3Face = (Real***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Real))) == NULL)
2055    goto on_error;
2056
2057  if ((U1d =      (Cons1D*)malloc(nmax*sizeof(Cons1D))) == NULL) goto on_error;
2058  if ((W   =      (Prim1D*)malloc(nmax*sizeof(Prim1D))) == NULL) goto on_error;
2059  if ((Wl  =      (Prim1D*)malloc(nmax*sizeof(Prim1D))) == NULL) goto on_error;
2060  if ((Wr  =      (Prim1D*)malloc(nmax*sizeof(Prim1D))) == NULL) goto on_error;
2061
2062  if ((Ul_x1Face = (Cons1D***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Cons1D)))
2063    == NULL) goto on_error;
2064  if ((Ur_x1Face = (Cons1D***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Cons1D)))
2065    == NULL) goto on_error;
2066  if ((Ul_x2Face = (Cons1D***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Cons1D)))
2067    == NULL) goto on_error;
2068  if ((Ur_x2Face = (Cons1D***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Cons1D)))
2069    == NULL) goto on_error;
2070  if ((Ul_x3Face = (Cons1D***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Cons1D))) 
2071    == NULL) goto on_error;
2072  if ((Ur_x3Face = (Cons1D***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Cons1D))) 
2073    == NULL) goto on_error;
2074  if ((x1Flux    = (Cons1D***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Cons1D))) 
2075    == NULL) goto on_error;
2076  if ((x2Flux    = (Cons1D***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Cons1D))) 
2077    == NULL) goto on_error;
2078  if ((x3Flux    = (Cons1D***)calloc_3d_array(Nx3,Nx2,Nx1, sizeof(Cons1D))) 
2079    == NULL) goto on_error;
2080
2081
2082#if defined MHD
2083  if ((dhalf = (Real***)calloc_3d_array(Nx3, Nx2, Nx1, sizeof(Real))) == NULL)
2084    goto on_error;
2085#else
2086  if(StaticGravPot != NULL){
2087    if ((dhalf = (Real***)calloc_3d_array(Nx3, Nx2, Nx1, sizeof(Real))) == NULL)
2088      goto on_error;
2089  }
2090#endif
2091
2092  return;
2093
2094  on_error:
2095  integrate_destruct();
2096  ath_error("[integrate_init]: malloc returned a NULL pointer\n");
2097}
2098
2099
2100/*=========================== PRIVATE FUNCTIONS ==============================*/
2101
2102/*----------------------------------------------------------------------------*/
2103/* integrate_emf1_corner
2104 * integrate_emf2_corner
2105 * integrate_emf3_corner
2106 *   Integrates face centered B-fluxes to compute corner EMFs.  Note:
2107 *   x1Flux.By = VxBy - BxVy = v1*b2-b1*v2 = -EMFZ
2108 *   x1Flux.Bz = VxBz - BxVz = v1*b3-b1*v3 = EMFY
2109 *   x2Flux.By = VxBy - BxVy = v2*b3-b2*v3 = -EMFX
2110 *   x2Flux.Bz = VxBz - BxVz = v2*b1-b2*v1 = EMFZ
2111 *   x3Flux.By = VxBy - BxVy = v3*b1-b3*v1 = -EMFY
2112 *   x3Flux.Bz = VxBz - BxVz = v3*b2-b3*v2 = EMFX
2113 */
2114
2115#ifdef MHD
2116static void integrate_emf1_corner(const Grid *pG)
2117{
2118  int i, is = pG->is, ie = pG->ie;
2119  int j, js = pG->js, je = pG->je;
2120  int k, ks = pG->ks, ke = pG->ke;
2121  Real de1_l2, de1_r2, de1_l3, de1_r3;
2122
2123  for (k=ks-1; k<=ke+2; k++) {
2124    for (j=js-1; j<=je+2; j++) {
2125      for (i=is-2; i<=ie+2; i++) {
2126/* NOTE: The x2-Flux of By is -E1. */
2127/*       The x3-Flux of Bz is +E1. */
2128        if (x2Flux[k-1][j][i].d > 0.0)
2129          de1_l3 = x3Flux[k][j-1][i].Bz - emf1_cc[k-1][j-1][i];
2130        else if (x2Flux[k-1][j][i].d < 0.0)
2131          de1_l3 = x3Flux[k][j][i].Bz - emf1_cc[k-1][j][i];
2132        else {
2133          de1_l3 = 0.5*(x3Flux[k][j-1][i].Bz - emf1_cc[k-1][j-1][i] +
2134                        x3Flux[k][][i].Bz - emf1_cc[k-1][][i] );
2135        }
2136
2137        if (x2Flux[k][j][i].d > 0.0)
2138          de1_r3 = x3Flux[k][j-1][i].Bz - emf1_cc[k][j-1][i];
2139        else if (x2Flux[k][j][i].d < 0.0)
2140          de1_r3 = x3Flux[k][j][i].Bz - emf1_cc[k][j][i];
2141        else {
2142          de1_r3 = 0.5*(x3Flux[k][j-1][i].Bz - emf1_cc[k][j-1][i] +
2143                        x3Flux[k][][i].Bz - emf1_cc[k][][i] );
2144        }
2145
2146        if (x3Flux[k][j-1][i].d > 0.0)
2147          de1_l2 = -x2Flux[k-1][j][i].By - emf1_cc[k-1][j-1][i];
2148        else if (x3Flux[k][j-1][i].d < 0.0)
2149          de1_l2 = -x2Flux[k][j][i].By - emf1_cc[k][j-1][i];
2150        else {
2151          de1_l2 = 0.5*(-x2Flux[k-1][j][i].By - emf1_cc[k-1][j-1][i]
2152                        -x2Flux[][j][i].By - emf1_cc[][j-1][i] );
2153        }
2154
2155        if (x3Flux[k][j][i].d > 0.0)
2156          de1_r2 = -x2Flux[k-1][j][i].By - emf1_cc[k-1][j][i];
2157        else if (x3Flux[k][j][i].d < 0.0)
2158          de1_r2 = -x2Flux[k][j][i].By - emf1_cc[k][j][i];
2159        else {
2160          de1_r2 = 0.5*(-x2Flux[k-1][j][i].By - emf1_cc[k-1][j][i]
2161                        -x2Flux[][j][i].By - emf1_cc[][j][i] );
2162        }
2163
2164        emf1[k][j][i] = 0.25*(  x3Flux[k][j][i].Bz + x3Flux[k][j-1][i].Bz
2165                              - x2Flux[k][j][i].By - x2Flux[k-1][j][i].By
2166                              + de1_l2 + de1_r2 + de1_l3 + de1_r3);
2167      }
2168    }
2169  }
2170
2171  return;
2172}
2173
2174
2175static void integrate_emf2_corner(const Grid *pG)
2176{
2177  int i, is = pG->is, ie = pG->ie;
2178  int j, js = pG->js, je = pG->je;
2179  int k, ks = pG->ks, ke = pG->ke;
2180  Real de2_l1, de2_r1, de2_l3, de2_r3;
2181
2182  for (k=ks-1; k<=ke+2; k++) {
2183    for (j=js-2; j<=je+2; j++) {
2184      for (i=is-1; i<=ie+2; i++) {
2185/* NOTE: The x1-Flux of Bz is +E2. */
2186/*       The x3-Flux of By is -E2. */
2187        if (x1Flux[k-1][j][i].d > 0.0)
2188          de2_l3 = -x3Flux[k][j][i-1].By - emf2_cc[k-1][j][i-1];
2189        else if (x1Flux[k-1][j][i].d < 0.0)
2190          de2_l3 = -x3Flux[k][j][i].By - emf2_cc[k-1][j][i];
2191        else {
2192          de2_l3 = 0.5*(-x3Flux[k][j][i-1].By - emf2_cc[k-1][j][i-1] 
2193                        -x3Flux[k][j][].By - emf2_cc[k-1][j][] );
2194        }
2195
2196        if (x1Flux[k][j][i].d > 0.0)
2197          de2_r3 = -x3Flux[k][j][i-1].By - emf2_cc[k][j][i-1];
2198        else if (x1Flux[k][j][i].d < 0.0)
2199          de2_r3 = -x3Flux[k][j][i].By - emf2_cc[k][j][i];
2200        else {
2201          de2_r3 = 0.5*(-x3Flux[k][j][i-1].By - emf2_cc[k][j][i-1] 
2202                        -x3Flux[k][j][].By - emf2_cc[k][j][] );
2203        }
2204
2205        if (x3Flux[k][j][i-1].d > 0.0)
2206          de2_l1 = x1Flux[k-1][j][i].Bz - emf2_cc[k-1][j][i-1];
2207        else if (x3Flux[k][j][i-1].d < 0.0)
2208          de2_l1 = x1Flux[k][j][i].Bz - emf2_cc[k][j][i-1];
2209        else {
2210          de2_l1 = 0.5*(x1Flux[k-1][j][i].Bz - emf2_cc[k-1][j][i-1] +
2211                        x1Flux[][j][i].Bz - emf2_cc[][j][i-1] );
2212        }
2213
2214        if (x3Flux[k][j][i].d > 0.0)
2215          de2_r1 = x1Flux[k-1][j][i].Bz - emf2_cc[k-1][j][i];
2216        else if (x3Flux[k][j][i].d < 0.0)
2217          de2_r1 = x1Flux[k][j][i].Bz - emf2_cc[k][j][i];
2218        else {
2219          de2_r1 = 0.5*(x1Flux[k-1][j][i].Bz - emf2_cc[k-1][j][i] +
2220                        x1Flux[][j][i].Bz - emf2_cc[k-1][j][i] );
2221        }
2222
2223        emf2[k][j][i] = 0.25*(  x1Flux[k][j][i].Bz + x1Flux[k-1][j][].Bz
2224                              - x3Flux[k][j][i].By - x3Flux[][j][i-1].By
2225                              + de2_l1 + de2_r1 + de2_l3 + de2_r3);
2226      }
2227    }
2228  }
2229
2230  return;
2231}
2232
2233static void integrate_emf3_corner(const Grid *pG)
2234{
2235  int i, is = pG->is, ie = pG->ie;
2236  int j, js = pG->js, je = pG->je;
2237  int k, ks = pG->ks, ke = pG->ke;
2238  Real de3_l1, de3_r1, de3_l2, de3_r2;
2239
2240  for (k=ks-2; k<=ke+2; k++) {
2241    for (j=js-1; j<=je+2; j++) {
2242      for (i=is-1; i<=ie+2; i++) {
2243/* NOTE: The x1-Flux of By is -E3. */
2244/*       The x2-Flux of Bx is +E3. */
2245        if (x1Flux[k][j-1][i].d > 0.0)
2246          de3_l2 = x2Flux[k][j][i-1].Bz - emf3_cc[k][j-1][i-1];
2247        else if (x1Flux[k][j-1][i].d < 0.0)
2248          de3_l2 = x2Flux[k][j][i].Bz - emf3_cc[k][j-1][i];
2249        else {
2250          de3_l2 = 0.5*(x2Flux[k][j][i-1].Bz - emf3_cc[k][j-1][i-1] + 
2251                        x2Flux[k][j][].Bz - emf3_cc[k][j-1][] );
2252        }
2253
2254        if (x1Flux[k][j][i].d > 0.0)
2255          de3_r2 = x2Flux[k][j][i-1].Bz - emf3_cc[k][j][i-1];
2256        else if (x1Flux[k][j][i].d < 0.0)
2257          de3_r2 = x2Flux[k][j][i].Bz - emf3_cc[k][j][i];
2258        else {
2259          de3_r2 = 0.5*(x2Flux[k][j][i-1].Bz - emf3_cc[k][j][i-1] + 
2260                        x2Flux[k][j][].Bz - emf3_cc[k][j][] );
2261        }
2262
2263        if (x2Flux[k][j][i-1].d > 0.0)
2264          de3_l1 = -x1Flux[k][j-1][i].By - emf3_cc[k][j-1][i-1];
2265        else if (x2Flux[k][j][i-1].d < 0.0)
2266          de3_l1 = -x1Flux[k][j][i].By - emf3_cc[k][j][i-1];
2267        else {
2268          de3_l1 = 0.5*(-x1Flux[k][j-1][i].By - emf3_cc[k][j-1][i-1]
2269                        -x1Flux[k][][i].By - emf3_cc[k][][i-1] );
2270        }
2271
2272        if (x2Flux[k][j][i].d > 0.0)
2273          de3_r1 = -x1Flux[k][j-1][i].By - emf3_cc[k][j-1][i];
2274        else if (x2Flux[k][j][i].d < 0.0)
2275          de3_r1 = -x1Flux[k][j][i].By - emf3_cc[k][j][i];
2276        else {
2277          de3_r1 = 0.5*(-x1Flux[k][j-1][i].By - emf3_cc[k][j-1][i]
2278                        -x1Flux[k][][i].By - emf3_cc[k][][i] );
2279        }
2280
2281        emf3[k][j][i] = 0.25*(  x2Flux[k][][i-1].Bz + x2Flux[k][j][i].Bz
2282                              - x1Flux[k][j-1][].By - x1Flux[k][j][i].By
2283                              + de3_l1 + de3_r1 + de3_l2 + de3_r2);
2284      }
2285    }
2286  }
2287
2288  return;
2289}
2290#endif /* MHD */
Note: See TracBrowser for help on using the repository browser.