FEAT: limit changes near the pole

src/solver_advection.cpp

@@ -318,7 +318,10 @@ void advect(Grid &grid,
 
   arma_mat area, xWidth, yWidth, geometry;
 
+  // These are all needed by the solver:
   neutrals.calc_mass_density();
+  neutrals.calc_mean_major_mass();
+  neutrals.calc_specific_heat();
 
   arma_mat t_to_e;
 
@@ -507,16 +510,37 @@ void advect(Grid &grid,
 
     neutrals.velocity_vcgc[0].slice(iAlt) = xVel;
     neutrals.velocity_vcgc[1].slice(iAlt) = yVel;
-    temp = totalE / rho - 0.5 * (xVel % xVel + yVel % yVel);
-
-    neutrals.temperature_scgc.slice(iAlt) = temp / t_to_e;
+    temp = (totalE / rho - 0.5 * (xVel % xVel + yVel % yVel)) / t_to_e;
+    precision_t fac, dm, dp;
+
+    for (int64_t j = nGCs; j < nY - nGCs; j++) {
+      for (int64_t i = nGCs; i < nX - nGCs; i++) {
+	fac = 1.0;
+	if (cos(grid.geoLat_scgc(i,j,iAlt)) < 0.2) {
+	  fac = fac * (0.2 - cos(grid.geoLat_scgc(i,j,iAlt)));
+	}
+	dm = (1.0 - fac) * neutrals.temperature_scgc(i,j,iAlt);
+	dp = (1.0 + fac) * neutrals.temperature_scgc(i,j,iAlt);
+	if (temp(i,j) < dm) temp(i,j) = dm;
+	if (temp(i,j) > dp) temp(i,j) = dp;
+	neutrals.temperature_scgc(i,j,iAlt) = temp(i,j);
+
+	dm = (1.0 - fac) * neutrals.rho_scgc(i,j,iAlt);
+	dp = (1.0 + fac) * neutrals.rho_scgc(i,j,iAlt);
+	if (rho(i,j) < dm) rho(i,j) = dm;
+	if (rho(i,j) > dp) rho(i,j) = dp;
+	neutrals.rho_scgc(i,j,iAlt) = rho(i,j);
+      }
+    }
     if (report.test_verbose(3) && iAlt == 8) {
       std::cout << "end t : " << neutrals.temperature_scgc.slice(iAlt).min() << " " << neutrals.temperature_scgc.slice(iAlt).max() << "\n";
       std::cout << "end temp : " << temp.min() << " " << temp.max() << "\n";
       std::cout << "end xVel : " << xVel.min() << " " << xVel.max() << "\n";
       std::cout << "end yVel : " << yVel.min() << " " << yVel.max() << "\n";
     }
   }
+  neutrals.calc_density_from_mass_concentration();
+
   report.exit(function);
   return;
-}
+}