Repository: CliMT/climt (language: python)

Issue:
Descriptive variable names in SCM example
In the example script [radiative_equilibrium_rrtmg.py](https://github.com/CliMT/climt-future/blob/master/examples/radiative_equilibrium_rrtmg.py), the names tp_profiles and mol_profiles could be more descriptive. Perhaps they should be named thermodynamic_profiles and mole_fraction_profiles? The arrays could also be saved in a npz file, which allows dictionary lookup on its contents, so that the standalone npz file describes what's inside it without needing to know beforehand.


Discussion / hints:
Fixed.

Fix patch:
diff --git a/AUTHORS.rst b/AUTHORS.rst
index 95021970..4a88f42e 100644
--- a/AUTHORS.rst
+++ b/AUTHORS.rst
@@ -22,3 +22,4 @@ Contributors
 * Dargan Frierson <dargan@atmos.washington.edu>
 * Jonathan Mitchell <jonmitch@ucla.edu>
 * Suhas D L <suhasd@iisc.ac.in>
+* Abel Shibu <abel.shibu@students.iiserpune.ac.in>
diff --git a/HISTORY.rst b/HISTORY.rst
index 19c2e660..80a85666 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -2,6 +2,13 @@
 History
 =======
 
+Latest
+--------
+
+* New component BucketHydrology that implements Manabe first generation land model
+* BucketHydrology calculates the sensible and latent heat flux within the component
+* Conservation test for the component also added
+
 v.0.16.4
 --------
 
@@ -44,7 +51,7 @@ v0.16.0
   so that partial cloud fractions can be used.
 * Add GMD example scripts to repository.
 * Fix docs to reflect API changes after refactor.
-* Fix wrong initialisation to use sigma values instead of pressure values 
+* Fix wrong initialisation to use sigma values instead of pressure values
   of optical depth for GrayLongwaveRadiation
 
 Breaking Changes
diff --git a/climt/__init__.py b/climt/__init__.py
index 1b35ebcf..4b34de91 100644
--- a/climt/__init__.py
+++ b/climt/__init__.py
@@ -12,7 +12,7 @@
     GridScaleCondensation, BergerSolarInsolation, SimplePhysics, RRTMGLongwave,
     RRTMGShortwave,
     EmanuelConvection, SlabSurface, GFSDynamicalCore,
-    DcmipInitialConditions, IceSheet, Instellation, DryConvectiveAdjustment)
+    DcmipInitialConditions, IceSheet, Instellation, DryConvectiveAdjustment, BucketHydrology)
 
 
 sympl.set_constant('top_of_model_pressure', 20., 'Pa')
@@ -26,6 +26,6 @@
     GridScaleCondensation, BergerSolarInsolation, SimplePhysics, RRTMGLongwave,
     RRTMGShortwave,
     EmanuelConvection, SlabSurface, GFSDynamicalCore, DcmipInitialConditions,
-    IceSheet, Instellation, DryConvectiveAdjustment)
+    IceSheet, Instellation, DryConvectiveAdjustment, BucketHydrology)
 
 __version__ = '0.16.3'
diff --git a/climt/_components/__init__.py b/climt/_components/__init__.py
index 40f29881..d69265d8 100644
--- a/climt/_components/__init__.py
+++ b/climt/_components/__init__.py
@@ -11,10 +11,11 @@
 from .dcmip import DcmipInitialConditions
 from .instellation import Instellation
 from .dry_convection import DryConvectiveAdjustment
