Update climlab/mcp_output/mcp_plugin/mcp_service.py

climlab/mcp_output/mcp_plugin/mcp_service.py

@@ -1,400 +1,739 @@
+import os
+import sys
+import numpy as np
+
+# Add the local source directory to sys.path
+source_path = os.path.join(os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))), "source")
+if source_path not in sys.path:
+    sys.path.insert(0, source_path)
+
 from fastmcp import FastMCP
 
-# 创建 FastMCP 服务应用
+# Import core modules from climlab using the correct API
+import climlab
+from climlab.domain.domain import single_column, zonal_mean_column, box_model_domain, zonal_mean_surface
+from climlab.domain.initial import column_state, surface_state
+from climlab.domain.field import Field, global_mean
+from climlab.model.ebm import EBM, EBM_annual, EBM_seasonal
+from climlab.model.column import GreyRadiationModel, RadiativeConvectiveModel, BandRCModel
+from climlab.radiation.insolation import P2Insolation, AnnualMeanInsolation, DailyInsolation, FixedInsolation
+from climlab.radiation.aplusbt import AplusBT
+from climlab.surface.albedo import ConstantAlbedo, P2Albedo, StepFunctionAlbedo
+from climlab.solar.insolation import daily_insolation, annual_insolation
+from climlab import constants as const
+
+# Create the FastMCP service application
 mcp = FastMCP("climlab_service")
 
-@mcp.tool(name="get_version", description="获取 climlab 的版本信息")
-def get_version() -> dict:
+
+# ==================== Domain Tools ====================
+
+@mcp.tool(name="create_column_state", description="Create a column state for climate modeling with atmospheric and surface temperatures")
+def create_column_state_tool(num_lev: int = 30, num_lat: int = 1, water_depth: float = 1.0) -> dict:
     """
-    获取 climlab 的版本信息。
+    Create a column state for climate modeling.
 
-    Returns:
-    - dict: 包含版本号的字典。
+    :param num_lev: The number of vertical levels (default: 30).
+    :param num_lat: The number of latitude points (default: 1).
+    :param water_depth: Depth of the slab ocean in meters (default: 1.0).
+    :return: A dictionary containing success status and the state object.
     """
     try:
-        from importlib import metadata
-        version = metadata.version("climlab")
-        return {"success": True, "version": version}
+        state = column_state(num_lev=num_lev, num_lat=num_lat, water_depth=water_depth)
+        result = {
+            "Ts": {
+                "shape": list(state['Ts'].shape),
+                "values": state['Ts'].tolist(),
+                "domain_type": state['Ts'].domain.domain_type,
+            },
+            "Tatm": {
+                "shape": list(state['Tatm'].shape),
+                "min": float(state['Tatm'].min()),
+                "max": float(state['Tatm'].max()),
+                "domain_type": state['Tatm'].domain.domain_type,
+            },
+            "num_lev": num_lev,
+            "num_lat": num_lat,
+        }
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
+
 
-@mcp.tool(name="create_ebm_model", description="创建一个 1D 能量平衡模型 (EBM)")
-def create_ebm_model() -> dict:
+@mcp.tool(name="create_surface_state", description="Create a surface state for EBM with latitude-dependent initial temperature")
+def create_surface_state_tool(num_lat: int = 90, water_depth: float = 10.0, T0: float = 12.0, T2: float = -40.0) -> dict:
     """
-    创建一个 1D 能量平衡模型。
+    Create a surface state for Energy Balance Models.
 
-    Returns:
-    - dict: 包含模型信息的字典。
+    :param num_lat: Number of latitude points (default: 90).
+    :param water_depth: Depth of the slab ocean in meters (default: 10.0).
+    :param T0: Global-mean initial temperature in °C (default: 12.0).
+    :param T2: 2nd Legendre coefficient for equator-to-pole gradient (default: -40.0).
+    :return: A dictionary containing success status and state information.
     """
     try:
-        from climlab import EBM
-        model = EBM()
-        return {
-            "success": True,
-            "model": str(model),
-            "state": model.state,
-            "diagnostics": list(model.diagnostics.keys())
+        state = surface_state(num_lat=num_lat, water_depth=water_depth, T0=T0, T2=T2)
+        Ts = state['Ts']
+        lat = Ts.domain.axes['lat'].points
+        result = {
+            "shape": list(Ts.shape),
+            "latitudes": lat.tolist(),
+            "temperatures": np.squeeze(Ts).tolist(),
+            "global_mean_temperature": float(global_mean(Ts)),
+            "equator_temperature": float(Ts[len(Ts)//2]) if len(Ts) > 1 else float(Ts[0]),
+            "pole_temperature": float(Ts[0]),
         }
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
 
-@mcp.tool(name="integrate_model", description="运行模型指定的时间步长")
-def integrate_model(model, years: float) -> dict:
-    """
-    运行模型指定的时间步长。
 
-    Parameters:
-    - model: climlab 模型实例。
-    - years: 运行的年数。
+# ==================== EBM Tools ====================
+
+@mcp.tool(name="create_ebm", description="Create an Energy Balance Model with diffusive heat transport")
+def create_ebm_tool(num_lat: int = 90, S0: float = 1365.2, A: float = 210.0, B: float = 2.0, 
+                   D: float = 0.555, a0: float = 0.3, a2: float = 0.078, ai: float = 0.62,
+                   Tf: float = -10.0, water_depth: float = 10.0) -> dict:
+    """
+    Create an Energy Balance Model (EBM) with specified parameters.
 
-    Returns:
-    - dict: 包含运行状态和结果的字典。
+    :param num_lat: Number of latitude points (default: 90).
+    :param S0: Solar constant in W/m² (default: 1365.2).
+    :param A: OLR parameter A in W/m² (default: 210.0).
+    :param B: OLR parameter B in W/m²/°C (default: 2.0).
+    :param D: Diffusion parameter in W/m²/°C (default: 0.555).
+    :param a0: Base albedo coefficient (default: 0.3).
+    :param a2: Second Legendre polynomial coefficient (default: 0.078).
+    :param ai: Ice albedo value (default: 0.62).
+    :param Tf: Freezing temperature in °C (default: -10.0).
+    :param water_depth: Slab ocean depth in meters (default: 10.0).
+    :return: A dictionary containing success status and model info.
     """
     try:
-        model.integrate_years(years)
-        return {
-            "success": True,
-            "state": model.state,
-            "diagnostics": model.diagnostics
+        model = EBM(num_lat=num_lat, S0=S0, A=A, B=B, D=D, 
+                   a0=a0, a2=a2, ai=ai, Tf=Tf, water_depth=water_depth)
+        lat = model.lat
+        result = {
+            "model_type": "EBM",
+            "parameters": {
+                "num_lat": num_lat,
+                "S0": S0,
+                "A": A,
+                "B": B,
+                "D": D,
+                "a0": a0,
+                "a2": a2,
+                "ai": ai,
+                "Tf": Tf,
+                "water_depth": water_depth,
+            },
+            "state_variables": list(model.state.keys()),
+            "subprocesses": list(model.subprocess.keys()),
+            "latitudes": lat.tolist(),
+            "initial_Ts": np.squeeze(model.Ts).tolist(),
+            "timestep_seconds": float(model.timestep),
         }
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
 
-@mcp.tool(name="list_available_models", description="列出 climlab 中的所有可用模型")
-def list_available_models() -> dict:
+
+@mcp.tool(name="run_ebm_to_equilibrium", description="Run an Energy Balance Model to equilibrium and return climate statistics")
+def run_ebm_to_equilibrium_tool(num_lat: int = 90, num_years: float = 5.0, 
+                                S0: float = 1365.2, A: float = 210.0, B: float = 2.0,
+                                D: float = 0.555, a0: float = 0.3, a2: float = 0.078, 
+                                ai: float = 0.62, Tf: float = -10.0) -> dict:
     """
-    列出 climlab 中的所有可用模型。
+    Run an Energy Balance Model to equilibrium.
 
-    Returns:
-    - dict: 包含模型名称的字典。
+    :param num_lat: Number of latitude points (default: 90).
+    :param num_years: Number of years to integrate (default: 5.0).
+    :param S0: Solar constant in W/m² (default: 1365.2).
+    :param A: OLR parameter A in W/m² (default: 210.0).
+    :param B: OLR parameter B in W/m²/°C (default: 2.0).
+    :param D: Diffusion parameter in W/m²/°C (default: 0.555).
+    :param a0: Base albedo coefficient (default: 0.3).
+    :param a2: Second Legendre polynomial coefficient (default: 0.078).
+    :param ai: Ice albedo value (default: 0.62).
+    :param Tf: Freezing temperature in °C (default: -10.0).
+    :return: A dictionary containing equilibrium results.
     """
     try:
-        models = ["EBM", "RCM", "GreyRadiationModel", "BandRCModel"]
-        return {"success": True, "models": models}
+        model = EBM(num_lat=num_lat, S0=S0, A=A, B=B, D=D, a0=a0, a2=a2, ai=ai, Tf=Tf)
+        model.integrate_years(num_years)
+        Ts = model.Ts
+        lat = model.lat
+        
+        # Get heat transport if available
+        heat_transport = None
+        if hasattr(model.subprocess.get('diffusion', None), 'heat_transport'):
+            heat_transport = model.subprocess['diffusion'].heat_transport.tolist()
+        
+        # Get OLR
+        OLR = None
+        if hasattr(model, 'OLR'):
+            OLR = np.squeeze(model.OLR).tolist()
+        
+        # Get ASR
+        ASR = None
+        if hasattr(model, 'ASR'):
+            ASR = np.squeeze(model.ASR).tolist()
+        
+        # Find ice edge latitude
+        ice_lat = None
+        if hasattr(model, 'icelat'):
+            ice_lat = float(model.icelat) if np.isfinite(model.icelat) else None
+        
+        result = {
+            "global_mean_temperature": float(global_mean(Ts)),
+            "max_temperature": float(Ts.max()),
+            "min_temperature": float(Ts.min()),
+            "equator_temperature": float(Ts[len(Ts)//2]) if len(Ts) > 1 else float(Ts[0]),
+            "pole_temperature": float(Ts[0]),
+            "latitudes": lat.tolist(),
+            "temperature_profile": np.squeeze(Ts).tolist(),
+            "ice_edge_latitude": ice_lat,
+            "OLR": OLR,
+            "ASR": ASR,
+            "heat_transport_PW": heat_transport,
+            "integrated_years": num_years,
+            "energy_balance": float(global_mean(model.ASR - model.OLR)) if hasattr(model, 'ASR') else None,
+        }
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
 
-@mcp.tool(name="create_model", description="根据名称创建 climlab 模型")
-def create_model(model_name: str) -> dict:
-    """
-    根据名称创建 climlab 模型。
 
-    Parameters:
-    - model_name: 模型名称 (例如 'EBM', 'RCM')。
+@mcp.tool(name="run_seasonal_ebm", description="Run a seasonal Energy Balance Model with time-varying insolation")
+def run_seasonal_ebm_tool(num_lat: int = 90, num_years: float = 5.0, 
+                          water_depth: float = 10.0) -> dict:
+    """
+    Run a seasonal Energy Balance Model with realistic seasonally varying insolation.
 
-    Returns:
-    - dict: 包含模型信息的字典。
+    :param num_lat: Number of latitude points (default: 90).
+    :param num_years: Number of years to integrate (default: 5.0).
+    :param water_depth: Slab ocean depth in meters (default: 10.0).
+    :return: A dictionary containing seasonal climate results.
     """
     try:
-        from climlab import EBM, RadiativeConvectiveModel, GreyRadiationModel, BandRCModel
-        models = {
-            "EBM": EBM,
-            "RCM": RadiativeConvectiveModel,
-            "GreyRadiationModel": GreyRadiationModel,
-            "BandRCModel": BandRCModel
-        }
-        if model_name not in models:
-            return {"success": False, "error": f"Model '{model_name}' not found."}
-
-        model = models[model_name]()
-        return {
-            "success": True,
-            "model": str(model),
-            "state": model.state,
-            "diagnostics": list(model.diagnostics.keys())
+        model = EBM_seasonal(num_lat=num_lat, water_depth=water_depth)
+        model.integrate_years(num_years)
+        Ts = model.Ts
+        lat = model.lat
+        
+        result = {
+            "global_mean_temperature": float(global_mean(Ts)),
+            "max_temperature": float(Ts.max()),
+            "min_temperature": float(Ts.min()),
+            "latitudes": lat.tolist(),
+            "temperature_profile": np.squeeze(Ts).tolist(),
+            "integrated_years": num_years,
+            "model_type": "EBM_seasonal",
         }
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
 
-@mcp.tool(name="create_convective_adjustment", description="创建一个对流调整模型")
-def create_convective_adjustment(adj_lapse_rate: float) -> dict:
-    """
-    创建一个对流调整模型。
 
-    Parameters:
-    - adj_lapse_rate: 调整的递减率 (单位: K/km)。
+# ==================== Insolation Tools ====================
 
-    Returns:
-    - dict: 包含模型信息的字典。
+@mcp.tool(name="calculate_daily_insolation", description="Calculate daily average insolation given latitude and day of year")
+def calculate_daily_insolation_tool(lat: float, day: int, S0: float = 1365.2,
+                                    ecc: float = 0.017236, obliquity: float = 23.446, 
+                                    long_peri: float = 281.37) -> dict:
+    """
+    Calculate daily average insolation at given latitude and day of year.
+
+    :param lat: Latitude in degrees (-90 to 90).
+    :param day: Calendar day (1-365), day 1 is January 1st.
+    :param S0: Solar constant in W/m² (default: 1365.2).
+    :param ecc: Orbital eccentricity (default: 0.017236 for present-day).
+    :param obliquity: Obliquity angle in degrees (default: 23.446).
+    :param long_peri: Longitude of perihelion in degrees (default: 281.37).
+    :return: A dictionary containing daily average insolation in W/m².
     """
     try:
-        from climlab.convection.convadj import ConvectiveAdjustment
-        model = ConvectiveAdjustment(adj_lapse_rate=adj_lapse_rate)
-        return {
-            "success": True,
-            "model": str(model),
-            "state": model.state
+        orb = {'ecc': ecc, 'obliquity': obliquity, 'long_peri': long_peri}
+        insolation = daily_insolation(lat=lat, day=day, orb=orb, S0=S0)
+        
+        result = {
+            "latitude": lat,
+            "day_of_year": day,
+            "insolation_W_m2": float(insolation),
+            "solar_constant": S0,
+            "orbital_parameters": orb,
         }
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
+
 
-@mcp.tool(name="create_emanuel_convection", description="创建一个Emanuel对流模型")
-def create_emanuel_convection() -> dict:
+@mcp.tool(name="calculate_annual_mean_insolation", description="Calculate annual mean insolation as a function of latitude")
+def calculate_annual_mean_insolation_tool(num_lat: int = 90, S0: float = 1365.2,
+                                          ecc: float = 0.017236, obliquity: float = 23.446,
+                                          long_peri: float = 281.37) -> dict:
     """
-    创建一个Emanuel对流模型。
+    Calculate annual mean insolation for a range of latitudes.
 
-    Returns:
-    - dict: 包含模型信息的字典。
+    :param num_lat: Number of latitude points (default: 90).
+    :param S0: Solar constant in W/m² (default: 1365.2).
+    :param ecc: Orbital eccentricity (default: 0.017236).
+    :param obliquity: Obliquity angle in degrees (default: 23.446).
+    :param long_peri: Longitude of perihelion in degrees (default: 281.37).
+    :return: A dictionary containing annual mean insolation profile.
     """
     try:
-        from climlab.convection.emanuel_convection import EmanuelConvection
-        model = EmanuelConvection()
-        return {
-            "success": True,
-            "model": str(model),
-            "state": model.state
+        lat = np.linspace(-90, 90, num_lat)
+        orb = {'ecc': ecc, 'obliquity': obliquity, 'long_peri': long_peri}
+        insolation = annual_insolation(lat=lat, orb=orb, S0=S0)
+        
+        result = {
+            "latitudes": lat.tolist(),
+            "insolation_W_m2": insolation.tolist() if hasattr(insolation, 'tolist') else [float(insolation)],
+            "global_mean_insolation": float(np.mean(insolation * np.cos(np.deg2rad(lat))) / np.mean(np.cos(np.deg2rad(lat)))),
+            "equator_insolation": float(insolation[num_lat//2]) if num_lat > 1 else float(insolation),
+            "pole_insolation": float(insolation[0]) if num_lat > 1 else float(insolation),
+            "solar_constant": S0,
         }
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
 
-@mcp.tool(name="create_grey_gas_model", description="创建一个灰气辐射模型")
-def create_grey_gas_model(absorptivity: float) -> dict:
-    """
-    创建一个灰气辐射模型。
 
-    Parameters:
-    - absorptivity: 吸收率。
+@mcp.tool(name="calculate_insolation_seasonal_cycle", description="Calculate the seasonal cycle of insolation at a single latitude")
+def calculate_insolation_seasonal_cycle_tool(lat: float, S0: float = 1365.2,
+                                             ecc: float = 0.017236, obliquity: float = 23.446,
+                                             long_peri: float = 281.37) -> dict:
+    """
+    Calculate the seasonal cycle of daily insolation at a given latitude.
 
-    Returns:
-    - dict: 包含模型信息的字典。
+    :param lat: Latitude in degrees (-90 to 90).
+    :param S0: Solar constant in W/m² (default: 1365.2).
+    :param ecc: Orbital eccentricity (default: 0.017236).
+    :param obliquity: Obliquity angle in degrees (default: 23.446).
+    :param long_peri: Longitude of perihelion in degrees (default: 281.37).
+    :return: A dictionary containing the seasonal cycle of insolation.
     """
     try:
-        from climlab.radiation.greygas import GreyGas
-        model = GreyGas(absorptivity=absorptivity)
-        return {
-            "success": True,
-            "model": str(model),
-            "diagnostics": model.diagnostics
+        days = np.arange(1, 366)
+        orb = {'ecc': ecc, 'obliquity': obliquity, 'long_peri': long_peri}
+        insolation = np.array([float(daily_insolation(lat=lat, day=d, orb=orb, S0=S0)) for d in days])
+        
+        result = {
+            "latitude": lat,
+            "days": days.tolist(),
+            "insolation_W_m2": insolation.tolist(),
+            "annual_mean": float(np.mean(insolation)),
+            "max_insolation": float(np.max(insolation)),
+            "min_insolation": float(np.min(insolation)),
+            "day_of_max": int(days[np.argmax(insolation)]),
+            "day_of_min": int(days[np.argmin(insolation)]),
         }
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
+
 
-@mcp.tool(name="create_grey_radiation_model", description="创建一个灰气辐射列模型")
-def create_grey_radiation_model() -> dict:
+# ==================== Radiation Tools ====================
+
+@mcp.tool(name="calculate_olr_aplusbt", description="Calculate Outgoing Longwave Radiation using A+BT parameterization")
+def calculate_olr_aplusbt_tool(temperature: float, A: float = 210.0, B: float = 2.0) -> dict:
     """
-    创建一个灰气辐射列模型。
+    Calculate Outgoing Longwave Radiation using the simple A+BT parameterization.
 
-    Returns:
-    - dict: 包含模型信息的字典。
+    :param temperature: Surface temperature in °C.
+    :param A: OLR parameter A in W/m² (default: 210.0).
+    :param B: OLR parameter B in W/m²/°C (default: 2.0).
+    :return: A dictionary containing the OLR value.
     """
     try:
-        from climlab.model.column import GreyRadiationModel
-        model = GreyRadiationModel()
-        return {
-            "success": True,
-            "model": str(model),
-            "state": model.state,
-            "diagnostics": list(model.diagnostics.keys())
+        OLR = A + B * temperature
+        result = {
+            "temperature_C": temperature,
+            "OLR_W_m2": float(OLR),
+            "A": A,
+            "B": B,
+            "description": f"OLR = {A} + {B} * T = {OLR:.2f} W/m²",
         }
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
+
 
-@mcp.tool(name="create_advection_diffusion", description="创建一个对流扩散模型")
-def create_advection_diffusion() -> dict:
+@mcp.tool(name="calculate_stefan_boltzmann_radiation", description="Calculate blackbody radiation using Stefan-Boltzmann law")
+def calculate_stefan_boltzmann_radiation_tool(temperature_K: float, emissivity: float = 1.0) -> dict:
     """
-    创建一个对流扩散模型。
+    Calculate blackbody (or greybody) radiation using Stefan-Boltzmann law.
 
-    Returns:
-    - dict: 包含模型信息的字典。
+    :param temperature_K: Temperature in Kelvin.
+    :param emissivity: Surface emissivity (0-1, default: 1.0 for blackbody).
+    :return: A dictionary containing the radiation flux.
     """
     try:
-        from climlab.dynamics.advection_diffusion import AdvectionDiffusion
-        model = AdvectionDiffusion()
-        return {
-            "success": True,
-            "model": str(model),
-            "diagnostics": model.diagnostics
+        sigma = const.sigma  # Stefan-Boltzmann constant
+        radiation = emissivity * sigma * temperature_K**4
+        
+        result = {
+            "temperature_K": temperature_K,
+            "temperature_C": temperature_K - 273.15,
+            "emissivity": emissivity,
+            "radiation_W_m2": float(radiation),
+            "stefan_boltzmann_constant": sigma,
         }
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
+
 
-@mcp.tool(name="create_surface_flux", description="创建一个地表湍流热量和水分通量模型")
-def create_surface_flux() -> dict:
+# ==================== Albedo Tools ====================
+
+@mcp.tool(name="calculate_p2_albedo", description="Calculate latitude-dependent albedo using P2 Legendre polynomial")
+def calculate_p2_albedo_tool(num_lat: int = 90, a0: float = 0.33, a2: float = 0.25) -> dict:
     """
-    创建一个地表湍流热量和水分通量模型。
+    Calculate latitude-dependent albedo using second-order Legendre polynomial.
+    α(φ) = a0 + a2 * P2(sin(φ)), where P2(x) = (3x² - 1)/2
 
-    Returns:
-    - dict: 包含模型信息的字典。
+    :param num_lat: Number of latitude points (default: 90).
+    :param a0: Base albedo parameter (default: 0.33).
+    :param a2: Second Legendre coefficient (default: 0.25).
+    :return: A dictionary containing the albedo profile.
     """
     try:
-        from climlab.surface.turbulent import SensibleHeatFlux, LatentHeatFlux
-        sensible_flux = SensibleHeatFlux()
-        latent_flux = LatentHeatFlux()
-        return {
-            "success": True,
-            "sensible_flux": str(sensible_flux),
-            "latent_flux": str(latent_flux)
+        from climlab.utils.legendre import P2
+        lat = np.linspace(-90, 90, num_lat)
+        phi = np.deg2rad(lat)
+        albedo = a0 + a2 * P2(np.sin(phi))
+        
+        result = {
+            "latitudes": lat.tolist(),
+            "albedo": albedo.tolist(),
+            "global_mean_albedo": float(np.average(albedo, weights=np.cos(phi))),
+            "equator_albedo": float(albedo[num_lat//2]),
+            "pole_albedo": float(albedo[0]),
+            "a0": a0,
+            "a2": a2,
         }
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
+
 
-@mcp.tool(name="create_time_dependent_process", description="创建一个时间相关的过程模型")
-def create_time_dependent_process() -> dict:
+@mcp.tool(name="calculate_ice_albedo_feedback", description="Calculate ice-albedo feedback with temperature-dependent ice line")
+def calculate_ice_albedo_feedback_tool(temperatures: list, Tf: float = -10.0,
+                                       a0: float = 0.3, a2: float = 0.078, ai: float = 0.62) -> dict:
     """
-    创建一个时间相关的过程模型。
+    Calculate albedo with ice-albedo feedback based on temperature.
 
-    Returns:
-    - dict: 包含模型信息的字典。
+    :param temperatures: List of temperatures in °C at each latitude.
+    :param Tf: Freezing temperature threshold in °C (default: -10.0).
+    :param a0: Unfrozen base albedo (default: 0.3).
+    :param a2: Second Legendre coefficient for unfrozen albedo (default: 0.078).
+    :param ai: Ice albedo value (default: 0.62).
+    :return: A dictionary containing the albedo with ice-albedo feedback.
     """
     try:
-        from climlab.process.process import TimeDependentProcess
-        model = TimeDependentProcess()
-        return {
-            "success": True,
-            "model": str(model),
-            "state": model.state
+        from climlab.utils.legendre import P2
+        T = np.array(temperatures)
+        num_lat = len(T)
+        lat = np.linspace(-90, 90, num_lat)
+        phi = np.deg2rad(lat)
+        
+        # Calculate unfrozen albedo using P2
+        albedo_unfrozen = a0 + a2 * P2(np.sin(phi))
+        
+        # Apply ice-albedo feedback: ice where T < Tf
+        albedo = np.where(T < Tf, ai, albedo_unfrozen)
+        
+        # Find ice edge latitude
+        ice_mask = T < Tf
+        if np.any(ice_mask) and not np.all(ice_mask):
+            # Find the latitude where ice starts
+            ice_lats = np.abs(lat[ice_mask])
+            ice_edge = float(np.min(ice_lats))
+        elif np.all(ice_mask):
+            ice_edge = 0.0  # Snowball Earth
+        else:
+            ice_edge = 90.0  # No ice
+        
+        result = {
+            "latitudes": lat.tolist(),
+            "temperatures": temperatures,
+            "albedo": albedo.tolist(),
+            "ice_mask": ice_mask.tolist(),
+            "ice_edge_latitude": ice_edge,
+            "global_mean_albedo": float(np.average(albedo, weights=np.cos(phi))),
+            "parameters": {"Tf": Tf, "a0": a0, "a2": a2, "ai": ai},
         }
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
+
 
-@mcp.tool(name="get_physical_constants", description="获取大气和海洋的物理常量")
-def get_physical_constants() -> dict:
+# ==================== Column Model Tools ====================
+
+@mcp.tool(name="run_grey_radiation_model", description="Run a grey-gas radiative column model")
+def run_grey_radiation_model_tool(num_lev: int = 30, num_years: float = 2.0,
+                                  albedo_sfc: float = 0.299, Q: float = 341.3) -> dict:
     """
-    获取大气和海洋的物理常量。
+    Run a grey-gas radiative column model to equilibrium.
 
-    Returns:
-    - dict: 包含常量的字典。
+    :param num_lev: Number of vertical levels (default: 30).
+    :param num_years: Number of years to integrate (default: 2.0).
+    :param albedo_sfc: Surface albedo (default: 0.299).
+    :param Q: Insolation in W/m² (default: 341.3).
+    :return: A dictionary containing the equilibrium temperature profile.
     """
     try:
-        from climlab.utils.constants import a, Lhvap, Lhsub, Lhfus, cp, Rd, kappa, Rv, cpv, eps, Omega, g, kBoltzmann, c_light, hPlanck, sigma, S0, ps, rho_w, cw
-        constants = {
-            "Earth_radius": a,
-            "Latent_heat_vaporization": Lhvap,
-            "Latent_heat_sublimation": Lhsub,
-            "Latent_heat_fusion": Lhfus,
-            "Specific_heat_dry_air": cp,
-            "Gas_constant_dry_air": Rd,
-            "Kappa": kappa,
-            "Gas_constant_water_vapor": Rv,
-            "Specific_heat_water_vapor": cpv,
-            "Epsilon": eps,
-            "Earth_rotation_rate": Omega,
-            "Gravitational_acceleration": g,
-            "Boltzmann_constant": kBoltzmann,
-            "Speed_of_light": c_light,
-            "Planck_constant": hPlanck,
-            "Stefan_Boltzmann_constant": sigma,
-            "Solar_constant": S0,
-            "Surface_pressure": ps,
-            "Density_of_water": rho_w,
-            "Specific_heat_liquid_water": cw
+        model = GreyRadiationModel(num_lev=num_lev, albedo_sfc=albedo_sfc, Q=Q)
+        model.integrate_years(num_years)
+        
+        Ts = float(np.squeeze(model.Ts))
+        Tatm = np.squeeze(model.Tatm).tolist()
+        lev = model.lev.tolist()
+        
+        result = {
+            "surface_temperature_K": Ts,
+            "surface_temperature_C": Ts - 273.15,
+            "atmospheric_temperature_K": Tatm,
+            "pressure_levels_hPa": lev,
+            "OLR_W_m2": float(np.squeeze(model.OLR)),
+            "ASR_W_m2": float(np.squeeze(model.ASR)),
+            "num_levels": num_lev,
+            "integrated_years": num_years,
         }
-        return {"success": True, "constants": constants}
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
 
-@mcp.tool(name="calculate_potential_temperature", description="计算大气的位温")
-def calculate_potential_temperature(T: float, p: float) -> dict:
-    """
-    计算大气的位温。
 
-    Parameters:
-    - T: 温度 (单位: K)。
-    - p: 压力 (单位: hPa)。
+@mcp.tool(name="run_radiative_convective_model", description="Run a radiative-convective column model with convective adjustment")
+def run_radiative_convective_model_tool(num_lev: int = 30, num_years: float = 5.0,
+                                        adj_lapse_rate: float = 6.5,
+                                        albedo_sfc: float = 0.299, Q: float = 341.3) -> dict:
+    """
+    Run a radiative-convective equilibrium model.
 
-    Returns:
-    - dict: 包含位温的字典。
+    :param num_lev: Number of vertical levels (default: 30).
+    :param num_years: Number of years to integrate (default: 5.0).
+    :param adj_lapse_rate: Convective adjustment lapse rate in K/km (default: 6.5).
+    :param albedo_sfc: Surface albedo (default: 0.299).
+    :param Q: Insolation in W/m² (default: 341.3).
+    :return: A dictionary containing the RCE temperature profile.
     """
     try:
-        from climlab.utils.thermo import potential_temperature
-        theta = potential_temperature(T, p)
-        return {"success": True, "potential_temperature": theta}
+        model = RadiativeConvectiveModel(num_lev=num_lev, adj_lapse_rate=adj_lapse_rate,
+                                         albedo_sfc=albedo_sfc, Q=Q)
+        model.integrate_years(num_years)
+        
+        Ts = float(np.squeeze(model.Ts))
+        Tatm = np.squeeze(model.Tatm).tolist()
+        lev = model.lev.tolist()
+        
+        result = {
+            "surface_temperature_K": Ts,
+            "surface_temperature_C": Ts - 273.15,
+            "atmospheric_temperature_K": Tatm,
+            "pressure_levels_hPa": lev,
+            "OLR_W_m2": float(np.squeeze(model.OLR)),
+            "ASR_W_m2": float(np.squeeze(model.ASR)),
+            "adj_lapse_rate_K_km": adj_lapse_rate,
+            "num_levels": num_lev,
+            "integrated_years": num_years,
+        }
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
+
+
+# ==================== Climate Sensitivity Tools ====================
 
-@mcp.tool(name="calculate_clausius_clapeyron", description="计算克劳修斯-克拉珀龙方程")
-def calculate_clausius_clapeyron(T: float) -> dict:
+@mcp.tool(name="calculate_climate_sensitivity", description="Calculate equilibrium climate sensitivity from EBM parameters")
+def calculate_climate_sensitivity_tool(B: float = 2.0, f: float = 0.0) -> dict:
     """
-    计算克劳修斯-克拉珀龙方程。
+    Calculate equilibrium climate sensitivity (ECS) from feedback parameters.
+    ECS = -ΔF / λ, where λ = B - f is the net feedback parameter.
 
-    Parameters:
-    - T: 温度 (单位: K)。
+    :param B: Planck feedback parameter in W/m²/°C (default: 2.0).
+    :param f: Sum of other feedback parameters in W/m²/°C (default: 0.0).
+    :return: A dictionary containing the climate sensitivity.
+    """
+    try:
+        # Net feedback parameter
+        lambda_net = B - f
+        
+        # Standard CO2 doubling forcing
+        delta_F_2xCO2 = 3.7  # W/m²
+        
+        # Equilibrium climate sensitivity
+        if lambda_net > 0:
+            ECS = delta_F_2xCO2 / lambda_net
+        else:
+            ECS = float('inf')  # Runaway climate
+        
+        result = {
+            "Planck_feedback_B": B,
+            "other_feedbacks_f": f,
+            "net_feedback_lambda": lambda_net,
+            "CO2_doubling_forcing_W_m2": delta_F_2xCO2,
+            "equilibrium_climate_sensitivity_C": float(ECS),
+            "description": f"For λ = {lambda_net:.2f} W/m²/°C, ECS = {ECS:.2f} °C per CO2 doubling",
+        }
+        return {"success": True, "result": result, "error": None}
+    except Exception as e:
+        return {"success": False, "result": None, "error": str(e)}
+
+
+@mcp.tool(name="run_ebm_co2_doubling", description="Run EBM to simulate global warming from CO2 doubling")
+def run_ebm_co2_doubling_tool(num_lat: int = 90, num_years: float = 50.0,
+                              A_control: float = 210.0, forcing: float = 3.7) -> dict:
+    """
+    Run EBM experiment comparing control run with CO2 doubling scenario.
 
-    Returns:
-    - dict: 包含计算结果的字典。
+    :param num_lat: Number of latitude points (default: 90).
+    :param num_years: Number of years to integrate (default: 50.0).
+    :param A_control: Control OLR parameter A in W/m² (default: 210.0).
+    :param forcing: Radiative forcing from CO2 doubling in W/m² (default: 3.7).
+    :return: A dictionary containing control and perturbed climate states.
     """
     try:
-        from climlab.utils.thermo import clausius_clapeyron
-        result = clausius_clapeyron(T)
-        return {"success": True, "result": result}
+        # Control run
+        model_ctrl = EBM(num_lat=num_lat, A=A_control)
+        model_ctrl.integrate_years(num_years)
+        T_ctrl = float(global_mean(model_ctrl.Ts))
+        
+        # Perturbed run (CO2 doubling reduces OLR, equivalent to reducing A)
+        model_2xCO2 = EBM(num_lat=num_lat, A=A_control - forcing)
+        model_2xCO2.integrate_years(num_years)
+        T_2xCO2 = float(global_mean(model_2xCO2.Ts))
+        
+        warming = T_2xCO2 - T_ctrl
+        
+        result = {
+            "control_global_mean_T_C": T_ctrl,
+            "2xCO2_global_mean_T_C": T_2xCO2,
+            "equilibrium_warming_C": warming,
+            "forcing_W_m2": forcing,
+            "effective_climate_sensitivity_C": warming,
+            "control_temperature_profile": np.squeeze(model_ctrl.Ts).tolist(),
+            "2xCO2_temperature_profile": np.squeeze(model_2xCO2.Ts).tolist(),
+            "latitudes": model_ctrl.lat.tolist(),
+        }
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
+
 
-@mcp.tool(name="test_ebm_model", description="测试 EBM 模型的创建和运行")
-def test_ebm_model() -> dict:
+# ==================== Physical Constants Tools ====================
+
+@mcp.tool(name="get_climate_constants", description="Get commonly used climate physical constants")
+def get_climate_constants_tool() -> dict:
     """
-    测试 EBM 模型的创建和运行。
+    Get commonly used climate physical constants from climlab.
 
-    Returns:
-    - dict: 测试结果。
+    :return: A dictionary containing physical constants.
     """
     try:
-        import numpy as np
-        import climlab
-        old_Jan1 = np.datetime64('2025-03-20T09:01') - np.timedelta64(80, 'D')
-        model = climlab.EBM_seasonal(initial_time=old_Jan1, water_depth=10.)
-        model.integrate_years(1)
-        return {
-            "success": True,
-            "state": model.state,
-            "diagnostics": model.diagnostics
+        result = {
+            "earth_radius_m": const.a,
+            "gravitational_acceleration_m_s2": const.g,
+            "solar_constant_W_m2": const.S0,
+            "stefan_boltzmann_constant_W_m2_K4": const.sigma,
+            "specific_heat_dry_air_J_kg_K": const.cp,
+            "gas_constant_dry_air_J_kg_K": const.Rd,
+            "latent_heat_vaporization_J_kg": const.Lhvap,
+            "latent_heat_fusion_J_kg": const.Lhfus,
+            "water_density_kg_m3": const.rho_w,
+            "specific_heat_water_J_kg_K": const.cw,
+            "seconds_per_day": const.seconds_per_day,
+            "days_per_year": const.days_per_year,
+            "earth_surface_area_m2": const.area_earth,
+            "present_day_orbital_parameters": const.orb_present,
         }
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
+
+
+# ==================== Utility Tools ====================
 
-@mcp.tool(name="test_rcm_model", description="测试 RCM 模型的创建和运行")
-def test_rcm_model() -> dict:
+@mcp.tool(name="calculate_global_mean", description="Calculate area-weighted global mean of a latitude-dependent field")
+def calculate_global_mean_tool(values: list, latitudes: list = None) -> dict:
     """
-    测试 RCM 模型的创建和运行。
+    Calculate the global mean of a field with proper area weighting.
 
-    Returns:
-    - dict: 测试结果。
+    :param values: List of field values at each latitude.
+    :param latitudes: List of latitudes in degrees (optional, defaults to even spacing).
+    :return: A dictionary containing the global mean.
     """
     try:
-        import numpy as np
-        import climlab
-        state = climlab.column_state(num_lev=40, num_lat=1, water_depth=5.)
-        h2o = climlab.radiation.ManabeWaterVapor(state=state, name='H2O')
-        convadj = climlab.convection.ConvectiveAdjustment(state=state, name='ConvectiveAdjustment', adj_lapse_rate=6.5)
-        rad = climlab.radiation.RRTMG(state=state, albedo=0.2, specific_humidity=h2o.q, name='Radiation')
-        rcm = climlab.couple([h2o, convadj, rad], name='RCM')
-        rcm.step_forward()
-        return {
-            "success": True,
-            "state": rcm.state,
-            "diagnostics": rcm.diagnostics
+        values_array = np.array(values)
+        if latitudes is None:
+            latitudes = np.linspace(-90, 90, len(values_array))
+        lat_array = np.array(latitudes)
+        
+        # Area weighting using cosine of latitude
+        weights = np.cos(np.deg2rad(lat_array))
+        global_mean_value = float(np.average(values_array, weights=weights))
+        
+        result = {
+            "global_mean": global_mean_value,
+            "simple_mean": float(np.mean(values_array)),
+            "max_value": float(np.max(values_array)),
+            "min_value": float(np.min(values_array)),
+            "num_points": len(values_array),
         }
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
 
-@mcp.tool(name="test_grey_radiation_model", description="测试灰气辐射模型的创建和运行")
-def test_grey_radiation_model() -> dict:
+
+@mcp.tool(name="calculate_meridional_heat_transport", description="Calculate implied meridional heat transport from energy imbalance")
+def calculate_meridional_heat_transport_tool(ASR: list, OLR: list, latitudes: list = None) -> dict:
     """
-    测试灰气辐射模型的创建和运行。
+    Calculate meridional heat transport from absorbed solar radiation and OLR.
 
-    Returns:
-    - dict: 测试结果。
+    :param ASR: List of absorbed shortwave radiation values (W/m²).
+    :param OLR: List of outgoing longwave radiation values (W/m²).
+    :param latitudes: List of latitudes in degrees (optional).
+    :return: A dictionary containing the heat transport profile.
     """
     try:
-        import numpy as np
-        import climlab
-        old_Jan1 = np.datetime64('2025-03-20T09:01') - np.timedelta64(80, 'D')
-        model = climlab.GreyRadiationModel(num_lev=30, num_lat=90)
-        model.integrate_years(1)
-        return {
-            "success": True,
-            "state": model.state,
-            "diagnostics": model.diagnostics
+        ASR_arr = np.array(ASR)
+        OLR_arr = np.array(OLR)
+        
+        if latitudes is None:
+            latitudes = np.linspace(-90, 90, len(ASR_arr))
+        lat = np.array(latitudes)
+        
+        # Net radiation at TOA
+        net_rad = ASR_arr - OLR_arr
+        
+        # Calculate heat transport by integrating from South Pole
+        phi = np.deg2rad(lat)
+        dlat = np.abs(lat[1] - lat[0]) if len(lat) > 1 else 1.0
+        dphi = np.deg2rad(dlat)
+        
+        # Integrate: H(φ) = 2πa² ∫ R(φ') cos(φ') dφ'
+        integrand = net_rad * np.cos(phi)
+        heat_transport = 2 * np.pi * const.a**2 * np.cumsum(integrand) * dphi
+        heat_transport_PW = heat_transport * 1e-15  # Convert to PW
+        
+        result = {
+            "latitudes": lat.tolist(),
+            "net_radiation_W_m2": net_rad.tolist(),
+            "heat_transport_PW": heat_transport_PW.tolist(),
+            "max_poleward_transport_PW": float(np.max(np.abs(heat_transport_PW))),
+            "latitude_of_max_transport": float(lat[np.argmax(np.abs(heat_transport_PW))]),
         }
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
+
 
-# 创建 FastMCP 应用实例
 def create_app() -> FastMCP:
     """
-    创建并返回 FastMCP 应用实例。
+    Create and return the FastMCP application instance.
 
-    Returns:
-    - FastMCP: FastMCP 应用实例。
+    :return: The FastMCP application instance.
     """
     return mcp