From b11f2fbc0c52f49dfa06340741f85ad63c03bca1 Mon Sep 17 00:00:00 2001 From: Abel Date: Sat, 23 Oct 2021 11:51:30 +0530 Subject: [PATCH 01/19] Component commit --- .../simple_boundary_layer/component.py | 305 ++++++++++++++++++ 1 file changed, 305 insertions(+) create mode 100644 climt/_components/simple_boundary_layer/component.py diff --git a/climt/_components/simple_boundary_layer/component.py b/climt/_components/simple_boundary_layer/component.py new file mode 100644 index 00000000..e18414e2 --- /dev/null +++ b/climt/_components/simple_boundary_layer/component.py @@ -0,0 +1,305 @@ +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_columns(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): + + 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): + + 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): + + 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_columns(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 From 88e8b6f3fdba2dcb768231f6a363ecc1e16e0a18 Mon Sep 17 00:00:00 2001 From: Abel Date: Sat, 23 Oct 2021 11:58:38 +0530 Subject: [PATCH 02/19] Calling the component at each level --- climt/__init__.py | 6 ++++-- climt/_components/__init__.py | 3 ++- climt/_components/simple_boundary_layer/__init__.py | 3 +++ 3 files changed, 9 insertions(+), 3 deletions(-) create mode 100644 climt/_components/simple_boundary_layer/__init__.py diff --git a/climt/__init__.py b/climt/__init__.py index 4b34de91..f9577e76 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.3' 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) From 006a21699f318bb083e18a5ac148e58b69bcf0c2 Mon Sep 17 00:00:00 2001 From: Abel Date: Sat, 23 Oct 2021 12:17:51 +0530 Subject: [PATCH 03/19] Changed history.rst --- HISTORY.rst | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/HISTORY.rst b/HISTORY.rst index 80a85666..0a03fe3e 100644 --- a/HISTORY.rst +++ b/HISTORY.rst @@ -5,6 +5,13 @@ History 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 From 775c71726581c05619d495fc433b051b94a8f146 Mon Sep 17 00:00:00 2001 From: Abel Date: Sat, 23 Oct 2021 12:46:48 +0530 Subject: [PATCH 04/19] Change version --- HISTORY.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/HISTORY.rst b/HISTORY.rst index eaf09aa3..6aaf1fb2 100644 --- a/HISTORY.rst +++ b/HISTORY.rst @@ -2,7 +2,7 @@ History ======= -v.0.16.11 +Latest --------- * New component SimpleBoundaryLayer that implements a boundary layer. From 8259683951a206cb340cc203e3dcc51b78fa5f12 Mon Sep 17 00:00:00 2001 From: Abel Date: Sat, 23 Oct 2021 12:53:18 +0530 Subject: [PATCH 05/19] Test_component added --- tests/test_components.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/tests/test_components.py b/tests/test_components.py index e733e9bc..94f17800 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, @@ -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): From f268f53fbb16d2860033a54746fd897e714bdd09 Mon Sep 17 00:00:00 2001 From: Abel Date: Sat, 23 Oct 2021 12:57:31 +0530 Subject: [PATCH 06/19] Add to test_initialization --- tests/test_initialization.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_initialization.py b/tests/test_initialization.py index 1bba7544..2a3558e3 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, get_grid ) import random from sympl import ( @@ -189,7 +189,7 @@ class ComponentQuantityInitializationTests(unittest.TestCase): RRTMGShortwave, EmanuelConvection, SlabSurface, DcmipInitialConditions, IceSheet, - Instellation + Instellation, SimpleBoundaryLayer ) pair_tests = 20 From c4c7d49d4282e4e589b6ce6351367ab67c55c679 Mon Sep 17 00:00:00 2001 From: Abel Date: Sat, 23 Oct 2021 13:07:26 +0530 Subject: [PATCH 07/19] Conservatio test added --- tests/test_conservation.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/tests/test_conservation.py b/tests/test_conservation.py index a0eff4f5..a02ba017 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.DryConvectiveAdjustment() + + 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): From ea325e9007ca5866da94ac0df263e44bbd11e025 Mon Sep 17 00:00:00 2001 From: Abel Date: Mon, 8 Nov 2021 23:38:30 +0530 Subject: [PATCH 08/19] Cache test and other changes --- .../TestSimpleBoundaryLayer-3d-0.cache | Bin 0 -> 12700 bytes .../TestSimpleBoundaryLayer-3d-1.cache | Bin 0 -> 459312 bytes .../TestSimpleBoundaryLayer-column-0.cache | Bin 0 -> 436 bytes .../TestSimpleBoundaryLayer-column-1.cache | Bin 0 -> 1520 bytes tests/test_components.py | 4 ++-- tests/test_initialization.py | 5 +++-- 6 files changed, 5 insertions(+), 4 deletions(-) create mode 100644 tests/cached_component_output/TestSimpleBoundaryLayer-3d-0.cache create mode 100644 tests/cached_component_output/TestSimpleBoundaryLayer-3d-1.cache create mode 100644 tests/cached_component_output/TestSimpleBoundaryLayer-column-0.cache create mode 100644 tests/cached_component_output/TestSimpleBoundaryLayer-column-1.cache 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 0000000000000000000000000000000000000000..e9082abb4e73a0b0a60bf7a95e968472c76421be GIT binary patch literal 12700 zcmeIyF;2rU6vpwS6%3^VNOXsrAb~nCF?F?!Vi8NOBRP($PQ3_6>(I3~=*W8xNDMGA zFo6D+@@%J89siVk_YccSHMW^_I-Ulb_0DA(cTvYTTHf1*^i#POu}it#*wmPfkB#Yb za(zE8`uNXUrgOP;el|p(Km9&=vSsniEPV(o8wOYA=Tht4sf@MATpvYS{VL_!**<@< z?JwWYa2gm_(fE*wP S$N_SI93ThC0dnB~=D-_=#N)yM literal 0 HcmV?d00001 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 0000000000000000000000000000000000000000..8498c793b136910b85072d34aca38f139af6173f GIT binary patch literal 459312 zcmeI$zfQtH0Kj2TVZi^piM)VO-+-uLad)(7q+(-1Qi{aM#3wU8m7Crb7PrR9nD3Hu zr5AeUv-}1wuZA0~)#!z3GaS81K1z!0IV6tjBhkx5BZX z=4Dc4lUY`#)uPPea=d38?Qq1scGh!eG0m&h7WMVcC>vk*!(te{WSAGlO+cc5SYSTQKg5rTu~z;obxpaVLf13I7sI-mnOpaVLf z13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7s zI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnO zpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf z13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7s zI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnO zpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf z13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7s zI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnO zpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf z13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7s zI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnO zpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mnOpaVLf13I7sI-mo8w*&Ei zV-O%ffB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N z0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+ z009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBly zK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF z5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk z1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs z0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZ zfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&U zAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C7 z2oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N z0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+ z009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBly zK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF z5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk z1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs z0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZ zfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&U zAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C7 z2oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N z0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+ z009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBly zK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF z5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk z1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs z0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZ zfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&U zAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C7 z2oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N z0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7e?-xK%-d@M6} literal 0 HcmV?d00001 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 0000000000000000000000000000000000000000..7f9d2f169cb3bd87363bf2d7b6828be4e7b57058 GIT binary patch literal 436 zcmZ>EabseD04^W}Vl(F?mViVU!R-7z7#l=`05=db1F>jceo;wAd16sYe0gSGN_=ri zQEG89NCm_k5DhYe2g+wH&C4u7l4lA?1Sw(Qig(M*$q7r$DNP0Q*?=Tt{SPo3M1c4l z5FXY^w7eK%bkU%W#j6@L)B_eH6nu3_PG9w$<4hXRoM${)cFZPrD z63^y~s$<~;xe$IVL#D)0sEmX@))hCHROp^2uRAiL0a1hKaH*LQ@j+;AaxGB9nF#Df zujZ4zJdws)P~8hfv8H9++QkRviQPvV+*jt=zC(4$qh7 zY{kOf`)IU}Y|j_oZU8OqB-Lhj;yPp}nS?Z?sBb+Evh{F6Jzx0Ny3o-QJTsrRg?Ig< dujfp4K#!g3!5#6=U-kcD`fPcW8~BqAJOHAmXmS7m literal 0 HcmV?d00001 diff --git a/tests/test_components.py b/tests/test_components.py index 94f17800..77a9dc92 100644 --- a/tests/test_components.py +++ b/tests/test_components.py @@ -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: diff --git a/tests/test_initialization.py b/tests/test_initialization.py index 2a3558e3..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, SimpleBoundaryLayer, 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, SimpleBoundaryLayer + Instellation, SimpleBoundaryLayer, + BucketHydrology ) pair_tests = 20 From f5f5ad90ae4ac8bd7c75a04fa904dc1cb908d2bd Mon Sep 17 00:00:00 2001 From: Abel Date: Mon, 8 Nov 2021 23:43:28 +0530 Subject: [PATCH 09/19] error fix --- tests/test_conservation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_conservation.py b/tests/test_conservation.py index a02ba017..a7249771 100644 --- a/tests/test_conservation.py +++ b/tests/test_conservation.py @@ -276,7 +276,7 @@ def get_quantity_forcing(self, state): class TestSimpleBoundaryLayerConservation(AtmosphereMoistEnthalpyConservation): def get_component_instance(self): - return climt.DryConvectiveAdjustment() + return climt.SimpleBoundaryLayer() def modify_state(self, state): state['eastward_wind'].values[:] = 3. From 4e11b8f6397f109c94df7730fd00b4c116119399 Mon Sep 17 00:00:00 2001 From: Abel Date: Mon, 8 Nov 2021 23:57:44 +0530 Subject: [PATCH 10/19] flake8 error1 fix --- climt/_components/simple_boundary_layer/component.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/climt/_components/simple_boundary_layer/component.py b/climt/_components/simple_boundary_layer/component.py index e18414e2..856fea01 100644 --- a/climt/_components/simple_boundary_layer/component.py +++ b/climt/_components/simple_boundary_layer/component.py @@ -112,8 +112,8 @@ def diffuse_profile(profile, p, p_int, rho, Diff, dt): 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]) + 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)) From 22484a7b577bb817d501ff1a3ac38ffcac35315c Mon Sep 17 00:00:00 2001 From: Abel Date: Tue, 9 Nov 2021 00:01:02 +0530 Subject: [PATCH 11/19] flake8 error1 try 2 --- climt/_components/simple_boundary_layer/component.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/climt/_components/simple_boundary_layer/component.py b/climt/_components/simple_boundary_layer/component.py index 856fea01..7636e68a 100644 --- a/climt/_components/simple_boundary_layer/component.py +++ b/climt/_components/simple_boundary_layer/component.py @@ -112,8 +112,8 @@ def diffuse_profile(profile, p, p_int, rho, Diff, dt): 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]) + 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)) From f99b539709a64f9ef02658e755060874e02f5d09 Mon Sep 17 00:00:00 2001 From: Abel Date: Tue, 9 Nov 2021 00:02:45 +0530 Subject: [PATCH 12/19] flake8 error1 try3 --- climt/_components/simple_boundary_layer/component.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/climt/_components/simple_boundary_layer/component.py b/climt/_components/simple_boundary_layer/component.py index 7636e68a..353d4457 100644 --- a/climt/_components/simple_boundary_layer/component.py +++ b/climt/_components/simple_boundary_layer/component.py @@ -112,8 +112,8 @@ def diffuse_profile(profile, p, p_int, rho, Diff, dt): 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]) + 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)) From 2ea7dc4836b87377b96bd67dd55c14a486d0d8c8 Mon Sep 17 00:00:00 2001 From: Abel Date: Tue, 9 Nov 2021 00:04:17 +0530 Subject: [PATCH 13/19] flake8 error 1 try4 --- climt/_components/simple_boundary_layer/component.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/climt/_components/simple_boundary_layer/component.py b/climt/_components/simple_boundary_layer/component.py index 353d4457..0f6f2ef0 100644 --- a/climt/_components/simple_boundary_layer/component.py +++ b/climt/_components/simple_boundary_layer/component.py @@ -112,8 +112,8 @@ def diffuse_profile(profile, p, p_int, rho, Diff, dt): 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]) + 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)) From 1c83cbf574e37738d97f44db97e3054a29fb93cd Mon Sep 17 00:00:00 2001 From: Abel <67835462+Ai33L@users.noreply.github.com> Date: Tue, 9 Nov 2021 12:01:49 +0530 Subject: [PATCH 14/19] Modify makefile --- climt/_lib/Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From fa11a0c4e661a5e74a212105021413651898e9aa Mon Sep 17 00:00:00 2001 From: Abel Date: Sun, 14 Nov 2021 18:58:09 +0530 Subject: [PATCH 15/19] Component documentation added --- .../simple_boundary_layer/component.py | 27 +++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/climt/_components/simple_boundary_layer/component.py b/climt/_components/simple_boundary_layer/component.py index 0f6f2ef0..968e76d6 100644 --- a/climt/_components/simple_boundary_layer/component.py +++ b/climt/_components/simple_boundary_layer/component.py @@ -176,6 +176,16 @@ def diffuse_profile(profile, p, p_int, rho, Diff, dt): 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': { @@ -253,6 +263,19 @@ class SimpleBoundaryLayer(Stepper): 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 @@ -272,6 +295,10 @@ def _update_constants(self): self._g = get_constant('gravitational_acceleration', 'm s^-2') def array_call(self, state, timestep): + """ + Takes temperature, humidty and wimd profiles for each column and + returns diffused temperature, humidty and wind profiles. + """ num_cols = state['air_temperature'].shape[1] From ca3774966563c777ea207de320ce8fe630a427c3 Mon Sep 17 00:00:00 2001 From: Abel <67835462+Ai33L@users.noreply.github.com> Date: Fri, 21 Jan 2022 11:16:52 +0530 Subject: [PATCH 16/19] Update component.py --- climt/_components/simple_boundary_layer/component.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/climt/_components/simple_boundary_layer/component.py b/climt/_components/simple_boundary_layer/component.py index 968e76d6..7eb4e04f 100644 --- a/climt/_components/simple_boundary_layer/component.py +++ b/climt/_components/simple_boundary_layer/component.py @@ -296,8 +296,8 @@ def _update_constants(self): def array_call(self, state, timestep): """ - Takes temperature, humidty and wimd profiles for each column and - returns diffused temperature, humidty and wind profiles. + Takes temperature, humidty and wind profiles for each column and + returns diffused temperature, humidity and wind profiles. """ num_cols = state['air_temperature'].shape[1] From 17fcdc43ec8ed9903dfb566bae99154f282ce85f Mon Sep 17 00:00:00 2001 From: Abel <67835462+Ai33L@users.noreply.github.com> Date: Fri, 21 Jan 2022 11:25:47 +0530 Subject: [PATCH 17/19] cosmetic change --- .../simple_boundary_layer/component.py | 52 +++++++++---------- 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/climt/_components/simple_boundary_layer/component.py b/climt/_components/simple_boundary_layer/component.py index 7eb4e04f..6820c00c 100644 --- a/climt/_components/simple_boundary_layer/component.py +++ b/climt/_components/simple_boundary_layer/component.py @@ -6,14 +6,14 @@ @jit(nopython=True, parallel=True) -def Parallel_columns(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): +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 @@ -310,23 +310,23 @@ def array_call(self, state, timestep): self.diagnostic_properties, state, self.input_properties ) - Parallel_columns(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()) + 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 From f5fd09a82616e9e654e92a45255eb53ef084d7f8 Mon Sep 17 00:00:00 2001 From: Abel <67835462+Ai33L@users.noreply.github.com> Date: Sat, 21 May 2022 10:25:04 +0530 Subject: [PATCH 18/19] Update _gfs_dynamics.pyx --- climt/_components/gfs/_gfs_dynamics.pyx | 74 +++++++++++++++++++++++++ 1 file changed, 74 insertions(+) diff --git a/climt/_components/gfs/_gfs_dynamics.pyx b/climt/_components/gfs/_gfs_dynamics.pyx index e668b457..a80b84ee 100644 --- a/climt/_components/gfs/_gfs_dynamics.pyx +++ b/climt/_components/gfs/_gfs_dynamics.pyx @@ -1,6 +1,7 @@ cimport numpy as cnp import numpy as np import cython +import copy # Typedef for function pointer returning void and taking no arguments (for now) @@ -400,6 +401,8 @@ def init_model( @cython.boundscheck(False) def init_spectral_arrays(spectral_dim, num_levs, num_tracers): + #print('Spectral arrays init PRINT') + global pyTracerSpec, pyTracerSpecTend, pyTopoSpec, \ pyLnPsSpec, pyLnPsSpecTend,\ pyDissSpec, pyDmpProf, pyDiffProf,\ @@ -455,11 +458,66 @@ def init_spectral_arrays(spectral_dim, num_levs, num_tracers): &pyTracerSpecTend[0,0,0], &pyLnPsSpecTend[0]) +##################################################################### +def reconvert(A,B): + A[:] = 0 + A[:] = (A + B)[:] + +def reinit_spectral_arrays(data): + + #print('Spectral arrays REinit PRINT') + + global pyTracerSpec, pyTracerSpecTend, pyTopoSpec, \ + pyLnPsSpec, pyLnPsSpecTend,\ + pyDissSpec, pyDmpProf, pyDiffProf,\ + pyVrtSpec, pyVrtSpecTend, pyDivSpec, pyDivSpecTend,\ + pyVirtTempSpec, pyVirtTempSpecTend + + print(np.array(pyVirtTempSpec).mean()) + print(pyVirtTempSpec.shape) + reconvert(pyVirtTempSpec,data[0]) + print(pyVirtTempSpec.shape) + print(np.array(pyVirtTempSpec).mean()) + reconvert(pyDivSpec,data[1]) + reconvert(pyVrtSpec,data[2]) + reconvert(pyTracerSpec,data[3]) + reconvert(pyLnPsSpec,data[4]) + reconvert(pyTopoSpec,data[5]) + reconvert(pyDissSpec,data[6]) + reconvert(pyDiffProf,data[7]) + reconvert(pyDmpProf,data[8]) + + #reconvert(pyVrtSpecTend,data[9]) + #reconvert(pyDivSpecTend,data[10]) + #reconvert(pyVirtTempSpecTend,data[11]) + #reconvert(pyTracerSpecTend,data[12]) + #reconvert(pyLnPsSpecTend,data[13]) + + + #gfs_initialise_spectral_arrays( + # &pyVrtSpec[0,0], + # &pyDivSpec[0,0], + # &pyVirtTempSpec[0,0], + # &pyTracerSpec[0,0,0], + # &pyTopoSpec[0], + # &pyLnPsSpec[0], + # &pyDissSpec[0], + # &pyDmpProf[0], + # &pyDiffProf[0]) + + #gfs_initialise_spectral_physics_arrays( + # &pyVrtSpecTend[0,0], + # &pyDivSpecTend[0,0], + # &pyVirtTempSpecTend[0,0], + # &pyTracerSpecTend[0,0,0], + # &pyLnPsSpecTend[0]) # Create the grid arrays (defined in grid_data.f90) def init_grid_arrays(num_lons, num_lats, num_levs, num_tracers): + #print('Grid arrays init PRINT') + global pyDlnpdtg, pyEtaDotg, pyPhis, pyDPhisdx, pyDPhisdy, \ pyDlnpsdt, pyPwat, tempVrtTend, tempDivTend @@ -561,13 +619,16 @@ def assign_grid_arrays( def update_spectral_arrays(): + print('Convert to spectral PRINT') gfs_convert_to_spectral() def take_one_step(): + #print('Step PRINT') gfs_take_one_step() def convert_to_grid(): + #print('Convert to grid PRINT') gfs_convert_to_grid() def calculate_pressure(): @@ -585,3 +646,16 @@ def shut_down_model(): #Remember to set time to zero! gfs_finalise() gfs_set_model_time(&zero_model_time) + + +def get_spectral_arrays(): + + return np.asarray(pyVirtTempSpec, dtype=np.cfloat),np.asarray(pyDivSpec, dtype=np.cfloat),np.asarray(pyVrtSpec, dtype=np.cfloat),np.asarray(pyTracerSpec, dtype=np.cfloat),\ + np.asarray(pyLnPsSpec, dtype=np.cfloat),np.asarray(pyTopoSpec, dtype=np.cfloat),np.asarray(pyDissSpec, dtype=np.cfloat),np.asarray(pyDiffProf, dtype=np.cfloat),np.asarray(pyDmpProf, dtype=np.cfloat),\ + np.asarray(pyVrtSpecTend, dtype=np.cfloat),np.asarray(pyDivSpecTend, dtype=np.cfloat),np.asarray(pyVirtTempSpecTend, dtype=np.cfloat),np.asarray(pyTracerSpecTend, dtype=np.cfloat),np.asarray(pyLnPsSpecTend, dtype=np.cfloat) + +def get_spectral_arrays2(): + + return pyVirtTempSpec,pyDivSpec,pyVrtSpec,pyTracerSpec,\ + pyLnPsSpec,pyTopoSpec,pyDissSpec,pyDiffProf,pyDmpProf,\ + pyVrtSpecTend,pyDivSpecTend,pyVirtTempSpecTend,pyTracerSpecTend,pyLnPsSpecTend From 252f578922007df65d4e07d81979641a8daa83ad Mon Sep 17 00:00:00 2001 From: Abel <67835462+Ai33L@users.noreply.github.com> Date: Tue, 24 May 2022 14:21:02 +0200 Subject: [PATCH 19/19] Update _gfs_dynamics.pyx --- climt/_components/gfs/_gfs_dynamics.pyx | 74 ------------------------- 1 file changed, 74 deletions(-) diff --git a/climt/_components/gfs/_gfs_dynamics.pyx b/climt/_components/gfs/_gfs_dynamics.pyx index a80b84ee..e668b457 100644 --- a/climt/_components/gfs/_gfs_dynamics.pyx +++ b/climt/_components/gfs/_gfs_dynamics.pyx @@ -1,7 +1,6 @@ cimport numpy as cnp import numpy as np import cython -import copy # Typedef for function pointer returning void and taking no arguments (for now) @@ -401,8 +400,6 @@ def init_model( @cython.boundscheck(False) def init_spectral_arrays(spectral_dim, num_levs, num_tracers): - #print('Spectral arrays init PRINT') - global pyTracerSpec, pyTracerSpecTend, pyTopoSpec, \ pyLnPsSpec, pyLnPsSpecTend,\ pyDissSpec, pyDmpProf, pyDiffProf,\ @@ -458,66 +455,11 @@ def init_spectral_arrays(spectral_dim, num_levs, num_tracers): &pyTracerSpecTend[0,0,0], &pyLnPsSpecTend[0]) -##################################################################### -def reconvert(A,B): - A[:] = 0 - A[:] = (A + B)[:] - -def reinit_spectral_arrays(data): - - #print('Spectral arrays REinit PRINT') - - global pyTracerSpec, pyTracerSpecTend, pyTopoSpec, \ - pyLnPsSpec, pyLnPsSpecTend,\ - pyDissSpec, pyDmpProf, pyDiffProf,\ - pyVrtSpec, pyVrtSpecTend, pyDivSpec, pyDivSpecTend,\ - pyVirtTempSpec, pyVirtTempSpecTend - - print(np.array(pyVirtTempSpec).mean()) - print(pyVirtTempSpec.shape) - reconvert(pyVirtTempSpec,data[0]) - print(pyVirtTempSpec.shape) - print(np.array(pyVirtTempSpec).mean()) - reconvert(pyDivSpec,data[1]) - reconvert(pyVrtSpec,data[2]) - reconvert(pyTracerSpec,data[3]) - reconvert(pyLnPsSpec,data[4]) - reconvert(pyTopoSpec,data[5]) - reconvert(pyDissSpec,data[6]) - reconvert(pyDiffProf,data[7]) - reconvert(pyDmpProf,data[8]) - - #reconvert(pyVrtSpecTend,data[9]) - #reconvert(pyDivSpecTend,data[10]) - #reconvert(pyVirtTempSpecTend,data[11]) - #reconvert(pyTracerSpecTend,data[12]) - #reconvert(pyLnPsSpecTend,data[13]) - - - #gfs_initialise_spectral_arrays( - # &pyVrtSpec[0,0], - # &pyDivSpec[0,0], - # &pyVirtTempSpec[0,0], - # &pyTracerSpec[0,0,0], - # &pyTopoSpec[0], - # &pyLnPsSpec[0], - # &pyDissSpec[0], - # &pyDmpProf[0], - # &pyDiffProf[0]) - - #gfs_initialise_spectral_physics_arrays( - # &pyVrtSpecTend[0,0], - # &pyDivSpecTend[0,0], - # &pyVirtTempSpecTend[0,0], - # &pyTracerSpecTend[0,0,0], - # &pyLnPsSpecTend[0]) # Create the grid arrays (defined in grid_data.f90) def init_grid_arrays(num_lons, num_lats, num_levs, num_tracers): - #print('Grid arrays init PRINT') - global pyDlnpdtg, pyEtaDotg, pyPhis, pyDPhisdx, pyDPhisdy, \ pyDlnpsdt, pyPwat, tempVrtTend, tempDivTend @@ -619,16 +561,13 @@ def assign_grid_arrays( def update_spectral_arrays(): - print('Convert to spectral PRINT') gfs_convert_to_spectral() def take_one_step(): - #print('Step PRINT') gfs_take_one_step() def convert_to_grid(): - #print('Convert to grid PRINT') gfs_convert_to_grid() def calculate_pressure(): @@ -646,16 +585,3 @@ def shut_down_model(): #Remember to set time to zero! gfs_finalise() gfs_set_model_time(&zero_model_time) - - -def get_spectral_arrays(): - - return np.asarray(pyVirtTempSpec, dtype=np.cfloat),np.asarray(pyDivSpec, dtype=np.cfloat),np.asarray(pyVrtSpec, dtype=np.cfloat),np.asarray(pyTracerSpec, dtype=np.cfloat),\ - np.asarray(pyLnPsSpec, dtype=np.cfloat),np.asarray(pyTopoSpec, dtype=np.cfloat),np.asarray(pyDissSpec, dtype=np.cfloat),np.asarray(pyDiffProf, dtype=np.cfloat),np.asarray(pyDmpProf, dtype=np.cfloat),\ - np.asarray(pyVrtSpecTend, dtype=np.cfloat),np.asarray(pyDivSpecTend, dtype=np.cfloat),np.asarray(pyVirtTempSpecTend, dtype=np.cfloat),np.asarray(pyTracerSpecTend, dtype=np.cfloat),np.asarray(pyLnPsSpecTend, dtype=np.cfloat) - -def get_spectral_arrays2(): - - return pyVirtTempSpec,pyDivSpec,pyVrtSpec,pyTracerSpec,\ - pyLnPsSpec,pyTopoSpec,pyDissSpec,pyDiffProf,pyDmpProf,\ - pyVrtSpecTend,pyDivSpecTend,pyVirtTempSpecTend,pyTracerSpecTend,pyLnPsSpecTend