+from .bucket_hydrology import BucketHydrology
 
 __all__ = (
     Frierson06LongwaveOpticalDepth, GrayLongwaveRadiation,
     HeldSuarez, GridScaleCondensation, BergerSolarInsolation, SimplePhysics,
     RRTMGLongwave, RRTMGShortwave, EmanuelConvection, SlabSurface,
     GFSDynamicalCore, DcmipInitialConditions, IceSheet,
-    Instellation, DryConvectiveAdjustment)
+    Instellation, DryConvectiveAdjustment, BucketHydrology)
diff --git a/climt/_components/bucket_hydrology/__init__.py b/climt/_components/bucket_hydrology/__init__.py
new file mode 100644
index 00000000..adaa293f
--- /dev/null
+++ b/climt/_components/bucket_hydrology/__init__.py
@@ -0,0 +1,3 @@
+from .component import BucketHydrology
+
+__all__ = (BucketHydrology)
diff --git a/climt/_components/bucket_hydrology/component.py b/climt/_components/bucket_hydrology/component.py
new file mode 100644
index 00000000..3efd43a5
--- /dev/null
+++ b/climt/_components/bucket_hydrology/component.py
@@ -0,0 +1,192 @@
+from math import sqrt
+from sympl import Stepper
+
+
+class BucketHydrology(Stepper):
+    """
+    Manages surface energy and moisture balance
+
+    This component assumes that the surface is a slab with some heat capacity and moisture holding capacity.
+    Calculates the sensible and latent heat flux, takes precipitation values as input.
+
+    """
+
+    input_properties = {
+        'downwelling_longwave_flux_in_air': {
+            'dims': ['*', 'interface_levels'],
+            'units': 'W m^-2',
+        },
+        'downwelling_shortwave_flux_in_air': {
+            'dims': ['*', 'interface_levels'],
+            'units': 'W m^-2',
+        },
+        'upwelling_longwave_flux_in_air': {
+            'dims': ['*', 'interface_levels'],
+            'units': 'W m^-2',
+        },
+        'upwelling_shortwave_flux_in_air': {
+            'dims': ['*', 'interface_levels'],
+            'units': 'W m^-2',
+        },
+        'surface_temperature': {
+            'dims': ['*'],
+            'units': 'degK',
+        },
+        'surface_material_density': {
+            'dims': ['*'],
+            'units': 'kg m^-3',
+        },
+        'soil_layer_thickness': {
+            'dims': ['*'],
+            'units': 'm',
+        },
+        'heat_capacity_of_soil': {
+            'dims': ['*'],
+            'units': 'J kg^-1 degK^-1',
+        },
+        'lwe_thickness_of_soil_moisture_content': {
+            'dims': ['*'],
+            'units': 'm',
+        },
+        'convective_precipitation_rate': {
+            'dims': ['*'],
+            'units': 'm s^-1 ',
+        },
+        'stratiform_precipitation_rate': {
+            'dims': ['*'],
+            'units': 'm s^-1 ',
+        },
+        'specific_humidity': {
+            'dims': ['mid_levels', '*'],
+            'units': 'kg/kg ',
+        },
+        'surface_specific_humidity': {
+            'dims': ['*'],
+            'units': 'kg/kg ',
+        },
+        'air_temperature': {
+            'dims': ['mid_levels', '*'],
+            'units': 'degK ',
+        },
+        'northward_wind': {
+            'dims': ['mid_levels', '*'],
+            'units': 'm s^-1',
+        },
+        'eastward_wind': {
+            'dims': ['mid_levels', '*'],
+            'units': 'm s^-1',
+        },
+    }
+
+    diagnostic_properties = {
+        'precipitation_rate': {
+            'dims': ['*'],
+            'units': 'm s^-1 ',
+        },
+        'surface_upward_latent_heat_flux': {
+            'dims': ['*'],
+            'units': 'W m^-2',
+        },
+        'surface_upward_sensible_heat_flux': {
+            'dims': ['*'],
+            'units': 'W m^-2',
+        },
+        'evaporation_rate': {
+            'dims': ['*'],
+            'units': 'm s^-1',
+        },
+    }
+
+    output_properties = {
+        'surface_temperature': {'units': 'degK'},
+        'lwe_thickness_of_soil_moisture_content': {'units': 'm'},
+    }
+
+    def __init__(self, soil_moisture_max=0.15, beta_parameter=0.75,
+                 specific_latent_heat_of_water=2260000, bulk_coefficient=0.0011, **kwargs):
+        """
+        Args:
+
+        soil_moisture_max:
+            The maximum moisture that can be held by the surface_temperature
+        beta_parameter:
+            A constant value that is used in the beta_factor calculation.
+        bulk_coefficient:
+            The bulk transfer coeffiecient that is used to calculate
+            maximum evaporation rate and sensible heat flux
+        """
+        self._smax = soil_moisture_max
+        self._g = beta_parameter
+        self._c = bulk_coefficient
+        self._l = specific_latent_heat_of_water
+        super(BucketHydrology, self).__init__(**kwargs)
+
+    def array_call(self, state, timestep):
+        '''
+        Calculates sensible and latent heat flux and returns
+        surface temperature and soil moisture after timestep.
+        '''
+
+        beta_factor = 0
+
+        north_wind_speed = state['northward_wind'][0][0]
+        east_wind_speed = state['eastward_wind'][0][0]
+
+        wind_speed = sqrt(pow(north_wind_speed, 2) +
+                          pow(east_wind_speed, 2))
+        potential_evaporation = self._c * wind_speed * (state['surface_specific_humidity'] - state['specific_humidity'][0])
+
+        precipitation_rate = state['convective_precipitation_rate'] + state['stratiform_precipitation_rate']
+
+        soil_moisture = state['lwe_thickness_of_soil_moisture_content'][0]
+
+        soil_moisture_tendency = 0
+
+        if soil_moisture >= self._g * self._smax:
+            beta_factor = 1
+        else:
+            beta_factor = soil_moisture/(self._g*self._smax)
+
+        evaporation_rate = beta_factor * potential_evaporation
+
+        if soil_moisture < self._smax or precipitation_rate <= evaporation_rate:
+            soil_moisture_tendency = precipitation_rate - evaporation_rate
+        else:
+            soil_moisture_tendency = 0
+
+        surface_upward_latent_heat_flux = self._l * evaporation_rate
+        surface_upward_sensible_heat_flux = self._c * wind_speed * (state['surface_temperature'] - state['air_temperature'][0])
+
+        net_heat_flux = (
+            state['downwelling_shortwave_flux_in_air'][:, 0] +
+            state['downwelling_longwave_flux_in_air'][:, 0] -
+            state['upwelling_shortwave_flux_in_air'][:, 0] -
+            state['upwelling_longwave_flux_in_air'][:, 0] -
+            surface_upward_sensible_heat_flux -
+            surface_upward_latent_heat_flux
+        )
+
+        mass_surface_slab = state['surface_material_density'] * \
+            state['soil_layer_thickness']
+        heat_capacity_surface = mass_surface_slab * state['heat_capacity_of_soil']
+
+        new_surface_temperature = state['surface_temperature'] + \
+            (net_heat_flux/heat_capacity_surface * timestep.total_seconds())
+        new_soil_moisture = state['lwe_thickness_of_soil_moisture_content'] + \
+            (soil_moisture_tendency * timestep.total_seconds())
+
+        new_soil_moisture[new_soil_moisture > 0.15] = 0.15
+
+        new_state = {
+            'surface_temperature': new_surface_temperature,
+            'lwe_thickness_of_soil_moisture_content': new_soil_moisture,
+        }
+
+        diagnostics = {
+            'precipitation_rate': precipitation_rate,
+            'surface_upward_sensible_heat_flux': surface_upward_sensible_heat_flux,
+            'surface_upward_latent_heat_flux': surface_upward_latent_heat_flux,
+            'evaporation_rate': evaporation_rate,
+        }
+
+        return diagnostics, new_state
diff --git a/climt/_core/initialization.py b/climt/_core/initialization.py
index 0655f6df..b7dc35b1 100644
--- a/climt/_core/initialization.py
+++ b/climt/_core/initialization.py
@@ -750,6 +750,9 @@ def get_exponent_for_sigma(b_half, num_sigma_levels):
 
     'solar_cycle_fraction': {'value': 0., 'units': 'dimensionless', 'domain': None},
     'flux_adjustment_for_earth_sun_distance': {'value': 1.0, 'units': 'dimensionless', 'domain': None},
