diff --git a/HISTORY.rst b/HISTORY.rst index 757d585e..6aaf1fb2 100644 --- a/HISTORY.rst +++ b/HISTORY.rst @@ -2,9 +2,16 @@ History ======= -v.0.16.11 +Latest --------- +* New component SimpleBoundaryLayer that implements a boundary layer. +* The boundary layer is created using the bulk Richardson number. +* Diffusion is implemented using the simplified Monin-Obukhov theory. + +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 diff --git a/climt/__init__.py b/climt/__init__.py index d9acb4e9..e4607a3a 100644 --- a/climt/__init__.py +++ b/climt/__init__.py @@ -12,7 +12,8 @@ GridScaleCondensation, BergerSolarInsolation, SimplePhysics, RRTMGLongwave, RRTMGShortwave, EmanuelConvection, SlabSurface, GFSDynamicalCore, - DcmipInitialConditions, IceSheet, Instellation, DryConvectiveAdjustment, BucketHydrology) + DcmipInitialConditions, IceSheet, Instellation, DryConvectiveAdjustment, + BucketHydrology, SimpleBoundaryLayer) sympl.set_constant('top_of_model_pressure', 20., 'Pa') @@ -26,6 +27,7 @@ GridScaleCondensation, BergerSolarInsolation, SimplePhysics, RRTMGLongwave, RRTMGShortwave, EmanuelConvection, SlabSurface, GFSDynamicalCore, DcmipInitialConditions, - IceSheet, Instellation, DryConvectiveAdjustment, BucketHydrology) + IceSheet, Instellation, DryConvectiveAdjustment, BucketHydrology, + SimpleBoundaryLayer) __version__ = '0.16.14' diff --git a/climt/_components/__init__.py b/climt/_components/__init__.py index d69265d8..20ac6bf6 100644 --- a/climt/_components/__init__.py +++ b/climt/_components/__init__.py @@ -12,10 +12,11 @@ from .instellation import Instellation from .dry_convection import DryConvectiveAdjustment from .bucket_hydrology import BucketHydrology +from .simple_boundary_layer import SimpleBoundaryLayer __all__ = ( Frierson06LongwaveOpticalDepth, GrayLongwaveRadiation, HeldSuarez, GridScaleCondensation, BergerSolarInsolation, SimplePhysics, RRTMGLongwave, RRTMGShortwave, EmanuelConvection, SlabSurface, GFSDynamicalCore, DcmipInitialConditions, IceSheet, - Instellation, DryConvectiveAdjustment, BucketHydrology) + Instellation, DryConvectiveAdjustment, BucketHydrology, SimpleBoundaryLayer) diff --git a/climt/_components/simple_boundary_layer/__init__.py b/climt/_components/simple_boundary_layer/__init__.py new file mode 100644 index 00000000..99fef5eb --- /dev/null +++ b/climt/_components/simple_boundary_layer/__init__.py @@ -0,0 +1,3 @@ +from .component import SimpleBoundaryLayer + +__all__ = (SimpleBoundaryLayer) diff --git a/climt/_components/simple_boundary_layer/component.py b/climt/_components/simple_boundary_layer/component.py new file mode 100644 index 00000000..6820c00c --- /dev/null +++ b/climt/_components/simple_boundary_layer/component.py @@ -0,0 +1,332 @@ +from sympl import initialize_numpy_arrays_with_properties, get_constant +from sympl import Stepper +import numpy as np +import numba +from numba import jit + + +@jit(nopython=True, parallel=True) +def Parallel_boundary(air_temperature, surface_temperature, air_pressure, + air_pressure_int, surface_pressure, specific_humidity, + surface_humidity, northward_wind, eastward_wind, + new_air_temperature, new_specific_humidity, + new_northward_wind, new_eastward_wind, + north_wind_stress, east_wind_stress, boundary_height, + Rd_val, Cp_val, g_val, k_val, z0_val, fb_val, P0_val, + Ric_val, num_cols, timestep): + + Rd = Rd_val + Cp_dry = Cp_val + g = g_val + k = k_val + z0 = z0_val + fb = fb_val + P0 = P0_val + Ri_c = Ric_val + + def K_b(Ri, Ri_a, u_fric, C, z): + + if Ri_a <= 0: + Kb = k*u_fric*np.sqrt(C)*z + else: + Kb = k*u_fric*np.sqrt(C)*z/(1+Ri/Ri_c*np.log(z/z0)/(1-Ri/Ri_c)) + + return(Kb) + + def TDMAsolver(a, b, c, d): + + n = len(d) + w = np.zeros(n-1) + g = np.zeros(n) + p = np.zeros(n) + + w[0] = c[0]/b[0] + g[0] = d[0]/b[0] + + for i in range(1, n-1): + w[i] = c[i]/(b[i] - a[i-1]*w[i-1]) + for i in range(1, n): + g[i] = (d[i] - a[i-1]*g[i-1])/(b[i] - a[i-1]*w[i-1]) + p[n-1] = g[n-1] + for i in range(n-1, 0, -1): + p[i-1] = g[i-1] - w[i-1]*p[i] + + return p + + def diffuse_profile(profile, p, p_int, rho, Diff, dt): + + num_layers = len(profile) + + diag_m = np.zeros(num_layers) + diag_p = np.zeros(num_layers) + diag = np.zeros(num_layers) + + for i in range(num_layers): + + if i != 0: + diag_m[i] = g*g*rho[i-1]*rho[i-1]*Diff[i-1]*dt/(p[i-1]-p[i])\ + * 1/(p_int[i]-p_int[i+1]) + + if i != num_layers-1: + diag_p[i] = g*g*rho[i]*rho[i]*Diff[i]*dt/(p[i]-p[i+1]) * \ + 1/(p_int[i]-p_int[i+1]) + + diag[i] = 1+diag_m[i]+diag_p[i] + + return TDMAsolver(-diag_m[1:], diag, -diag_p[:-1], profile) + + for col in numba.prange(num_cols): + + col_air_temperature = air_temperature[:, col] + col_surface_temperature = surface_temperature[col] + col_air_pressure = air_pressure[:, col] + col_air_pressure_int = air_pressure_int[:, col] + col_surface_pressure = surface_pressure[col] + col_specific_humidity = specific_humidity[:, col] + col_surface_humidity = surface_humidity[col] + col_north_wind = northward_wind[:, col] + col_east_wind = eastward_wind[:, col] + + col_north_wind_int = 0.5*(col_north_wind[1:]+col_north_wind[:-1]) + col_east_wind_int = 0.5*(col_east_wind[1:]+col_east_wind[:-1]) + col_air_temperature_int = 0.5*(col_air_temperature[1:] + + col_air_temperature[:-1]) + col_specific_humidity_int = 0.5*(col_specific_humidity[1:] + + col_specific_humidity[:-1]) + col_rho = col_air_pressure_int[1:-1]/(Rd * (1+0.608 * + col_specific_humidity_int) * + col_air_temperature_int) + + n = len(col_air_temperature_int) + + pot_virt_temp = col_air_temperature_int *\ + np.power((P0/col_air_pressure_int[1:-1]), Rd/Cp_dry) * \ + (1+0.608*col_specific_humidity_int) + pot_virt_temp_surf = col_surface_temperature *\ + np.power((P0/col_surface_pressure), Rd/Cp_dry) *\ + (1+0.608*col_surface_humidity) + + z = np.zeros(n) + z[0] = (Rd*(1+0.608*col_specific_humidity[0])*col_air_temperature[0] / + g) * np.log(col_surface_pressure/col_air_pressure_int[1:-1][0]) + for i in range(1, n): + z[i] = z[i-1]+(Rd*(1+0.608*col_specific_humidity[i]) * + col_air_temperature[i]/g) *\ + np.log(col_air_pressure_int[1:-1][i-1] / + col_air_pressure_int[1:-1][i]) + + wind_int = np.sqrt(np.power(col_north_wind_int, 2) + + np.power(col_east_wind_int, 2)) + for i in range(len(wind_int)): + if wind_int[i] < 1: + wind_int[i] = 1 + + Ri_a = g*z[0]*(pot_virt_temp[0]-pot_virt_temp_surf)/( + pot_virt_temp_surf*wind_int[0]*wind_int[0]) + if Ri_a < 0: + C = k*k*np.power(np.log(z[0]/z0), -2) + elif Ri_a < Ri_c: + C = k*k*np.power(np.log(z[0]/z0), -2)*np.power((1-Ri_a/Ri_c), 2) + else: + C = 0 + + count = 0 + Rich = np.zeros(n) + for i in range(n): + Rich[i] = g*z[i]*(pot_virt_temp[i]-pot_virt_temp[0])/( + pot_virt_temp[0]*wind_int[i]*wind_int[i]) + if Rich[i] > Ri_c: + count = i+1 + break + h = z[count-1] + boundary_height[col] = h + + north_wind_stress[col] = col_rho[0]*C*wind_int[0]*col_north_wind_int[0] + east_wind_stress[col] = col_rho[0]*C*wind_int[0]*col_east_wind_int[0] + + u_fric = wind_int[0] + + diff = np.zeros(n) + + for i in range(count): + if z[i] < fb*h: + diff[i] = K_b(Rich[i], Ri_a, u_fric, C, z[i]) + + else: + diff[i] = K_b(Rich[i], Ri_a, u_fric, C, fb*h)*z[i]/(h*fb) *\ + np.power(1-(z[i]-fb*h)/((1-fb)*h), 2) + + new_temp = diffuse_profile(col_air_temperature, col_air_pressure, + col_air_pressure_int, col_rho, diff, + timestep) + new_humidity = diffuse_profile(col_specific_humidity, col_air_pressure, + col_air_pressure_int, col_rho, diff, + timestep) + new_north_wind = diffuse_profile(col_north_wind, col_air_pressure, + col_air_pressure_int, col_rho, diff, + timestep) + new_east_wind = diffuse_profile(col_east_wind, col_air_pressure, + col_air_pressure_int, col_rho, diff, + timestep) + + new_air_temperature[:, col] = new_temp + new_specific_humidity[:, col] = new_humidity + new_northward_wind[:, col] = new_north_wind + new_eastward_wind[:, col] = new_east_wind + + +class SimpleBoundaryLayer(Stepper): + """ + This is a simple boundary layer component that diffuses heat, humidity and + momemtum upwards from the lowest model level. + + This component assumes that a surface flux component has been already run, + which has made the changes due to surface fluxes at the lowest model + level. This component then diffuses heat, humidity and momentum using + diffusion coefficients calculated using the simplified Monin-Obukhov + theory. + """ + + input_properties = { + 'air_temperature': { + 'dims': ['mid_levels', '*'], + 'units': 'degK ', + }, + 'specific_humidity': { + 'dims': ['mid_levels', '*'], + 'units': 'kg/kg', + }, + 'air_pressure': { + 'dims': ['mid_levels', '*'], + 'units': 'Pa', + }, + 'air_pressure_on_interface_levels': { + 'dims': ['interface_levels', '*'], + 'units': 'Pa', + }, + 'northward_wind': { + 'dims': ['mid_levels', '*'], + 'units': 'm s^-1', + }, + 'eastward_wind': { + 'dims': ['mid_levels', '*'], + 'units': 'm s^-1', + }, + 'surface_air_pressure': { + 'dims': ['*'], + 'units': 'Pa', + }, + 'surface_temperature': { + 'dims': ['*'], + 'units': 'degK', + }, + 'surface_specific_humidity': { + 'dims': ['*'], + 'units': 'kg/kg', + }, + } + + output_properties = { + 'air_temperature': { + 'dims': ['mid_levels', '*'], + 'units': 'degK ', + }, + 'specific_humidity': { + 'dims': ['mid_levels', '*'], + 'units': 'kg/kg', + }, + 'northward_wind': { + 'dims': ['mid_levels', '*'], + 'units': 'm s^-1', + }, + 'eastward_wind': { + 'dims': ['mid_levels', '*'], + 'units': 'm s^-1', + }, + } + + diagnostic_properties = { + 'northward_wind_stress': { + 'dims': ['*'], + 'units': 'Pa', + }, + 'eastward_wind_stress': { + 'dims': ['*'], + 'units': 'Pa', + }, + 'boundary_layer_height': { + 'dims': ['*'], + 'units': 'm', + }, + } + + def __init__(self, von_karman_constant=0.4, roughness_length=0.0000321, + specific_fraction=0.1, reference_pressure=100000, + critical_richardson_number=1, **kwargs): + """ + Args: + roughness_length: + A measure of the surface roughness. + specific_fraction: + A parameter used in the calculation of diffusion coefficients. + reference_pressure: + The reference pressure used in the potential temperature + calculations. + critical_richardson_number: + A set threshold value which determines the diffusion coefficients + and the height of the boundary layer. + """ + + self._k = von_karman_constant + self._z0 = roughness_length + self._fb = specific_fraction + self._P0 = reference_pressure + self._Ric = critical_richardson_number + self._update_constants() + + super(SimpleBoundaryLayer, self).__init__(**kwargs) + + def _update_constants(self): + + self._Rd = get_constant('gas_constant_of_dry_air', 'J kg^-1 K^-1') + self._Cp =\ + get_constant('heat_capacity_of_dry_air_at_constant_pressure', + 'J kg^-1 K^-1') + self._g = get_constant('gravitational_acceleration', 'm s^-2') + + def array_call(self, state, timestep): + """ + Takes temperature, humidty and wind profiles for each column and + returns diffused temperature, humidity and wind profiles. + """ + + num_cols = state['air_temperature'].shape[1] + + new_state = initialize_numpy_arrays_with_properties( + self.output_properties, state, self.input_properties + ) + + diagnostics = initialize_numpy_arrays_with_properties( + self.diagnostic_properties, state, self.input_properties + ) + + Parallel_boundary(state['air_temperature'], + state['surface_temperature'], + state['air_pressure'], + state['air_pressure_on_interface_levels'], + state['surface_air_pressure'], + state['specific_humidity'], + state['surface_specific_humidity'], + state['northward_wind'], state['eastward_wind'], + new_state['air_temperature'], + new_state['specific_humidity'], + new_state['northward_wind'], + new_state['eastward_wind'], + diagnostics['northward_wind_stress'], + diagnostics['eastward_wind_stress'], + diagnostics['boundary_layer_height'], + self._Rd, self._Cp, self._g, self._k, self._z0, + self._fb, self._P0, self._Ric, num_cols, + timestep.total_seconds()) + + return diagnostics, new_state diff --git a/climt/_lib/Makefile b/climt/_lib/Makefile index 94706c65..6b5711c6 100644 --- a/climt/_lib/Makefile +++ b/climt/_lib/Makefile @@ -61,7 +61,7 @@ sht_lib: $(CLIMT_ARCH)/libshtns_omp.a blas_lib: $(CLIMT_ARCH)/libopenblas.a $(CLIMT_ARCH)/libopenblas.a: - if [ ! -d $(BLAS_DIR) ]; then tar -xvzf $(BLAS_GZ) > log 2>&1; fi; cd $(BLAS_DIR); CFLAGS=$(CLIMT_CFLAGS) make NO_SHARED=1 NO_LAPACK=0 > log 2>&1 ; make PREFIX=$(BASE_DIR) NO_SHARED=1 NO_LAPACK=0 install > log.install.blas 2>&1 ; cp ../lib/libopenblas.a $(LIB_DIR) + if [ ! -d $(BLAS_DIR) ]; then tar -xvzf $(BLAS_GZ) > log 2>&1; fi; cd $(BLAS_DIR); CFLAGS=$(CLIMT_CFLAGS) make NO_SHARED=1 NO_LAPACK=0 > log 2>&1 ; cat log; make PREFIX=$(BASE_DIR) NO_SHARED=1 NO_LAPACK=0 install > log.install.blas 2>&1 ; cp ../lib/libopenblas.a $(LIB_DIR) # GFS Configuration diff --git a/tests/cached_component_output/TestSimpleBoundaryLayer-3d-0.cache b/tests/cached_component_output/TestSimpleBoundaryLayer-3d-0.cache new file mode 100644 index 00000000..e9082abb Binary files /dev/null and b/tests/cached_component_output/TestSimpleBoundaryLayer-3d-0.cache differ diff --git a/tests/cached_component_output/TestSimpleBoundaryLayer-3d-1.cache b/tests/cached_component_output/TestSimpleBoundaryLayer-3d-1.cache new file mode 100644 index 00000000..8498c793 Binary files /dev/null and b/tests/cached_component_output/TestSimpleBoundaryLayer-3d-1.cache differ diff --git a/tests/cached_component_output/TestSimpleBoundaryLayer-column-0.cache b/tests/cached_component_output/TestSimpleBoundaryLayer-column-0.cache new file mode 100644 index 00000000..7f9d2f16 Binary files /dev/null and b/tests/cached_component_output/TestSimpleBoundaryLayer-column-0.cache differ diff --git a/tests/cached_component_output/TestSimpleBoundaryLayer-column-1.cache b/tests/cached_component_output/TestSimpleBoundaryLayer-column-1.cache new file mode 100644 index 00000000..9f9d77b0 Binary files /dev/null and b/tests/cached_component_output/TestSimpleBoundaryLayer-column-1.cache differ diff --git a/tests/test_components.py b/tests/test_components.py index e733e9bc..77a9dc92 100644 --- a/tests/test_components.py +++ b/tests/test_components.py @@ -13,7 +13,7 @@ RRTMGShortwave, SlabSurface, EmanuelConvection, DcmipInitialConditions, GFSDynamicalCore, BucketHydrology, IceSheet, Instellation, DryConvectiveAdjustment, - get_grid) + SimpleBoundaryLayer, get_grid) import climt from sympl import ( Stepper, TendencyStepper, TimeDifferencingWrapper, @@ -46,7 +46,7 @@ def load_dictionary(filename): def state_3d_to_1d(state): return_state = {} for name, value in state.items(): - if name is 'time': + if name == 'time': return_state[name] = value else: dim_list = [] @@ -62,7 +62,7 @@ def state_3d_to_1d(state): def transpose_state(state, dims=None): return_state = {} for name, value in state.items(): - if name is 'time': + if name == 'time': return_state[name] = state[name] else: if dims is None: @@ -333,6 +333,12 @@ def get_component_instance(self): return SimplePhysics() +class TestSimpleBoundaryLayer(ComponentBaseColumn, ComponentBase3D): + + def get_component_instance(self): + return SimpleBoundaryLayer() + + class TestSimplePhysicsImplicitPrognostic(ComponentBaseColumn, ComponentBase3D): def get_component_instance(self): diff --git a/tests/test_conservation.py b/tests/test_conservation.py index a0eff4f5..a7249771 100644 --- a/tests/test_conservation.py +++ b/tests/test_conservation.py @@ -273,6 +273,19 @@ def get_quantity_forcing(self, state): return 0 +class TestSimpleBoundaryLayerConservation(AtmosphereMoistEnthalpyConservation): + + def get_component_instance(self): + return climt.SimpleBoundaryLayer() + + def modify_state(self, state): + state['eastward_wind'].values[:] = 3. + unstable_level = 5 + state['air_temperature'][:unstable_level] += 10 + state['specific_humidity'][:unstable_level] = 0.05 + return state + + class TestDryConvectionCondensibleConservation(AtmosphereTracerConservation): def get_component_instance(self): diff --git a/tests/test_initialization.py b/tests/test_initialization.py index 1bba7544..65e90fdb 100644 --- a/tests/test_initialization.py +++ b/tests/test_initialization.py @@ -2,7 +2,7 @@ get_default_state, Frierson06LongwaveOpticalDepth, GrayLongwaveRadiation, HeldSuarez, GridScaleCondensation, BergerSolarInsolation, SimplePhysics, RRTMGLongwave, RRTMGShortwave, EmanuelConvection, SlabSurface, GFSDynamicalCore, DcmipInitialConditions, IceSheet, - Instellation, get_grid + Instellation, SimpleBoundaryLayer, BucketHydrology, get_grid ) import random from sympl import ( @@ -189,7 +189,8 @@ class ComponentQuantityInitializationTests(unittest.TestCase): RRTMGShortwave, EmanuelConvection, SlabSurface, DcmipInitialConditions, IceSheet, - Instellation + Instellation, SimpleBoundaryLayer, + BucketHydrology ) pair_tests = 20