Cosipy educational NBs and interface

.gitignore

@@ -1,3 +1,4 @@
+.ipynb_checkpoints
 .idea/
 venv/
 htmlcov/

cosipy/cpkernel/cosipy_core.py

@@ -19,9 +19,10 @@
 
 from cosipy.cpkernel.init import init_snowpack, load_snowpack
 from cosipy.cpkernel.io import IOClass
+from cosipy.utils.options import read_opt
 
 
-def cosipy_core(DATA, indY, indX, GRID_RESTART=None, stake_names=None, stake_data=None):
+def cosipy_core(DATA, indY, indX, GRID_RESTART=None, stake_names=None, stake_data=None, opt_dict=None):
     """ Cosipy core function, which perform the calculations on one core.
 
     Params
@@ -49,6 +50,9 @@ def cosipy_core(DATA, indY, indX, GRID_RESTART=None, stake_names=None, stake_dat
         max_layers = int(DATA.max_layers.values)
         z = float(DATA.ZLVL.values)
 
+    # Replace imported variables with content of the opt_dict. If it's empty
+    # nothing happens.
+    read_opt(opt_dict, globals())
     # Local variables
     nt = len(DATA.time.values)         #accessing DATA is expensive 		
     _RRR = np.full(nt, np.nan)
@@ -96,14 +100,14 @@ def cosipy_core(DATA, indY, indX, GRID_RESTART=None, stake_names=None, stake_dat
     # Initialize snowpack or load restart grid
     #--------------------------------------------
     if GRID_RESTART is None:
-        GRID = init_snowpack(DATA)
+        GRID = init_snowpack(DATA, opt_dict)
     else:
         GRID = load_snowpack(GRID_RESTART)
 
     # Create the local output datasets if not coupled
     RESTART = None
     if not WRF_X_CSPY:
-        IO = IOClass(DATA)
+        IO = IOClass(DATA, opt_dict)
         RESTART = IO.create_local_restart_dataset()
 
     # hours since the last snowfall (albedo module)
@@ -222,12 +226,12 @@ def cosipy_core(DATA, indY, indX, GRID_RESTART=None, stake_names=None, stake_dat
         #--------------------------------------------
         # Calculate albedo and roughness length changes if first layer is snow
         #--------------------------------------------
-        alpha = updateAlbedo(GRID)
+        alpha = updateAlbedo(GRID, opt_dict)
 
         #--------------------------------------------
         # Update roughness length
         #--------------------------------------------
-        z0 = updateRoughness(GRID)
+        z0 = updateRoughness(GRID, opt_dict)
 
         #--------------------------------------------
         # Surface Energy Balance
@@ -237,7 +241,7 @@ def cosipy_core(DATA, indY, indX, GRID_RESTART=None, stake_names=None, stake_dat
 
         # Penetrating SW radiation and subsurface melt
         if SWnet > 0.0:
-            subsurface_melt, G_penetrating = penetrating_radiation(GRID, SWnet, dt)
+            subsurface_melt, G_penetrating = penetrating_radiation(GRID, SWnet, dt, opt_dict)
         else:
             subsurface_melt = 0.0
             G_penetrating = 0.0
@@ -249,12 +253,14 @@ def cosipy_core(DATA, indY, indX, GRID_RESTART=None, stake_names=None, stake_dat
             # Find new surface temperature (LW is used from the input file)
             fun, surface_temperature, lw_radiation_in, lw_radiation_out, sensible_heat_flux, latent_heat_flux, \
                 ground_heat_flux, rain_heat_flux, rho, Lv, MOL, Cs_t, Cs_q, q0, q2 \
-                = update_surface_temperature(GRID, dt, z, z0, T2[t], RH2[t], PRES[t], sw_radiation_net, U2[t], RAIN, SLOPE, LWin=LWin[t])
+                = update_surface_temperature(GRID, dt, z, z0, T2[t], RH2[t], PRES[t],
+                                             sw_radiation_net, U2[t], RAIN, SLOPE, LWin=LWin[t], opt_dict=opt_dict)
         else:
             # Find new surface temperature (LW is parametrized using cloud fraction)
             fun, surface_temperature, lw_radiation_in, lw_radiation_out, sensible_heat_flux, latent_heat_flux, \
                 ground_heat_flux, rain_heat_flux, rho, Lv, MOL, Cs_t, Cs_q, q0, q2 \
-                = update_surface_temperature(GRID, dt, z, z0, T2[t], RH2[t], PRES[t], sw_radiation_net, U2[t], RAIN, SLOPE, N=N[t])
+                = update_surface_temperature(GRID, dt, z, z0, T2[t], RH2[t], PRES[t],
+                                             sw_radiation_net, U2[t], RAIN, SLOPE, N=N[t], opt_dict=opt_dict)
 
         #--------------------------------------------
         # Surface mass fluxes [m w.e.q.]
@@ -286,22 +292,22 @@ def cosipy_core(DATA, indY, indX, GRID_RESTART=None, stake_names=None, stake_dat
         #--------------------------------------------
         # Percolation
         #--------------------------------------------
-        Q  = percolation(GRID, melt + condensation + RAIN/1000.0 + lwc_from_melted_layers, dt)
+        Q  = percolation(GRID, melt + condensation + RAIN/1000.0 + lwc_from_melted_layers, dt, opt_dict)
 
         #--------------------------------------------
         # Refreezing
         #--------------------------------------------
-        water_refreezed = refreezing(GRID)
+        water_refreezed = refreezing(GRID, opt_dict)
 
         #--------------------------------------------
         # Solve the heat equation
         #--------------------------------------------
-        solveHeatEquation(GRID, dt)
+        solveHeatEquation(GRID, dt, opt_dict)
 
         #--------------------------------------------
         # Calculate new density to densification
         #--------------------------------------------
-        densification(GRID, SLOPE, dt)
+        densification(GRID, SLOPE, dt, opt_dict)
 
         #--------------------------------------------
         # Calculate mass balance

cosipy/cpkernel/init.py

@@ -3,10 +3,14 @@
 from constants import *
 from config import *
 from cosipy.cpkernel.grid import *
+from cosipy.utils.options import read_opt
 
-def init_snowpack(DATA):
+def init_snowpack(DATA, opt_dict):
     ''' INITIALIZATION '''
 
+    # Read and set options
+    read_opt(opt_dict, globals())
+
     ##--------------------------------------
     ## Check for WRF data
     ##--------------------------------------

cosipy/cpkernel/io.py

@@ -13,12 +13,15 @@
 from constants import *
 from config import * 
 import configparser
+from cosipy.utils.options import read_opt
 
 class IOClass:
 
-    def __init__(self, DATA=None):
+    def __init__(self, DATA=None, opt_dict=None):
         """ Init IO Class"""
 
+        # Read and set options
+        read_opt(opt_dict, globals())
         # Read variable list from file
         config = configparser.ConfigParser()
         config.read('./cosipy/output')

cosipy/modules/albedo.py

@@ -1,9 +1,12 @@
 import numpy as np
 from constants import albedo_method, albedo_fresh_snow, albedo_firn, albedo_ice, \
                       albedo_mod_snow_aging, albedo_mod_snow_depth, snow_ice_threshold
+from cosipy.utils.options import read_opt
 
-def updateAlbedo(GRID):
+def updateAlbedo(GRID, opt_dict):
     """ This methods updates the albedo """
+    # Read and set options
+    read_opt(opt_dict, globals())
     albedo_allowed = ['Oerlemans98']
     if albedo_method == 'Oerlemans98':
         alphaMod = method_Oerlemans(GRID)

cosipy/modules/densification.py

@@ -2,13 +2,16 @@
 from constants import densification_method, snow_ice_threshold, minimum_snow_layer_height, \
                       zero_temperature, ice_density
 from numba import njit
+from cosipy.utils.options import read_opt
 
-def densification(GRID,SLOPE,dt):
+def densification(GRID,SLOPE,dt, opt_dict=None):
     """ Densification of the snowpack
     Args:
         GRID    ::  GRID-Structure
 	dt      ::  integration time
     """
+    # Read and set options
+    read_opt(opt_dict, globals())
 
     densification_allowed = ['Boone', 'Vionnet', 'empirical', 'constant']
     if densification_method == 'Boone':

cosipy/modules/evaluation.py

@@ -1,11 +1,14 @@
 import numpy as np
 from config import eval_method, obs_type
+from cosipy.utils.options import read_opt
 
 
-def evaluate(stake_names, stake_data, df_):
+def evaluate(stake_names, stake_data, df_, opt_dict=None):
     """ This methods evaluates the simulation with the stake measurements
         stake_name  ::  """
-
+
+    # Read and set options
+    read_opt(opt_dict, globals())
     if eval_method == 'rmse':
         stat = rmse(stake_names, stake_data, df_)
     else:

cosipy/modules/heatEquation.py

@@ -1,8 +1,22 @@
 import numpy as np
 from numba import njit
+from cosipy.utils.options import read_opt
+
+def solveHeatEquation(GRID, dt, opt_dict=None):
+
+    # Read and set options
+    read_opt(opt_dict, globals())
+
+    heatEquation_method = 'default'
+    heatEquation_allowed = ['default']
+
+    if heatEquation_method == 'default':
+        heatEquation_default(GRID,dt)
+    else:
+        raise ValueError("Heat equation = \"{:s}\" is not allowed, must be one of {:s}".format(heatEquation_method, ", ".join(heatEquation_allowed)))
 
 @njit
-def solveHeatEquation(GRID, dt):
+def heatEquation_default(GRID, dt):
     """ Solves the heat equation on a non-uniform grid
 
     dt  ::  integration time