+    'lwe_thickness_of_soil_moisture_content': {'value': 0, 'units': 'm', 'domain': 'surface'},
+    'convective_precipitation_rate': {'value': 0., 'units': 'mm day^-1', 'domain': 'surface'},
+    'stratiform_precipitation_rate': {'value': 0., 'units': 'm s^-1', 'domain': 'surface'},
 }
 for d in default_values.values():
     if 'grid' not in d.keys():
diff --git a/docs/components.rst b/docs/components.rst
index cd3686b8..1a40ee18 100644
--- a/docs/components.rst
+++ b/docs/components.rst
@@ -62,6 +62,9 @@ Surface Processes
     SlabSurface
     SlabSurface.__call__
 
+    BucketHydrology
+    BucketHydrology.__call__
+
 Ice and Snow
 ------------
 


Test patch:
diff --git a/tests/cached_component_output/TestBucketHydrology-3d-0.cache b/tests/cached_component_output/TestBucketHydrology-3d-0.cache
new file mode 100644
index 00000000..3b498648
Binary files /dev/null and b/tests/cached_component_output/TestBucketHydrology-3d-0.cache differ
diff --git a/tests/cached_component_output/TestBucketHydrology-3d-1.cache b/tests/cached_component_output/TestBucketHydrology-3d-1.cache
new file mode 100644
index 00000000..ae35b820
Binary files /dev/null and b/tests/cached_component_output/TestBucketHydrology-3d-1.cache differ
diff --git a/tests/cached_component_output/TestBucketHydrology-column-0.cache b/tests/cached_component_output/TestBucketHydrology-column-0.cache
new file mode 100644
index 00000000..21ba230c
Binary files /dev/null and b/tests/cached_component_output/TestBucketHydrology-column-0.cache differ
diff --git a/tests/cached_component_output/TestBucketHydrology-column-1.cache b/tests/cached_component_output/TestBucketHydrology-column-1.cache
new file mode 100644
index 00000000..9159dfdf
Binary files /dev/null and b/tests/cached_component_output/TestBucketHydrology-column-1.cache differ
diff --git a/tests/test_bucket_hydrology.py b/tests/test_bucket_hydrology.py
new file mode 100644
index 00000000..507ac1a9
--- /dev/null
+++ b/tests/test_bucket_hydrology.py
@@ -0,0 +1,82 @@
+from datetime import timedelta
+import numpy as np
+from climt import (BucketHydrology, get_default_state)
+
+def get_quantities(state):
+    heat = state['surface_material_density'].values*state['soil_layer_thickness'].values*\
+           state['surface_temperature'].values*state['heat_capacity_of_soil'].values
+    moisture = state['lwe_thickness_of_soil_moisture_content'].values
+
+    return heat, moisture
+
+time_step = timedelta(seconds=1)
+state = get_default_state([BucketHydrology()])
+
+state['upwelling_shortwave_flux_in_air'].values[:] = 40.
+state['upwelling_longwave_flux_in_air'].values[:] = 40.
+state['downwelling_shortwave_flux_in_air'].values[:] = 40.
+state['downwelling_longwave_flux_in_air'].values[:] = 40.
+state['soil_layer_thickness'].values[:] = 1.
+state['surface_specific_humidity'].values[:] = 1
+state['northward_wind'].values[:] = 3
+state['eastward_wind'].values[:] = 4
+state['surface_temperature'].values[:] = 300
+state['stratiform_precipitation_rate'].values[:] = 0.005
+state['convective_precipitation_rate'].values[:] = 0.5
+state['lwe_thickness_of_soil_moisture_content'].values[:] = 0.07
+
+
+diag, new_state = BucketHydrology()(state, time_step)
+
+
+surf_forcing = 0
+
+if 'upwelling_shortwave_flux_in_air' in state:
+    surf_forcing -= state['upwelling_shortwave_flux_in_air'].to_units('W/m^2')[0].values
+
+if 'upwelling_longwave_flux_in_air' in state:
+    surf_forcing -= state['upwelling_longwave_flux_in_air'].to_units('W/m^2')[0].values
+
+if 'downwelling_shortwave_flux_in_air' in state:
+    surf_forcing += state['downwelling_shortwave_flux_in_air'].to_units('W/m^2')[0].values
+
+if 'downwelling_longwave_flux_in_air' in state:
+    surf_forcing += state['downwelling_longwave_flux_in_air'].to_units('W/m^2')[0].values
+
+if 'surface_upward_sensible_heat_flux' in diag:
+    surf_forcing -= diag['surface_upward_sensible_heat_flux'].to_units('W/m^2').values
+
+if 'surface_upward_latent_heat_flux' in diag:
+    surf_forcing -= diag['surface_upward_latent_heat_flux'].to_units('W/m^2').values
+
+
+mois_forcing = 0
+
+if 'convective_precipitation_rate' in state:
+    mois_forcing += state['convective_precipitation_rate'].to_units('m s^-1').values
+
+if 'stratiform_precipitation_rate' in state:
+    mois_forcing += state['stratiform_precipitation_rate'].to_units('m s^-1').values
+
+if 'evaporation_rate' in diag:
+    mois_forcing -= diag['evaporation_rate'].to_units('m s^-1').values
+
+old_heat, old_moisture = get_quantities(state)
+
+state.update(diag)
+state.update(new_state)
+
+new_heat, new_moisture = get_quantities(state)
+
+
+#Test starts
+
+def test_bucket_hydrology_heat_conservation():
+    forcing_amount = surf_forcing * time_step.total_seconds()
+    assert np.isclose(new_heat - old_heat, forcing_amount,
+                      rtol=0, atol=1e-3)
+
+def test_bucket_hydrology_moisture_conservation():
+    forcing_amount = mois_forcing * time_step.total_seconds()
+    assert np.isclose(new_moisture - old_moisture, forcing_amount,
+                      rtol=0, atol=1e-3)
diff --git a/tests/test_components.py b/tests/test_components.py
index 85da42be..e733e9bc 100644
--- a/tests/test_components.py
+++ b/tests/test_components.py
@@ -11,7 +11,7 @@
     Frierson06LongwaveOpticalDepth, GridScaleCondensation,
     BergerSolarInsolation, SimplePhysics, RRTMGLongwave,
     RRTMGShortwave, SlabSurface, EmanuelConvection,
-    DcmipInitialConditions, GFSDynamicalCore,
+    DcmipInitialConditions, GFSDynamicalCore, BucketHydrology,
     IceSheet, Instellation, DryConvectiveAdjustment,
     get_grid)
 import climt
@@ -437,6 +437,11 @@ def get_component_instance(self):
         return SlabSurface()
 
 
+class TestBucketHydrology(ComponentBaseColumn, ComponentBase3D):
+    def get_component_instance(self):
+        return BucketHydrology()
+
+
 class TestEmanuel(ComponentBaseColumn, ComponentBase3D):
     def get_component_instance(self):
         emanuel = EmanuelConvection()