Merge branch 'master' into calibration_updates

updating from master

Author: Benjamin-Maglio

README.md

@@ -29,6 +29,13 @@ command line.
 
 More details are available in the [User Guide - Basic Model Setup and Run](https://uaf-arctic-eco-modeling.github.io/dvm-dos-tem/examples_and_tutorials/basic_model_run.html#basic-model-setup-and-run).
 
+> [!TIP]
+> These instructions should work for users with [Podman](https://podman.io)
+> rather than Docker, but several additional steps may be necessary to ensure
+> that Podman can run with non-root access. See this
+> [issue](https://github.com/uaf-arctic-eco-modeling/dvm-dos-tem/issues/776)
+> for more information.
+
 The settings shown here will use the input data that is shipped with the code in the
 repository's ``demo-data/`` directory. The run will output a single variable (GPP), 
 and will run for 2 pixels.

config/output_spec.csv

@@ -18,14 +18,14 @@ DEADC,Standing dead C,g/m2,,,invalid,invalid,invalid,invalid,double,
 DEADN,Standing dead N,g/m2/time,,,invalid,invalid,invalid,invalid,double,
 DEEPC,Amorphous SOM C,g/m2,,,invalid,invalid,invalid,invalid,double,
 DEEPDZ,Amorphous SOM horizon thickness,m,,invalid,invalid,invalid,invalid,invalid,double,
-DRIVINGNIRR,Input driving NIRR,W/m2,invalid,invalid,,invalid,invalid,invalid,
+DRIVINGNIRR,Input driving NIRR,W/m2,invalid,invalid,,invalid,invalid,invalid,float,
 DRIVINGRAINFALL,Input driving precip data after being split from snowfall,mm,,,,invalid,invalid,invalid,float,
 DRIVINGSNOWFALL,Input driving precip data after being split from rainfall,mm,,,,invalid,invalid,invalid,float,
 DRIVINGTAIR,Input driving air temperature data,degree_C,invalid,invalid,,invalid,invalid,invalid,float,
 DRIVINGVAPO,Input driving vapor pressure,hPa,invalid,invalid,,invalid,invalid,invalid,float,
 DWDC,Dead woody debris C,g/m2/time,,,invalid,invalid,invalid,invalid,double,
 DWDN,Dead woody debris N,g/m2/time,,,invalid,invalid,invalid,invalid,double,
-EET,Actual ET,mm/m2/time,,,,,invalid,invalid,
+EET,Actual ET,mm/m2/time,,,,invalid,invalid,invalid,
 FRONTSDEPTH,Depth of fronts,mm,,,,invalid,invalid,,double,
 FRONTSTYPE,Front types,,,,,invalid,invalid,,int,
 GPP,GPP,g/m2/time,y,,invalid,,,invalid,double,
@@ -45,12 +45,12 @@ LAI,LAI,m2/m2,,,invalid,,invalid,invalid,double,
 LAYERDEPTH,Layer depth from the surface,m,,,invalid,invalid,invalid,forced,double,
 LAYERDZ,Thickness of layer,m,,,invalid,invalid,invalid,forced,double,
 LAYERTYPE,0:moss 1:shlw 2:deep 3:mineral,,,,invalid,invalid,invalid,forced,int,
-LTRFALC,Litterfall C,g/m2/time,,,invalid,,,invalid,double,
-LTRFALN,Litterfall N,g/m2/time,,,invalid,,,invalid,double,
+LFNVC,Litterfall C from non-vascular PFTs,g/m2/time,,,invalid,invalid,invalid,invalid,double,
+LFNVN,Litterfall N from non-vascular PFTs,g/m2/time,,,invalid,invalid,invalid,invalid,double,
+LFVC,Litterfall C from vascular PFTs,g/m2/time,,,invalid,,,invalid,double,
+LFVN,Litterfall N from vascular PFTs,g/m2/time,,,invalid,,,invalid,double,
 LWCLAYER,LWC by layer,m3/m3,,,invalid,invalid,invalid,forced,double,
 MINEC,Mineral SOM C,g/m2,,,invalid,invalid,invalid,invalid,double,
-MOSSDEATHC,Moss death carbon,g/m2,,,invalid,invalid,invalid,invalid,double,
-MOSSDEATHN,Moss death nitrogen,g/m2,,,invalid,invalid,invalid,invalid,double,
 MOSSDZ,Moss thickness,m,,invalid,invalid,invalid,invalid,invalid,double,
 NDRAIN,N losses from drainage (AVLN),g/m2/time,invalid,invalid,invalid,invalid,invalid,invalid,double,
 NETNMIN,Soil net N mineralization,g/m2/time,,,invalid,invalid,invalid,,double,
@@ -64,7 +64,7 @@ NUPTAKEST,Structural N uptake by plants,g/m2/time,,,invalid,,,invalid,double,
 ORGN,Total soil organic N,g/m2,,,invalid,invalid,invalid,,double,
 PERCOLATION,Percolation by layer,,invalid,invalid,,invalid,invalid,forced,double,
 PERMAFROST,Permafrost (1 or 0),,,invalid,invalid,invalid,invalid,invalid,int,
-PET,Potential ET,mm/m2/time,,,,,invalid,invalid,double,
+PET,Potential ET,mm/m2/time,,,,invalid,invalid,invalid,double,
 QDRAINAGE,Water drainage quotient (~ratio),,,,invalid,invalid,invalid,invalid,double,
 QDRAINLAYER,Qdrain per layer,mm,invalid,invalid,,invalid,invalid,forced,double,
 QINFILTRATION,Water infiltration quotient (~ratio),,,,invalid,invalid,invalid,invalid,double,
@@ -107,7 +107,9 @@ TMINEC,Temperature bottom mineral,degree_C,,,,invalid,invalid,invalid,double,
 TRANSPIRATION,Transpiration,mm/day,,,,invalid,invalid,invalid,double,
 TSHLW,Temperature in fibrous horizon,degree_C,,,,invalid,invalid,invalid,double,
 VEGC,Total veg. biomass C,g/m2,,,invalid,,,invalid,double,
-VEGN,Total veg. biomass N,g/m2,,,invalid,,,invalid,double,
+VEGNTOT,Total veg. biomass N (structural and labile),g/m2,,,invalid,invalid,invalid,invalid,double,
+VEGNLAB,Veg. biomass labile N,g/m2,,,invalid,invalid,invalid,invalid,double,
+VEGNSTR,Veg. biomass structural N,g/m2,,,invalid,,,invalid,double,
 VWCDEEP,VWC in amorphous horizon,m3/m3,,,,invalid,invalid,invalid,double,
 VWCLAYER,VWC by layer,m3/m3,,,invalid,invalid,invalid,forced,double,
 VWCMINEA,VWC top mineral,m3/m3,,,,invalid,invalid,invalid,double,

docs_src/sphinx/source/examples_and_tutorials/basic_model_run.rst

@@ -20,6 +20,15 @@ and visualization, see the :ref:`Plotting Example
   For users who have Docker and Git installed and some familiarity with the 
   command line.
 
+  .. note::
+
+    These instructions should work for users with `Podman <https://podman.io>`_
+    rather than Docker, but several additional steps may be necessary to ensure
+    that Podman can run with non-root access. See this
+    `issue <https://github.com/uaf-arctic-eco-modeling/dvm-dos-tem/issues/776>`_
+    for more information.
+
+
   The settings shown here will use the input data that is shipped with the code
   in the repository's ``demo-data/`` directory. The run will output a single
   variable (GPP), and will run for 2 pixels.
@@ -148,6 +157,14 @@ The dependencies for a piece of software are isolated inside the container. This
 isolation allows software with conflicting dependencies to run on the host
 system.
 
+.. hint::
+
+  These instructions should work for users with `Podman <https://podman.io>`_
+  rather than Docker, but several additional steps may be necessary to ensure
+  that Podman can run with non-root access. See this
+  `issue <https://github.com/uaf-arctic-eco-modeling/dvm-dos-tem/issues/776>`_
+  for more information.
+
 There are other ways you can get the environment necessary for running,
 compiling and developing ``dvmdostem``, such as a native installation, or using
 a virtual machine (e.g. with Vagrant and Virtual Box or VMWare). Each path to a

include/BgcData.h

@@ -72,6 +72,7 @@ class BgcData {
   void soil_beginOfMonth();
   void soil_endOfMonth(const int currmind);
 
+  void soil_updatePools(const int currmind);
 private:
 
 

include/EnvData.h

@@ -144,7 +144,7 @@ class EnvData {
   void grnd_beginOfDay();
 
   // accumulating/averaging monthly variables at the end of day
-  void atm_endOfDay(const int & dinm);
+  void atm_endOfDay(const int & dinm, const int & dayidx);
   void veg_endOfDay(const int & dinm);
   void grnd_endOfDay(const int & dinm, const int & doy);
 

include/OutputHolder.h

@@ -39,8 +39,8 @@ class OutputHolder{
   std::vector<double> deepc_for_output;
   std::vector<double> dwdc_for_output;
   std::vector<double> minec_for_output;
-  std::vector<double> mossdeathc_for_output;
-  std::vector<double> mossdeathn_for_output;
+  std::vector<double> lfnvc_for_output; //non-vascular litterfall C
+  std::vector<double> lfnvn_for_output; //non-vascular litterfall N
   std::vector<double> mossdz_for_output;
   std::vector<double> qdrain_for_output;
   std::vector<double> qinfil_for_output;
@@ -105,10 +105,10 @@ class OutputHolder{
   std::vector<std::array<double, NUM_PFT>> gpp_pft_for_output;
   std::vector<std::array<std::array<double, NUM_PFT>, NUM_PFT_PART>> gpp_for_output; 
 
-  std::vector<double> ltrfalc_tot_for_output;
-  std::vector<std::array<double, NUM_PFT_PART>> ltrfalc_part_for_output;
-  std::vector<std::array<double, NUM_PFT>> ltrfalc_pft_for_output;
-  std::vector<std::array<std::array<double, NUM_PFT>, NUM_PFT_PART>> ltrfalc_for_output; 
+  std::vector<double> lfvc_tot_for_output;
+  std::vector<std::array<double, NUM_PFT_PART>> lfvc_part_for_output;
+  std::vector<std::array<double, NUM_PFT>> lfvc_pft_for_output;
+  std::vector<std::array<std::array<double, NUM_PFT>, NUM_PFT_PART>> lfvc_for_output; 
 
   std::vector<double> ingpp_tot_for_output;
   std::vector<std::array<double, NUM_PFT_PART>> ingpp_part_for_output;

scripts/simulation_comparison_report.py

@@ -737,43 +737,48 @@ def vegdynamic(simpath,simlist,vlist,sclist,oname,stage):
 
 os.mkdir(os.path.join(POD,'results'))
 
+stage = 'eq'
+
 print("Generating timeseries flux plots...")
-VARlist=['RHSOM','GPP','NPP','LTRFALC','RG','RM']
-ts_flux(POD,PODlist,VARlist,scenariolist,colorlist,widthlist,'carbon','Yearly Carbon Flux Time Series','eq')
+VARlist=['RHSOM','GPP','NPP','LFVC','RG','RM']
+ts_flux(POD,PODlist,VARlist,scenariolist,colorlist,widthlist,'carbon','Yearly Carbon Flux Time Series',stage)
 
-VARlist=['NUPTAKEST','NUPTAKELAB','NRESORB','NIMMOB','NETNMIN','LTRFALN']
-ts_flux(POD,PODlist,VARlist,scenariolist,colorlist,widthlist,'nitrogen','Yearly Nitrogen Flux Time Series', 'eq')
+VARlist=['NUPTAKEST','NUPTAKELAB','NRESORB','NIMMOB','NETNMIN','LFVN']
+ts_flux(POD,PODlist,VARlist,scenariolist,colorlist,widthlist,'nitrogen','Yearly Nitrogen Flux Time Series', stage)
 
 VARlist=['BURNSOIL2AIRC','BURNVEG2AIRC','BURNVEG2DEADC','RHDWD', 'BURNAIR2SOILN','BURNSOIL2AIRN','BURNVEG2AIRN','BURNVEG2DEADN']
-ts_flux(POD,PODlist,VARlist,scenariolist,colorlist,widthlist, 'fire', 'Yearly Wildfire Flux Time Series', 'eq')
+ts_flux(POD,PODlist,VARlist,scenariolist,colorlist,widthlist, 'fire', 'Yearly Wildfire Flux Time Series', stage)
 
 print("Generating ts stock plots...")
 VARlist=['SHLWC','DEEPC','MINEC','SOMRAWC','SOMA','SOMPR','SOMCR','VEGC']
-ts_stock(POD,PODlist,VARlist,scenariolist,colorlist,widthlist, 'carbon','Yearly Carbon Stock Time Series', 'eq')
+ts_stock(POD,PODlist,VARlist,scenariolist,colorlist,widthlist, 'carbon','Yearly Carbon Stock Time Series', stage)
 
-VARlist=['AVLN','ORGN','VEGN','NETNMIN','LTRFALN']
-ts_stock(POD,PODlist,VARlist,scenariolist,colorlist,widthlist,'nitrogen','Yearly Nitrogen Stock Time Series','eq')
+VARlist=['AVLN','ORGN','VEGNSTR','NETNMIN','LFVN']
+ts_stock(POD,PODlist,VARlist,scenariolist,colorlist,widthlist,'nitrogen','Yearly Nitrogen Stock Time Series',stage)
 
 VARlist=['DEADC','DWDC','DEADN','DWDN']
-ts_stock(POD,PODlist,VARlist,scenariolist,colorlist,widthlist,'wildfire','Yearly Burned C, N Stock in Time series', 'eq')
+ts_stock(POD,PODlist,VARlist,scenariolist,colorlist,widthlist,'wildfire','Yearly Burned C, N Stock in Time series', stage)
 
 print("Generating seasonality plots...")
 VARlist=['GPP','RHSOM','LAI']
-seasonality(POD,PODlist,VARlist,scenariolist,colorlist,'Seas_Bio','eq')
+seasonality(POD,PODlist,VARlist,scenariolist,colorlist,'Seas_Bio',stage)
 
-VARlist=['SNOWTHICK','EET','PET','TRANSPIRATION','WATERTAB']
-seasonality(POD,PODlist,VARlist,scenariolist,colorlist,'Seas_Env', 'eq')
+#Commented out when EET by PFT was disabled. This will need to
+# be modified to work with ecosystem EET or wait until EET by
+# PFT is fixed.
+#VARlist=['SNOWTHICK','EET','PET','TRANSPIRATION','WATERTAB']
+#seasonality(POD,PODlist,VARlist,scenariolist,colorlist,'Seas_Env', stage)
 
 print("Generating soil CN plots...")
-soilcnprofile(POD,PODlist,scenariolist,'eq')
+soilcnprofile(POD,PODlist,scenariolist,stage)
 
 print("Generating soil env plots...")
 VARlist=['TLAYER','VWCLAYER']
-soilenvprofile(POD,PODlist,VARlist,scenariolist,'Profile','eq')
+soilenvprofile(POD,PODlist,VARlist,scenariolist,'Profile',stage)
 
 print("Generating vegetation dynamic plots...")
-VARlist=['VEGC','VEGN']
-vegdynamic(POD,PODlist,VARlist,scenariolist,'Veg_Dyn','eq')
+VARlist=['VEGC','VEGNSTR']
+vegdynamic(POD,PODlist,VARlist,scenariolist,'Veg_Dyn',stage)
 
 
 

scripts/tests/m2y_analysis_hg.py

@@ -0,0 +1,187 @@
+#!/usr/bin/env python
+
+#### Run command
+### 1. Prepare the output_spec.csv files for four scenarios (yearly/ecosystem, monthly/ecosystem, monthly/pft/layer, monthly/pft/part)
+### 2. cd /data/workflows/TILE_testing/
+### 3. cp ./run-mask.nc ./run-mask_original.nc
+### 4. ncap2 -O -h -s'run(:,:) = 0' ./run-mask.nc ./run-mask.nc
+### 5. ncap2 -O -h -s'run(0,0) = 1' ./run-mask.nc ./run-mask.nc
+### 6. cp output_spec_scenario.csv output_spec.csv
+### 7. /work/dvmdostem -l err -f ./config/config.js -p 100 -e 200 -s 10 -t 3
+### 8. mv -r ./output ./output_scenario
+### 9. mkdir ./output
+
+#### Analysis of the results
+
+import os
+import xarray as xr
+import pandas as pd
+
+tolerance = 0.000001
+outdir = 'temp_outputs/reco_output_debugging/3_outputs_for_analysis/'
+sc_list = ('y', 'm', 'm_pft_layer', 'm_pft_part')
+
+
+
+
+
+
+### Month to year test
+
+
+varlist = ('AVLN','DEEPC','EET','GPP','HKDEEP','LAI','LFVC','LFVN','MINEC','NETNMIN','NIMMOB','NPP','NRESORB','NUPTAKELAB','NUPTAKEST','ORGN','PET','QDRAINAGE','QINFILTRATION','QRUNOFF','RAINFALL','RG','RHSOM','RM','SHLWC','SNOWFALL','SNOWTHICK','SOMA','SOMCR','SOMPR','SOMRAWC','SWE','TCDEEP','TDEEP','TRANSPIRATION','TSHLW','VEGC','VEGNSTR','VWCDEEP','WATERTAB')
+statlist = ('stock','stock','flux','flux','mean','mean','flux','flux','stock','flux','flux','flux','flux','flux','flux','stock','flux','flux','flux','flux','flux','flux','flux','flux','stock','flux','stock','stock','stock','stock','stock','mean','mean','mean','flux','mean','stock','stock','mean','mean')
+
+sc = 'y'
+yearly = pd.DataFrame(columns = ['year'])
+for i, var in enumerate(varlist):
+  print(var)
+  dd = xr.open_dataset(os.path.join(outdir,'output_'+str(sc), var + '_yearly_tr.nc')).isel(x=0, y=0)
+  dd['time'] = dd.indexes['time'].to_datetimeindex()
+  dd = dd.to_dataframe()
+  dd.reset_index(inplace=True)
+  dd['year'] = dd['time'].dt.year
+  dd = dd.drop(columns=['time'])
+  dd = dd.drop(columns=['albers_conical_equal_area'])
+  yearly = pd.merge(yearly, dd, on=['year'], how='outer')
+
+
+sc = 'm'
+monthly = pd.DataFrame(columns = ['year','month'])
+m2y = pd.DataFrame(columns = ['year'])
+for i, var in enumerate(varlist):
+  dd = xr.open_dataset(os.path.join(outdir,'output_'+str(sc), var + '_monthly_tr.nc')).isel(x=0, y=0)
+  dd['time'] = dd.indexes['time'].to_datetimeindex()
+  dd = dd.to_dataframe()
+  dd.reset_index(inplace=True)
+  dd['year'] = dd['time'].dt.year
+  dd['month'] = dd['time'].dt.month
+  dd = dd.drop(columns=['time'])
+  dd = dd.drop(columns=['albers_conical_equal_area'])
+  monthly = pd.merge(monthly, dd, on=['year','month'], how='outer')
+  if statlist[i] == 'flux':
+    dy = pd.DataFrame(dd.groupby('year')[var].sum())
+  if statlist[i] == 'mean':
+    dy = pd.DataFrame(dd.groupby('year')[var].mean())
+  if statlist[i] == 'stock':
+    dy = dd[dd['month'] == 12]
+    dy = dy.drop(columns=['month'])
+  dy.reset_index(inplace=True)
+  if 'index' in dy.columns:
+    dy = dy.drop(columns=['index'])
+#  dy = dy.rename(columns={var: var + '_m2y'}) 
+  m2y = pd.merge(m2y, dy, on=['year'], how='outer')
+
+
+test = pd.DataFrame((m2y-yearly).mean())
+test.reset_index(inplace=True)
+test = test.rename(columns={'index': 'varname', 0:'diffmean'}) 
+
+test[abs(test['diffmean']) > tolerance]
+
+print("Monthly to yearly")
+print(varlist)
+print(test)
+
+
+###  pft to ecosystem test
+
+varlist = ('EET','GPP','LAI','LFVC','LFVN','NPP','NRESORB','NUPTAKELAB','NUPTAKEST','PET','RG','RM','VEGC','VEGNSTR')
+
+sc = 'm_pft_layer'
+pft = pd.DataFrame(columns = ['year','month','pft'])
+for i, var in enumerate(varlist):
+  dd = xr.open_dataset(os.path.join(outdir,'output_'+str(sc), var + '_monthly_tr.nc')).isel(x=0, y=0)
+  dd['time'] = dd.indexes['time'].to_datetimeindex()
+  dd = dd.to_dataframe()
+  dd.reset_index(inplace=True)
+  dd['year'] = dd['time'].dt.year
+  dd['month'] = dd['time'].dt.month
+  dd = dd.drop(columns=['time'])
+  dd = dd.drop(columns=['albers_conical_equal_area'])
+  pft = pd.merge(pft, dd, on=['year','month','pft'], how='outer')
+
+pft2eco = pft.groupby(['year','month']).sum()
+pft2eco.reset_index(inplace=True)
+pft2eco = pft2eco.drop(columns=['pft'])
+
+test = pd.DataFrame((pft2eco-monthly).mean())
+test.reset_index(inplace=True)
+test = test.rename(columns={'index': 'varname', 0:'diffmean'}) 
+
+test[abs(test['diffmean']) > tolerance]
+
+#breakpoint()
+
+print("PFT to Ecosystem")
+print(varlist)
+print(test)
+
+
+
+
+###  part to pft test
+
+varlist = ('GPP','LFVC','LFVN','NPP','NRESORB','NUPTAKEST','RG','RM','VEGC','VEGNSTR')
+
+sc = 'm_pft_part'
+part = pd.DataFrame(columns = ['year','month','pft','pftpart'])
+for i, var in enumerate(varlist):
+  dd = xr.open_dataset(os.path.join(outdir,'output_'+str(sc), var + '_monthly_tr.nc')).isel(x=0, y=0)
+  dd['time'] = dd.indexes['time'].to_datetimeindex()
+  dd = dd.to_dataframe()
+  dd.reset_index(inplace=True)
+  dd['year'] = dd['time'].dt.year
+  dd['month'] = dd['time'].dt.month
+  dd = dd.drop(columns=['time'])
+  dd = dd.drop(columns=['albers_conical_equal_area'])
+  part = pd.merge(part, dd, on=['year','month','pft','pftpart'], how='outer')
+
+part2pft = part.groupby(['year','month','pft']).sum()
+part2pft.reset_index(inplace=True)
+part2pft = part2pft.drop(columns=['pftpart'])
+
+test = pd.DataFrame((part2pft-pft).mean())
+test.reset_index(inplace=True)
+test = test.rename(columns={'index': 'varname', 0:'diffmean'}) 
+
+test[abs(test['diffmean']) > tolerance]
+
+#breakpoint()
+print("Compartment to PFT")
+print(varlist)
+print(test)
+
+
+
+###  layer to column test
+
+varlist = ('AVLN','NETNMIN','NIMMOB','ORGN','RHSOM','SOMA','SOMCR','SOMPR','SOMRAWC')
+
+sc = 'm_pft_layer'
+layer = pd.DataFrame(columns = ['year','month','layer'])
+for i, var in enumerate(varlist):
+  dd = xr.open_dataset(os.path.join(outdir,'output_'+str(sc), var + '_monthly_tr.nc')).isel(x=0, y=0)
+  dd['time'] = dd.indexes['time'].to_datetimeindex()
+  dd = dd.to_dataframe()
+  dd.reset_index(inplace=True)
+  dd['year'] = dd['time'].dt.year
+  dd['month'] = dd['time'].dt.month
+  dd = dd.drop(columns=['time'])
+  dd = dd.drop(columns=['albers_conical_equal_area'])
+  layer = pd.merge(layer, dd, on=['year','month','layer'], how='outer')
+
+layer2column = layer.groupby(['year','month']).sum()
+layer2column.reset_index(inplace=True)
+layer2column = layer2column.drop(columns=['layer'])
+
+test = pd.DataFrame((layer2column-monthly).mean())
+test.reset_index(inplace=True)
+test = test.rename(columns={'index': 'varname', 0:'diffmean'}) 
+
+test[abs(test['diffmean']) > tolerance]
+
+#breakpoint()
+print("Layer to Ecosystem")
+print(varlist)
+print(test)