Merge pull request microhh#324 from bartvstratum/main_rho_restart

Write base state `rho` to its own restart file.

Author: Chiil

include/fields.h

@@ -117,6 +117,8 @@ class Fields
 
         void save(int);
         void load(int);
+        void save_rhoref();
+        void load_rhoref();
 
         TF check_momentum();
         TF check_tke();

python/microhh_tools.py

@@ -495,6 +495,7 @@ def restart_pre(origin, timestr):
     fnames += glob.glob('../' + origin + '/grid.0000000')
     fnames += glob.glob('../' + origin + '/fftwplan.0000000')
     fnames += glob.glob('../' + origin + '/thermo_basestate.0000000')
+    fnames += glob.glob('../' + origin + '/rhoref.0000000')
     fnames += glob.glob('../' + origin + '/*.' + timestr)
     for file in fnames:
         shutil.copy(file, '.')

src/fields.cxx

@@ -1208,6 +1208,10 @@ void Fields<TF>::create_stats(Stats<TF>& stats)
         // (Turbulence) Kinetic Energy
         stats.add_prof("ke" , "Kinetic energy" , "m2 s-2", "z", group_name);
         stats.add_prof("tke", "Turbulent kinetic energy" , "m2 s-2", "z", group_name);
+
+        // Base state density.
+        stats.add_fixed_prof("rhoref",  "Full level density dynamic core", "kg m-3", "z" , group_name, this->rhoref);
+        stats.add_fixed_prof("rhorefh", "Half level density dynamic core", "kg m-3", "zh", group_name, this->rhorefh);
     }
 
     // Add time series of scalar surface values
@@ -1245,6 +1249,7 @@ void Fields<TF>::save(int n)
     auto tmp1 = get_tmp();
     auto tmp2 = get_tmp();
 
+    // Save all prognostic fields.
     int nerror = 0;
     for (auto& f : ap)
     {
@@ -1336,6 +1341,81 @@ void Fields<TF>::load(int n)
         throw std::runtime_error("Error loading fields");
 }
 
+template<typename TF>
+void Fields<TF>::save_rhoref()
+{
+    auto& gd = grid.get_grid_data();
+
+    char filename[256];
+    int nerror = 0;
+    std::sprintf(filename, "%s.%07d", "rhoref", 0);
+    master.print_message("Saving \"%s\" ... ", filename);
+
+    if (master.get_mpiid() == 0)
+    {
+        FILE *pFile;
+        pFile = fopen(filename, "ab");
+
+        if (pFile == NULL)
+        {
+            master.print_message("FAILED\n");
+            nerror++;
+        }
+        else
+        {
+            master.print_message("OK\n");
+
+            fwrite(&rhoref [gd.kstart], sizeof(TF), gd.kmax, pFile);
+            fwrite(&rhorefh[gd.kstart], sizeof(TF), gd.kmax+1, pFile);
+
+            fclose(pFile);
+        }
+    }
+
+    master.sum(&nerror, 1);
+    if (nerror)
+        throw std::runtime_error("Error in writing basestate density.");
+}
+
+
+template<typename TF>
+void Fields<TF>::load_rhoref()
+{
+    auto& gd = grid.get_grid_data();
+
+    char filename[256];
+    int nerror = 0;
+    std::sprintf(filename, "%s.%07d", "rhoref", 0);
+    master.print_message("Loading \"%s\" ... ", filename);
+
+    if (master.get_mpiid() == 0)
+    {
+        FILE* pFile;
+        pFile = fopen(filename, "rb");
+        if (pFile == NULL)
+        {
+            master.print_message("FAILED\n");
+            ++nerror;
+        }
+        else
+        {
+            master.print_message("OK\n");
+
+            if (fread(&rhoref[gd.kstart], sizeof(TF), gd.ktot  , pFile) != (unsigned)gd.ktot )
+                ++nerror;
+            if (fread(&rhorefh[gd.kstart], sizeof(TF), gd.ktot+1, pFile) != (unsigned)gd.ktot+1)
+                ++nerror;
+            fclose(pFile);
+        }
+    }
+
+    if (nerror)
+        throw std::runtime_error("Error in loading basestate density.");
+
+    master.broadcast(&rhoref[gd.kstart], gd.ktot  );
+    master.broadcast(&rhorefh[gd.kstart], gd.ktot+1);
+}
+
 #ifndef USECUDA
 template<typename TF>
 TF Fields<TF>::check_momentum()

