Paper deep dive
Systematic LLM Translation of Legacy Scientific Code to Differentiable Frameworks: Application to a Land Surface Model
Aya Lahlou, Linnia Hawkins, Pierre Gentine
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 94%
Last extracted: 7/9/2026, 2:23:10 AM
Summary
The paper introduces a five-phase LLM-based agentic pipeline that systematically translates legacy Fortran scientific code into the JAX differentiable programming framework. Applied to the CLM-ml-v2 land surface model, the pipeline leverages static dependency analysis, state documentation, iterative compile-repair loops (Ralph loops), and a Fortran reference oracle to ensure numerical parity. The resulting differentiable model, clm-ml-jax, computes full Jacobians in a single backward pass, achieves a 24x wall-clock speedup over sequential Fortran, and enables efficient gradient-based parameter estimation and sensitivity analysis.
Entities (12)
Relation Signals (11)
LLM → translates → Fortran
confidence 97% · We present a five-phase LLM-based agentic pipeline that translates legacy Fortran into JAX
LLM → translatesto → JAX
confidence 97% · translates legacy Fortran into JAX
Pipeline → appliedto → CLM-ml-v2
confidence 96% · We instantiate and evaluate the pipeline on CLM-ml-v2, a 19,000-line Fortran land surface model
clm-ml-jax → achieves → 24x speedup
confidence 95% · achieves a 24 times wall-clock speedup over sequential Fortran at ensemble size N=2,048
Pipeline → produces → clm-ml-jax
confidence 95% · producing clm-ml-jax: a fully differentiable and GPU enabled validated re-implementation
Pipeline → uses → Ralph Loop
confidence 94% · Both cycles are implemented as Ralph loops, an orchestration pattern that intercepts the agent’s exit signal
Pipeline → uses → plan.md
confidence 93% · plan.md tracks per-module, per-phase completion flags
Pipeline → uses → CLAUDE.md
confidence 93% · CLAUDE.md encodes project conventions, requirements, and Fortran-to-JAX mapping rules
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Differentiable programming offers transformative capabilities for scientific modeling, enabling gradient-based parameter estimation, sensitivity analysis, and data assimilation. Yet, migrating legacy codebases into differentiable frameworks remains a challenge. We present a five-phase LLM-based agentic pipeline that translates legacy Fortran into JAX: static dependency analysis determines module translation order from the full call graph; iterative compile-repair loops correct errors autonomously; and a Fortran reference oracle enforces numerical parity at the module level before integration and gradient verification. We instantiate and evaluate the pipeline on CLM-ml-v2, a 19,000-line Fortran land surface model, and analyze agent behavior across 73 module translation tasks. The resulting differentiable model computes the complete Jacobian in a single backward pass, recovers physical parameters in eight times fewer steps than gradient-free optimization, and achieves a 24 times wall-clock speedup over sequential Fortran at ensemble size N=2,048. Both the translated model and pipeline infrastructure are released as a reusable framework for differentiating other Earth system model components.
Tags
Links
- Source: https://arxiv.org/abs/2606.07681v1
- Canonical: https://arxiv.org/abs/2606.07681v1
Trouble viewing inline? Open PDF directly →
Full Text
59,693 characters extracted from source content.
Expand or collapse full text
Systematic LLM Translation of Legacy Scientific Code to Differentiable Frameworks: Application to a Land Surface Model Aya Lahlou 1 Linnia Hawkins 1 Pierre Gentine 1 Abstract Differentiable programming offers transforma- tive capabilities for scientific modeling, enabling gradient-based parameter estimation, sensitivity analysis, and data assimilation. Yet, migrating legacy codebases into differentiable frameworks remains a challenge. We present a five-phase LLM-based agentic pipeline that translates legacy Fortran into JAX: static dependency analysis de- termines module translation order from the full call graph; iterative compile-repair loops correct errors autonomously; and a Fortran reference ora- cle enforces numerical parity at the module level before integration and gradient verification. We instantiate and evaluate the pipeline on CLM-ml- v2, a 19,000-line Fortran land surface model, and analyze agent behavior across 73 module transla- tion tasks. The resulting differentiable model com- putes the complete Jacobian in a single backward pass, recovers physical parameters in eight times fewer steps than gradient-free optimization, and achieves a 24 times wall-clock speedup over se- quential Fortran at ensemble size N=2,048. Both the translated model and pipeline infrastructure are released as a reusable framework for differen- tiating other Earth system model components. 1. Introduction Legacy Scientific Models, such as Earth system models (ESMs), are among the most consequential yet most out- dated large software systems in active use (Neumann et al., 2019). Popular ESMs comprise 500,000 to 1.3 million lines of Fortran code (M ́ endez et al., 2014; M ́ endez & Tinetti, 2016), which support international climate policy and fore- cast weather for billions of people. Yet, surveys show per- vasive deprecated components, poor modularity, missing test coverage (M ́ endez et al., 2014), and a steadily shrinking 1 Department of Earth and Environmental Engineering, Columbia University, New York, United States. Correspondence to: Aya Lahlou <al4385@columbia.edu>. Preprint. Submitted to ICML 2026 DL4C workshop. May, 2026. developer workforce (ECMWF, 2024). Despite active community efforts to modernize these code- bases, no comprehensive, fully differentiable Earth system model exists to date (Gelbrecht et al., 2023). End-to-end differentiable models enable gradient-based joint calibra- tion of physical parameters (Raoult et al., 2025; Qu et al., 2024), automatic adjoint generation for variational data as- similation, replacing decade-long hand-coded adjoint main- tenance projects (Le Dimet & Talagrand, 1986; Talagrand & Courtier, 1987; Giering & Kaminski, 1998), and online train- ing of neural network closures inside the dynamical core, the only configuration shown to yield stable hybrid climate integrations (Kochkov et al., 2024; Rasp et al., 2018). Re- verse automatic differentiation computes the gradient of any scalar loss with respect to all p parameters in a single back- ward pass (Baydin et al., 2018), collapsing theO(10 4 –10 5 ) model evaluations per loss required by ensemble Kalman filtering (Evensen, 2003) or MCMC (Vrugt et al., 2009) to a single forward-backward pair. LLM-based code translation is promising but brittle. Pan et al. (2024) showed that single-pass LLMs achieve low functional correctness even across high-resource language pairs, and Fortran’s under-representation in pre-training cor- pora compounds the difficulty further (Jin et al., 2024; Ka- dosh et al., 2025). Iterative compile-execute-repair cycles are critical for functionally correct output (Nichols et al., 2025; Jin et al., 2024; IBM Research, 2023). Thus, the open problem is not whether LLMs can translate Fortran, but how to produce a numerically equivalent, fully differentiable translation of a scientific codebase. We present a five-phase agentic pipeline that addresses this problem through three design principles missing from prior work: (i) static dependency analysis ensuring topologicaly ordered context-aware translation; (i) modular numerical parity testing using Fortran reference inputs/outputs; and (i) comprehensive gradient verification. We apply the pipeline to CLM-ml-v2 (Bonan et al., 2021), a 19,000-line, 102-module multilayer canopy land surface model coupling leaf photosynthesis (FvCB (Farquhar et al., 1980)), stomatal conductance (Medlyn (Medlyn et al., 2011)), and Harman- Finnigan roughness sublayer turbulence (Harman & Finni- gan, 2007; 2008) with a Runge-Kutta solver (Butcher, 1964), 1 arXiv:2606.07681v1 [cs.SE] 4 Jun 2026 Systematic LLM Translation of Legacy Scientific Code to Differentiable Frameworks producingclm-ml-jax: a fully differentiable and GPU enabled validated re-implementation. 2. Related Work 2.1. Differentiable Programming for Earth System Models Code modernization in Earth system modeling has taken two distinct paths. The first is clean rewriting in modern frameworks, most of which aim for end-to-end differentia- bility (Wang & Frankenberg, 2023; H ̈ afner et al., 2018), but few models to date have achieved it (Kochkov et al., 2024). (Kochkov et al., 2024; Meunier et al., 2025) rebuilds atmospheric dynamics in JAX and is fully differentiable end- to-end. The second path is faithful translation: preserving the physics and parameterizations of an existing operational model. NoahPy (Tian et al., 2026) translates the Noah LSM to PyTorch, (Davenport et al., 2026) implements a differen- tiable version of SPEEDY-based physics (Kucharski et al., 2013) parametrization in JAX. JAX is a predominant framework for differentiable physics- based simulations (Baydin et al., 2018; Gelbrecht et al., 2023; Shen et al., 2023; Campagne et al., 2023; Bezgin et al., 2023; Jiang et al., 2025) enabling composable autodiff, vectorization, and GPU/TPU simulations without manual porting. Gelbrecht et al. (2023) identifies automated code translation as a promising path toward making operational ESMs differentiable; this work investigates that direction using a legacy land surface model. 2.2. Parameter Estimation in Land Surface Models Parameter estimation for Earth System Models, and LSMs in particular, has traditionally relied on sampling methods such as Monte Carlo Markov Chain (MCMC) or ensemble Kalman filtering, requiringO(10 3 –10 5 )model evaluations (Vrugt et al., 2009; Evensen, 2003). Gradient-based alternatives, identified as a key emerging avenue by Raoult et al. (2025), reduce per-step cost to a sin- gle forward-backward pass, independent of parameter count. AdJULES (Raoult et al., 2016) demonstrated gradient-based LSM calibration via a proprietary adjoint compiler (TAF (Giering & Kaminski, 1998)). Still, TAF carries a com- mercial license (Fastopt GmbH), operates via a Fortran- only source-to-source transformation with no GPU support, and cannot differentiate using iterative solvers (e.g., bisec- tion root-finding, Monin–Obukhov stability iteration) with- out manually reformulating each solver.clm-ml-jax removes all three constraints: the translation pipeline is open, and LLM-driven, JAX/XLA compiles gradients na- tively to GPU, and iterative solvers are differentiated via implicit function theorem adjoints (Blondel et al., 2022b), a capability absent from prior adjoint-based LSM work. 2.3. LLM-Assisted Scientific Code Translation LLM-assisted translation of legacy Fortran has been demon- strated for particle-physics codes (CodeScribe (Dhruv & Dubey, 2024)), atmospheric HPC kernels (Gupta et al., 2025), and general Fortran-to-C++ tasks (Ranasinghe et al., 2025; Li et al., 2024).Translating entire scientific codebases remains difficult to scale; performance degrades when mov- ing beyond isolated functions, largely due to cross-file de- pendencies and build-system complexity. A further chal- lenge is numerical fidelity. Beyond syntax, translated code must reproduce the original physics to high precision, re- quiring domain-specific testing infrastructure (Pietrini et al., 2024; Ranasinghe et al., 2025). Recent benchmarks reveal a significant deficit in evaluating LLMs for scientific code generation: no standard benchmark exists that measures the scientific validity of LLM-generated code. For Earth system models in particular, Zhou et al. (2024) shows that ChatGPT can translate isolated land-surface routines into JAX. Still, their method operates on individual functions without de- pendency analysis, provides no systematic numerical parity validation against the Fortran reference, and does not verify the correctness of gradients required for scientific deploy- ment. Unlike general software correctness benchmarks such as SWE-bench (Jimenez et al., 2024), which define correct- ness via test-suite execution, scientific translation demands a domain-specific oracle that enforces numerical parity with the original physical model. 3. Methodology We describe our five-phase workflow for producing clm-ml-jax. The workflow is instantiated here for CLM- ml-v2 but is designed to be transferred to other ESM parts in the future as part of a large refactoring project. A model is suitable for this workflow if it satisfies three conditions: (a) deterministic single-column execution under fixed inputs and physical constants, (b) accessible Fortran source with a standard batch build system, and (c) external linkage against production object files. 3.1. Phase 1: Scoping and Dependency Analysis Naive file-by-file agent prompting fails on coupled code- bases because LLMs hallucinate data types for arguments whose types are implicit inUSE-statement chains spanning multiple files (Pietrini et al., 2024; Li et al., 2024). We therefore apply static analysis to map the full codebase de- pendency structure before any Translation begins. We subject the CLM-ml-v2 module dependency graph to a topological sort (Kahn, 1962) to yield a dependency- respecting translation order: When moduleAis translated, all modules it imports have already been constructed, elimi- nating the hallucinated-type failures that arise in single-file 2 Systematic LLM Translation of Legacy Scientific Code to Differentiable Frameworks prompting. The translation order and other static analysis findings are saved in state documents (Section 3.2), serving as the agent’s initialization context and task list. This structured decompo- sition reduces the translation of a 19,000-line, 102-module codebase to a sequence of bounded, context-complete tasks, making autonomous agentic translation tractable at this scale. Static Analyzer ESM Architecture & Data Flow Topologically sorted Module List Plan.md Claude.md Human Developer CHANGELOG.md Fortran Codebase Figure 1. Flowchart of Static Dependency Analysis and State Documentation Setup. 3.2. Phase 2: State Documentation Setup We establish three persistent state documents checked into the repository (Figure 1):CLAUDE.mdencodes project conventions, requirements, and Fortran-to-JAX mapping rules;plan.mdtracks per-module, per-phase completion flags; andCHANGELOG.mdwhere agents record current session status and failed attempts. Logging failed attempts is critical for preventing successive sessions from revisiting the same dead-ends. Together, these documents function as external agent memory, preserving project state across context-window boundaries and enabling coherent session resumption without reliance on prior conversation history. 3.3. Phase 3: Fortran Oracle Construction and Functional Testing Testing each translated module against its Fortran refer- ence ensures that deviations from scientific ground truth are caught early and attributed to a specific component rather than obscured by the full model’s complexity. We task a team of agents to develop a functional test suite for CLM-ml-v2 following the workflow described in fig- ure 2, classifying subroutines into three tiers by testability. Tier 1 subroutines have pure scalar interfaces and are di- rectly callable. Tier 2 subroutines contain isolable physics inside derived-type loops, requiring minor non-breaking visibility changes or scalar helper extraction before test- ing. Tier 3 subroutines are fully coupled to model state and cannot be tested in isolation; they are instead validated indirectly through full-column parity. This phase produced 26 executables covering tests for 32 subroutines across 10 modules. Plan.md Claude.md CHANGELOG.md Functional Testing Session ESM scanner Agent Test Generator Agent Fortran .exe i Generate Reference I/O module i reference I/O json Fortran Codebase Fortran module i + dependencies Pytest script i RL exit If complete RL start Figure 2. Fortran Oracle Construction and Functional Testing Flowchart. RL refers to Ralph Loop. Dotted arrows refer to updates to the state documents. The Fortran tests are executed once under fixed inputs, sav- ing the reference input-output pairs at each checkpoint for validation during the translation phase. 3.4. Iterative Translation, Testing, and Repair Following the bottom-up dependency order from Phase 1, each module is processed by two nested autonomous loops (Figure 3). A translator agent, guided byCLAUDE.md, produces a JAX equivalent of the current module. A test- ing agent then constructs a numerical parity test, feeding recorded Fortran inputs from the Phase 3 golden I/O files into the JAX module and asserting that all outputs agree with the Fortran reference to within a relative tolerance of 10 −4 . A passing test advances the module’splan.mdflag and the outer loop moves to the next module. A failing test engages the inner repair loop: a repair agent reads the test output, diagnoses the discrepancy, edits the module, and re-runs the parity test, iterating until the test passes or the module is flagged inCHANGELOG.mdfor human inspection. Both cycles are implemented as Ralph loops (Huntley, 2026), an orchestration pattern that intercepts the agent’s exit signal and re-injects the driving prompt iteratively un- til a verifiable goal condition is met, enabling unattended multi-hour sessions: while true; do cat PROMPT_translate.md | claude --dangerously-skip-permissions done This design addresses two failure modes of current LLMs in long running tasks. First, agentic laziness(Kwa et al., 3 Systematic LLM Translation of Legacy Scientific Code to Differentiable Frameworks 2025): agents tend to signal completion prematurely; the loop forces re-evaluation of allplan.mdtasks before exit. Second, statelessness:plan.mdandCHANGELOG.md provide the external memory required for coherent resump- tion across sessions (Park et al., 2023). Human over- sight is maintained asynchronously viaCHANGELOG.md; the--dangerously-skip-permissionsflag sup- presses interactive prompts for headless operation and should be used only within an appropriate sandbox. A module passes when its output relative RMSE against the Fortran oracle is below 1% across held-out inputs. Claude.md CHANGELOG.md Translation & Validation Session Translator Agent JAX module i Tester Agent module i tests FailPass Repair Agent fixed JAX module i Plan.md Fortran module i T. RL exit If complete T. RL start R. RL start R. RL exit Figure 3. Autonomous Agentic Translation, Testing, and Repair Flowchart. 3.5. Phase 5: Integration and Differentiability Validation Following module-level translation, an Integration Agent assembles the full column pipeline and validates end-to-end behavior against the Full column Fortran oracle as seen in figure 4. This phase revealed errors not exposed by module- level parity tests, such as mismatched array shapes and indexing conventions across module boundaries. These errors were resolved using the Full Column Repair Agent with shape annotations and systematic indexing checks. The integrated model then achieved full-column parity with the Fortran reference. Then, a Differentiability Agent is deployed to iteratively identify and resolve JAX-specific issues required for correct and efficient differentiation. The differentiability validation runs iteratively until all active (output, parameter) gradient pairs fall within 1% of the finite-difference estimates. Teams targeting other frameworks may encounter analogous issues (Table 1 lists the general method and JAX-specific fixes). Lessons Learned from validating the differentiability of clm-ml-jax are as follow: Static loop unrolling The target framework must trace through the loop body, not execute it iteratively at the Python Claude.md CHANGELOG.md Integration & Differentiability Session Integration Diags. Agent Full column Parity Tests FailPass Full Column Repair Agent fixed JAX codebase Plan.md I. RL start I RL exit Differentiability Diags. Agent Diff. Repair Agent Diff. Validation Tests FailPass fixed JAX codebase D RL exit D. RL start Figure 4. Integration and Differentiability Validation Flowchart. level. FortranDOloops over Runge-Kutta sub-steps are re- placed withjax.lax.scan, tracing the loop body once at compile time, resulting in a single XLA kernel spanning all sub-steps. This transformation is necessary for efficient gradient computation — JAX can differentiate Pythonfor loops by unrolling them into the computation graph, but do- ing so forO(10 2 )Runge-Kutta sub-steps creates an unman- ageably large computation graph and incurs per-iteration Python–XLA dispatch overhead — and yields a∼200× reduction in per-step wall-clock time. Static vs. Dynamic Typing: Any value that can change between calls must not enter the compile-time cache key. Closures that capture Python scalars extracted from JAX ar- rays introduce a new value into the XLA cache key on every call, triggering full recompilation (∼290 s per run). We elim- inate this by converting all such constants to JAX scalars before closing over them injit-compiled functions and caching kernel factories withfunctools.lrucache. After this fix, subsequent runs reuse the compiled kernel (∼0.3 s cached lookup). Branch-safe Guard JAX evaluates both branches of a jnp.whereexpression during the backward pass. Branch- local operations of the formx ** nor1/xwithx = 0pro- duceinfgradients; the product0×∞ = NaNpropagates through the computation graph. We eliminate this by ap- plying safe lower bounds (jnp.maximum(x, 1e-30)) before fractional powers and divisions in five identified mod- ules. Implicit function theorem (IFT) adjoints for iterative solvers.Two physics routines implement iterative root- finding viajax.lax.foriloop: the WUE Cowan– Farquhar stomatal conductance bisection (Cowan & Far- quhar, 1977) and the Monin–Obukhov length solver (Monin & Obukhov, 1954).Differentiating through the loop directly yields incorrect gradients becausejnp.where propagates gradients through both branches simultane- 4 Systematic LLM Translation of Legacy Scientific Code to Differentiable Frameworks ously at every iteration. We resolve this using the well- established IFT adjoint (Bai et al., 2019; Blondel et al., 2022a): given the converged rootx ∗ , we formx IFT = x ∗ − F(x ∗ ;θ) / stopgradient(∂F/∂x), so that in the forward passx IFT ≈ x ∗ , and in the backward pass ∂x IFT /∂θ = (∂F/∂θ)/(∂F/∂x). 3.6. Agentic Pipeline Characterisation To characterize the process for teams applying this method- ology to other ESMs, we tracked all debugging sessions, bug instances, and resolution attempts of Phase 5’s differentia- bility validation via structured logging inCHANGELOG.md. Achieving correct end-to-end gradients required fixing 48 documented bugs across 9 mechanistic types over an esti- mated 46 agentic sessions (Appendix C, Table 2). The domi- nant class, T1 (NaN gradients fromjnp.where), accounts for 35.4% of bugs (17/48) and follows one mechanical fix pattern:jnp.maximum(x,1e-30)applied at 35+ sites across 9 modules; all 17 were resolved on the agent’s first attempt. The hardest class, T5 (gradient explosion from iterative solvers), required an average of 2.0 attempts per bug and the IFT insight above; without IFT, the Obukhov secant-solver gradient reached9.95×10 144 . Overall, 81.2% of bugs were resolved autonomously on the first attempt; the remaining 18.8% required human-directed re-specification (principally T5 and T2 parameter-injection failures, where the failure mode requires detailed knowledge of JAX’s trac- ing model). Estimated human oversight:∼10–15 hours across 46 agentic sessions. Full taxonomy, timeline, module density map, and failure analysis are in Appendix C. 4. Experiments We evaluate clm-ml-jax along two axes: (i) numerical cor- rectness versus the Fortran oracle, and (i) gradient correct- ness using central finite differences. We also run proof-of- concept experiments demonstrating the merits of differen- tiability for parameter tuning in the multilayer canopy land model. All experiments use the CHATS7 walnut orchard AmeriFlux site (Patton et al., 2011) with May 2007 mete- orological forcing data (30-min resolution). Hardware and software specifications are included in the appendix. 4.1. Validation of Numerical Equivalence We compareclm-ml-jaxoutputs against the Fortran reference across a 31-day simulation in May 2007 at the CHATS7 site. Figure 5 shows time series and scatter com- parisons across the full month; canopy profile comparisons are shown in Figure 10. 0200400600 0 200 400 600 JAX CLM-ml Bias : -0.004 Corr : 1.000 NSE : 1.000 (a) Rn (W m −2 ) −100−50050100150 −100 −50 0 50 100 150 Bias : 0.000 Corr : 1.000 NSE : 1.000 (b) H (W m −2 ) −200−1000100200 −200 −100 0 100 200 Bias : -0.028 Corr : 1.000 NSE : 1.000 (c) LE (W m −2 ) 0100200300400500 0 100 200 300 400 500 Bias : 0.025 Corr : 1.000 NSE : 1.000 (d) Rabs (W m −2 ) 051015202530 0 5 10 15 20 25 30 JAX CLM-ml Bias : 0.001 Corr : 1.000 NSE : 1.000 (e) LWup (W m −2 ) 0.000.250.500.751.001.251.501.75 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 Bias : 0.000 Corr : 1.000 NSE : 1.000 (f) ET (m s −1 ) 0255075100125150175 0 25 50 75 100 125 150 175 Bias : 0.000 Corr : 1.000 NSE : 1.000 (g) GPP (μmol m −2 s −1 ) 340360380400420440460480500 340 360 380 400 420 440 460 480 500 Bias : 0.004 Corr : 1.000 NSE : 1.000 (h) SWdn (W m −2 ) 280285290295300305 280 285 290 295 300 305 JAX CLM-ml Bias : 0.000 Corr : 1.000 NSE : 1.000 (i) Tair (K) −75−50−250255075100125 −75 −50 −25 0 25 50 75 100 125 Bias : -0.000 Corr : 1.000 NSE : 1.000 (j) Tveg (K) −250255075100125150 −25 0 25 50 75 100 125 150 Bias : 0.002 Corr : 1.000 NSE : 1.000 (k) Tsoil (K) −5051015202530 −5 0 5 10 15 20 25 30 Bias : 0.001 Corr : 1.000 NSE : 1.000 (l) u* (m s −1 ) 05101520253035 Fortran CLM-ml 0 5 10 15 20 25 30 35 JAX CLM-ml Bias : 0.001 Corr : 1.000 NSE : 1.000 (m) CO 2 (ppm) 0100200300400500 Fortran CLM-ml 0 100 200 300 400 500 Bias : 0.024 Corr : 1.000 NSE : 1.000 (n) Rnet_veg (W m −2 ) 010203040506070 Fortran CLM-ml 0 10 20 30 40 50 60 70 Bias : -0.000 Corr : 1.000 NSE : 1.000 (o) Rnet_soil (W m −2 ) −3−2−101234 Fortran CLM-ml −3 −2 −1 0 1 2 3 4 Bias : 0.000 Corr : 1.000 NSE : 1.000 (p) Obu (m) Figure 5. Oracle validation: time series (left) and scatter compari- son (right) of JAX vs. Fortran outputs for sensible heat (H), latent heat (LE), net radiation (Rn), and GPP across 1488 half-hourly timesteps (May 2007, CHATS7). 4.2. Validation of Backpropagation Capability We validatejax.gradover the CLM-ml-jax column us- ing central finite differences for four key parameters (leaf absorptivityα sw ,V c,max 25 , stomatal slope, and canopy con- ductance) and three outputs (GPP, LE, H), using leaf-level fluxes as differentiable proxies (column-level aggregates such as total ET involve non-differentiable accumulation paths in the test harness; leaf-level fluxes isolate the differ- entiable physics kernel). All four active parameters achieve relative errors less than1× 10 −4 against central FD for all three output fluxes, well within the 1% acceptance thresh- old. A stage-by-stage isolation experiment confirms the gradient path forα sw :d(apar)/d(α sw )matches FD to 1.2 × 10 −9 (solar radiation is exactly differentiable via the Norman two-stream scheme (Norman, 1979) ), and d(agross)/d(α sw )matches to1.8× 10 −7 after the im- plicit function theorem (IFT) fix applied to the stomatal con- ductance solver. The stomatal conductanceg s is determined implicitly using a nonlinear equation; direct differentiation through the iterative solver introduces truncation error, so we apply the IFT to obtain exact gradients at the converged fixed point. The initial four-parameter validation uses the WUE stomatal model; extending to the Medlyn stomatal model (Medlyn et al., 2011) with a 10-parameter sweep, 7 of the 10 parameters are structurally active (the remaining 3 are inactive at this stomatal configuration and produce identically zero gradients). All 7 active parameters achieve relative errors less than 1× 10 −4 . 4.3. Jacobian-Based Sensitivity Analysis We compute∂(GPP,H, LE)/∂θwhereθcomprises five scale parameters on:V c,max 25 , air temperature, shortwave radiation, specific humidityq, and plant area index per canopy layer (dpai). GPP is the dpai-weighted sum of gross photosynthesis over canopy layers; H and LE are 5 Systematic LLM Translation of Legacy Scientific Code to Differentiable Frameworks Vcmax25T_airSW_radq_refdpai GPP H LE 1.051.220.80-1.280.80 1.973.032.281.621.13 2.001.512.101.651.85 (a) Log scale magnitude: log| output/scale| Vcmax25T_airSW_radq_refdpai GPP H LE +0.68-1.00+0.37-0.00+0.37 -0.09-1.00+0.18+0.04+0.01 +0.80+0.26+1.00-0.35+0.56 (b) Normalised sensitivity: output/scale 1 0 1 2 3 log |J| 1.0 0.5 0.0 0.5 1.0 Figure 6.Jacobian-based sensitivity analysis: heatmap of ∂(GPP, H, LE)/∂θto forcing parameters (V c,max 25 , air temper- ature, shortwave radiation, specific humidityq, dpai (plant area index per canopy layer)), computed viajax.jacrevin one backward pass. Stomatal model: WUE. GPP, H, and LE are dpai- weighted canopy sums of leaf-level fluxes. a) Log-scale Jacobian magnitude. b) Row-normalized relative sensitivity (each row di- vided by its maximum absolute value). Tair dominates H and LE; SWrad and V c,max 25 are the leading GPP drivers. dpai-weighted sums of sensible heat and latent heat respec- tively.V c,max 25 is injected via the per-PFT carboxylation capacity dynamic tensor. Scalar-loss Calibration: Measured AD vs. FD In the scalar-loss calibration regime (n = 1output), the crossover betweenjax.gradand central FD is governed by the ratioT b /T f , the cost of one backward pass relative to one forward pass. We measure this directly forp∈1, 2, 3, 5 parameters: the backward pass costsT b ≈ 4.4T f (median across p). Figure 7. AD vs. finite-difference timing as a function of number of parametersp(CHATS7, GPU, scalar GPP loss). Blue bars: T AD (jax.grad, constant inp). Orange bars:T FD (2pforward evaluations, linear inp). AD becomes cheaper than FD atp≥ 3, confirming the theoretical crossover atp cross = T b /(2T f )≈ 2.2. Atp = 5, AD is2.2×faster than FD; this advantage grows linearly with p. Parameter Calibration - Proof of Concept As a proof- of-concept, we calibrate three parameters: stomatal effi- ciencyι, maximum carboxylation rateV c,max 25 , andT ref , against a synthetic target at a single timestep and site (CHATS7, May 2007), bounding the scope of this proof-of- concept to synthetic parameter recovery under controlled conditions. Parameters are initialized from perturbed val- ues; the objective is to recover the known ground-truthθ ∗ from the loss landscape. We compare three optimizers: L- BFGS-B+AD (Zhu et al., 1997) , Adam+AD (Kingma & Ba, 2014) (100 cosine-annealed steps), and gradient-free Nelder-Mead (Nelder & Mead, 1965) (Figure 8). L-BFGS-B+AD reaches machine-precision loss (∼ 10 −19 ) in fewer than 50 evaluations and recovers all three pa- rameters exactly ( ˆ θ/θ ∗ = 1forι,V c,max 25 , andT ref ). Nelder-Mead requires approximately8×as many evalu- ations to achieve the same recovery as a gradient-informed quasi-Newton search compared to a gradient-free simplex. Adam+AD stalls at∼ 10 −2 after 100 cosine-annealed steps, reflecting insufficient iterations for a first-order method on this landscape rather than a failure of differentiability; L- BFGS-B exploits exact curvature information and converges ∼ 17orders of magnitude lower within the same evaluation budget. This efficient gradient-based calibration — previ- ously unavailable in CLM-ml-v2 — is the direct scientific payoff of the five-phase agentic translation pipeline. 050100150200 Evaluations (Gradient and Forward-pass) 10 17 10 14 10 11 10 8 10 5 10 2 10 1 10 4 Loss (a) Convergence Adam+AD L-BFGS-B+AD Nelder-Mead Vcmax 25 T ref 0.8 1.0 1.2 1.4 1.6 Parameter scale factor (b) Parameter recovery = 1 0 (init) (Adam) (L-BFGS-B+AD) (Nelder-Mead) Figure 8. Proof-of-concept parameter recovery forι,V c,max 25 , andT ref (CHATS7, single timestep, May 2007). (a) Loss vs. evalu- ations (gradient and forward-pass counts combined) for Adam+AD (blue solid), L-BFGS-B+AD (orange dash-dot), and Nelder-Mead (green dotted). L-BFGS-B+AD reaches∼ 10 −19 in< 50eval- uations; Nelder-Mead requires≈ 8×more. Adam+AD stalls at ∼ 10 −2 within 100 steps. (b) Parameter recovery ratios ˆ θ/θ ∗ for all three methods; dashed line at1.0indicates exact recovery (θ ∗ ); grey bars show initial perturbation θ 0 . 4.4. Performance Benchmark We characterize the wall-clock cost ofjax.vmapover pa- rameter ensembles relative to sequential Fortran execution, measuring throughput as a function of ensemble sizeN on an Quadro RTX 8000 GPU and sequential CPU, with hardware and software specifications in Appendix A. Throughput scaling. Figure 9 shows amortized cost per sample as a function ofN. GPU cost falls from 24.9ms/sample atN = 1to11.4ms atN = 2,048, plateauing at4.7×below the Fortran sequential cost of ≈54ms/sample. AtN = 2,048, Fortran sequential requires 553s while JAX on GPU requires23.4s, a24×total wall- clock reduction. This makes multi-start ensemble calibra- tion and parameter uncertainty quantification tractable at the single-site scale. 6 Systematic LLM Translation of Legacy Scientific Code to Differentiable Frameworks Scan dispatch overhead. Replacing the Python-level Runge-Kutta loop withjax.lax.scantraces the loop body once at compile time, eliminating per-iteration Python– XLA dispatch overhead. This yields a164×speedup for Euler timestepping and a3,100×speedup for fourth-order RK4. Scalar-loss calibration: AD vs. finite differences. In the scalar-loss regime (n = 1output), the crossover be- tweenjax.gradand central finite differences is gov- erned byT b /T f , the cost of one backward pass relative to one forward pass. We measureT b ≈ 4.4T f (median acrossp; Figure 7). AD becomes cheaper than finite dif- ferences atp ≥ 3, confirming the theoretical crossover at p cross = T b /(2T f )≈ 2.2. Atp = 5, AD is2.2×faster than finite differences; this advantage grows linearly with p. Numerical precision. The model is dominated by memory-bandwidth-bound operations (exp/log/sqrtin stomatal and radiative transfer kernels);float32provides no throughput benefit overfloat64on modern GPUs (≤5% difference across all N ; Table 3). Compile-time scaling.Applying JIT-compiledvmapexe- cution on CPU requires XLA to materialise a flatO(N×M) LLVM compilation graph at trace time, whereMis the num- ber of model operations; atN = 128this exhausts virtual address space regardless of available physical RAM, as the latter affects only the time to failure, not the outcome. On GPU,vmaptiles across the batch dimension rather than un- rolling it, keeping compile timeO(M)independently ofN. GPU execution is therefore a hard requirement for ensemble sizes N > 1, not merely a performance optimization. 183212851210242048 Ensemble size N 2 × 10 1 3 × 10 1 4 × 10 1 6 × 10 1 ms / sample 24.9 13.1 11.8 11.7 11.5 11.6 11.4 24.3 19.8 18.5 18.2 58.5 51.4 52.2 54.054.054.0 54.0 (a) Amortized cost per sample JAX GPU (Quadro RTX 8000) JAX CPU Fortran sequential (N runs) 183212851210242048 Ensemble size N 1 2 3 4 5 Reference ms/sample / GPU ms/sample 0.98× 1.55× 2.4× 4.7× (b) GPU speedup over reference (per sample) JAX GPU speedup vs JAX CPU JAX GPU speedup vs Fortran seq Figure 9. Throughput scaling of CLM-ML-JAX with ensemble size N(CHATS7 walnut orchard, 46 layers, Quadro RTX 8000).Nis the number of independent forward passes—each with a distinct parameter vector—run simultaneously, as needed for ensemble calibration or uncertainty quantification. (a) GPU amortized cost falls from 24.9 to 11.4 ms/sample asNgrows from 1 to 2048; Fortran executes serially and stays flat at≈54 ms/sample. (b) Per- sample GPU speedup over Fortran reaches4.7×atN = 2048and continues to grow with N . 5. Discussion 5.1. Scientific Implications The oracle validation establishes thatclm-ml-jaxis a scientifically equivalent re-implementation of the Fortran model, preserving all PFT parameterizations without re- calibration. Verified reverse-mode gradients through the full 46-layer coupled physics column—encompassing IFT solvers, a Runge-Kutta sub-step loop, and tightly coupled leaf-canopy-turbulence physics—confirm that autodiff is accurate through non-trivial numerical structures. The Ja- cobian analysis yields a concrete scientific result: air tem- perature and shortwave radiation dominate sensitivities of GPP and energy fluxes at the CHATS site, withV c,max 25 ex- erting moderate but distinct control over GPP. Full-column parameter sensitivities, previously intractable in a single backward pass, directly inform which observations most constrain model predictions at a given site. End-to-end differentiability further positionsclm-ml-jaxfor hybrid physics-ML architectures in which neural parameterizations of uncertain processes are jointly trained with the physical column (Aboelyazeed et al., 2023; Shen et al., 2023). 5.2. Limitations and Future Work Headless operation requires skipping permissions; opera- tors should apply sandboxing and reviewCHANGELOG.md before deployment. The calibration experiment remains a proof-of-concept: 3 parameters, a single site, and a single timestep. Future work includes couplingclm-ml-jaxto JAX-based atmospheric models (H ̈ afner et al., 2018) for fully differentiable land–atmosphere simulation, applying the translation methodology to CLM5 and the CLUBB con- vection scheme in CAM, and training neural parameteriza- tions end-to-end within the differentiable column. 6. Conclusion We have proven three things. First, that a 19,000-line val- idated Fortran land surface model can be translated to a numerically equivalent, fully differentiable JAX implemen- tation via a five-phase agentic pipeline without manual reim- plementation and with confirmed numerical equivalence. Second, we validated the model’s end-to-end differentia- bility: gradients pass finite-difference validation through the full coupled physics column, including IFT-corrected iterative solvers, and the Jacobian recovers physically in- terpretable sensitivities in a single backward pass. Third, gradient-based parameter optimization (L-BFGS-B) recov- ers three land surface parameters to machine precision in fewer than 50 evaluations 8×faster than gradient-free search on the same problem. Beyond the model itself, this paper is a case study in AI as a tool for scientific code modernization. We present a 7 Systematic LLM Translation of Legacy Scientific Code to Differentiable Frameworks five-step code translation workflow with autonomous multi- day repair loops, domain-specific oracle validation, and a version-controlled audit trail. The 81.2% first-attempt au- tonomous resolution rate and∼10–15 hours of total human oversight across 46 sessions quantify what co-authorship looks like in practice for legacy scientific code. We release the model and the methodology infrastructure to support the community in applying this approach to other ESM components. Acknowledgements We acknowledge funding from NSF through the Learning the Earth with Artificial intelligence and Physics (LEAP) Science and Technology Center (STC) (Award #2019625). Impact Statement Differentiable land surface models enable gradient-based parameter estimation and sensitivity analysis that were pre- viously intractable, with direct relevance to reducing uncer- tainty in land-atmosphere carbon and energy flux projec- tions. The calibration experiment presented here is a proof- of-concept at single-site, single-timestep scale; extension to multi-site, long-period calibration could meaningfully reduce projection uncertainty in carbon cycle feedbacks. The agentic pipeline introduces governance considerations specific to autonomous scientific software development. As agentic pipelines assume roles traditionally held by domain scientists, the field requires audit trails and validation stan- dards commensurate with the scientific trust placed in the resulting artifacts. The state-document infrastructure pre- sented here represents a step toward reproducible and au- ditable agentic workflows in scientific computing. References Aboelyazeed, D., Xu, C., Hoffman, F. M., Liu, J., Jones, A. W., Rackauckas, C., Vangala, A., Sanyal, S., et al. A differentiable, physics-informed ecosystem modeling and learning framework for large-scale inverse problems. Biogeosciences, 20:2671–2692, 2023. doi: 10.5194/ bg-20-2671-2023. Bai, S., Kolter, J. Z., and Koltun, V. Deep equilibrium models. In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019. Baydin, A. G., Pearlmutter, B. A., Radul, A. A., and Siskind, J. M. Automatic differentiation in machine learning: a survey. Journal of Machine Learning Research, 18(153): 1–43, 2018. Bezgin, D. A., Buhendwa, A. B., and Adams, N. A. Jax- fluids: A fully-differentiable high-order computational fluid dynamics solver for compressible two-phase flows. Computer Physics Communications, 282:108527, 2023. ISSN 0010-4655.doi: 10.1016/j.cpc.2022.108527. URLhttps://w.sciencedirect.com/ science/article/pii/S0010465522002466. Blondel, M., Berthet, Q., Cuturi, M., Frostig, R., Hoyer, S., Llinares-L ́ opez, F., Pedregosa, F., and Vert, J.-P. Efficient and modular implicit differentiation. In Advances in Neu- ral Information Processing Systems, volume 35, 2022a. URL https://arxiv.org/abs/2105.15183. Blondel, M., Berthet, Q., Cuturi, M., Frostig, R., Hoyer, S., Llinares-L ́ opez, F., Pedregosa, F., and Vert, J.-P. Efficient and modular implicit differentiation. Advances in Neural Information Processing Systems, 35, 2022b. Bonan, G. B., Lombardozzi, D. L., Wieder, W. R., Oleson, K. W., Lawrence, D. M., Hoffman, F. M., and Collier, N. Moving beyond the incorrect but useful paradigm: reevaluating big-leaf and multilayer plant canopies to model biosphere–atmosphere fluxes – a review. Agricul- tural and Forest Meteorology, 306:108435, 2021. doi: 10.1016/j.agrformet.2021.108435. Butcher, J. C. On runge-kutta processes of high order. Jour- nal of the Australian Mathematical Society, 4(2):179–194, 1964. doi: 10.1017/S1446788700023387. Campagne, J.-E., Lanusse, F., Zuntz, J., Boucaud, A., Casas, S., Karamanis, M., Kirkby, D., Lanzieri, D., Li, Y., and Peel, A. Jax-cosmo: An end-to-end differentiable and gpu accelerated cosmology library. The Open Journal of Astrophysics, 6, 2023. doi: 10.21105/astro.2302.05163. URL https://arxiv.org/abs/2302.05163. Cowan, I. R. and Farquhar, G. D. Stomatal function in relation to leaf metabolism and environment. Symposia of the Society for Experimental Biology, 31:471–505, 1977. Davenport, E. H., Madan, J. V., Gjini, R., Brzenski, J., Ho, N., Hsu, T., Liang, Y., Liu, Z., Manivannan, V., Pham, E., Vutukuru, R., Williams, A. I. L., Yang, Z., Yu, R., Lutsko, N. J., Hoyer, S., and Watson-Parris, D. Jcm v1.0: A differentiable, intermediate-complexity at- mospheric model. EGUsphere, 2026. doi: 10.5194/ egusphere-2025-6266. URLhttps://doi.org/10. 5194/egusphere-2025-6266. Preprint. Dhruv, A. and Dubey, A.Leveraging large language models for code translation and software development in scientific computing.In Platform for Advanced Scientific Computing Conference (PASC), 2024. doi: 10.1145/3732775.3733572. arXiv:2410.24119. ECMWF. Modernisation of the integrated forecasting sys- tem. ECMWF Newsletter, 182, 2024. 8 Systematic LLM Translation of Legacy Scientific Code to Differentiable Frameworks Evensen, G. The ensemble Kalman filter: theoretical formu- lation and practical implementation. Ocean Dynamics, 53:343–367, 2003. doi: 10.1007/s10236-003-0036-9. Farquhar, G. D., von Caemmerer, S., and Berry, J. A. A biochemical model of photosynthetic CO 2 assimilation in leaves of C3 species. Planta, 149(1):78–90, 1980. doi: 10.1007/BF00386231. Gelbrecht, M., White, A., Bathiany, S., and Boers, N. Differ- entiable programming for Earth system modeling. Geo- scientific Model Development, 16(11):3123–3135, 2023. doi: 10.5194/gmd-16-3123-2023. Giering, R. and Kaminski, T. Recipes for adjoint code con- struction. ACM Transactions on Mathematical Software, 24(4):437–474, 1998. doi: 10.1145/293686.293695. Gupta, A. et al.From legacy Fortran to portable Kokkos: An autonomous agentic AI workflow, 2025. arXiv:2509.12443; poster SC25. H ̈ afner, D., Jacobsen, R. L., Eden, C., Kristensen, M. R. B., Jochum, M., Nuterman, R., and Vinter, B. Veros v0.1 – a fast and versatile ocean simulator in pure Python. Geo- scientific Model Development, 11(8):3299–3312, 2018. doi: 10.5194/gmd-11-3299-2018. JAX backend released 2021; original model 2018. Harman, I. N. and Finnigan, J. J. A simple unified theory for flow in the canopy and roughness sublayer. Boundary- Layer Meteorology, 123(2):339–363, 2007. doi: 10.1007/ s10546-006-9145-6. Harman, I. N. and Finnigan, J. J. Scalar concentration profiles in the canopy and roughness sublayer. Boundary- Layer Meteorology, 129:323–351, 2008. doi: 10.1007/ s10546-008-9328-4. Huntley, G.Everything is a ralph loop.https:// ghuntley.com/loop/ , January 2026. Accessed: 2026. IBM Research. IBM unveils watsonx generative AI capabil- ities to accelerate mainframe application modernization, 2023. Jiang, P., Shen, C., et al. JAX-CanVeg: A differentiable land surface model. Water Resources Research, 61, 2025. doi: 10.1029/2024WR038X. Jimenez, C. E., Yang, J., Wettig, A., Yao, S., Pei, K., Press, O., and Narasimhan, K. SWE-bench: Can language mod- els resolve real-world GitHub issues? In International Conference on Learning Representations (ICLR), 2024. arXiv:2310.06770. Jin, M. et al. Fortran2CPP: Automating Fortran-to-C++ migration using LLMs via multi-turn dialogue and dual- agent integration. arXiv preprint arXiv:2412.19770, 2024. URL https://arxiv.org/abs/2412.19770. Kadosh, T. et al. LLM-assisted translation of legacy FOR- TRAN codes to C++: A cross-platform study. arXiv preprint arXiv:2504.15424, 2025.URLhttps:// arxiv.org/abs/2504.15424. Kahn, A. B. Topological sorting of large networks. Commun. ACM, 5(11):558–562, November 1962. ISSN 0001-0782. doi: 10.1145/368996.369025. URLhttps://doi. org/10.1145/368996.369025. Kingma, D. P. and Ba, J. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014. Kochkov, D., Yuval, J., Langmore, I., Norgaard, P., Smith, J., Mooers, G., Kl ̈ ower, M., Lottes, J., Rasp, S., D ̈ uben, P., Hatfield, S., Battaglia, P., Sanchez-Gonzalez, A., Will- son, M., Brenner, M. P., and Hoyer, S. Neural general circulation models for weather and climate. Nature, 632: 1060–1066, 2024. doi: 10.1038/s41586-024-07744-y. Kucharski, F., Molteni, F., King, M. P., Farneti, R., Kang, I.-S., and Feudale, L. On the need of intermediate com- plexity general circulation models: a ‘SPEEDY’ exam- ple. Bulletin of the American Meteorological Society, 94(1):25–30, 2013. doi: 10.1175/BAMS-D-11-00238. 1. URLhttp://journals.ametsoc.org/doi/ pdf/10.1175/BAMS-D-11-00238.1. Kwa, A., Clymer, J., McDougall, C., Greenblatt, R., Deni- son, C., Roger, G., and Bowman, S. R. Measuring the persuasiveness of language model generated free-text ar- guments. In arXiv preprint, 2025. The term “agentic laziness” referring to premature task completion by LLM agents is attributed to this work. Le Dimet, F.-X. and Talagrand, O. Variational algorithms for analysis and assimilation of meteorological obser- vations: theoretical aspects. Tellus A: Dynamic Mete- orology and Oceanography, 38(2):97–110, 1986. doi: 10.1111/j.1600-0870.1986.tb00459.x. Li, Z. et al. Fortran2CPP: Automating Fortran-to-C++ trans- lation using LLMs via multi-turn dialogue and dual-agent integration. arXiv preprint, 2024. arXiv:2412.19770. Medlyn, B. E., Duursma, R. A., Eamus, D., Ellsworth, D. S., Prentice, I. C., Barton, C. V. M., Crous, K. Y., De An- gelis, P., Freeman, M., and Wingate, L. Reconciling the optimal and empirical approaches to modelling stomatal conductance. Global Change Biology, 17(6):2134–2144, 2011. doi: 10.1111/j.1365-2486.2010.02375.x. 9 Systematic LLM Translation of Legacy Scientific Code to Differentiable Frameworks M ́ endez, M. and Tinetti, F. G. Legacy Fortran software: Ap- plying syntactic metrics to global climate models. Journal of Computer Science & Technology, 2016. M ́ endez, M., Tinetti, F. G., and Overbey, J. L. Climate models: Challenges for Fortran development tools. In 2nd International Workshop on Software Engineering for High Performance Computing in Computational Science and Engineering (SE-HPCCSE). IEEE, 2014. Meunier, E., Ouala, S., Frezat, H., Sommer, J. L., and Fa- blet, R. Towards fully differentiable neural ocean model with veros, 2025. URLhttps://arxiv.org/abs/ 2511.17427. Monin, A. S. and Obukhov, A. M. Basic laws of turbulent mixing in the surface layer of the atmosphere. Trudy Geofiz. Inst. Akad. Nauk SSSR, 24(151):163–187, 1954. Nelder, J. A. and Mead, R. A simplex method for func- tion minimization. The Computer Journal, 7(4):308–313, January 1965. ISSN 1460-2067. doi: 10.1093/comjnl/ 7.4.308.URLhttp://dx.doi.org/10.1093/ comjnl/7.4.308. Neumann, P., D ̈ uben, P., Adamidis, P., Bauer, P., Br ̈ uck, M., Kornblueh, L., Klocke, D., Stevens, B., Wedi, N., and Biercamp, J. Assessing the scales in numerical weather and climate predictions: will exascale be the rescue? Philosophical Transactions of the Royal Society A, 377 (2142):20180148, 2019. doi: 10.1098/rsta.2018.0148. Nichols, J. A. et al.From legacy Fortran to portable Kokkos: An autonomous agentic AI workflow. arXiv preprint arXiv:2509.12443, 2025.URLhttps:// arxiv.org/abs/2509.12443. Norman, J. M. Modeling the complete canopy. In Barfield, B. J. and Gerber, J. F. (eds.), Modification of the Aerial Environment of Crops, p. 249–277. American Society of Agricultural Engineers, St. Joseph, Michigan, 1979. Pan, R., Ibrahimzada, A. R., Krishna, R., Sankar, D., Wassi, L. P., Merler, M., Sobolev, A., Pavuluri, M., Sinha, S., and Jabbarvand, R. Lost in translation: A study of bugs introduced by large language models while translating code. In Proceedings of the 46th IEEE/ACM International Conference on Software Engineering (ICSE), 2024. URL https://arxiv.org/abs/2308.03109. Park, J. S., O’Brien, J. C., Cai, C. J., Morris, M. R., Liang, P., and Bernstein, M. S. Generative agents: Interactive simulacra of human behavior. In Proceedings of the 36th Annual ACM Symposium on User Interface Software and Technology, UIST ’23. ACM, 2023. doi: 10.1145/ 3586183.3606763. Patton, E. G., Horst, T. W., Sullivan, P. P., Lenschow, D. H., Oncley, S. P., Brown, W. O. J., Burns, S. P., Guenther, A. B., Held, A., Karl, T., Mayor, S. D., Rizzo, L. V., Spuler, S. M., Sun, J., Turnipseed, A. A., Allwine, E. J., Edburg, S. L., Lamb, B. K., Avissar, R., Calhoun, R. J., Kleissl, J., Massman, W. J., Paw U, K. T., and Weil, J. C. The canopy horizontal array turbulence study. Bulletin of the American Meteorological Society, 92(5):593–611, 2011. doi: 10.1175/2010BAMS2614.1. Pietrini, M. et al. Bridging the gap: LLM-assisted trans- lation of legacy Fortran code to modern C++. In SC Workshops, 2024. Qu, Y., Bhouri, M. A., and Gentine, P. Joint parameter and parameterization inference with uncertainty quan- tification through differentiable programming. arXiv preprint arXiv:2403.02215, 2024.URLhttps:// arxiv.org/abs/2403.02215. Ranasinghe, K. et al. LLM-assisted modernization of legacy Fortran scientific code. arXiv preprint arXiv:2501.x, 2025. Raoult, N. M., Jupp, T. E., Cox, P. M., and Luke, C. M. ad- JULES: a variational data assimilation system for JULES. Geoscientific Model Development, 9(7):2683–2696, 2016. doi: 10.5194/gmd-9-2683-2016. Raoult, N. M. et al. Parameter estimation in land surface models: Challenges and opportunities with data assimila- tion and machine learning. Journal of Advances in Model- ing Earth Systems, 2025. doi: 10.1029/2024MS004733. Rasp, S., Pritchard, M. S., and Gentine, P. Deep learning to represent subgrid processes in climate models. Pro- ceedings of the National Academy of Sciences, 115(39): 9684–9689, 2018. doi: 10.1073/pnas.1810286115. Shen, C., Appling, A. P., Gentine, P., Bandai, T., Gupta, H., Tartakovsky, A., Baity-Jesi, M., Fenicia, F., Kifer, D., Li, L., et al. Differentiable modelling to unify machine learning and physical models for geosciences. Nature Reviews Earth & Environment, 4:552–567, 2023. doi: 10.1038/s43017-023-00450-9. Talagrand, O. and Courtier, P. Variational assimilation of meteorological observations with the adjoint vorticity equation. I: Theory. Quarterly Journal of the Royal Me- teorological Society, 113(478):1311–1328, 1987. doi: 10.1002/qj.49711347812. Tian, W., Yu, H., Zhao, S., Cao, Y., Yi, W., Xu, J., and Nan, Z. NoahPy: a differentiable Noah land surface model for simulating permafrost thermo-hydrology. Geo- scientific Model Development, 19:57–72, 2026. doi: 10.5194/gmd-19-57-2026. 10 Systematic LLM Translation of Legacy Scientific Code to Differentiable Frameworks Vrugt, J. A., ter Braak, C. J. F., Diks, C. G. H., Robinson, B. A., Hyman, J. M., and Higdon, D. DREAM(D): an adaptive Markov chain Monte Carlo simulation algorithm to solve discrete, noncontinuous, and combinatorial pos- terior parameter estimation problems. Hydrology and Earth System Sciences, 13(12):1913–1940, 2009. doi: 10.5194/hess-13-1913-2009. Wang, Y. and Frankenberg, C. CliMA Land: a differentiable land surface model for the CliMA earth system model. Journal of Advances in Modeling Earth Systems, 15, 2023. doi: 10.1029/2023MS003583. Zhou, A., Hawkins, L., and Gentine, P.A proof of concept: Using ChatGPT to translate and modernize an earth system model from Fortran to Python/JAX. arXiv preprint arXiv:2405.00018, 2024. URLhttps: //arxiv.org/abs/2405.00018. Zhu, C., Byrd, R. H., Lu, P., and Nocedal, J. Algorithm 778: L-bfgs-b: Fortran subroutines for large-scale bound- constrained optimization. ACM Transactions on mathe- matical software (TOMS), 23(4):550–560, 1997. 11 Systematic LLM Translation of Legacy Scientific Code to Differentiable Frameworks A. Hardware and Software Specifications Performance Benchmarking experiments were run on NVIDIA A40 GPU and CPU nodes. Fortran reference timing measured on NCAR Derecho (Intel Ice Lake CPUs). All other experiments were run on an NVIDIA A100-PCIE-40GB GPU—software: JAX 0.9.2, Python 3.11, CUDA 12.8, gfortran 12. Table 1. Common Fortran patterns encountered in CLM-ml-v2, the general differentiable-programming concept each requires, and the JAX-specific equivalent. The General concept column is framework-agnostic; teams targeting PyTorch, Julia Flux, or other frameworks can apply the same concept using their framework’s primitives. Fortran PatternGeneral conceptJAX EquivalentNotes DO loop over sub-stepsStatic loop unrolling for gradi- ent tracing jax.lax.scanPythonforloops produce unrolled graphs;lax.scan traces once In-place array mutationFunctional (copy-on-write) up- date x.at[i].set(v)JAX arrays are immutable; all target frameworks require ex- plicit copy COMMONblocks / module state Explicit argument threadingExplicit argument passing; NamedTuple state No global mutable state in any differentiable framework IF/ELSEwith loop-carried deps Differentiable conditional jax.lax.cond / jax.lax.switch Avoids Python-level branch- ing in JIT-compiled code Implicit typingExplicit dtype annotation jnp.float64 throughoutType consistency required for correct gradient accumulation Conditionalbranchwith x ** n or 1/x, x = 0 Branch-safe guard jnp.maximum(x,1e-30) before power/division Both branches evaluated in backward pass;0×∞ = NaN Iterative root-finding (DOcon- vergence loop) Custom adjoint via implicit differentiation Newton-refinement IFT wrap- per Loop-through-convergence gradients incorrect; IFT gives exact adjoint Module-level mutable arraysImmutable state containerNamedTuplewith .replace(field=value) Enables functional differentia- tion; PyTorch:dataclass; Julia: struct B. Static Analysis: CLM-ml-v2 B.1. Dependency Graph Analysis The inter-moduleUSEdependency structure resolves to a directed graphG M of 76 nodes and 315 directed edges. The graph density isρ = 315/(76× 75) = 0.0553, indicating moderately sparse but non-trivial inter-module coupling, with an average of 8.29USEreferences per module.G M is acyclic: the absence of circular dependencies yields a strict directed acyclic graph (DAG), which is a necessary precondition for safe topological translation ordering. 295.30295.35295.40295.45295.50 0 5 10 15 20 Height (m) (a) Air temperature (K) Fortran JAX 0.20.40.60.81.01.2 (b) Wind speed (m s −1 ) 4.85.05.25.45.6 (c) Specific humidity (g kg −1 ) 0.00.51.01.52.0 (d) CO 2 (ppm) Figure 10. CHATS7 May 1, 2007: Canopy profiles at noon (timestep 24): JAX vs. Fortran reference across all 46 canopy layers. Variables shown include air temperature, wind speed, specific humidity, CO 2 concentration, and leaf-level photosynthesis and stomatal conductance. 12 Systematic LLM Translation of Legacy Scientific Code to Differentiable Frameworks B.2. Implications for Translation Ordering and LLM Context Budgeting The strict DAG structure directly enables a topologically safe translation order. Modules are processed in increasing depth order, ensuring every dependency of a given module is translated, with its numerical parity verified against Fortran golden data, before that module’s own LLM translation prompt is constructed. The 12 depth levels partition all 73 modules into parallel batches. The maximum achievable parallelism is 15 concurrent translations at depth 2. Context budgeting. Module size varies considerably across the codebase, producing a wide range of per-call token demands. Deep-stack modules require loading the signatures of all previously translated dependencies alongside the source itself.MLCanopyFluxesMod.F90is representative of this challenge: at dependency depth 7 with 40 declared internal dependencies, a single translation call must incorporate up to 40 JAX module signatures in addition to a 1,511-line source file. Static analysis identifies high-effort modules which receive targeted treatment: reference I/O data are generated across physically diverse atmospheric forcing conditions for robust numerical validation, and a higher-capacity model (Claude Opus 4.6) is allocated for agentic translation and repair. C. Characterization of Agentic Differentiability Repair We report statistics of the Differentiability Optimization phase (Phase 5, Part 2), derived fromCHANGELOG.md. These statistics characterize the repair effort a team may encounter when applying this methodology to a comparably scaled ESM. Bug taxonomy. Table 2 classify 48 documented bugs into 9 mechanistic types. Types T1 and T2 jointly account for 54.2T1 bugs follow a single recurring template: insertingjnp.maximum(x,1e-30)prior to any division or fractional power within a guarded branch. This fix was applied at 35 or more sites across 9 modules. T2 bugs decompose into four sub-patterns. (a) Pythonfloat()casts silently break the JAX autodiff tape by detaching values from the trace. (b) Physics state is overwritten byMLCanopyFluxes.initbefore the traced computation executes. (c) Functions decorated with@jax.jitcapture module-global parameter values as XLA compile-time constants, rendering subsequent module-level mutations invisible to JAX. (d)from module import Xcreates a local binding to the original object that is not updated when the module variable is reassigned. Sub-patterns (c) and (d) were particularly costly to diagnose: together they required three sessions and two failed repair attempts before the root cause was identified (bugs B36–B37, B39). Module bug density.MLLeafPhotosynthesisModandMLCanopyTurbulenceModeach contain 8 bugs, jointly accounting for 33% of the total (Figure 2). Both modules contain iterative solvers—the photosynthesisc i scan and the Obukhov length secant method—that generate T5 gradient explosions, and both produce the highest per-sub-step arithmetic density driving T1 patterns.MLLeafPhotosynthesisModadditionally couples parameter tables viaMLpftcon injection, exhibiting T2, T1, T3, and T5 bugs;MLCanopyTurbulenceModexhibits T1, T3, T5, and T6. Teams translating comparable coupled land surface models should pre-audit these module classes prior to invoking jax.grad. 13 Systematic LLM Translation of Legacy Scientific Code to Differentiable Frameworks Table 2. Differentiability bug taxonomy. Avg. att.: mean agent debugging iterations per bug instance. n: total bugs of this type. All statistics from CHANGELOG.md sessions 1–46 (April 1 – May 8 2026). CodeTypen%Avg. att.Root mechanism T1NaN Gradient1735.41.0JAXevaluatesbothjnp.where branches; masked branch withx=0gives 0×∞=NaN T2Zero/Wrong Gradient918.81.3 Parameter not reaching JAX trace (Python cast, state overwrite, JIT constant, local binding) T3XLA Recompilation24.21.0JIT cache miss per call (unstable closure, missing lrucache) T4Memory / OOM36.31.7 Trace graph too large for device (gradient unrolling, CPU vmap unroll, large tensor) T5Gradient Explosion36.32.0JacobianaccumulationthroughN- iteration solver: |J| N →∞ T6Device–Host Sync36.31.0 np.asarray(jaxarr)orfloat() inside hot loops forces GPU→CPU copies T7Optimization Algorithm48.31.5Adamβ 2 stall, step-index arithmetic, un- derdetermined system T8Crash / Compile36.31.3XLA backend bug, GPU contention, par- allel agent race condition T9Diagnostic Reliability36.31.3 FDεinstability,spvalcontamination in loss, wrong timing barrier All481001.3381.2% fixed on first attempt; 18.8% re- quired human re-specification Table 3. Float32 vs. float64 throughput on the NVIDIA A40 GPU. Values are amortized cost per ensemble member (ms/sample), median of 5 repeats. Float32 offers no measurable advantage; the model is memory-bandwidth limited. Nf64 (ms/sample)f32 (ms/sample)f32/f64 ratio 129.227.60.95 813.513.51.00 3211.511.91.04 12811.211.51.03 51211.411.61.02 102411.811.40.97 204810.911.51.05 14 Systematic LLM Translation of Legacy Scientific Code to Differentiable Frameworks Apr 6 Apr 13Apr 20Apr 27 May 4 0 4 8 12 Bugs fixed per session Phase 5b differentiability repair campaign (1 Apr 8 May 2026, 46 sessions) NaN elim. (T1) Zero-grad root cause IFT + param injection Calibration (T7) Bugs fixed (left axis) Cumulative (right axis) 0 10 20 30 40 50 Cumulative bugs fixed First jax.grad 12 NaN bugs 5 params GPU-verified 7-param Jacobian (all non-zero) p=10 calibration + Tikhonov Figure 11. Differentiability repair campaign timeline (April 1 – May 8, 2026, 46 sessions). Stems show bugs fixed per session. The dashed line shows the cumulative count (right axis). Background shading marks the Four dominant phases: NaN elimination (orange), zero-gradient root cause (salmon), IFT and parameter-injection fixes (blue), calibration and optimization (green). Full 7-parameter Jacobian with all non-zero columns was confirmed on April 27, 27 days after the first jax.grad call. 15