Issue425 Example runs should run

core/src/RectGridIO.cpp

@@ -21,6 +21,7 @@
 #include <ncFile.h>
 #include <ncVar.h>
 
+#include <algorithm>
 #include <list>
 #include <map>
 #include <string>
@@ -37,6 +38,8 @@ void dimensionSetter(
     for (size_t d = 0; d < nDims; ++d) {
         dims[d] = dataGroup.getVar(fieldName).getDim(d).getSize();
     }
+    // The dimensions in the netCDF are in the reverse order compared to ModelArray
+    std::reverse(dims.begin(), dims.end());
     // A special case for Type::Z: use NZLevels for the third dimension
     if (type == ModelArray::Type::Z)
         dims[2] = NZLevels::get();
@@ -72,6 +75,7 @@ ModelState RectGridIO::getModelState(const std::string& filePath)
     state.data[ticeName] = ModelArray::ZField();
     std::vector<size_t> startVector = { 0, 0, 0 };
     std::vector<size_t> zArrayDims = ModelArray::dimensions(ModelArray::Type::Z);
+    std::reverse(zArrayDims.begin(), zArrayDims.end());
     dataGroup.getVar(ticeName).getVar(startVector, zArrayDims, &state.data[ticeName][0]);
 
     ncFile.close();
@@ -100,9 +104,9 @@ void RectGridIO::dumpModelState(const ModelState& state, const ModelMetadata& me
     // data or the parent group
     netCDF::NcDim xDim = dataGroup.addDim(dimensionNames[0], nx);
     netCDF::NcDim yDim = dataGroup.addDim(dimensionNames[1], ny);
-    std::vector<netCDF::NcDim> dims2 = { xDim, yDim };
+    std::vector<netCDF::NcDim> dims2 = { yDim, xDim };
     netCDF::NcDim zDim = dataGroup.addDim(dimensionNames[2], nz);
-    std::vector<netCDF::NcDim> dims3 = { xDim, yDim, zDim };
+    std::vector<netCDF::NcDim> dims3 = { zDim, yDim, xDim };
 
     for (const auto entry : state.data) {
         const std::string& name = entry.first;

physics/src/IceGrowth.cpp

@@ -155,6 +155,7 @@ void IceGrowth::initializeThicknesses()
     overElements(std::bind(&IceGrowth::initializeThicknessesElement, this, std::placeholders::_1,
                      std::placeholders::_2),
         TimestepTime());
+    iVertical->initialiseTice();
 }
 
 // Divide by ice concentration to go from cell-averaged to ice-averaged values,

physics/src/modules/ConstantAtmosphereBoundary.cpp

@@ -26,6 +26,7 @@ void ConstantAtmosphereBoundary::setData(const ModelState::DataMap& ms)
     evap = 0; // somehow...
     uwind = 0;
     vwind = 0;
+    penSW = 0.;
 }
 
 void ConstantAtmosphereBoundary::update(const TimestepTime& tst)

physics/src/modules/include/IAtmosphereBoundary.hpp

@@ -28,6 +28,7 @@ class IAtmosphereBoundary : public ModelComponent {
     , emp(ModelArray::Type::H)
     , uwind(ModelArray::Type::U)
     , vwind(ModelArray::Type::V)
+    , penSW(ModelArray::Type::H)
     {
         m_couplingArrays.resize(static_cast<size_t>(CouplingFields::COUNT));
         m_couplingArrays[static_cast<size_t>(CouplingFields::SUBL)] = &subl;

physics/src/modules/include/IIceThermodynamics.hpp

@@ -40,6 +40,8 @@ class IIceThermodynamics : public ModelComponent {
      */
     virtual void update(const TimestepTime& tsTime) = 0;
 
+    virtual void initialiseTice() { tice = tice0; }
+
     inline static std::string getKappaSConfigKey() { return "thermo.ks"; }
 
     virtual size_t getNZLevels() const = 0;

run/era5_topaz4_maker.py

@@ -185,7 +185,7 @@ def rotate_velocities(u, v, angle):
     if args.prefix is not None:
         filepfx = args.prefix + "."
     else:
-        filepfx = "."
+        filepfx = ""
 
     # read a grid spec (from a restart file)
     root = netCDF4.Dataset(args.file, "r", format = "NETCDF4")
@@ -231,6 +231,7 @@ def rotate_velocities(u, v, angle):
     # ERA5 data
 
     era5_out_file = f"{filepfx}ERA5_{args.start}_{args.stop}.nc"
+    print(f"Writing ERA5 data to {era5_out_file}")
     era_root = netCDF4.Dataset(era5_out_file, "w", format="NETCDF4")
     structgrp = era_root.createGroup("structure")
     structgrp.type = target_structure
@@ -288,8 +289,8 @@ def rotate_velocities(u, v, angle):
                 data[target_t_index, :, :] = time_data
         else:
             # Also handle the wind components along with the wind speed
-            u_var = datagrp.createVariable("u", "f8", ("time", "x", "y"))
-            v_var = datagrp.createVariable("v", "f8", ("time", "x", "y"))
+            u_var = datagrp.createVariable("u", "f8", timefield_dims)
+            v_var = datagrp.createVariable("v", "f8", timefield_dims)
             for target_t_index in range(len(unix_times_e)):
                 # get the source data
                 u_file = netCDF4.Dataset(era5_source_file_name("u10", unix_times_e[target_t_index]), "r")
@@ -337,6 +338,7 @@ def rotate_velocities(u, v, angle):
 
     topaz_out_file = f"{filepfx}TOPAZ4_{args.start}_{args.stop}.nc"
     topaz_root = netCDF4.Dataset(topaz_out_file, "w", format="NETCDF4")
+    print(f"Writing TOPAZ4 data to {topaz_out_file}")
     structgrp = topaz_root.createGroup("structure")
     structgrp.type = target_structure
 

run/init_topaz128x128.py

@@ -207,10 +207,10 @@ def topaz4_interpolate(target_lon_deg, target_lat_deg, data, lat_array):
     sst[:, :] = sst_data * noice + mu * sss_data * isice
 
     # Ice starts at rest
-    u = datagrp.createVariable("u", "f8", ("x", "y"))
+    u = datagrp.createVariable("u", "f8", hfield_dims)
     u[:, :] = 0
 
-    v = datagrp.createVariable("v", "f8", ("x", "y"))
+    v = datagrp.createVariable("v", "f8", hfield_dims)
     v[:, :] = 0
 
     root.close()