src/model.cxx

@@ -252,6 +252,8 @@ void Model<TF>::load()
 
     // Load the fields, and create the field statistics
     fields->load(timeloop->get_iotime());
+    fields->load_rhoref();
+
     fields->create_stats(*stats);
     fields->create_column(*column);
 
@@ -315,6 +317,7 @@ void Model<TF>::save()
 
     thermo->create_basestate(*input, *input_nc, *timeloop);
     thermo->save(timeloop->get_iotime());
+    fields->save_rhoref();
 
     boundary->create_cold_start(*input_nc);
     boundary->save(timeloop->get_iotime(), *thermo);

src/thermo_dry.cxx

@@ -662,11 +662,11 @@ void Thermo_dry<TF>::create_stats(Stats<TF>& stats)
     if (stats.get_switch())
     {
         bs_stats = bs;
+
         // Add base state profiles to statistics
-        stats.add_fixed_prof("rhoref",  "Full level basic state density",  "kg m-3", "z",  group_name, fields.rhoref );
-        stats.add_fixed_prof("rhorefh", "Half level basic state density",  "kg m-3", "zh", group_name, fields.rhorefh);
         stats.add_fixed_prof("thref",   "Full level basic state potential temperature", "K", "z" , group_name, bs_stats.thref);
         stats.add_fixed_prof("threfh",  "Half level basic state potential temperature", "K", "zh", group_name, bs_stats.thref);
+
         if (bs_stats.swbasestate == Basestate_type::anelastic)
         {
             stats.add_fixed_prof("phydro",  "Full level hydrostatic pressure", "Pa", "z" , group_name, bs_stats.pref );

src/thermo_moist.cxx

@@ -1916,24 +1916,23 @@ void Thermo_moist<TF>::create_stats(Stats<TF>& stats)
     // Add variables to the statistics
     if (stats.get_switch())
     {
-        /* Add fixed base-state density and temperature profiles. Density should probably be in fields (?), but
-           there the statistics are initialized before thermo->create() is called */
-        stats.add_fixed_prof("rhoref",  "Full level basic state density", "kg m-3", "z" , group_name, bs.rhoref );
-        stats.add_fixed_prof("rhorefh", "Half level basic state density", "kg m-3", "zh", group_name, bs.rhorefh);
+        // Do we also want this time dependent with `swupdatebasestate`?
         stats.add_fixed_prof("thvref", "Full level basic state virtual potential temperature", "K", "z" , group_name, bs.thvref);
         stats.add_fixed_prof("thvrefh", "Half level basic state virtual potential temperature", "K", "zh", group_name, bs.thvrefh);
 
         if (bs_stats.swupdatebasestate)
         {
             stats.add_prof("phydro", "Full level hydrostatic pressure", "Pa", "z" , group_name);
             stats.add_prof("phydroh", "Half level hydrostatic pressure", "Pa", "zh", group_name);
-            stats.add_prof("rho",  "Full level density", "kg m-3", "z" , group_name);
-            stats.add_prof("rhoh", "Half level density", "kg m-3", "zh", group_name);
+            stats.add_prof("rho",  "Full level thermodynamic density", "kg m-3", "z" , group_name);
+            stats.add_prof("rhoh", "Half level thermodynamic density", "kg m-3", "zh", group_name);
         }
         else
         {
             stats.add_fixed_prof("pydroh",  "Full level hydrostatic pressure", "Pa", "z" , group_name, bs.pref);
             stats.add_fixed_prof("phydroh", "Half level hydrostatic pressure", "Pa", "zh", group_name, bs.prefh);
+            stats.add_fixed_prof("rho",  "Full level thermodynamic density", "kg m-3", "z" , group_name, bs.rhoref);
+            stats.add_fixed_prof("rhoh", "Half level thermodynamic density", "kg m-3", "zh", group_name, bs.rhorefh);
         }
 
         auto thv = fields.get_tmp();