@@ -19,62 +19,63 @@ void Neutrals::add_sources(Times time, Planets planet, Grid grid) {
19
19
precision_t dt = time.get_dt ();
20
20
21
21
heating_sources_total = heating_euv_scgc
22
- + heating_chemical_scgc
23
- + heating_ion_friction_scgc
24
- // + heating_ion_heat_transfer_scgc
25
- - O_cool_scgc
26
- - NO_cool_scgc;
27
-
22
+ + heating_chemical_scgc
23
+ + heating_ion_friction_scgc
24
+ // + heating_ion_heat_transfer_scgc
25
+ - O_cool_scgc
26
+ - NO_cool_scgc;
27
+
28
28
// Solve the laplace equations using the source terms,
29
29
// updating the neutral temperature:
30
30
update_temperature (grid, time);
31
31
32
32
int64_t iDir, iSpec, iSpecies;
33
33
double tSim = time.get_simulation_time ();
34
-
34
+
35
35
// Horizontal winds use bulk winds:
36
- if (input.get_use_coriolis ())
36
+ if (input.get_use_coriolis ())
37
37
acc_coriolis = coriolis (velocity_vcgc, planet.get_omega (), grid.geoLat_scgc );
38
38
39
- /*
40
- // Vertical winds use species winds:
41
- for (iSpec = 0; iSpec < nSpeciesAdvect; iSpec++) {
42
- // Pick out the advected neutral species:
43
- species_chars & advected_neutral = species[species_to_advect[iSpec]];
44
-
45
- iDir = 2;
46
- // update velocities based on acceleration:
47
- // reduce neutral friction until solver is added
48
- advected_neutral.velocity_vcgc[iDir] =
49
- advected_neutral.velocity_vcgc[iDir] +
50
- dt * (ramp * grid.cent_acc_vcgc[iDir] +
51
- ramp * acc_coriolis[iDir] +
52
- advected_neutral.acc_neutral_friction[iDir] / 4.0 +
53
- advected_neutral.acc_ion_drag[iDir] +
54
- advected_neutral.acc_eddy);
55
- }
39
+ /*
40
+ // Vertical winds use species winds:
41
+ for (iSpec = 0; iSpec < nSpeciesAdvect; iSpec++) {
42
+ // Pick out the advected neutral species:
43
+ species_chars & advected_neutral = species[species_to_advect[iSpec]];
56
44
57
- calc_mass_density();
58
- // Calculate bulk vertical winds:
59
- velocity_vcgc[2].zeros();
60
- for (iSpecies = 0; iSpecies < nSpecies; iSpecies++)
61
- if (species[iSpecies].DoAdvect) {
62
- velocity_vcgc[2] = velocity_vcgc[2] +
63
- species[iSpecies].mass * species[iSpecies].density_scgc %
64
- species[iSpecies].velocity_vcgc[2] / rho_scgc;
45
+ iDir = 2;
46
+ // update velocities based on acceleration:
47
+ // reduce neutral friction until solver is added
48
+ advected_neutral.velocity_vcgc[iDir] =
49
+ advected_neutral.velocity_vcgc[iDir] +
50
+ dt * (ramp * grid.cent_acc_vcgc[iDir] +
51
+ ramp * acc_coriolis[iDir] +
52
+ advected_neutral.acc_neutral_friction[iDir] / 4.0 +
53
+ advected_neutral.acc_ion_drag[iDir] +
54
+ advected_neutral.acc_eddy);
65
55
}
66
-
67
- */
56
+
57
+ calc_mass_density();
58
+ // Calculate bulk vertical winds:
59
+ velocity_vcgc[2].zeros();
60
+ for (iSpecies = 0; iSpecies < nSpecies; iSpecies++)
61
+ if (species[iSpecies].DoAdvect) {
62
+ velocity_vcgc[2] = velocity_vcgc[2] +
63
+ species[iSpecies].mass * species[iSpecies].density_scgc %
64
+ species[iSpecies].velocity_vcgc[2] / rho_scgc;
65
+ }
66
+
67
+ */
68
68
69
69
// Add Velocity sources to bulk winds:
70
70
for (iDir = 0 ; iDir < 2 ; iDir++) {
71
71
velocity_vcgc[iDir] =
72
72
velocity_vcgc[iDir] + dt * (
73
- grid.cent_acc_vcgc [iDir] +
74
- acc_coriolis[iDir] +
75
- acc_ion_collisions[iDir]);
73
+ grid.cent_acc_vcgc [iDir] +
74
+ acc_coriolis[iDir] +
75
+ acc_ion_collisions[iDir]);
76
76
acc_sources_total[iDir].zeros ();
77
77
}
78
+
78
79
// Apply Viscosity:
79
80
update_horizontal_velocity (grid, time);
80
81
@@ -83,17 +84,17 @@ void Neutrals::add_sources(Times time, Planets planet, Grid grid) {
83
84
for (iDir = 0 ; iDir < 2 ; iDir++)
84
85
species[iSpecies].velocity_vcgc [iDir] = velocity_vcgc[iDir];
85
86
86
- /*
87
+ /*
87
88
// If we only consider the bulk winds in the horizontal direction:
88
89
if (input.get_advection_neutrals_bulkwinds()) {
89
90
// Calculate Coriolis:
90
- if (input.get_use_coriolis())
91
+ if (input.get_use_coriolis())
91
92
acc_coriolis = coriolis(velocity_vcgc, planet.get_omega(), grid.geoLat_scgc);
92
93
// Add Velocity sources to bulk winds:
93
94
for (int iDir = 0; iDir < 3; iDir++) {
94
95
velocity_vcgc[iDir] = velocity_vcgc[iDir] + dt * (
95
96
grid.cent_acc_vcgc[iDir] +
96
- acc_coriolis[iDir] +
97
+ acc_coriolis[iDir] +
97
98
acc_ion_collisions[iDir]);
98
99
acc_sources_total[iDir].zeros();
99
100
}
@@ -104,18 +105,18 @@ void Neutrals::add_sources(Times time, Planets planet, Grid grid) {
104
105
// Pick out the advected neutral species:
105
106
species_chars & advected_neutral = species[species_to_advect[iSpec]];
106
107
// Calculate Coriolis:
107
- if (input.get_use_coriolis())
108
- acc_coriolis = coriolis(advected_neutral.velocity_vcgc,
109
- planet.get_omega(),
108
+ if (input.get_use_coriolis())
109
+ acc_coriolis = coriolis(advected_neutral.velocity_vcgc,
110
+ planet.get_omega(),
110
111
grid.geoLat_scgc);
111
112
112
113
for (int iDir = 0; iDir < 2; iDir++) {
113
114
// update velocities based on acceleration:
114
115
// reduce neutral friction until solver is added
115
116
advected_neutral.velocity_vcgc[iDir] =
116
117
advected_neutral.velocity_vcgc[iDir] +
117
- dt * (grid.cent_acc_vcgc[iDir] +
118
- acc_coriolis[iDir] +
118
+ dt * (grid.cent_acc_vcgc[iDir] +
119
+ acc_coriolis[iDir] +
119
120
advected_neutral.acc_neutral_friction[iDir] / 4.0 +
120
121
advected_neutral.acc_ion_drag[iDir]);
121
122
// eddy acceleration is only in the vertical direction:
0 commit comments