cosipy/modules/penetratingRadiation.py

@@ -2,8 +2,13 @@
 from constants import penetrating_method, snow_ice_threshold, spec_heat_ice, \
                       zero_temperature, ice_density, water_density, lat_heat_melting
 from numba import njit
+from cosipy.utils.options import read_opt
+
+def penetrating_radiation(GRID, SWnet, dt, opt_dict=None):
+
+    # Read and set options
+    read_opt(opt_dict, globals())
 
-def penetrating_radiation(GRID, SWnet, dt):
     penetrating_allowed = ['Bintanja95']
     if penetrating_method == 'Bintanja95':
         subsurface_melt, Si = method_Bintanja(GRID, SWnet, dt)

cosipy/modules/percolation.py

@@ -1,8 +1,26 @@
 import numpy as np
 from numba import njit
+from cosipy.utils.options import read_opt
+
+
+def percolation(GRID, water, dt, opt_dict=None):
+
+    # Read and set options
+    read_opt(opt_dict, globals())
+
+    percolation_method = 'bucket'
+    percolation_allowed = ['bucket']
+
+    if percolation_method == 'bucket':
+        Q = bucket(GRID,water,dt)
+    else:
+        raise ValueError("Percolation method = \"{:s}\" is not allowed, must be one of {:s}".format(percolation_method, ", ".join(percolation_allowed)))
+
+    return Q
+
 
 @njit
