Cleanup to ascii mesh converter

src/core_landice/Registry.xml

@@ -1110,11 +1110,11 @@ is the value of that variable from the *previous* time level!
 	       persistence="scratch"
                />
 	  <var name="slopeCellX" type="real" dimensions="nCells" units="none"
-	       descriptiom="x-component of slope on cell centers"
+	       description="x-component of slope on cell centers"
 	       persistence="scratch"
 	       />
 	  <var name="slopeCellY" type="real" dimensions="nCells" units="none"
-	       descriptiom="y-component of slope on cell centers"
+	       description="y-component of slope on cell centers"
 	       persistence="scratch"
 	       />
 	  <var name="vertexIndices" type="integer" dimensions="nVertices" units="none"

src/core_landice/mode_forward/Interface_velocity_solver.cpp

@@ -2197,6 +2197,7 @@ int prismType(long long int const* prismVertexMpasIds, int& minIndex)
     write_ascii_mesh_field(thicknessUncertaintyData, "thickness_uncertainty");
 
     write_ascii_mesh_field(elevationData, "surface_height");
+    write_ascii_mesh_field(bedTopographyData, "bed_topography");
 
     write_ascii_mesh_field(betaData, "basal_friction");
 
@@ -2233,10 +2234,10 @@ int prismType(long long int const* prismVertexMpasIds, int& minIndex)
     outfile.open ("surface_velocity.ascii", std::ios::out | std::ios::trunc);
     if (outfile.is_open()) {
        outfile << nVertices << " " << 2 << "\n";  //number of vertices, number of components per vertex
-       for(int i = 0; i<nVertices; ++i) {
+       for(int i = 0; i<nVertices; ++i)
          outfile << observedVeloXData[i] << "\n";
+       for(int i = 0; i<nVertices; ++i)
          outfile << observedVeloYData[i] << "\n";
-       }
        outfile.close();
     }
     else {

testing_and_setup/compass/landice/EISMINT2/visualize_output_EISMINT2.py

@@ -76,7 +76,7 @@ def contourMPAS(field, contour_levs=None):
   if contour_levs == None:
      im = plt.contour(xi, yi, zi)
   else:
-     im = plt.contour(xi, yi, zi, contour_levs)
+     im = plt.contour(xi, yi, zi, contour_levs, cmap=plt.cm.jet)
 
   #plt.scatter(xCell, yCell, c=temperature[timelev,:,-1], s=100, vmin=zi.min(), vmax=zi.max())  # to see the raw data on top
   plt.colorbar(im)
@@ -149,11 +149,15 @@ def contourMPAS(field, contour_levs=None):
 #markershape='h'
 
 # print ice locations with gray hexagons
-iceIndices = np.where(thickness[timelev,:]>0.0)[0]
+iceIndices = np.where(thickness[timelev,:]>10.0)[0]
 plt.scatter(xCell[iceIndices], yCell[iceIndices], markersize, (0.8, 0.8, 0.8), marker=markershape, edgecolors='none') # print ice locations with gray hexagons
 
 # add contours of ice temperature over the top
-contourMPAS(basalTemperature[timelev,:], np.linspace(240.0, 275.0, 8))
+basalTemp = basalTemperature[timelev,:]
+# fill places below dynamic limit with non-ice value of 273.15
+basalTemp[np.where(thickness[timelev,:]<10.0)] = 273.15
+#plt.scatter(xCell, yCell, markersize, basalTemp[:], marker=markershape, edgecolors='none') # print ice locations with gray hexagons; plt.plot()
+contourMPAS(basalTemp, np.linspace(240.0, 275.0, 8))
 
 plt.axis('equal')
 plt.title('Modeled basal temperature (K) \n at time ' + netCDF4.chartostring(xtime)[timelev].strip() )

testing_and_setup/compass/landice/hydro-radial/plot_hydro-radial_profile.py

@@ -215,6 +215,9 @@
 plt.legend()
 plt.xlabel('delta x (m)')
 plt.ylabel('error in W (m)')
+print "avg W err=", Werr[hasice].mean()
+print "max W err=", Werr[hasice].max()
+
 
 ax = fig.add_subplot(2,1,2)
 plt.plot(dx, Perr[hasice].mean()/1.0e5, 's', label='avg P err')
@@ -224,6 +227,8 @@
 plt.legend()
 plt.xlabel('delta x (m)')
 plt.ylabel('error in P (bar)')
+print "avg P err=", Perr[hasice].mean()/1.0e5
+print "max P err=", Perr[hasice].max()/1.0e5
 
 
 print "plotting complete"