-def percolation(GRID, water, dt):
+def bucket(GRID, water, dt):
     """ Percolation and refreezing of melt water through the snow- and firn pack
 
     Args:

cosipy/modules/refreezing.py

@@ -2,9 +2,25 @@
 from constants import zero_temperature, spec_heat_ice, ice_density, \
                       water_density, lat_heat_melting
 from numba import njit
+from cosipy.utils.options import read_opt
 
-@njit
-def refreezing(GRID):
+def refreezing(GRID, opt_dict=None):
+
+    # Read and set options
+    read_opt(opt_dict, globals())
+
+    refreezing_method = 'default'
+    refreezing_allowed = ['default']
+
+    if refreezing_method == 'default':
+        water_refreezed = refreezing_default(GRID)
+    else:
+        raise ValueError("Refreezing method = \"{:s}\" is not allowed, must be one of {:s}".format(refreezing_method, ", ".join(refreezing_allowed)))
+
+    return water_refreezed
+
+#@njit
+def refreezing_default(GRID):
 
     # water refreezed
     water_refreezed = 0.0

cosipy/modules/roughness.py

@@ -1,8 +1,12 @@
 from constants import roughness_method, roughness_fresh_snow, \
                       roughness_firn, roughness_ice, snow_ice_threshold, \
                       aging_factor_roughness
+from cosipy.utils.options import read_opt
 
-def updateRoughness(GRID):
+def updateRoughness(GRID, opt_dict=None):
+
+    # Read and set options
+    read_opt(opt_dict, globals())
 
     roughness_allowed = ['Moelg12']
     if roughness_method == 'Moelg12':

cosipy/modules/surfaceTemperature.py

@@ -5,9 +5,11 @@
 from scipy.optimize import minimize, newton
 from numba import njit
 from types import SimpleNamespace
+from cosipy.utils.options import read_opt
 
 
-def update_surface_temperature(GRID, dt, z, z0, T2, rH2, p, SWnet, u2, RAIN, SLOPE, LWin=None, N=None):
+def update_surface_temperature(GRID, dt, z, z0, T2, rH2, p, SWnet, u2, RAIN, SLOPE, LWin=None,
+                               N=None, opt_dict=None):
     """ This methods updates the surface temperature and returns the surface fluxes
 
     Given:
@@ -47,6 +49,9 @@ def update_surface_temperature(GRID, dt, z, z0, T2, rH2, p, SWnet, u2, RAIN, SLO
         phi     ::  Stability correction term [-]
     """
 
+    # Read and set options
+    read_opt(opt_dict, globals())
+
     #Interpolate subsurface temperatures to selected subsurface depths for GHF computation
     B_Ts = interp_subT(GRID)