CAdenv1.f

C=====================================================================
C CA-DENDRITIC MICROSTRUCTURES
C Two-Dimensional Cellular Automaton Model for Primary Dendritic
C Solidification in Binary Alloys
C
C Copyright (c) 2026 David R. Johnson
C Licensed under the MIT License (see LICENSE file for details)
C
C Compile:
C    gfortran -O3 -Wall -Wextra -o CAsim <source>.f
C
C Note:
C    This code is intended for teaching, experimentation, and
C    preservation of a legacy cellular automaton solidification model.
C    It is not a validated solidification simulation package.
C=====================================================================
C OVERVIEW
C --------
C This program implements a 2-D cellular automaton (CA) model for
C primary dendritic solidification in binary alloys.
C
C Scope:
C   Primary solidification using an effective binary alloy.
C   Multicomponent alloys may be approximated using an effective
C   binary description of the primary solidification path.
C
C   Full mushy-zone evolution and multiphase reactions are outside the
C   intended model scope. The model is intended to generate
C   mass-balanced primary microstructure states that may be used as
C   starting points for subsequent higher-fidelity simulations.
C
C Main features:
C   - Explicit liquid-phase solute diffusion on the CA grid
C   - Local equilibrium partitioning during interface solidification
C   - 4-fold anisotropic growth through orientation-based weighting
C   - Analytical tip-scale guidance from KGT/marginal stability models
C   - Macro-scale thermal histories coupled to local thermal analysis
C
C APPROACH
C --------
C The CA grid resolves evolving morphology and liquid-phase transport,
C while unresolved dendrite-tip physics is supplied analytically
C through the KGT/marginal-stability framework.
C
C Interface growth couples these two descriptions through:
C   - Local equilibrium partitioning
C   - Orientation-dependent capture rules
C   - Staged solute release from unresolved local reservoirs
C
C INPUT & OUTPUT WORKFLOW
C -----------------------
C Input files:
C   xparam.txt   : simulation parameters
C   xgrid.txt    : starting CA grid
C   xthermal.csv : imposed macro thermal history
C Output files:
C   data.csv     : runtime diagnostics
C   grid.txt     : final CA grid state
C   outN.ppm     : optional PPM movie frames
C
C BOUNDARY CONDITIONS
C -------------------
C   IBCTYPE = 0 : closed   X, closed   Y
C   IBCTYPE = 1 : periodic X, closed   Y (default)
C   IBCTYPE = 2 : closed   X, periodic Y
C   IBCTYPE = 3 : periodic X, periodic Y
C
C THERMAL INPUT
C -------------
C The macro-thermal field is read from THFIL, normally xthermal.csv,
C by LOADTHERMCSV().
C
C IEMBED controls the initial thermal alignment:
C  IEMBED = 1: use imposed thermal history directly
C  IEMBED = 0: shift imposed thermal history so the initial seed/front
C               aligns with TL0 - DTSEED (standalone, default)
C
C NUMERICAL NOTES & CONSTRAINTS
C -----------------------------
C  - ZUNITS controls subcell resolution within each CA cell
C  - Explicit solute diffusion uses ALPHA = D*DT/DX^2
C  - Stability and resolution depend on DCOEFF, VEL, DX, timestep size,
C    thermal history, and interface update frequency
C  - DFRZ is the local freezing/capture scale used by the interface
C    update; it is evaluated internally from the current KGT-based
C    tip-scale estimate.
C
C DIAGNOSTICS
C -----------
C   - Tip velocity
C   - Tip temperature
C   - KGT curvature/radius
C   - Peclet number
C   - Supersaturation comparison (OmegaCA vs OmegaKGT)
C
C REFERENCES & MODEL BASIS
C ------------------------
C Cellular automaton solidification framework
C -------------------------------------------
C  Rappaz, M. and Gandin, C.-A.,
C  "Probabilistic modelling of microstructure formation in
C  solidification processes,"
C  Acta Metallurgica et Materialia, 41 (1993) 345-360.
C
C Dendrite growth theory, marginal stability, and KGT routines
C ------------------------------------------------------------
C  Kurz, W., Giovanola, B., and Trivedi, R.,
C  "Theory of microstructural development during rapid solidification,"
C  Acta Metallurgica, 34 (1986) 823-830.
C
C  W. Kurz, D.J. Fisher, "Fundamentals of Solidification, 4th Ed.,"
C  Trans Tech Publications, 1998, pp. 242-244.
C
C Ivantsov diffusion solution and supersaturation relationships
C --------------------------------------------------------------
C  Ivantsov, G.P.,
C  "Temperature field around a spherical, cylindrical, and
C   needle-shaped crystal growing in a supercooled melt,"
C  Doklady Akademii Nauk SSSR, 58 (1947) 567-569.
C
C  Walter Gautschi and William F. Cahill, "Exponential Integral and
C  Related Functions," in Handbook of Mathematical Functions, 
C  Milton Abramowitz and I.A. Stegun, eds., Dover Publications
C  New York, 1965, pp. 227-254.
C
C Solute trapping & velocity-dependent partition coefficient
C ----------------------------------------------------------
C  Aziz, M.J.,
C  "Model for solute redistribution during rapid solidification,"
C  Journal of Applied Physics, 53 (1982) 1158-1168.
C
C Gibbs-Thomson capillarity relation
C ----------------------------------
C  Standard Gibbs-Thomson curvature undercooling relation used in
C  dendritic solidification theory and KGT marginal-stability models.
C
C Notes
C -----
C  Specific routine references are noted directly in the in-line
C  comments. Otherwise, the present implementation combines analytical
C  dendrite-tip physics, explicit liquid-phase solute diffusion, and
C  cellular automaton interface capture/growth logic.
C
C  The staged-solute transport and DFRZ update rule used here are part
C  of the present CA model.
C
C  Fixed code parameters:
C    GRID slot names and the CA subcell scale are defined consistently
C    throughout the code using PARAMETER constants:
C
C      INTEGER LCELLS, LIQUID, SOLID, PHASE, TCELL, CELLDIR
C      PARAMETER (LCELLS=1, LIQUID=2, SOLID=3, PHASE=4)
C      PARAMETER (TCELL=5, CELLDIR=6)
C
C      REAL ZUNITS
C      PARAMETER (ZUNITS=1000.0)
C=====================================================================
      PROGRAM CADEN
      IMPLICIT NONE
C Size limits
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
C Sparse-list capacities
C   These scale with maximum compiled grid size.
C   Active use is normally much smaller and is guarded at runtime.
      INTEGER MAXUPD, MAXREJ, MAXTOUCH
      PARAMETER (MAXUPD = NXMAX*NYMAX/5)
      PARAMETER (MAXREJ = NXMAX*NYMAX/50)
      PARAMETER (MAXTOUCH = NXMAX*NYMAX/5)
C Grid slot tags
C Fixed CA & geometry constants
      INTEGER BOX2, TOPEDGE
      REAL TOL
      PARAMETER (BOX2=3, TOPEDGE=60)
      PARAMETER (TOL=1.0E-6)
C Boundary-condition type
C   0 = closed   X, closed   Y
C   1 = periodic X, closed   Y   current default
C   2 = closed   X, periodic Y
C   3 = periodic X, periodic Y
      INTEGER IBCTYPE
      DATA IBCTYPE/1/
C Main arrays
      REAL GRID(NXMAX, NYMAX, 6, 2)
      REAL XADD(NXMAX, NYMAX)
      REAL XRES(NXMAX, NYMAX)
      REAL XRESADD(NXMAX, NYMAX)
      REAL XRND(NXMAX, NYMAX)
C Domain size
      INTEGER XLEN, YLEN
      DATA XLEN/200/, YLEN/400/
C Screen & loop control
      INTEGER I, J
      INTEGER SCRN, SCUR, SNEW
      INTEGER IEMBED
      INTEGER DONE
      INTEGER NTOK, NPRINTS
      DATA IEMBED/0/
      DATA DONE/0/
C Boundary & seeding controls
      INTEGER WTOP, WBOT, WLEFT, WRIGHT
      DATA WTOP/1/, WBOT/1/, WLEFT/1/, WRIGHT/1/
C I/O control
      INTEGER IUNIT, DSPNUM, ISTOPMODE
      CHARACTER*128 PARFIL, GRIDFIL
      PARAMETER (IUNIT=21, DSPNUM=200)
      LOGICAL MOVIE
      DATA MOVIE /.TRUE./
C Counters
      INTEGER GCOUNT1, GCOUNT2
      INTEGER GINCFRZ, GINCPRI, GINCEUT
      INTEGER ECOUNT, NGROW
C Tip & live diagnostic variables
      INTEGER TIPX, TIPY, OTIPY
      REAL VELMEAS, CURVKGT
      REAL OTIME, TTIP
      REAL TSTART
      REAL GTIME, DEUTE
      DATA OTIPY/0/
      DATA GTIME/0.0/
      DATA OTIME/0.25/
C Runtime control & closure parameters
      REAL DFRZ, DTSEED
      REAL DTAD
      REAL FIN, FIN2
      LOGICAL CASTOP
      DATA DFRZ/150.0/, DTSEED/1.0/
      DATA DTAD/10.0/
      DATA FIN/10.0/
C Transport & growth parameters
      REAL DCOEFF, DT, DX
      REAL GRAD, VEL
      REAL DMASK
      REAL DEFF, DALPHA
      DATA DCOEFF/1.0E-9/, DT/1.0E-4/, DX/1.0E-6/
      DATA GRAD/0.015/, VEL/200.0/
      DATA DMASK/0.60/
C Phase diagram & alloy system
      REAL TMP, MLIQ, KALPHA
      REAL C0, CEUT, TL0, TEUT
      REAL TLIQBIN
      EXTERNAL TLIQBIN
      DATA TMP/327.502/, MLIQ/-2.326/, KALPHA/0.296/
      DATA C0/8.0/, CEUT/61.9/
C Capillarity / anisotropy / orientation
      REAL THETA1, THETA2, FR, GKGT0, VV
      REAL GIBT, DD, SPARE
      DATA THETA1/1.04720/, THETA2/1.39626/
      DATA GIBT/2.0E-7/, DD/0.05/, SPARE/1.0/
C Sparse work arrays
      INTEGER IUPD(MAXUPD), JUPD(MAXUPD)
      INTEGER IREJ(MAXREJ), JREJ(MAXREJ)
      INTEGER ITCH(MAXTOUCH), JITCH(MAXTOUCH)
C Thermal terms
      REAL QLAT(NXMAX, NYMAX)
      REAL ATHERM, LHEAT, CPHEAT
      DATA ATHERM/2.5E-9/, LHEAT/1.0E4/, CPHEAT/5.0E2/    
C Imposed macro thermal-history table
      INTEGER MAXTHERM, NTHERM
      PARAMETER (MAXTHERM=20000)
      REAL TTAB(MAXTHERM), TREFTAB(MAXTHERM), GTAB(MAXTHERM)
      REAL RTAB(MAXTHERM), VTAB(MAXTHERM)
      CHARACTER*128 THFIL
      INTEGER JREFTH
      REAL TSHIFT
C -------------------------------------------------------------------
C Batch input file and I/O preference
C -------------------------------------------------------------------
      PARFIL = 'xparam.txt'
      GRIDFIL = 'xgrid.txt'
      THFIL = 'xthermal.csv'
      WRITE(*,*) 'Dendritic CA simulation'
      WRITE(*,*) 
      WRITE(*,*) 'Enter solidification parameter file:'
      READ(*,*) PARFIL
      WRITE(*,*)
      WRITE(*,*) 'Stop control:'
      WRITE(*,*) '1=time, 2=dendrite at top, 3=whichever occurs first'
      READ (*,*) ISTOPMODE
      WRITE(*,*)
      WRITE(*,*) 'Make movie?  1=yes, 0=no'
      READ (*,*) I
      MOVIE = .FALSE.
      IF (I.EQ.1) MOVIE = .TRUE.
C-------------------------------------------------------------------
C Initialize CA update screens
C-------------------------------------------------------------------
      SCRN = 1
      SNEW = 3 - SCRN
      SCUR = SCRN
C -------------------------------------------------------------------
C Open CSV output and write header
C--------------------------------------------------------------------
      OPEN(UNIT=IUNIT, FILE='data.csv', STATUS='UNKNOWN',
     &     FORM='FORMATTED')
      CALL CSVHEAD(IUNIT)
C -------------------------------------------------------------------
C Load parameters from CAparam.txt
C -------------------------------------------------------------------
      CALL LOADCA(PARFIL, PARFIL, GRIDFIL, THFIL,
     &     TMP, MLIQ, KALPHA, C0, CEUT, GIBT,
     &     DD, DCOEFF, DT, DX, FIN, GRAD, VEL, SPARE,
     &     THETA1, THETA2, DTSEED, DMASK, DTAD,
     &     XLEN, YLEN, ATHERM, LHEAT, IEMBED, IBCTYPE, CPHEAT)
      IF (XLEN .GT. NXMAX .OR. YLEN .GT. NYMAX) THEN
        WRITE(*,*) 'ERROR: Requested grid size exceeds limits!'
        STOP
      ENDIF
C -------------------------------------------------------------------
C Update derived quantities after file load
C -------------------------------------------------------------------
      NTHERM = 0
      JREFTH = 1
      TSHIFT = 0.0
      CALL LOADTHERMCSV(THFIL,TTAB,TREFTAB,GTAB,RTAB,VTAB,
     &     MAXTHERM,NTHERM)
      TL0 = TLIQBIN(C0,TMP,MLIQ,0.0)
      TEUT = TLIQBIN(CEUT,TMP,MLIQ,0.0)
      TSTART = TL0 - DTSEED
      VV = VEL*DX
      GKGT0 = GRAD/DX
      FR = TL0 - TEUT
      DEFF = DCOEFF/2.0
      DALPHA = DCOEFF*DT/(DX*DX)
      FIN2 = FIN    
C --------------------------------------------------------------------
C Stage A: load prepared starting GRID from xgrid.txt
C --------------------------------------------------------------------
      CALL GETGRID(GRID, XLEN, YLEN, GRIDFIL)
C The loaded grid may contain an old thermal field from setup.
C Reset TCELL to the imposed macro field at t=0 so that the initial
C perturbation is zero. Otherwise the first HEATPERT2D call would
C interpret the old TCELL field as a huge latent-heat perturbation.
      CALL INITTHERMCSV(GRID, XLEN, YLEN, SCRN, SNEW, DX,
     &   GRAD, VEL, TTAB, TREFTAB, GTAB, RTAB, VTAB, NTHERM,
     &   IEMBED, TL0, DTSEED, JREFTH, TSHIFT)
C --------------------------------------------------------------------
C Initialize counters and runtime variables
C --------------------------------------------------------------------
      CASTOP = .FALSE.
      GTIME = 0.0
      ECOUNT = 0
      GCOUNT1 = 0
      GCOUNT2 = 0
      NGROW = 0
      DONE = 0
C --------------------------------------------------------------------
C Initialize diagnostic counters and accumulators
C --------------------------------------------------------------------
      NTOK = 0
      NPRINTS = 0
      GINCFRZ = 0
      GINCPRI = 0
      GINCEUT = 0
      DEUTE = 0.0
C --------------------------------------------------------------------
C Initialize thermal analysis variables
C --------------------------------------------------------------------
      DO J=1,YLEN
      DO I=1,XLEN
         XRES(I,J) = 0.0
         XADD(I,J) = 0.0
         XRESADD(I,J) = 0.0
         XRND(I,J) = 0.0
         QLAT(I,J) = 0.0
      ENDDO
      ENDDO       
C --------------------------------------------------------------------
C Main CA simulation loop
C --------------------------------------------------------------------
      IF (.NOT. CASTOP) THEN
         CALL CALOOP(
C State arrays
     &    GRID, XADD, XRES, XRESADD, XRND, QLAT,     
C Domain size
     &    XLEN, YLEN,
C Screen state
     &    SCRN,
C Time and grid scales
     &    GTIME, FIN, DT, DX, TOL,
C Runtime control and output
     &    IUNIT, DSPNUM,
C Transport and thermal driving
     &    VEL, GRAD,
     &    DCOEFF, DMASK, DALPHA,
     &    ATHERM, LHEAT, CPHEAT,
     &    TTAB, TREFTAB, GTAB, RTAB, VTAB, NTHERM,
     &    JREFTH, TSHIFT,     
C Interface / capillarity / closure
     &    DD, GIBT, DFRZ, DTSEED,
C Phase diagram / alloy parameters
     &    TEUT, TL0, MLIQ, KALPHA, C0, CEUT,
C KGT and box-scale diagnostics
     &    CURVKGT,
     &    BOX2, TOPEDGE,
C Tip tracking / time history
     &    OTIPY, OTIME, TTIP,
C Counters and run statistics
     &    GCOUNT1, GCOUNT2, ECOUNT, NGROW,
     &    GINCFRZ, GINCPRI, GINCEUT,
     &    DEUTE, DONE,
C Measured outputs
     &    VELMEAS, TIPX, TIPY,
C Sparse work lists
     &    IUPD, JUPD, MAXUPD,
     &    IREJ, JREJ, MAXREJ,
     &    ITCH, JITCH, MAXTOUCH, IBCTYPE,
C Stop flag
     &    MOVIE, CASTOP, ISTOPMODE)
      ENDIF
C Save at least one ppm-image of the CA run.
      CALL MOVPPM(GRID, XLEN, YLEN, SCRN, C0, CEUT)
C --------------------------------------------------------------------
C Finalize run and write output files
C --------------------------------------------------------------------
      CLOSE(IUNIT)
      CALL DOTABLE(GRID, SCRN, XLEN, YLEN)
      WRITE(*,*)
      WRITE(*,*) '   --Finished--'
      WRITE(*,*) ' Parameters used: ', PARFIL
      WRITE(*,*) ' Starting grid:   ', GRIDFIL
      WRITE(*,*) ' Runtime data:    data.csv'
      WRITE(*,*) ' Final grid:      grid.txt'
      WRITE(*,*) ' Bitmap image:    out#.ppm'
      WRITE(*,*)
      IF (MOVIE) THEN
        WRITE(*,*) 'Possible command for movie scripting:'
        WRITE(*,*) 'ffmpeg -framerate 20 -i out%d.ppm movie.mp4'
      ENDIF
      END
C
C=====================================================================
C  CALOOP == Main CA time-stepping loop
C=====================================================================
      SUBROUTINE CALOOP(
     &    GRID, XADD, XRES, XRESADD, XRND, QLAT,
     &    XLEN, YLEN,
     &    SCRN,
     &    GTIME, FIN, DT, DX, TOL,
     &    IUNIT, DSPNUM,
     &    VEL, GRAD,
     &    DCOEFF, DMASK, DALPHA,
     &    ATHERM, LHEAT, CPHEAT,
     &    TTAB, TREFTAB, GTAB, RTAB, VTAB, NTHERM,
     &    JREFTH, TSHIFT,     
     &    DD, GIBT, DFRZ, DTSEED,
     &    TEUT, TL0, MLIQ, KALPHA, C0, CEUT,
     &    CURVKGT,
     &    BOX2, TOPEDGE,
     &    OTIPY, OTIME, TTIP,
     &    GCOUNT1, GCOUNT2, ECOUNT, NGROW,
     &    GINCFRZ, GINCPRI, GINCEUT,
     &    DEUTE, DONE,
     &    VELMEAS, TIPX, TIPY,
     &    IUPD, JUPD, MAXUPD,
     &    IREJ, JREJ, MAXREJ,
     &    ITCH, JITCH, MAXTOUCH, IBCTYPE,
     &    MOVIE, CASTOP, ISTOPMODE)
      IMPLICIT NONE
C Parameters
      INTEGER NXMAX, NYMAX, LENBOX, IBCTYPE
      PARAMETER (NXMAX=1000, NYMAX=1000)
      PARAMETER (LENBOX=7)
C Arrays
      REAL GRID(NXMAX,NYMAX,6,2)
      REAL XADD(NXMAX,NYMAX)
      REAL XRES(NXMAX,NYMAX)
      REAL XRESADD(NXMAX,NYMAX)
      REAL XRND(NXMAX,NYMAX)
      REAL QLAT(NXMAX,NYMAX)      
C Scalars
      INTEGER XLEN,YLEN
      INTEGER SCRN,IUNIT,DSPNUM
      REAL GTIME,FIN,DT,DX,TOL
      REAL VEL,GRAD,DCOEFF,DMASK,DALPHA
      REAL ATHERM, LHEAT, CPHEAT
      INTEGER NTHERM, JREFTH
      REAL TTAB(*), TREFTAB(*), GTAB(*), RTAB(*), VTAB(*)
      REAL TSHIFT
      REAL DD,GIBT,DFRZ,DTSEED
      REAL TEUT,TL0,MLIQ,KALPHA,C0, CEUT
      REAL CURVKGT
      INTEGER BOX2,TOPEDGE
      INTEGER OTIPY
      REAL OTIME,TTIP
      INTEGER GCOUNT1,GCOUNT2,ECOUNT,NGROW
      INTEGER GINCFRZ,GINCPRI,GINCEUT
      REAL DEUTE
      INTEGER DONE, ISTOPMODE
      REAL VELMEAS
      INTEGER TIPX,TIPY
      INTEGER MAXUPD,MAXREJ,MAXTOUCH
      INTEGER IUPD(MAXUPD),JUPD(MAXUPD)
      INTEGER IREJ(MAXREJ),JREJ(MAXREJ)
      INTEGER ITCH(MAXTOUCH),JITCH(MAXTOUCH)
      LOGICAL MOVIE, CASTOP
C Locals
      INTEGER I,J,NREJ,NTOUCH,NSTAGEBOX
      INTEGER OTIPX,ITIPLOAD
      INTEGER ISTAGEBOX(64),JSTAGEBOX(64),IFRZSTAGE(64)
      LOGICAL HAVETIP,VINIT,EINIT,CINIT,UINIT
      DOUBLE PRECISION MASSBAL,M0
      REAL OEMATIME,CTAVG,UNCLAVG
      REAL DTBEST,TLBEST
      REAL VELAVG,VKGTMS,GGKGTUSE
      REAL GX,GY
      REAL GTXTIP,GTYTIP,GGTIP
      REAL DT2KGT,VV,GKGT,RTIP,PE,IVR,DELT0
      REAL DFRZKGT,DFRZMS,PHIKGT,RHOTIP
      REAL VDTRAP,KEQ,KPART
      REAL TLX,TLY,TLGX,TLGY,TLCLTIP,TLLDIF,TLWID,TLLEN
      REAL TOTDEUTE,MAXFRATE,MTDIFF
      INTEGER TOTFRZ,TOTPRI,TOTEUT
C External
      REAL KPARTBIN
      EXTERNAL MASSBAL,KPARTBIN
C Initialization
      M0 = MASSBAL(GRID, XADD,XRES,XRND,XLEN,YLEN,SCRN,CASTOP)
      DO J=1,YLEN
      DO I=1,XLEN
         XADD(I,J)=0.0
         XRESADD(I,J)=0.0
      ENDDO
      ENDDO
      TIPX=1
      TIPY=1
      OTIPX=TIPX
      OTIPY=TIPY
      OTIME=GTIME
      TTIP=TL0-DTSEED
      VELMEAS=0.0
      VELAVG=0.0
      VKGTMS=0.0
      GGKGTUSE=0.0
      GX=0.0
      GY=1.0
      VINIT=.FALSE.
      HAVETIP=.FALSE.
      OEMATIME=0.0
      CTAVG=0.0
      UNCLAVG=0.0
      EINIT=.FALSE.
      CINIT=.FALSE.
      UINIT=.FALSE.
      DTBEST=1.0E30
      TLBEST=-1.0E30
      DT2KGT=DTSEED
      VV=0.0
      GKGT=0.0
      RTIP=0.0
      PE=0.0
      IVR=0.0
      DELT0=0.0
      DFRZKGT=0.0
      DFRZMS=DFRZ
      PHIKGT=1.0
      RHOTIP=0.0
      VDTRAP=1.0E-2
      KEQ=KALPHA
      TLX=0.0
      TLY=0.0
      TLGX=0.0
      TLGY=1.0
      GTXTIP=0.0
      GTYTIP=1.0
      GGTIP=ABS(GRAD)/DX
      TLCLTIP=C0
      TLLDIF=0.0
      TLWID=0.0
      TLLEN=0.0
      ITIPLOAD=0
      DEUTE=0.0
      MTDIFF=0.0
      IF (DALPHA > 0.25) THEN
         WRITE(*,*) 'Error: DALPHA > 0.25; stability issues!'
         CASTOP = .TRUE.
      ENDIF
      CALL FNDTDR(GRID,XLEN,YLEN,SCRN,GX,GY,TIPX,TIPY,HAVETIP)
      IF (.NOT. HAVETIP) THEN
         TIPX = XLEN/2
         TIPY = 1
      ENDIF
      CALL DIAGINIT(TOTFRZ,TOTPRI,TOTEUT,TOTDEUTE,MAXFRATE)
C=====================================================================
C Runtime loop
C=====================================================================
      IF (.NOT. CASTOP) THEN
      DO WHILE (.NOT. CASTOP)
C Runtime advance
C -------------------------
         CALL STEPINC(GTIME,DT,ECOUNT,GCOUNT1,GCOUNT2,NREJ)
C -------------------------
C Block A: transport
C   Energy transport is updated first.  
C -------------------------
         CALL HEATPERT2D(GRID,QLAT,XLEN,YLEN,SCRN,
     &      DT,DX,ATHERM,GTIME,VEL,GRAD,
     &      TTAB,TREFTAB,GTAB,RTAB,VTAB,NTHERM,
     &      JREFTH,TSHIFT)
C   Solute transport follows the thermal update.
         CALL DIFFUSE(GRID,XLEN,YLEN,SCRN,DALPHA,DMASK, IBCTYPE)
C ----------------------------------
C Block B: runtime interface state
C ----------------------------------
         CALL XRUNSTATE(VEL,GRAD,DX,DCOEFF,DT,
     &      GIBT,MLIQ,KALPHA,C0,
     &      DT2KGT,VV,GKGT,RTIP,CURVKGT,PE,IVR,DELT0,
     &      DFRZKGT,DFRZMS,PHIKGT,RHOTIP,DFRZ)
         KEQ = KPARTBIN(TLCLTIP,TTIP,TL0,MLIQ,KALPHA,C0)
         CALL XKPARTRUN(VV,KEQ,VDTRAP,KPART)
C Select scan-band mode from the current geometric tip direction.
C FNDTDR locates the active tip; the tip's stored CELLDIR selects
C vertical, horizontal, or rectangular banding for XSTAGE.
         CALL DTINFC(GRID,XLEN,YLEN,SCRN)
C -----------------------------------------------------------
C Block C: dendrite growth.  Note the two separate variables:
C   - DFRZKGT controls how much solid is staged at the tip
C   - DFRZ controls when remaining liquid is flushed outward
C-------------------------------------------------------------
         CALL XTIPBLD(GRID,XLEN,YLEN,SCRN,BOX2,
     &      TIPX,TIPY,TTIP,
     &      MLIQ,TL0,C0, IBCTYPE,
     &      DFRZKGT,64,NSTAGEBOX,ISTAGEBOX,JSTAGEBOX,
     &      IFRZSTAGE,3)    
C --
C Current leading-tip thermal-gradient state for local KGT use.
C This is kept separate from the diagnostic XDIAGTIP cadence.
         CALL TIPGRAD(GRID,XLEN,YLEN,SCRN,TIPX,TIPY,DX,
     &      GTXTIP,GTYTIP,GGTIP)
C --         
         CALL XGROW(GRID,XADD,XRES,XRND,QLAT,XLEN,YLEN,
     &    SCRN,BOX2,TOL,TEUT,TL0,MLIQ,KPART,C0,
     &    GIBT,DD,DFRZ,CURVKGT,LHEAT,CPHEAT,
     &    VEL,GRAD,DX,DCOEFF,DT,
     &    GINCFRZ,GINCPRI,GINCEUT,
     &    DTBEST,TLBEST,TIPX,TIPY,HAVETIP,TTIP,
     &    TLX,TLY,TLGX,TLGY,
     &    GTXTIP,GTYTIP,GGTIP,
     &    TLCLTIP,TLLDIF,TLWID,TLLEN,
     &    ITIPLOAD,
     &    NSTAGEBOX,ISTAGEBOX,JSTAGEBOX,IFRZSTAGE,
     &    IUPD,JUPD,MAXUPD,NGROW,LENBOX, IBCTYPE, CASTOP) 
C --------------------------------
C Block D0: staged solute transfer and mass balance
C --------------------------------
         CALL CATRANSFER(GRID,XLEN,YLEN,SCRN,
     &     XADD,XRES,XRESADD,XRND,
     &     NREJ,IREJ,JREJ,MAXREJ,
     &     NTOUCH,ITCH,JITCH,MAXTOUCH,
     &     MLIQ,TL0,C0,
     &     GCOUNT1,DSPNUM,M0,MTDIFF, IBCTYPE, CASTOP)
C --------------------------------
C Block D: diagnostics and output
C Diagnostic KGT state only
C --------------------------------
         IF (MOD(GCOUNT1,DSPNUM).EQ.0) THEN
            CALL XDIAGTIP(GRID,XRES,XLEN,YLEN,SCRN,
     &         GTIME,DT,DX,DSPNUM,
     &         IUNIT,
     &         TEUT,TL0,DFRZ,
     &         OTIPX,OTIPY,OTIME,OEMATIME,EINIT,TTIP,
     &         VELAVG,VINIT,VKGTMS,GGKGTUSE,
     &         TIPX,TIPY,DONE,
     &         CTAVG,UNCLAVG,CINIT,UINIT, IBCTYPE,
     &         C0,MLIQ,KALPHA,DD,GIBT,DCOEFF,GX,GY)
            CALL XUPDATE(GCOUNT1,GINCFRZ,GINCPRI,GINCEUT,
     &         DEUTE,XLEN,YLEN,MTDIFF,
     &         TOTFRZ,TOTPRI,TOTEUT,TOTDEUTE,MAXFRATE)
            IF (MOVIE) CALL MOVPPM(GRID, XLEN, YLEN, SCRN, C0, CEUT)
         ENDIF
C Loop tail
         CALL XRESET(GINCFRZ,GINCPRI,GINCEUT,DEUTE)
C Stop control:
C   ISTOPMODE = 1 : stop at FIN only
C   ISTOPMODE = 2 : stop when dendrite reaches TOPEDGE
C   ISTOPMODE = 3 : stop at FIN or TOPEDGE, whichever comes first
         IF (ISTOPMODE .EQ. 1) THEN
            IF (GTIME .GE. FIN) CASTOP = .TRUE.
         ELSE IF (ISTOPMODE .EQ. 2) THEN
            CALL XFINISH(TIPY,YLEN,TOPEDGE,GTIME,FIN)
         ELSE
            IF (GTIME .GE. FIN) CASTOP = .TRUE.
            CALL XFINISH(TIPY,YLEN,TOPEDGE,GTIME,FIN)
         ENDIF
         IF (CASTOP) GTIME = FIN
      ENDDO
      ENDIF
      RETURN
      END
C
C=====================================================================
C Binary primary-solidification thermodynamic wrappers
C
C Default model: linear liquidus with constant partition coefficient.
C The CA growth logic should call these routines rather than embedding
C liquidus/solidus equations directly. 
C=====================================================================
C Liquidus temperature
      REAL FUNCTION TLIQBIN(CL,TREF,MLIQ,CREF)
      IMPLICIT NONE
      REAL CL,TREF,MLIQ,CREF
      TLIQBIN = TREF + MLIQ*(CL - CREF)
      RETURN
      END
C Liquid composition
      REAL FUNCTION CLIQBIN(T,TREF,MLIQ,CREF)
      IMPLICIT NONE
      REAL T,TREF,MLIQ,CREF
      CLIQBIN = CREF
      IF (ABS(MLIQ) .GT. 1.0E-12) THEN
         CLIQBIN = CREF + (T - TREF)/MLIQ
      ENDIF
      RETURN
      END
C Partition coefficient
      REAL FUNCTION KPARTBIN(CL,T,TREF,MLIQ,KALPHA,CREF)
      IMPLICIT NONE
      REAL CL,T,TREF,MLIQ,KALPHA,CREF
C Current equilibrium model is constant k.  The extra arguments are
C retained for later T- or C-dependent thermodynamic partitioning.
      KPARTBIN = KALPHA + 0.0*(CL+T+TREF+MLIQ+CREF)
      RETURN
      END
C Solid composition
      REAL FUNCTION CSOLBIN(CL,T,TREF,MLIQ,KALPHA,CREF)
      IMPLICIT NONE
      REAL CL,T,TREF,MLIQ,KALPHA,CREF
      REAL KPARTBIN
      EXTERNAL KPARTBIN
      CSOLBIN = KPARTBIN(CL,T,TREF,MLIQ,KALPHA,CREF)*CL
      RETURN
      END
C
C ===============================================================
C LOADCA -- read simulation starting data as named value pairs
C Numeric entries: NAME  VALUE
C String/file entries:
C      PARFIL   filename
C      GRIDFIL  filename
C      THFIL    filename
C ===============================================================
      SUBROUTINE LOADCA(FNAME, PARFIL, GRIDFIL, THFIL,
     &     TMP, MLIQ, KALPHA, C0, CEUT, GIBT,
     &     DD, DCOEFF, DT, DX, FIN, GRAD, VEL, SPARE,
     &     THETA1, THETA2, DTSEED, DMASK, DTAD, XLEN, YLEN,
     &     ATHERM, LHEAT, IEMBED, IBCTYPE, CPHEAT)
      IMPLICIT NONE
      CHARACTER*(*) FNAME, PARFIL, GRIDFIL, THFIL
      CHARACTER*128 LINE, NAME, STRVAL
      INTEGER IOS, JSTAT, UNIT, XLEN, YLEN, IFOUND, NFOUND
      INTEGER IEMBED, IBCTYPE
      REAL TMP, MLIQ, KALPHA, C0, CEUT, GIBT
      REAL DD, DCOEFF, DT, DX, FIN, GRAD, VEL
      REAL THETA1, THETA2, VALUE, DTSEED, SPARE
      REAL DMASK, DTAD, ATHERM, LHEAT, CPHEAT
      UNIT   = 20
      NFOUND = 0
      OPEN(UNIT=UNIT, FILE=FNAME, STATUS='OLD', IOSTAT=IOS)
      IF (IOS .EQ. 0) THEN
         READ(UNIT,'(A)',IOSTAT=IOS) LINE
         DO WHILE (IOS .EQ. 0)
          IF (LINE .NE. ' ') THEN
           IF (LINE(1:1) .NE. '#') THEN
              IFOUND = 0
C First try numeric named-value input.
              READ(LINE,*,IOSTAT=JSTAT) NAME, VALUE
              IF (JSTAT .EQ. 0) THEN
                CALL SETPARNUM(NAME, VALUE,
     &               TMP, MLIQ, KALPHA, C0, CEUT, GIBT,
     &               DD, DCOEFF, DT, DX, FIN,
     &               GRAD, VEL, THETA1, THETA2, DTSEED, DMASK,
     &               DTAD, XLEN, YLEN, SPARE,
     &               ATHERM, LHEAT, CPHEAT,
     &               IEMBED, IBCTYPE, IFOUND)
             ELSE
C If numeric read fails, try string/file input.
                READ(LINE,*,IOSTAT=JSTAT) NAME, STRVAL
                IF (JSTAT .EQ. 0) THEN
                   CALL SETPARSTR(NAME, STRVAL,
     &                  PARFIL, GRIDFIL, THFIL, IFOUND)
                ENDIF
              ENDIF
              IF (IFOUND .EQ. 0) THEN
                WRITE(*,*) 'LOADCA: unknown name ignored: ',NAME
              ELSE
                NFOUND = NFOUND + 1
              ENDIF
           ENDIF
          ENDIF
          READ(UNIT,'(A)',IOSTAT=IOS) LINE
         END DO
         CLOSE(UNIT)
         WRITE(*,*) 'LOADCA: loaded ', NFOUND,
     &              ' parameters from ', FNAME
      ELSE
         WRITE(*,*) 'LOADCA: could not open file ', FNAME
      ENDIF
      RETURN
      END
C
C=====================================================================
C SETPARNUM -- assign one numeric named parameter.
C=====================================================================
      SUBROUTINE SETPARNUM(NAME, VALUE,
     &     TMP, MLIQ, KALPHA, C0, CEUT, GIBT,
     &     DD, DCOEFF, DT, DX, FIN, GRAD, VEL,
     &     THETA1, THETA2, DTSEED, DMASK,
     &     DTAD, XLEN, YLEN, SPARE,
     &     ATHERM, LHEAT, CPHEAT, IEMBED, IBCTYPE, IFOUND)
      IMPLICIT NONE
      CHARACTER*(*) NAME
      INTEGER XLEN, YLEN, IEMBED, IBCTYPE, IFOUND
      REAL VALUE
      REAL TMP, MLIQ, KALPHA, C0, CEUT, GIBT
      REAL DD, DCOEFF, DT, DX, FIN, GRAD, VEL
      REAL THETA1, THETA2, DTSEED, DMASK, DTAD, SPARE
      REAL ATHERM, LHEAT, CPHEAT
      IFOUND = 1
      IF (NAME .EQ. 'TMP') THEN
         TMP = VALUE
      ELSE IF (NAME .EQ. 'MLIQ') THEN
         MLIQ = VALUE
      ELSE IF (NAME .EQ. 'KALPHA') THEN
         KALPHA = VALUE
      ELSE IF (NAME .EQ. 'C0') THEN
         C0 = VALUE
      ELSE IF (NAME .EQ. 'CEUT') THEN
         CEUT = VALUE
      ELSE IF (NAME .EQ. 'GIBT') THEN
         GIBT = VALUE
      ELSE IF (NAME .EQ. 'DD') THEN
         DD = VALUE
      ELSE IF (NAME .EQ. 'DCOEFF') THEN
         DCOEFF = VALUE
      ELSE IF (NAME .EQ. 'DT') THEN
         DT = VALUE
      ELSE IF (NAME .EQ. 'DX') THEN
         DX = VALUE
      ELSE IF (NAME .EQ. 'FIN') THEN
         FIN = VALUE
      ELSE IF (NAME .EQ. 'GRAD') THEN
         GRAD = VALUE
      ELSE IF (NAME .EQ. 'VEL') THEN
         VEL = VALUE
      ELSE IF (NAME .EQ. 'THETA1') THEN
         THETA1 = VALUE
      ELSE IF (NAME .EQ. 'THETA2') THEN
         THETA2 = VALUE
      ELSE IF (NAME .EQ. 'DTSEED') THEN
         DTSEED = VALUE
      ELSE IF (NAME .EQ. 'DMASK') THEN
         DMASK = VALUE
      ELSE IF (NAME .EQ. 'DTAD') THEN
         DTAD = VALUE
      ELSE IF (NAME .EQ. 'XLEN') THEN
         XLEN = INT(VALUE)
      ELSE IF (NAME .EQ. 'YLEN') THEN
         YLEN = INT(VALUE)
      ELSE IF (NAME .EQ. 'SPARE') THEN
         SPARE = VALUE
      ELSE IF (NAME .EQ. 'ATHERM') THEN
         ATHERM = VALUE
      ELSE IF (NAME .EQ. 'LHEAT') THEN
         LHEAT = VALUE
      ELSE IF (NAME .EQ. 'CPHEAT') THEN
         CPHEAT = VALUE
      ELSE IF (NAME .EQ. 'IEMBED') THEN
         IEMBED = INT(VALUE)
      ELSE IF (NAME .EQ. 'IBCTYPE') THEN
         IBCTYPE = INT(VALUE)
      ELSE
         IFOUND = 0
      ENDIF
      RETURN
      END
C
C=====================================================================
C SETPARSTR -- assign one string/file named parameter.
C=====================================================================
      SUBROUTINE SETPARSTR(NAME,STRVAL,PARFIL, GRIDFIL, THFIL, IFOUND)
      IMPLICIT NONE
      CHARACTER*(*) NAME, STRVAL
      CHARACTER*(*) PARFIL, GRIDFIL, THFIL
      INTEGER IFOUND
      IFOUND = 1
      IF (NAME .EQ. 'PARFIL') THEN
         PARFIL = STRVAL
      ELSE IF (NAME .EQ. 'GRIDFIL') THEN
         GRIDFIL = STRVAL
      ELSE IF (NAME .EQ. 'THFIL') THEN
         THFIL = STRVAL
      ELSE
         IFOUND = 0
      ENDIF
      RETURN
      END
C   
C=====================================================================
C CSVHEAD -- write CSV header
C=====================================================================
      SUBROUTINE CSVHEAD(IUNIT)
      IMPLICIT NONE
      INTEGER IUNIT
      WRITE(IUNIT,1000)
1000  FORMAT('time(s),vel(um/s),tipX,tipY,tipT(degC),',
     &     'tipCavg(wt%),undercoolavg(degC),',
     &     'radiusKGT(um),curvKGT(1/m),theta(rad),DFRZ,',
     &     'DTKGT(degC),peclet,Clstar(wt%),OmegaCA,OmegaKGT')
      RETURN
      END
C
C ====================================================================
C GETGRID -- read the CA grid from a text file
C            and load into BOTH screens
C ====================================================================
      SUBROUTINE GETGRID(GRID, XLEN, YLEN, FILNAM)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN
      REAL GRID(NXMAX, NYMAX, 6, 2)
      CHARACTER FILNAM*128
      INTEGER UNIT, IOS, I, J, K, IX, IY
      INTEGER LCELLS, LIQUID, SOLID, PHASE, TCELL, CELLDIR
      REAL LQCELLS, LCOMP, SCOMP, PHASEV, TEMP, ANGLE
      PARAMETER (LCELLS=1, LIQUID=2, SOLID=3, PHASE=4, TCELL=5)
      PARAMETER (CELLDIR=6)
      CHARACTER LINE*120
      LOGICAL ERROR
      UNIT = 20
      ERROR = .FALSE.
      OPEN (UNIT=UNIT, FILE=FILNAM, STATUS='OLD', IOSTAT=IOS)
      IF (IOS .EQ. 0) THEN
C Skip the header line
        READ (UNIT, '(A)', IOSTAT=IOS) LINE
        IF (IOS .NE. 0) THEN
           WRITE (*,*) ' *** ERROR reading header in ', FILNAM
           ERROR = .TRUE.
        ENDIF
C Read the data into BOTH screens (K=1,2)
        J = 1
        DO WHILE (J .LE. YLEN .AND. .NOT. ERROR)
         I = 1
         DO WHILE (I .LE. XLEN .AND. .NOT. ERROR)
           READ (UNIT, *, IOSTAT=IOS) IX, IY,
     &           LQCELLS, LCOMP, SCOMP, PHASEV, TEMP, ANGLE
           IF (IOS .NE. 0) THEN
             WRITE (*,*) ' *** ERROR reading data at I,J =',I,J
             ERROR = .TRUE.
           ELSE
             IF (IX .NE. I .OR. IY .NE. J) THEN
                WRITE(*,*) ' *** GETGRID order mismatch: I,J=',I,J,
     &                     ' IX,IY=',IX,IY
             ENDIF
C Success: Fill both time-steps (SCUR/SNEW)
              DO K = 1, 2
               GRID(I,J,LCELLS,K) = LQCELLS
               GRID(I,J,LIQUID,K) = LCOMP
               GRID(I,J,SOLID,K) = SCOMP
               GRID(I,J,PHASE,K) = PHASEV
               GRID(I,J,TCELL,K) = TEMP
               GRID(I,J,CELLDIR,K) = ANGLE
              END DO
           ENDIF
           I = I + 1
         END DO
         J = J + 1
        END DO
        IF (.NOT. ERROR) THEN
           WRITE(*,*) 'Successfully loaded grid from ', FILNAM
        ENDIF
        CLOSE (UNIT)
      ELSE
         WRITE (*,*) ' *** ERROR opening file ', FILNAM
      ENDIF
      RETURN
      END
C
C ====================================================================
C DOTABLE - Save the full GRID as a text file (grid.txt)
C ====================================================================
      SUBROUTINE DOTABLE (GRID,SCRNUM,XLEN,YLEN)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER SCRNUM, XLEN, YLEN
      REAL GRID(NXMAX, NYMAX, 6, 2)
      INTEGER LCELLS, LIQUID, SOLID, PHASE, TCELL, CELLDIR
      INTEGER I,J,UNIT
      CHARACTER*20 FILENAME
      PARAMETER (LCELLS=1, LIQUID=2, SOLID=3, PHASE=4, TCELL=5)
      PARAMETER (CELLDIR=6)
      DATA FILENAME /'grid.txt'/
      UNIT = 10
      OPEN(UNIT, FILE=FILENAME, STATUS='UNKNOWN')
      WRITE(UNIT,'(A)') ' X Y #liqCells Lcomp Scomp Phase Temp Orient'
      DO J = 1, YLEN
         DO I = 1, XLEN
            WRITE(UNIT, '(2I4,3F10.3,2F10.2,F10.4)') I, J,
     &           GRID(I, J, LCELLS, SCRNUM),  
     &           GRID(I, J, LIQUID, SCRNUM),  
     &           GRID(I, J, SOLID, SCRNUM),  
     &           GRID(I, J, PHASE, SCRNUM),  
     &           GRID(I, J, TCELL, SCRNUM),
     &           GRID(I, J, CELLDIR, SCRNUM)   
         END DO
      END DO
      CLOSE(UNIT)
      PRINT *, 'Full grid data written to ', FILENAME
      RETURN
      END
C
C=====================================================================
C DIAGINIT -- Initialize diagnostic accumulators (caller-owned state)
C=====================================================================
      SUBROUTINE DIAGINIT(TOTFRZ, TOTPRI, TOTEUT, TOTDEUTE, MAXFRATE)
      IMPLICIT NONE
      INTEGER TOTFRZ, TOTPRI, TOTEUT
      REAL    TOTDEUTE, MAXFRATE
      TOTFRZ = 0
      TOTPRI = 0
      TOTEUT = 0
      TOTDEUTE = 0.0
      MAXFRATE = 0.0
      RETURN
      END
C
C=====================================================================
C XRESET -- reset per-step counters
C=====================================================================
      SUBROUTINE XRESET(GINCFRZ,GINCPRI,GINCEUT,DEUTE)
      IMPLICIT NONE
      INTEGER GINCFRZ, GINCPRI, GINCEUT
      REAL DEUTE
      GINCFRZ = 0
      GINCPRI = 0
      GINCEUT = 0
      DEUTE = 0.0
      RETURN
      END
C
C=====================================================================
C XFINISH -- finish check when the tip reaches the top buffer
C=====================================================================
      SUBROUTINE XFINISH(TIPY,YLEN,TOPEDGE,GTIME,FIN)
      IMPLICIT NONE
      INTEGER TIPY, YLEN, TOPEDGE
      REAL GTIME, FIN
      IF (TIPY .GT. (YLEN-TOPEDGE)) FIN = GTIME
      RETURN
      END
C
C=====================================================================
C STEPINC -- advance time and per-step counters
C=====================================================================
      SUBROUTINE STEPINC(GTIME,DT,ECOUNT,GCOUNT1,GCOUNT2,NREJ)
      IMPLICIT NONE
      REAL GTIME, DT
      INTEGER ECOUNT, GCOUNT1, GCOUNT2, NREJ
      GTIME = GTIME + DT
      ECOUNT = ECOUNT + 1
      GCOUNT1 = GCOUNT1 + 1
      GCOUNT2 = GCOUNT2 + 1
      NREJ = 0
      RETURN
      END
C
C=====================================================================
C XRUNSTATE -- build runtime interface state for the current step
C   Runtime:
C     - imposed VEL, GRAD
C     - runtime KGT state
C     - runtime DFRZ
C   This routine does NOT use measured tip velocity or diagnostics.
C=====================================================================
      SUBROUTINE XRUNSTATE(VEL,GRAD,DX,DCOEFF,DT,
     &    GIBT,MLIQ,KALPHA,C0,
     &    DT2KGT,VV,GKGT,RTIP,CURVKGT,PE,IVR,DELT0,
     &    DFRZKGT,DFRZMS,PHIKGT,RHOTIP,DFRZ)
      IMPLICIT NONE
      REAL VEL,GRAD,DX,DCOEFF,DT,ZUNITS
      REAL GIBT,MLIQ,KALPHA,C0
      REAL DT2KGT,VV,GKGT,RTIP,CURVKGT,PE,IVR,DELT0
      REAL DFRZKGT,DFRZMS,PHIKGT,RHOTIP,DFRZ
      REAL CURV2KGT, SOLNEW, AA0
      REAL LDIFF, LTRANS, WSTART
      PARAMETER (ZUNITS=1000.0)
C Constant setting the strength of the runtime DFRZ update.
      DATA AA0 /0.12/  
C KGT analysis
      VV = VEL*DX
      GKGT = GRAD/DX
      CALL KGT(VV,GKGT,GIBT,DCOEFF,MLIQ,KALPHA,C0,
     &         DT2KGT,RTIP,CURV2KGT,PE,IVR,DELT0)
      CURVKGT = CURV2KGT
      IF (VV.GT.0.0 .AND. RTIP.GT.0.0 .AND.
     &    DT2KGT.GT.0.0 .AND. PE.GT.0.0) THEN
         CALL XKGTFRZCAP(VV,DT,DX,RTIP,DFRZKGT,DFRZMS,PHIKGT,RHOTIP)     
C Startup weighting based on diffusion length
C diffusion length with a few cells threshold (LTRANS)
      LDIFF = DCOEFF / MAX(VV,1.0E-12)     
      LTRANS = 3.0*DX                      
      WSTART = LDIFF/(LDIFF + LTRANS)
      IF (WSTART .LT. 0.0) WSTART = 0.0
      IF (WSTART .GT. 1.0) WSTART = 1.0     
C Runtime DFRZ threshold for final capture/ejection.
C This is different than DFRZKGT, which is the staged tip-freezing
C amount from the KGT dendrite cap geometry.
         IF (RTIP .GT. 0.0) THEN
          SOLNEW = ZUNITS*AA0*DX/RTIP  
          DFRZ = WSTART*SOLNEW + (1.0 - WSTART)*50.0
          IF (DFRZ .LT. 10.0) DFRZ = 10.0
          IF (DFRZ .GT. 0.8*ZUNITS) DFRZ = 0.8*ZUNITS
         ENDIF
      ENDIF
      RETURN
      END
C
C=====================================================================
C XKPARTRUN -- runtime Aziz partition coefficient.  VV [m/s]
C=====================================================================
      SUBROUTINE XKPARTRUN(VV,KALPHA,VDTRAP,KPART)
      IMPLICIT NONE
      REAL VV,KALPHA,VDTRAP,KPART
      REAL KVEFF
      EXTERNAL KVEFF
      KPART = KVEFF(VV,KALPHA,VDTRAP)
      RETURN
      END
C
C=====================================================================
C XTIPBLD -- build the local staged-freezing allocation around
C the current leading tip.
C
C Role of the model
C -----------------
C This is a key hybrid-interface routine.  It converts the analytical
C tip-scale freezing amount DFRZKGT into a set of integer local
C freezing increments on nearby CA cells.
C
C In this sense, XTIPBLD replaces a purely grid-rule interface
C advance with a discrete, analytically guided allocation step:
C  - KGT provides the target local freezing amount at the active tip
C  - nearby eligible CA cells are identified around (TIPX,TIPY)
C  - a weighted share of that freezing budget is assigned to them
C  - the output is an integer staged-freezing list for later execution
C
C This routine does not update GRID directly.  It only builds the
C sparse staged list:
C
C   ISTG, JSTG   : cell locations to be staged
C   IFRZSTG      : integer subcell freezing increments for each cell
C   NSTG         : number of active staged cells
C
C Selection logic
C ---------------
C A candidate cell must:
C   - lie in the local 5x5 box around the leading tip
C   - still contain liquid/subcells (S1 > 0)
C   - have at least one neighboring solid cell (NNN test)
C   - be locally undercooled relative to the liquidus
C
C Weighting logic
C ---------------
C Eligible cells are weighted using:
C   - local orientation alignment from BOX versus stored crystal angle
C   - a thermal backoff factor based on TTIP - TL
C   - a minimum floor FMIN so off-axis cells can still participate
C
C Allocation logic
C ----------------
C The raw weights are normalized to match the requested freezing target
C TARG = DFRZKGT. The resulting real-valued freezing amounts are then:
C
C   - clipped by the local liquid capacity of each cell
C   - converted to integer subcell increments
C   - corrected by largest-fraction remainder so the final integer sum
C     tracks the target as closely as possible
C
C Interpretation
C --------------
C XTIPBLD is the local analytical-to-discrete handoff step for
C interface advance.  
C=====================================================================
      SUBROUTINE XTIPBLD(GRID,XLEN,YLEN,SCRN,BOX2,
     &    TIPX,TIPY,TTIP,
     &    MLIQ,TL0,C0, IBCTYPE,
     &    DFRZKGT,MAXSTG,NSTG,ISTG,JSTG,IFRZSTG, LENBOX)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN, SCRN, BOX2, MAXSTG
      INTEGER LCELLS, LIQUID, TCELL, CELLDIR
      INTEGER TIPX, TIPY, NSTG, IBCTYPE
      INTEGER ISTG(MAXSTG), JSTG(MAXSTG), IFRZSTG(MAXSTG)
      REAL GRID(NXMAX,NYMAX,6,2)
      REAL TTIP
      REAL MLIQ, TL0, C0
      REAL TLIQBIN
      PARAMETER (LCELLS=1, LIQUID=2)
      PARAMETER (TCELL=5, CELLDIR=6)
      LOGICAL DONE
      EXTERNAL TLIQBIN
      REAL DFRZKGT
C Working arrays for the sparse staged set:
C   WSTG    : raw candidate weight
C   DFRZSTG : real-valued freezing amount before integer rounding
C   CAPSTG  : local liquid/subcell capacity available to freeze
C   FRCSTG  : fractional remainder used in final cleanup
      REAL WSTG(64), DFRZSTG(64), CAPSTG(64), FRCSTG(64)
      INTEGER I, J, I1, I2, J1, J2, NC, IP, JP, K, KBEST
      INTEGER NTARG, NSUM, NLEFT
      REAL S1, TL, COMP, LNN, THETA, PTHETA
      REAL PI, HALFPI, FMIN, DT0, DT1
      REAL ANGDIFF, CANG, WANG, WTHERM, DTBACK, FRAC
      REAL WRAW, WSUM, SCALE, TARG, BESTF
      INTEGER LENBOX
      PARAMETER (PI=3.1415927)
      PARAMETER (HALFPI=1.5707963)
C Minimum participation weight for eligible but poorly aligned cells.
      PARAMETER (FMIN=0.15)
C Thermal backoff window:
C   DTBACK <= DT0  -> full thermal weight
C   DTBACK >= DT1  -> zero thermal weight
C   DT0 < DTBACK < DT1 -> linear ramp down
      PARAMETER (DT0=0.04)
      PARAMETER (DT1=0.16)
C Initialize empty staged list and requested freeze target.
      NSTG = 0
      WSUM = 0.0
C Target local freezing budget from analytical tip-scale closure.
      TARG = DFRZKGT
C Basic guards.
      IF (TIPX .LT. 1 .OR. TIPX .GT. XLEN) RETURN
      IF (TIPY .LT. 1 .OR. TIPY .GT. YLEN) RETURN
      IF (TARG .LE. 0.0) RETURN
C Search only in a local 5x5 box around the current tip.
C Y limits also respect the thermal buffer BOX2.
      I1 = MAX(1, TIPX-2)
      I2 = MIN(XLEN, TIPX+2)
      J1 = MAX(BOX2+1, TIPY-2)
      J2 = MIN(YLEN-BOX2-1, TIPY+2)
      DO J = J1, J2
      DO I = I1, I2
C Candidate must still contain liquid/subcells.
         S1 = GRID(I,J,LCELLS,SCRN)
         IF (S1 .GT. 1.0E-6) THEN
C Candidate must lie on/near the active interface:
C NNN returns the count of solid neighbors and nearest solid direction.
          CALL NNN(GRID,I,J,XLEN,YLEN,SCRN,NC,LNN,IP,JP, IBCTYPE)
          IF (NC .GT. 0) THEN
             TL   = GRID(I,J,TCELL,SCRN)
             COMP = GRID(I,J,LIQUID,SCRN)
C Candidate must be locally undercooled relative to the liquidus.
          IF (TLIQBIN(COMP,TL0,MLIQ,C0) .GT. TL) THEN
C Estimate local interface direction from BOX and compare with the
C stored crystal direction.  Better alignment gives larger weight.
             CALL BOX(GRID,XLEN,YLEN,SCRN,I,J,LENBOX,THETA, IBCTYPE)
             PTHETA = GRID(I,J,CELLDIR,SCRN)
             ANGDIFF = ABS(THETA - PTHETA)
             IF (ANGDIFF .GT. HALFPI) ANGDIFF = PI - ANGDIFF
             CANG = COS(ANGDIFF)
             WANG = CANG*CANG
C Thermal weighting: cells far behind the current tip temperature are
C downweighted.  This concentrates staged freezing near the active tip.
             DTBACK = TTIP - TL
             IF (DTBACK .LE. DT0) THEN
                WTHERM = 1.0
             ELSE IF (DTBACK .GE. DT1) THEN
                WTHERM = 0.0
             ELSE
                FRAC = (DTBACK - DT0)/(DT1 - DT0)
                WTHERM = 1.0 - FRAC
             ENDIF
C Final raw weight = thermal/orientation preference plus a floor.
             WRAW = FMIN + (1.0-FMIN)*WTHERM*WANG
C Store accepted candidate in sparse list.
             IF (WRAW .GT. 0.0) THEN
             IF (NSTG .LT. MAXSTG) THEN
               NSTG = NSTG + 1
               ISTG(NSTG) = I
               JSTG(NSTG) = J
               WSTG(NSTG) = WRAW
               CAPSTG(NSTG) = MAX(0.0, S1)
               IFRZSTG(NSTG) = 0
               WSUM = WSUM + WRAW
             ENDIF
             ENDIF
          ENDIF
          ENDIF
         ENDIF
      ENDDO
      ENDDO
C No accepted cells or zero total weight -> nothing to stage.
      IF (NSTG .LE. 0) RETURN
      IF (WSUM .LE. 0.0) RETURN
C Normalize weights so the total real-valued freezing matches TARG.
      SCALE = TARG/WSUM
      NSUM = 0
      DO K = 1, NSTG
        DFRZSTG(K) = WSTG(K) * SCALE
C Do not exceed the local liquid capacity of the cell.
        IF (DFRZSTG(K) .GT. CAPSTG(K)) DFRZSTG(K) = CAPSTG(K)
C Integer staging amount from the truncated real value.
        IFRZSTG(K) = INT(DFRZSTG(K))
        IF (IFRZSTG(K) .GT. INT(CAPSTG(K))) IFRZSTG(K)=INT(CAPSTG(K))
C Save fractional remainder for final leftover assignment.
        FRCSTG(K) = DFRZSTG(K) - FLOAT(IFRZSTG(K))
        NSUM = NSUM + IFRZSTG(K)
      ENDDO
C Assign any leftover integer subcells by largest fractional remainder.
      NTARG = NINT(TARG)
      NLEFT = NTARG - NSUM
      DONE = .FALSE.
      DO WHILE ((NLEFT .GT. 0).AND.(.NOT. DONE))
        BESTF = -1.0
        KBEST = 0
        DO K = 1, NSTG
         IF (IFRZSTG(K) .LT. INT(CAPSTG(K))) THEN
           IF (FRCSTG(K) .GT. BESTF) THEN
            BESTF = FRCSTG(K)
            KBEST = K
           ENDIF
         ENDIF
        ENDDO
        IF (KBEST .LE. 0) THEN
          DONE = .TRUE.
        ELSE
          IFRZSTG(KBEST) = IFRZSTG(KBEST) + 1
          FRCSTG(KBEST) = -1.0
          NLEFT = NLEFT - 1
        ENDIF
      ENDDO      
      RETURN
      END
C
C=====================================================================
C XGROWGATE -- Decide whether XSTAGE should run this timestep.
C If a staged tip list exists, use it as a quick undercooling screen.
C If no staged list exists, allow XSTAGE to scan normally. This avoids
C startup deadlock before a valid tip/stage list has been built.
C=====================================================================
      SUBROUTINE XGROWGATE(GRID,XLEN,YLEN,SCUR,
     &     TL0,MLIQ,C0,
     &     NSTAGEBOX,ISTAGEBOX,JSTAGEBOX,
     &     NGROW,DOGROW)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      REAL GRID(NXMAX, NYMAX, 6, 2)
      INTEGER XLEN, YLEN
      INTEGER LIQUID, TCELL, SCUR
      REAL TL0, MLIQ, C0
      INTEGER NSTAGEBOX
      INTEGER ISTAGEBOX(64), JSTAGEBOX(64)
      INTEGER NGROW
      LOGICAL DOGROW
      INTEGER K, II, JJ
      REAL CL, TLOC, TLIQ, DELT, DELTMAX
      REAL TLIQBIN
      PARAMETER (LIQUID=2, TCELL=5)
      EXTERNAL TLIQBIN
      DOGROW = .TRUE.
      IF (NSTAGEBOX .LE. 0) RETURN
      DELTMAX = -1.0E12
      DO K = 1, NSTAGEBOX
         II = ISTAGEBOX(K)
         JJ = JSTAGEBOX(K)
         IF (II.GE.1 .AND. II.LE.XLEN .AND.
     &       JJ.GE.1 .AND. JJ.LE.YLEN) THEN
            CL   = GRID(II,JJ,LIQUID,SCUR)
            TLOC = GRID(II,JJ,TCELL,SCUR)
            TLIQ = TLIQBIN(CL,TL0,MLIQ,C0)
            DELT = TLIQ - TLOC
            IF (DELT .GT. DELTMAX) DELTMAX = DELT
         ENDIF
      ENDDO
      IF (DELTMAX .LE. 0.0) THEN
         NGROW = 0
         DOGROW = .FALSE.
      ENDIF
      RETURN
      END
C
C=====================================================================
C XGROW -- Coordinate one growth/update pass.
C   1) decide whether growth should be attempted,
C   2) stage interface updates on SNEW using XSTAGE,
C   3) stop cleanly if the sparse update list overflows,
C   4) commit staged updates back to the active screen.
C=====================================================================
      SUBROUTINE XGROW(GRID, XADD, XRES, XRND, QLAT, XLEN, YLEN,
     &    SCRN, BOX2, TOL, TEUT, TL0, MLIQ, KPART, C0,
     &    GIBT, DD, DFRZ, CURVKGT, LHEAT, CPHEAT,
     &    VEL, GRAD, DX, DCOEFF, DT,
     &    GINCFRZ, GINCPRI, GINCEUT,
     &    DTBEST, TLBEST, TIPX, TIPY, HAVETIP, TTIP,
     &    TLX, TLY, TLGX, TLGY,
     &    GTXTIP, GTYTIP, GGTIP,
     &    TLCLTIP, TLLDIF, TLWID, TLLEN,
     &    ITIPLOAD,
     &    NSTAGEBOX, ISTAGEBOX, JSTAGEBOX, IFRZSTAGE,
     &    IUPD, JUPD, MAXUPD, NGROW, LENBOX, IBCTYPE, 
     &    CASTOP)
      IMPLICIT NONE
C Size limits
      INTEGER NXMAX, NYMAX, XLEN, YLEN
      PARAMETER (NXMAX=1000, NYMAX=1000)
C Main state arrays
      REAL GRID(NXMAX, NYMAX, 6, 2)
      REAL XADD(NXMAX, NYMAX)
      REAL XRES(NXMAX, NYMAX)
      REAL XRND(NXMAX, NYMAX)
      REAL QLAT(NXMAX, NYMAX)
      REAL LHEAT, CPHEAT
C Grid slot tags
C Active screen and local stencil controls
      INTEGER SCRN, BOX2, LENBOX, IBCTYPE
C Numerical controls
      REAL TOL
C Alloy and local equilibrium inputs
      REAL TEUT, TL0, MLIQ, KPART, C0
C Interface closure, capillarity, and anisotropy
      REAL GIBT, DD, DFRZ, CURVKGT
      REAL VEL, GRAD, DX, DCOEFF, DT
C Growth counters
      INTEGER GINCFRZ, GINCPRI, GINCEUT
      INTEGER NGROW
C Tip-tracking and local growth diagnostics
      REAL DTBEST, TLBEST
      REAL TTIP, TTIPREF
      INTEGER TIPX, TIPY
      LOGICAL HAVETIP
      REAL TLX, TLY, TLGX, TLGY
      REAL GTXTIP, GTYTIP, GGTIP
      REAL TLCLTIP, TLLDIF, TLWID, TLLEN
      INTEGER ITIPLOAD
C Sparse staged-box list
      INTEGER NSTAGEBOX
      INTEGER ISTAGEBOX(64), JSTAGEBOX(64), IFRZSTAGE(64)
C Sparse updated-cell list
      INTEGER MAXUPD
      INTEGER IUPD(MAXUPD), JUPD(MAXUPD)
C Local control
      INTEGER SCUR, SNEW, NUPD
      LOGICAL UPDOVER, CASTOP, DOGROW
      SCUR = SCRN
      SNEW = 3 - SCUR
C Use current tip temperature as the local growth reference.
      TTIPREF = TTIP
C Decide whether this timestep has any staged undercooling to grow.
      CALL XGROWGATE(GRID,XLEN,YLEN,SCUR,
     &     TL0,MLIQ,C0,
     &     NSTAGEBOX,ISTAGEBOX,JSTAGEBOX,
     &     NGROW,DOGROW)
      IF (DOGROW) THEN
C Build staged interface updates on SNEW and collect touched cells.
         CALL XSTAGE(GRID,XADD,XRES,XRND,QLAT,XLEN,YLEN,SCUR,SNEW,
     &      BOX2,TOL,
     &      TEUT,TL0,MLIQ,KPART,C0,
     &      GIBT,DD,DFRZ,CURVKGT,LHEAT,CPHEAT,
     &      VEL,GRAD,DX,DCOEFF,DT,
     &      GINCFRZ,GINCPRI,GINCEUT,
     &      DTBEST,TLBEST,TIPX,TIPY,TTIPREF,HAVETIP,
     &      TLX, TLY, TLGX, TLGY,
     &      GTXTIP, GTYTIP, GGTIP,
     &      TLCLTIP, TLLDIF, TLWID, TLLEN,
     &      ITIPLOAD,
     &      NSTAGEBOX, ISTAGEBOX, JSTAGEBOX, IFRZSTAGE,
     &      IUPD,JUPD,MAXUPD,NUPD,NGROW,UPDOVER,LENBOX,IBCTYPE)
C Stop cleanly if the sparse update list overflows.
         IF (UPDOVER) THEN
            WRITE(*,*) 'XGROW: MAXUPD overflow'
            CASTOP = .TRUE.
         ENDIF
C Commit updated cells from SNEW back to the active screen.
         CALL XCOMMIT(GRID,SCUR,SNEW, IUPD,JUPD,NUPD)
      ENDIF
      SCRN = SCUR
      RETURN
      END
C
C=====================================================================
C CATRANSFER -- apply staged solute redistribution and monitor mass
C   - applies queued solute deposits stored in XADD
C   - drains and redistributes any local residual solute in XRES
C   - evaluates mass-balance drift at output cadence
C=====================================================================
      SUBROUTINE CATRANSFER(GRID,XLEN,YLEN,SCRN,
     &     XADD,XRES,XRESADD,XRND,
     &     NREJ,IREJ,JREJ,MAXREJ,
     &     NTOUCH,ITCH,JITCH,MAXTOUCH,
     &     MLIQ,TL0,C0,
     &     GCOUNT1,DSPNUM,M0,MTDIFF, IBCTYPE, CASTOP)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN, SCRN
      REAL GRID(NXMAX,NYMAX,6,2)
      REAL XADD(NXMAX,NYMAX), XRES(NXMAX,NYMAX)
      REAL XRESADD(NXMAX,NYMAX)
      REAL XRND(NXMAX,NYMAX)
      INTEGER MAXREJ, MAXTOUCH, IBCTYPE
      INTEGER NREJ, NTOUCH
      INTEGER IREJ(MAXREJ), JREJ(MAXREJ)
      INTEGER ITCH(MAXTOUCH), JITCH(MAXTOUCH)
      REAL MLIQ, TL0, C0
      INTEGER GCOUNT1, DSPNUM
      DOUBLE PRECISION M0, CURMB
      DOUBLE PRECISION MASSBAL
      EXTERNAL MASSBAL
      REAL MTDIFF
      LOGICAL CASTOP
C Apply staged solute deposits recorded during interface ejection
      CALL DO1XADD(GRID,XLEN,YLEN,SCRN,XADD,XRES,
     &                NREJ, IREJ, JREJ, MAXREJ, CASTOP)
C Drain and redistribute any local residual solute reservoir
      CALL DO2XADD(GRID,XLEN,YLEN,SCRN,XRES,XRESADD, NREJ, IREJ,JREJ,
     &      NTOUCH, ITCH,JITCH, MAXTOUCH, MLIQ,TL0,C0, IBCTYPE,CASTOP)
C Update mass-balance drift only at output cadence
      MTDIFF = 0.0
      IF (MOD(GCOUNT1,DSPNUM) .EQ. 0) THEN
          CURMB = MASSBAL(GRID, XADD,XRES,XRND,XLEN,YLEN,
     &                      SCRN, CASTOP)
          MTDIFF = REAL((CURMB - M0)/(DBLE(XLEN)*DBLE(YLEN)))
      ENDIF
      RETURN
      END
C
C===================================================================
C TIPGRAD -- Current leading-tip thermal-gradient direction and mag.
C
C Computes a centered-difference thermal-gradient vector at the current
C leading tip.  The returned vector components are normalized direction
C cosines, while GGTIP is the physical gradient magnitude in K/m.
C
C This routine is intended for runtime KGT/projection use and is called
C every timestep before XGROW.  It is separate from XDIAGTIP, which is
C only called at the diagnostic/output cadence.
C===================================================================
      SUBROUTINE TIPGRAD(GRID,XLEN,YLEN,SCRN,TIPX,TIPY,DX,
     &     GTXTIP,GTYTIP,GGTIP)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      REAL GRID(NXMAX,NYMAX,6,2)
      INTEGER XLEN,YLEN,SCRN,TIPX,TIPY,TCELL
      REAL DX,GTXTIP,GTYTIP,GGTIP
      INTEGER IXM,IXP,IYM,IYP
      REAL DTDX,DTDY,GMAG,TOL
        PARAMETER (TCELL=5, TOL=1.0E-20)      
      IXM = MAX(1,TIPX-1)
      IXP = MIN(XLEN,TIPX+1)
      IYM = MAX(1,TIPY-1)
      IYP = MIN(YLEN,TIPY+1)
      DTDX = 0.0
      DTDY = 0.0
      IF (IXP .NE. IXM) THEN
        DTDX = (GRID(IXP,TIPY,TCELL,SCRN)-GRID(IXM,TIPY,TCELL,SCRN))
     &           /(FLOAT(IXP-IXM)*DX)
      ENDIF
      IF (IYP .NE. IYM) THEN
        DTDY = (GRID(TIPX,IYP,TCELL,SCRN)-GRID(TIPX,IYM,TCELL,SCRN))
     &           /(FLOAT(IYP-IYM)*DX)
      ENDIF
      GMAG = SQRT(DTDX*DTDX + DTDY*DTDY)
      GGTIP = GMAG
      IF (GMAG .GT. TOL) THEN
         GTXTIP = DTDX/GMAG
         GTYTIP = DTDY/GMAG
      ELSE
         GTXTIP = 0.0
         GTYTIP = 1.0
         GGTIP = 0.0
      ENDIF
      RETURN
      END
C
C=====================================================================
C XDIAGTIP -- tip diagnostics, EMA updates, live output, CSV write
C
C Computes the leading dendrite tip geometrically with FNDTDR() and
C evaluates measured tip velocity at the output cadence without
C modifying the CA state. For the KGT comparison, the local thermal-
C gradient direction near the tracked tip defines the analytical
C growth axis. The measured tip speed is projected onto that axis
C before calling KGT.
C=====================================================================
      SUBROUTINE XDIAGTIP(GRID,XRES,XLEN,YLEN,SCRN,
     &    GTIME,DT,DX,DSPNUM,
     &    IUNIT,
     &    TEUT,TL0,DFRZ,
     &    OTIPX,OTIPY,OTIME,OEMATIME,EINIT,TTIP,
     &    VELAVG,VINIT,VKGTMS,GGKGTUSE,
     &    TIPX,TIPY,DONE,
     &    CTAVG,UNCLAVG,CINIT,UINIT,IBCTYPE,
     &    C0,MLIQ,KALPHA,DD,GIBT,DCOEFF,
     &    GX,GY)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN, SCRN
      REAL GRID(NXMAX,NYMAX,6,2)
      REAL XRES(NXMAX,NYMAX)
      REAL XRMASS, XRSUMAVG
      REAL GTIME, DT, DX
      INTEGER DSPNUM, IUNIT, IBCTYPE
      REAL TEUT, TL0, DFRZ
      INTEGER TCELL, LIQUID
      INTEGER TIPX, TIPY, DONE
      INTEGER OTIPX, OTIPY
      REAL OTIME, OEMATIME, TTIP
      LOGICAL EINIT
      REAL VELAVG, VKGTMS, GGKGTUSE
      LOGICAL VINIT
      REAL CTAVG, UNCLAVG
      LOGICAL CINIT, UINIT
      REAL C0, MLIQ, KALPHA, DD, GIBT, DCOEFF
      REAL GX, GY
      INTEGER CELLDIR
      REAL TAU, VELUM, CTIP, UNCL
      REAL THETA, PTHETA
      REAL CURVKGT, DTKGT, PE
      REAL VVMEAS, GGKGT, RTIPKGT, IVKGT, DELT0KGT
      REAL DTCURVTIP
      REAL TCURV
      PARAMETER (LIQUID=2, TCELL=5, CELLDIR=6)
      EXTERNAL TCURV, XRSUMAVG
      INTEGER LENBOX
      PARAMETER (LENBOX = 7) 
      LOGICAL FOUNDTIP
      INTEGER IXF, IYF, IXP, IXM, IYP, IYM
      REAL AGX, AGY
      REAL GTXRAW, GTYRAW, GTMAG, GTXUSE, GTYUSE, PROJTG
      REAL GTXEMA, GTYEMA, EPSG, ALPHAG
      LOGICAL GINIT
      PARAMETER (EPSG=1.0E-12, ALPHAG=0.25)
      SAVE GTXEMA, GTYEMA, GINIT
      DATA GINIT /.FALSE./
      VELUM = 0.0
      CURVKGT = 0.0
      DTKGT = 0.0
      PE = 0.0
      THETA = 0.0
      DTCURVTIP = 0.0
      CALL FNDTDR(GRID,XLEN,YLEN,SCRN,GX,GY,TIPX,TIPY,FOUNDTIP)
      IF (.NOT. FOUNDTIP) THEN
         DONE = 0
         TIPX = OTIPX
         TIPY = OTIPY
         IF (TIPX .LT. 1) TIPX = 1
         IF (TIPX .GT. XLEN) TIPX = XLEN
         IF (TIPY .LT. 1) TIPY = 1
         IF (TIPY .GT. YLEN) TIPY = YLEN
      ELSE
         CALL SPLOT(GRID,SCRN,XLEN,YLEN,TIPX,TIPY,GTIME,TEUT,DONE)
      ENDIF
C Velocity update uses OTIME only
      TAU = 6.0E-2
      CALL VELTIPDIR(TIPX,TIPY,GTIME,OTIPX,OTIPY,OTIME,
     &               DX,TAU,GX,GY,VELAVG,VINIT,VELUM)
      TTIP = GRID(TIPX,TIPY,TCELL,SCRN)
C Sample liquid composition one cell ahead along growth direction
      AGX = ABS(GX)
      AGY = ABS(GY)
      IXF = TIPX
      IYF = TIPY
      IF (AGY .GE. AGX) THEN
         IF (GY .GE. 0.0) THEN
            IYF = TIPY + 1
         ELSE
            IYF = TIPY - 1
         ENDIF
      ELSE
         IF (GX .GE. 0.0) THEN
            IXF = TIPX + 1
         ELSE
            IXF = TIPX - 1
         ENDIF
      ENDIF
      IF (IXF .LT. 1) IXF = 1
      IF (IXF .GT. XLEN) IXF = XLEN
      IF (IYF .LT. 1) IYF = 1
      IF (IYF .GT. YLEN) IYF = YLEN
      CTIP = GRID(IXF,IYF,LIQUID,SCRN)
      UNCL = TL0 - TTIP
C EMA update uses independent EMA time
      IF (.NOT. EINIT) THEN
         OEMATIME = GTIME
         EINIT = .TRUE.
      ENDIF
      CALL DOEMA(DT, DSPNUM, GTIME, OEMATIME,
     &           CTIP, UNCL, CTAVG, CINIT, UNCLAVG, UINIT)
      OEMATIME = GTIME
C Estimate local thermal-gradient direction near the tracked tip.
C This defines the analytical growth axis used in the KGT comparison.
      IXP = TIPX + 1
      IF (IXP .GT. XLEN) IXP = 1
      IXM = TIPX - 1
      IF (IXM .LT. 1) IXM = XLEN
      IYP = TIPY + 1
      IF (IYP .GT. YLEN) IYP = YLEN
      IYM = TIPY - 1
      IF (IYM .LT. 1) IYM = 1
      GTXRAW=0.5*(GRID(IXP,TIPY,TCELL,SCRN)-GRID(IXM,TIPY,TCELL,SCRN))
      GTYRAW=0.5*(GRID(TIPX,IYP,TCELL,SCRN)-GRID(TIPX,IYM,TCELL,SCRN))
      GTMAG = SQRT(GTXRAW*GTXRAW + GTYRAW*GTYRAW)
      IF (GTMAG .GT. EPSG) THEN
         GTXRAW = GTXRAW/GTMAG
         GTYRAW = GTYRAW/GTMAG
         IF (.NOT. GINIT) THEN
            GTXEMA = GTXRAW
            GTYEMA = GTYRAW
            GINIT = .TRUE.
         ELSE
            GTXEMA = (1.0-ALPHAG)*GTXEMA + ALPHAG*GTXRAW
            GTYEMA = (1.0-ALPHAG)*GTYEMA + ALPHAG*GTYRAW
         ENDIF
         GTMAG = SQRT(GTXEMA*GTXEMA + GTYEMA*GTYEMA)
         IF (GTMAG .GT. EPSG) THEN
            GTXUSE = GTXEMA/GTMAG
            GTYUSE = GTYEMA/GTMAG
         ELSE
            GTXUSE = 0.0
            GTYUSE = 1.0
         ENDIF
         GGKGTUSE = SQRT((0.5*(GRID(IXP,TIPY,TCELL,SCRN)
     &                -GRID(IXM,TIPY,TCELL,SCRN))/DX)**2 +
     &                (0.5*(GRID(TIPX,IYP,TCELL,SCRN)
     &                -GRID(TIPX,IYM,TCELL,SCRN))/DX)**2)
      ELSE
         IF (GINIT) THEN
            GTMAG = SQRT(GTXEMA*GTXEMA + GTYEMA*GTYEMA)
            IF (GTMAG .GT. EPSG) THEN
               GTXUSE = GTXEMA/GTMAG
               GTYUSE = GTYEMA/GTMAG
            ELSE
               GTXUSE = 0.0
               GTYUSE = 1.0
            ENDIF
         ELSE
            GTXUSE = 0.0
            GTYUSE = 1.0
         ENDIF
         GGKGTUSE = 0.0
      ENDIF
C Project measured tip speed into the thermal-gradient frame for KGT.
      PROJTG = ABS(GX*GTXUSE + GY*GTYUSE)
      IF (PROJTG .GT. 1.0) PROJTG = 1.0
      VVMEAS = VELUM*1.0E-6
      VKGTMS = VVMEAS*PROJTG
      GGKGT = GGKGTUSE
      CALL KGT(VKGTMS,GGKGT,GIBT,DCOEFF,MLIQ,KALPHA,C0,
     &         DTKGT,RTIPKGT,CURVKGT,PE,IVKGT,DELT0KGT)
C Orientation and curvaure
      CALL BOX(GRID,XLEN,YLEN,SCRN,TIPX,TIPY,LENBOX, THETA,IBCTYPE)
      PTHETA = GRID(TIPX,TIPY,CELLDIR,SCRN)
      DTCURVTIP = TCURV(GIBT,THETA,PTHETA,CURVKGT,DD)
C save and show output
      XRMASS = XRSUMAVG(XRES,XLEN,YLEN)
      CALL CSVOUT(IUNIT, GTIME, VELUM, TIPX, TIPY, TTIP,
     &     CTAVG, UNCLAVG, DTKGT, CURVKGT, THETA, DFRZ,
     &     DCOEFF, TL0, MLIQ, C0, KALPHA, DTCURVTIP, XRMASS)
      FLUSH(IUNIT)
      CALL LIVEOUT(VELUM,TL0,TTIP,CURVKGT,DTKGT, PE,DFRZ,
     &                   MLIQ,C0,KALPHA,DTCURVTIP,XRMASS)
      RETURN
      END
C
C=====================================================================
C XUPDATE -- Accumulate compact run metrics and diagnostics
C Inputs:
C GCOUNT1  - global step counter
C GINCFRZ  - sub-cells frozen this step
C GINCPRI  - primary cells captured this step
C GINCEUT  - eutectic cells formed this step
C DEUTE    - solute mass to eutectic this step
C MTDIFF   - mass delta wt% from initial
C Persistent:
C TOTFRZ, TOTPRI, TOTEUT, TOTDEUTE, MAXFRATE
C=====================================================================
      SUBROUTINE XUPDATE(GCOUNT1, GINCFRZ, GINCPRI, GINCEUT,
     &                   DEUTE, XLEN, YLEN,  MTDIFF,
     &                   TOTFRZ, TOTPRI, TOTEUT, TOTDEUTE, MAXFRATE)
      IMPLICIT NONE
      REAL    DEUTE, ZUNITS, MTDIFF
      INTEGER GCOUNT1, GINCFRZ, GINCPRI, GINCEUT, XLEN, YLEN
      PARAMETER (ZUNITS=1000.0)
C Persistent state (caller-owned)
      INTEGER TOTFRZ, TOTPRI, TOTEUT
      REAL    TOTDEUTE, MAXFRATE
C Locals
      REAL    SCRATCH, TOTCELLS, FRZSTEP, FRZCUM
C Accumulate
      TOTFRZ = TOTFRZ   + GINCFRZ
      TOTPRI = TOTPRI   + GINCPRI
      TOTEUT = TOTEUT   + GINCEUT
      TOTDEUTE = TOTDEUTE + DEUTE
C Compute fractions
      SCRATCH = REAL(XLEN*YLEN)
      TOTCELLS = SCRATCH*ZUNITS
      FRZSTEP = REAL(GINCFRZ)/TOTCELLS
      FRZCUM = REAL(TOTFRZ)/TOTCELLS
C Update max rate
      IF (FRZSTEP .GT. MAXFRATE) MAXFRATE = FRZSTEP
      WRITE(*,9002) GINCFRZ, TOTFRZ
      WRITE(*,9003) GINCEUT, TOTEUT
      WRITE(*,9004) GCOUNT1, MTDIFF
9002  FORMAT(' Frozen sub-cells (step,total):',I8,1X,I10)
9003  FORMAT(' Eutectic cells  (step,total):',I8,1X,I10)
9004  FORMAT(' Step=',I7,'   Total mass delta(wt%)=',E15.7)
      RETURN
      END
C
C=====================================================================
C KGT -- Analytical dendritic tip kinetics (Kurz-Giovanola-Trivedi)
C Computes steady-state tip undercooling and radius using Ivantsov
C transport with marginal-stability (MSC) tip selection.
C Formulation:
C   Omega = Iv(Pe) = Pe*exp(Pe)*E1(Pe)     (solute transport)
C   R set by Sigma* ~ 0.025                (MSC selection)
C   DT includes thermal-gradient correction (G)
C Inputs:
C   V       (m/s)   Growth velocity
C   G       (K/m)   Thermal gradient
C   GIBT    (K-m)   Gibbs-Thomson coefficient
C   DCOEFF  (m2/s) Liquid diffusivity
C   MLIQ    (K/wt%) Liquidus slope
C   KALPHA  (-)     Partition coefficient
C   C0      (wt%)   Nominal composition
C Outputs:
C   DTKGT   (K)     Tip undercooling
C   RTIP    (m)     Tip radius
C   CURVKGT (1/m)   Tip curvature
C   PE      (-)     Peclet number (V*R/2D)
C   IVR     (-)     Iv(Pe)/Pe
C   DELT0   (K)     Freezing range (-m*C0*(1-k)/k)
C References:
C   Kurz et al., Acta Metall. 34 (1986)
C   Langer & Muller-Krumbhaar, Acta Metall. 26 (1978)
C   Ivantsov, Dokl. Akad. Nauk SSSR 58 (1947)
C=====================================================================
      SUBROUTINE KGT(V,G,GIBT,DCOEFF,MLIQ,KALPHA,C0,
     &                DTKGT,RTIP,CURVKGT,PE,IVR,DELT0)
      IMPLICIT NONE
      REAL V,G,GIBT,DCOEFF,MLIQ,KALPHA,C0
      REAL DTKGT,RTIP,CURVKGT,PE,IVR,DELT0
      REAL TOL, IV, XDELT0, XRTIP, IVANTSOV
      PARAMETER (TOL = 1.0E-12)
C Set Baseline Defaults
      DTKGT = 0.0
      RTIP = 1.0E10
      CURVKGT = 0.0
      PE = 0.0
      IVR = 1.0
      DELT0 = 0.0
      IF (V .GT. 0.0 .AND. DCOEFF .GT. TOL .AND. GIBT .GT. 0.0) THEN
       DELT0 = XDELT0(MLIQ, C0, KALPHA)
       IF (DELT0 .GT. TOL) THEN
         RTIP = XRTIP(V, G, GIBT, DCOEFF, DELT0)         
         IF (RTIP .GT. TOL) CURVKGT = 1.0/RTIP
         PE = V*RTIP/(2.0*DCOEFF)
         IV = IVANTSOV(PE)
         IF (PE .GT. TOL) THEN
            IVR = IV/PE
         ELSE
            IVR = 1.0
         ENDIF
C Final undercooling corrected by a thermal gradient factor 
         DTKGT = DELT0*IV/(1.0+(G*GIBT/(V*DELT0*DELT0)))
        ENDIF
      ENDIF
      RETURN
      END
C
C=====================================================================
C FNDTDR -- Geometric leading interface tip along direction (GX,GY)
C   Scans the grid and identifies the interface cell that is most
C   advanced along the specified direction vector (GX,GY).
C   The projected coordinate is   S = X*GX + Y*GY
C   and the cell with the largest S is selected as the geometric tip.
C   This makes the tip detection independent of the grid orientation.
C   Tie-break rule (canidates are interface cells):
C     the cell with the higher temperature TCELL is chosen.
C=====================================================================
      SUBROUTINE FNDTDR(GRID,XLEN,YLEN,SCRN,
     &    GX,GY,TIPX,TIPY,FOUND)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX, LCELLS, TCELL
      PARAMETER (NXMAX=1000, NYMAX=1000, LCELLS=1, TCELL=5)
      INTEGER XLEN, YLEN, SCRN
      REAL GRID(NXMAX,NYMAX,6,2)
      REAL GX, GY
      INTEGER TIPX, TIPY
      LOGICAL FOUND
      INTEGER I, J
      REAL SBEST, TLBEST, STOL, ZUNITS
      REAL S, TL, S1
      REAL XI, YJ, YTERM
        PARAMETER (STOL = 1.0E-6, ZUNITS = 1000.0)
      FOUND = .FALSE.
      TIPX = 1
      TIPY = 1
      SBEST = -1.0E30
      TLBEST = -1.0E30
      DO J = 1, YLEN
      YJ = FLOAT(J)
      YTERM = YJ*GY
       DO I = 1, XLEN
        S1 = GRID(I,J,LCELLS,SCRN)
        IF (S1 .GT. 0.0 .AND. S1 .LT. ZUNITS) THEN
          XI = FLOAT(I)
          S = XI*GX + YTERM
          IF (S .GT. SBEST + STOL) THEN
            TL = GRID(I,J,TCELL,SCRN)
            SBEST = S
            TLBEST = TL
            TIPX = I
            TIPY = J
            FOUND = .TRUE.
          ELSE IF (ABS(S-SBEST) .LE. STOL) THEN
            TL = GRID(I,J,TCELL,SCRN)
            IF (TL .GT. TLBEST) THEN
              TLBEST = TL
              TIPX = I
              TIPY = J
              FOUND = .TRUE.
            ENDIF
          ENDIF
        ENDIF
       ENDDO
      ENDDO
      RETURN
      END
C
C=====================================================================
C VELTIPDIR -- EMA-smoothed dendrite tip velocity along direction
C   Computes the dendrite tip velocity from the change in projected
C   tip position along a specified direction vector:
C        S = X*GX + Y*GY
C   The raw speed is DS/DT (cells/s) where DS is the change in the
C   projected tip coordinate.  An exponential moving average (EMA)
C   with time scale TAU is applied to stabilize the reported speed.
C   Output: VELAVG = EMA velocity (cells/s)
C           VELUM  = velocity converted to um/s
C   OTIPX, OTIPY, OTIME store the previous reference state.
C=====================================================================
      SUBROUTINE VELTIPDIR(TIPX,TIPY,GTIME,
     &   OTIPX,OTIPY,OTIME,
     &   DX,TAU,GX,GY,
     &   VELAVG,VINIT,VELUM)
      IMPLICIT NONE
      INTEGER TIPX, TIPY, OTIPX, OTIPY
      REAL GTIME, OTIME, DX, TAU
      REAL GX, GY
      REAL VELAVG, VELUM
      LOGICAL VINIT
      REAL DT, ALPHA
      REAL SNEW, SOLD, DS, VRAW
      REAL TAUEFF
C Define reference state so first report is V=0
      IF (.NOT. VINIT) THEN
         OTIPX = TIPX
         OTIPY = TIPY
         OTIME = GTIME
         VELAVG = 0.0
         VELUM = 0.0
         VINIT = .TRUE.
      ELSE
         DT = GTIME - OTIME
         IF (DT .GT. 0.0) THEN
            SOLD = FLOAT(OTIPX)*GX + FLOAT(OTIPY)*GY
            SNEW = FLOAT(TIPX )*GX + FLOAT(TIPY )*GY
            DS = SNEW - SOLD
C Use only forward advance along the chosen direction
            IF (DS .LT. 0.0) DS = 0.0
            VRAW = DS/DT
            TAUEFF = AMAX1(TAU, 0.0)
            ALPHA = DT/(TAUEFF + DT)
            VELAVG = VELAVG + ALPHA*(VRAW - VELAVG)
         ENDIF
C Convert cells/s to um/s
         VELUM = VELAVG*DX*1.0E6
C Update reference state
         OTIPX = TIPX
         OTIPY = TIPY
         OTIME = GTIME
      ENDIF
      RETURN
      END
C  
C=====================================================================
C DOEMA -- EMA-smoothed diagnostics for composition and undercooling
C Only smooths:
C   - tip liquid composition (CTIP to CTAVG)
C   - tip undercooling    (UNCL to UNCLAVG)
C Default time constants tied to output cadence:
C   DTOUT = DSPNUM*DT
C   TAUC  = 3.0*DTOUT    (composition)
C   TAUU  = 2.0*DTOUT    (undercooling)
C Uses DTEMA = GTIME - OTIME (actual time between output calls)
C=====================================================================
      SUBROUTINE DOEMA(DT, DSPNUM, GTIME, OTIME,
     &                 CTIP, UNCL, CTAVG, CINIT, UNCLAVG, UINIT)
      IMPLICIT NONE
      REAL DT, GTIME, OTIME
      INTEGER DSPNUM
      REAL CTIP, UNCL
      REAL CTAVG, UNCLAVG
      LOGICAL CINIT, UINIT
      REAL DTEMA, DTOUT, TAUC, TAUU
C composition & undercooling smoothing
      DTOUT = FLOAT(DSPNUM)*DT
      TAUC = 3.0*DTOUT           
      TAUU = 2.0*DTOUT           
      DTEMA = GTIME - OTIME
      CALL EMASCL(CTIP, TAUC, DTEMA, CTAVG, CINIT)
      CALL EMASCL(UNCL, TAUU, DTEMA, UNCLAVG, UINIT)
      RETURN
      END
C    
C=====================================================================
C  BOX  –  Crystallographic orientation via second-moment tensor
C  Estimates the preferred growth direction (principal axis angle θ)
C  from the spatial distribution of solid phase in a local square
C  neighborhood. The local weight is based on solid fraction:
C        w = 1 - fL   where fL is the liquid fraction.
C  The weighted second-moment tensor about the box center is
C        M = [ Sxx  Sxy ]    with  Sxx = sum(w*dx*dx)
C            [ Sxy  Syy ]          Syy = sum(w*dy*dy)
C                                  Sxy = sum(w*dx*dy)
C  The eigenvector associated with the largest eigenvalue gives the
C  direction of maximum elongation. The corresponding orientation is
C      theta = 0.5*atan2(2*Sxy, Sxx - Syy)   on [0, pi) radians
C  Reference:
C  Bigun, J., Granlund, G.H., Wiklund, J.,
C  "Multidimensional orientation estimation with applications
C  to texture analysis and optical flow",
C  IEEE Trans. Pattern Anal. Mach. Intell. 13(8) (1991) 775–790
C  doi: 10.1109/34.85668
C=====================================================================
      SUBROUTINE BOX(GRID, XLEN, YLEN, SCRN, X, Y, LEN,THETA, IBCTYPE)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN, X, Y, LEN, SCRN
      INTEGER XI, YJ, II, JJ, NS, IBCTYPE
      REAL GRID(NXMAX, NYMAX, 6, 2), THETA, ZUNITS
      REAL    SXX, SYY, SXY
      REAL    LFRAC, WGT, TOL, PI
      INTEGER LCELLS
      PARAMETER (LCELLS = 1)
      PARAMETER (TOL =1.0E-6, PI=3.14159265, ZUNITS=1000.0)
        SXX = 0.0
        SYY = 0.0
        SXY = 0.0
        NS = 0
      DO JJ = -LEN, LEN
      DO II = -LEN, LEN
C Use direct indexing for interior stencil points.  Call GETBCIJ only
C when the requested stencil point falls outside the active domain.
         XI = X + II
         YJ = Y + JJ
         IF (XI .LT. 1 .OR. XI .GT. XLEN .OR.
     &       YJ .LT. 1 .OR. YJ .GT. YLEN) THEN
            CALL GETBCIJ(XI,YJ,XLEN,YLEN,IBCTYPE,XI,YJ)
         ENDIF
C Liquid fraction and corresponding solid-based weight
         LFRAC = GRID(XI, YJ, LCELLS, SCRN)/ZUNITS
         WGT = 1.0 - AMIN1(1.0, AMAX1(0.0, LFRAC))
C Accumulate weighted second moments
         IF (WGT .GT. TOL) THEN
            NS = NS + 1
            SXX = SXX + WGT*FLOAT(II*II)
            SYY = SYY + WGT*FLOAT(JJ*JJ)
            SXY = SXY + WGT*FLOAT(II*JJ)
         ENDIF
      END DO
      END DO
C Principal-axis orientation from the second-moment tensor
      IF (NS .GT. 0) THEN
         THETA = 0.5*ATAN2(2.0*SXY, SXX-SYY)
      ELSE
         THETA = 0.0
      ENDIF
      IF (THETA .LT. 0.0) THETA = THETA + PI
      RETURN
      END
C
C =============================================================
C Tcurv - curvature undercooling (optimized version)
C Thermodynamic: less undercooling in preferred direction 
C        TEMP = 1.0 - 15.0*DD*COS(4.0*(THETA-PTHETA))
C        TCURV = GIBT*CURV*TEMP
C Optimization:
C   Use COS(4D)= 1-2*SIN^2(2D), only one SIN() call is needed.
C =============================================================
      REAL FUNCTION TCURV(GIBT, THETA, PTHETA, CURV, DD)
      IMPLICIT NONE
      REAL GIBT, THETA, PTHETA, CURV, DD
      REAL PI2, PI4, D, AF, S2, C4, TEMP
      PARAMETER (PI2 = 2.0*ATAN(1.0))
      PARAMETER (PI4 = ATAN(1.0))
C D = THETA - PTHETA, folded into [-pi/4, +pi/4)
      D = THETA - PTHETA
      C4 = MOD(D + PI4, PI2)
      IF (C4 .LT. 0.0) C4 = C4 + PI2
      D = C4 - PI4
C SIMPLE ANISOTROPY CORRECTION FACTORS (by 20%)
C   FTH4(D): 1 + 0.2*sin^2(2D)
      S2 = SIN(2.0*D)
      AF = 1.0 + 0.2*S2*S2
C COS(4D) via identity: COS(4D) = 1 - 2*SIN^2(2D)
      C4 = 1.0 - 2.0*S2*S2
      TEMP = 1.0 - 15.0*DD*C4*AF
      TCURV = GIBT*CURV*TEMP
      RETURN
      END
C
C=====================================================================
C CSVOUT -- lean diagnostics
C=====================================================================
      SUBROUTINE CSVOUT(IUNIT, GTIME, VELUM, TIPX, TIPY, TTIP,
     &     CTAVG, UNCLAVG, DTKGT, CURVKGT, THETA, DFRZ,
     &     DCOEFF, TL0, MLIQ, C0, KALPHA, DTCURVTIP, XRMASS)
      IMPLICIT NONE
      INTEGER IUNIT, TIPX, TIPY
      REAL GTIME, VELUM, TTIP, CTAVG, UNCLAVG, DTKGT, CURVKGT
      REAL THETA, DFRZ, DCOEFF
      REAL TL0, MLIQ, C0, KALPHA, DTCURVTIP
      REAL VELM, PECLET, RADKGT, XRMASS
      REAL CLSTAR, OMEGACA, OMEGAKGT, DEN, SP
      REAL CLIQBIN
      EXTERNAL CLIQBIN
C VELUM is in um/s, convert to m/s
      VELM = VELUM*1.0E-6
C Peclet number: P = V*R/(2D) = V/(2Dk), where k = CURVKGT
      PECLET = 0.0
      IF (DCOEFF .GT. 0.0 .AND. CURVKGT .GT. 1.0E-12) THEN
         PECLET = VELM/(2.0*DCOEFF*CURVKGT)
      ENDIF
C KGT radius in um
      RADKGT = 0.0
      IF (CURVKGT .GT. 1.0E-12) THEN
         RADKGT = 1.0E6/CURVKGT
      ENDIF
C Tip liquid composition corrected for curvature undercooling
      IF (ABS(MLIQ) .GT. 1.0E-12) THEN
         CLSTAR = CLIQBIN(TTIP + DTCURVTIP,TL0,MLIQ,C0)
      ELSE
         CLSTAR = C0
      ENDIF
C Instantaneous CA supersaturation
      DEN = (1.0 - KALPHA)*CLSTAR
      IF (ABS(DEN) .GT. 1.0E-12) THEN
         OMEGACA = (CLSTAR - C0)/DEN
      ELSE
         OMEGACA = 0.0
      ENDIF
C 2-D Ivantsov supersaturation from Peclet number
      OMEGAKGT = 0.0
      IF (PECLET .GT. 0.0) THEN
         SP = SQRT(PECLET)
         OMEGAKGT = SQRT(3.14159265358979*PECLET)
     &            * EXP(PECLET) * ERFC(SP)
      ENDIF
      WRITE(IUNIT,9200) GTIME, VELUM, TIPX, TIPY, TTIP,
     &     CTAVG, UNCLAVG, RADKGT, CURVKGT, THETA, DFRZ,
     &     DTKGT, PECLET, CLSTAR, OMEGACA, OMEGAKGT, XRMASS

9200  FORMAT(E13.6,',',E13.6,',',I6,',',I6,',',E13.6,',',E13.6,
     &       ',',E13.6,',',E13.6,',',E13.6,',',E13.6,',',E13.6,
     &       ',',E13.6,',',E13.6,',',E13.6,',',E13.6,',',E13.6, E13.6)
      RETURN
      END
C
C=====================================================================
C LIVEOUT -- print compact live tip diagnostics.
C=====================================================================
      SUBROUTINE LIVEOUT(VELUM,TL0,TTIP,CURVKGT,DTKGT, PE,DFRZ,
     &                   MLIQ,C0,KALPHA,DTCURVTIP, XRMASS)
      IMPLICIT NONE
      REAL VELUM, TL0, TTIP, CURVKGT, PE, DFRZ
      REAL MLIQ, C0, KALPHA, DTCURVTIP, DTKGT, XRMASS
      REAL RTIPUM, DTCA, CLSTAR, DEN, OMEGACA, OMEGAKGT
      REAL CLIQBIN
      EXTERNAL CLIQBIN
      REAL SP
      DTCA = TL0 - TTIP
      RTIPUM = 0.0
      IF (CURVKGT .GT. 1.0E-6) RTIPUM = 1.0E6/CURVKGT
C Snapshot CA supersaturation using curvature-corrected tip liquid.
      IF (ABS(MLIQ) .GT. 1.0E-12) THEN
         CLSTAR = CLIQBIN(TTIP + DTCURVTIP,TL0,MLIQ,C0)
      ELSE
         CLSTAR = C0
      ENDIF
      DEN = (1.0 - KALPHA)*CLSTAR
      IF (ABS(DEN) .GT. 1.0E-12) THEN
         OMEGACA = (CLSTAR - C0)/DEN
      ELSE
         OMEGACA = 0.0
      ENDIF
C KGT supersaturation for LIVEOUT comparison only.
C Use the 2-D Ivantsov form here so the live output matches
C the CA/Python comparison plot.  
      OMEGAKGT = 0.0
      IF (PE .GT. 0.0) THEN
         SP = SQRT(PE)
         OMEGAKGT = SQRT(3.14159265358979*PE)*EXP(PE)*ERFC(SP)
      ENDIF
      WRITE(*,*) 'VEL[um/s]=',VELUM,'  Peclet=',PE
      WRITE(*,*) 'omegaKGT =',OMEGAKGT,' omegaCA=',OMEGACA
      WRITE(*,*) 'curvKGT(1/m)=',CURVKGT,' RtipKGT[um]=',RTIPUM
      WRITE(*,*) 'delT:KGT[K]=',DTKGT,' delT:CA[K]=',DTCA
      WRITE(*,*) 'DFRZ=',DFRZ,' XRMASS=',XRMASS
      RETURN
      END
C
C ======================================================== 
C Live diagonostic output using terminal text screen
C ========================================================
      SUBROUTINE SPLOT(GRID,SCRN,XLEN,YLEN,
     &                 TIPX,TIPY,GTIME,TEUT,DONE)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN, SCRN
      REAL GRID(NXMAX,NYMAX,6,2), GTIME, TLOCAL, TEUT
      INTEGER*1 SCREEN(70,20)
      INTEGER I, J, IMAX, JMAX, TIPX, TIPY, DONE, TOPDONE
      INTEGER LCELLS, TCELL
      PARAMETER (LCELLS=1, TCELL=5)
C Blank the ascii output screen
      DONE = 0
      IMAX = 0
      JMAX = 0
      DO J = 1, 20
      DO I = 1, 70
         SCREEN(I,J) = 32
      ENDDO
      ENDDO
C Set plot bounds around supplied tip 
      IMAX = TIPX - 35
      IF (IMAX .LT. 0) IMAX = 0
      IF (XLEN .LE. 70) THEN
         IMAX = 0
      ELSE IF ((IMAX+70) .GT. XLEN) THEN
         IMAX = XLEN - 70
      ENDIF
      JMAX = TIPY - 19
      IF (JMAX .LT. 0) JMAX = 0
      IF (YLEN .LE. 20) THEN
         JMAX = 0
      ELSE IF ((JMAX+20) .GT. YLEN) THEN
         JMAX = YLEN - 20
      ENDIF
C Plot CA state into ascii screen
      DO J = 1, 20
      DO I = 1, 70
       IF (GRID(I+IMAX,J+JMAX,LCELLS,SCRN) .LT. 1.0) SCREEN(I,J) = 42
      ENDDO
      ENDDO
C Print ascii plot
      DO J = 20, 1, -1
         WRITE(*,*) (CHAR(SCREEN(I,J)), I=1,70)
      ENDDO
      TLOCAL = GRID(TIPX,TIPY,TCELL,SCRN)
      TOPDONE = 60
      IF (YLEN .LT. TOPDONE) TOPDONE = YLEN - 1
      IF (TOPDONE .LT. 1) TOPDONE = 1
      IF ((TIPY .GT. YLEN-TOPDONE) .OR. (TLOCAL .LT. TEUT)) DONE = 1
      WRITE(*,1300) GTIME, TIPX, TIPY, TLOCAL
1300  FORMAT(' TIME=',F9.4,'  TIPX=',I5,'  TIPY=',I5,'  TIPTEMP=',F7.2)
      RETURN
      END
C
C=====================================================================
C DO1XADD -- apply staged local solute and build residual re-eject list
C - releases a fraction of fresh staged solute from XADD each step
C - adds available solute directly to local liquid when liquid remains
C - retains unresolved local solute in XRES
C - records fully solid residual-source cells for later re-ejection
C=====================================================================
      SUBROUTINE DO1XADD(GRID,XLEN,YLEN,ISCR,XADD,XRES,
     &                NREJ, IREJ, JREJ, MAXREJ, CASTOP)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN, ISCR
      INTEGER I, J
      REAL GRID(NXMAX,NYMAX,6,2)
      REAL XADD(NXMAX,NYMAX), XRES(NXMAX,NYMAX)
      INTEGER NREJ, MAXREJ
      INTEGER IREJ(MAXREJ), JREJ(MAXREJ)
      INTEGER LCELLS, LIQUID
      PARAMETER (LCELLS=1, LIQUID=2)
      REAL MINS1
      PARAMETER (MINS1=0.5)
C Fraction of staged solute released each step (XADD to liquid)
C (1-BETAADD) remains temporarily in XADD, providing smoothing of
C solute injection and reducing startup/transient oscillations
      REAL BETAADD
      PARAMETER (BETAADD=0.8) 
      REAL S1, XNEW, XHOLD, PEND
      LOGICAL CASTOP
      NREJ = 0
      DO J = 1, YLEN
      DO I = 1, XLEN
C Split fresh staged solute into released and held parts
         XNEW = BETAADD*XADD(I,J)
         XHOLD = XADD(I,J) - XNEW
         IF (XHOLD .LT. 0.0) XHOLD = 0.0
C Solute made available during this step
         PEND = XNEW + XRES(I,J)
         IF (PEND .GT. 0.0 .OR. XHOLD .GT. 0.0) THEN
           S1 = GRID(I,J,LCELLS,ISCR)
           IF (S1 .GT. MINS1) THEN
C Add available solute directly to the local liquid
             IF (PEND .GT. 0.0) THEN
              GRID(I,J,LIQUID,ISCR)=GRID(I,J,LIQUID,ISCR)+PEND/S1
             ENDIF
C Local residual is drained; delayed fraction remains in XADD
            XRES(I,J) = 0.0
            XADD(I,J) = XHOLD
           ELSE
C Keep available solute local; only fully solid cells enter
C the residual re-eject list
            XRES(I,J) = PEND
            XADD(I,J) = XHOLD
            IF (S1 .LE. 0.0 .AND. XRES(I,J) .GT. 0.0) THEN
              IF (NREJ .LT. MAXREJ) THEN
                 NREJ = NREJ + 1
                 IREJ(NREJ) = I
                 JREJ(NREJ) = J
              ELSE
                 WRITE(*,*) 'DO1XADD: MAXREJ overflow'
                 CASTOP = .TRUE.
              ENDIF
            ENDIF
           ENDIF
         ENDIF
      END DO
      END DO
      RETURN
      END
C
C=====================================================================
C DO2XADD -- eject local residual solute from listed source cells.
C   - Processes only the cells listed in (IREJ,JREJ).
C   - Stages accepted neighbor deposits in XRESADD.
C   - Records flagged neighbors in (ITCH,JITCH) for sparse apply/clear.
C
C This routine is often inactive:
C   - It preserves mass if staged rejected solute cannot be deposited
C     locally because nearby liquid capacity vanishes.
C   - The routine is retained for large domains, high undercooling,
C     and fast growth cases, such as additive manufacturing.
C=====================================================================
      SUBROUTINE DO2XADD(GRID,XLEN,YLEN,ISCR,XRES,
     &           XRESADD, NREJ, IREJ, JREJ, NTOUCH,
     &           ITCH, JITCH, MAXTOUCH,MLIQ,TL0,C0, IBCTYPE, CASTOP)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN, ISCR
      REAL GRID(NXMAX,NYMAX,6,2)
      REAL XRES(NXMAX,NYMAX)
      REAL XRESADD(NXMAX,NYMAX)
      REAL MLIQ, TL0, C0
      INTEGER NREJ
      INTEGER IREJ(*), JREJ(*)
      INTEGER NTOUCH, MAXTOUCH, IBCTYPE
      INTEGER ITCH(MAXTOUCH), JITCH(MAXTOUCH)
      INTEGER K, I, J, II, JJ, I2, J2
      INTEGER LCELLS, LIQUID, TCELL
      PARAMETER (LCELLS=1, LIQUID=2, TCELL=5)
      REAL DIFFMASK, TOL
      PARAMETER (DIFFMASK=0.0, TOL=1.0E-8)
      REAL S1, PEND, CL, TL, CLSTAR, DELCL
      REAL CLIQBIN
      EXTERNAL CLIQBIN
      REAL CAP, CAPSUM, W, PUT, LEFT
      LOGICAL CASTOP
      NTOUCH = 0
      IF (NREJ .LE. 0) RETURN
      IF (ABS(MLIQ) .LE. TOL) RETURN
C Build staged neighbor deposits only from listed residual-source cells
      DO K = 1, NREJ
         I = IREJ(K)
         J = JREJ(K)
         PEND = XRES(I,J)
         IF (PEND .GT. 0.0) THEN
C Re-eject only if the source cell remains fully solid
           S1 = GRID(I,J,LCELLS,ISCR)
           IF (S1 .LE. 0.0) THEN
C Sum admissible receiver capacity in the outward neighborhood
              CAPSUM = 0.0
              DO JJ = -1, 1
              DO II = -1, 1
                IF (.NOT.(II.EQ.0 .AND. JJ.EQ.0)) THEN
                 CALL GETBCIJ(I+II,J+JJ,XLEN,YLEN,IBCTYPE,I2,J2)
                 S1 = GRID(I2,J2,LCELLS,ISCR)
                 IF (S1 .GT. DIFFMASK) THEN
                    CL = GRID(I2,J2,LIQUID,ISCR)
                    TL = GRID(I2,J2,TCELL ,ISCR)
                    CLSTAR = CLIQBIN(TL,TL0,MLIQ,C0)
                    DELCL = CLSTAR - CL
                    IF (DELCL .GT. 0.0) THEN
                      CAP = S1*DELCL
                      CAPSUM = CAPSUM + CAP
                    ENDIF
                 ENDIF
                ENDIF
              END DO
              END DO
C No admissible local receiver capacity; keep residual in place
              IF (CAPSUM .LE. TOL) THEN
                XRES(I,J) = PEND
              ELSE
C Accepted portion of the residual source mass
                PUT = AMIN1(PEND, CAPSUM)
                LEFT = PEND - PUT
C Stage deposits into XRESADD in proportion to local capacity
                DO JJ = -1, 1
                DO II = -1, 1
                 IF (.NOT.(II.EQ.0 .AND. JJ.EQ.0)) THEN
                   CALL GETBCIJ(I+II,J+JJ,XLEN,YLEN,IBCTYPE,I2,J2)
                   S1 = GRID(I2,J2,LCELLS,ISCR)
                   IF (S1 .GT. DIFFMASK) THEN
                     CL = GRID(I2,J2,LIQUID,ISCR)
                     TL = GRID(I2,J2,TCELL ,ISCR)
                     CLSTAR = CLIQBIN(TL,TL0,MLIQ,C0)
                     DELCL = CLSTAR - CL
                     IF (DELCL .GT. 0.0) THEN
                        CAP = S1*DELCL
                        W = CAP/CAPSUM
                        XRESADD(I2,J2) = XRESADD(I2,J2) + PUT*W
C Record touched neighbors for sparse apply/clear
                        IF (NTOUCH .LT. MAXTOUCH) THEN
                           NTOUCH = NTOUCH + 1
                           ITCH(NTOUCH) = I2
                           JITCH(NTOUCH) = J2
                        ELSE
                           WRITE(*,*) 'DO2XADD: MAXTOUCH overflow'
                           CASTOP = .TRUE.
                        ENDIF
                     ENDIF
                   ENDIF
                 ENDIF
                END DO
                END DO
C Keep any unresolved remainder in the local residual reservoir
                  XRES(I,J) = LEFT
              ENDIF
            ENDIF
         ENDIF
      END DO
C---------------------------------------------------------------------
C Sparse apply/clear of staged deposits (duplicate touches allowed)
C---------------------------------------------------------------------
      DO K = 1, NTOUCH
         I2 = ITCH(K)
         J2 = JITCH(K)
         IF (XRESADD(I2,J2) .GT. TOL) THEN
            XRES(I2,J2) = XRES(I2,J2) + XRESADD(I2,J2)
            XRESADD(I2,J2) = 0.0
         ENDIF
      END DO
      RETURN
      END
C
C=====================================================================
C FUNCTION MASSBAL -- Total solute mass in GRID plus staged reservoirs
C  GRID contribution:
C       FL*CL + (1-FL)*CS    wt% per cell
C  XADD and XRES are stored as wt% * subcells so are divided
C  by ZUNITS to convert to wt% per cell before summing.
C=====================================================================
      DOUBLE PRECISION FUNCTION MASSBAL(GRID, XADD, XRES, XRND,
     &                                  XLEN, YLEN, SCRN, CASTOP)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN, SCRN
      INTEGER I, J
      REAL GRID(NXMAX, NYMAX, 6, 2)
      REAL XADD(NXMAX, NYMAX)
      REAL XRES(NXMAX, NYMAX)
      REAL XRND(NXMAX, NYMAX)
      INTEGER LCELLS, LIQUID, SOLID
      LOGICAL CASTOP
      REAL ZUNITS
      PARAMETER (LCELLS=1, LIQUID=2, SOLID=3)
      PARAMETER (ZUNITS=1000.0)
      DOUBLE PRECISION AMTOT, FL, CL, CS
      DOUBLE PRECISION XR, XA, XQ
      DOUBLE PRECISION XRSUM, XASUM, XQSUM, RESALL
      DOUBLE PRECISION RATIO, TINY
      PARAMETER (TINY = 1.0D-30)
C Mass conservation thresholds
      DOUBLE PRECISION XRELMX, XABSMX
      PARAMETER (XRELMX = 1.0D-3)
      PARAMETER (XABSMX = 1.0D-4)
      AMTOT = 0.0D0
      XRSUM = 0.0D0
      XASUM = 0.0D0
      XQSUM = 0.0D0
      DO J = 1, YLEN
      DO I = 1, XLEN
         FL = DBLE(GRID(I,J,LCELLS,SCRN))/DBLE(ZUNITS)
         IF (FL .LT. 0.0D0) FL = 0.0D0
         IF (FL .GT. 1.0D0) FL = 1.0D0
         CL = DBLE(GRID(I,J,LIQUID,SCRN))
         CS = DBLE(GRID(I,J,SOLID ,SCRN))
C Persistent reservoir and staged-release reservoir.
C Both are stored as wt%*subcells.
         XR = DBLE(XRES(I,J))/DBLE(ZUNITS)
         XA = DBLE(XADD(I,J))/DBLE(ZUNITS)
         XQ = DBLE(XRND(I,J))/DBLE(ZUNITS)
         AMTOT = AMTOT + FL*CL + (1.0D0-FL)*CS
         XRSUM = XRSUM + XR
         XASUM = XASUM + XA
         XQSUM = XQSUM + XQ
      END DO
      END DO
      RESALL = XRSUM + XASUM
C Reservoir guard.
      RATIO = RESALL/(DABS(AMTOT) + TINY)
      IF ( (DABS(RESALL) .GT. XABSMX) .AND.
     &     (DABS(RATIO) .GT. XRELMX) ) THEN
         WRITE(*,*) '  '
         WRITE(*,*) ' Reservoir inventory non-negligible: STOP'
         WRITE(*,*) '   Interpretation: solute is being retained in'
         WRITE(*,*) '   staged/rejected reservoirs instead of GRID.'
         WRITE(*,*) '  '
         WRITE(*,*) '  SCRN        = ', SCRN
         WRITE(*,*) '  XRES sum    = ', XRSUM
         WRITE(*,*) '  XADD sum    = ', XASUM
         WRITE(*,*) '  XRND sum    = ', XQSUM
         WRITE(*,*) '  AMTOT       = ', AMTOT
         WRITE(*,*) '  RES/AMTOT   = ', RATIO
         WRITE(*,*) '  Limits: rel=', XRELMX,' abs=', XABSMX
         CASTOP = .TRUE.
      ENDIF
      MASSBAL = AMTOT + XRSUM + XASUM + XQSUM
      RETURN
      END
C
C=====================================================================
C XCOMMIT -- commit staged updates from SNEW back to SCUR
C   Refactor Pass-3: extracted from CALOOP for readability.
C   Does NOT overwrite TCELL.
C=====================================================================
      SUBROUTINE XCOMMIT(GRID,SCUR,SNEW, IUPD,JUPD,NUPD)
      IMPLICIT NONE
      INTEGER SCUR, SNEW
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
C     INTEGER XLEN, YLEN   ! removed from args
      REAL GRID(NXMAX,NYMAX,6,2)
      INTEGER LCELLS, LIQUID, SOLID, PHASE, CELLDIR
      PARAMETER (LCELLS=1, LIQUID=2, SOLID=3, PHASE=4, CELLDIR=6)
      INTEGER NUPD, IUPD(NUPD), JUPD(NUPD)
      INTEGER K, I, J
      DO 10 K = 1, NUPD
         I = IUPD(K)
         J = JUPD(K)
         GRID(I,J,LCELLS ,SCUR) = GRID(I,J,LCELLS ,SNEW)
         GRID(I,J,LIQUID ,SCUR) = GRID(I,J,LIQUID ,SNEW)
         GRID(I,J,SOLID  ,SCUR) = GRID(I,J,SOLID  ,SNEW)
         GRID(I,J,PHASE  ,SCUR) = GRID(I,J,PHASE  ,SNEW)
         GRID(I,J,CELLDIR,SCUR) = GRID(I,J,CELLDIR,SNEW)
   10 CONTINUE
      RETURN
      END
C
C=====================================================================
C XKGTFRZCAP -- KGT-based interface closure (1-cell spherical cap)
C Computes interface advance and residual liquid from KGT tip
C kinematics using a one-cell spherical-cap geometry.
C Inputs:
C   VELMS   (m/s)   Tip velocity
C   DT      (s)     Timestep
C   DX      (m)     Cell size
C   RTIP    (m)     KGT tip radius
C   ZUNITS  (-)     Subcells per cell
C Outputs:
C   DFRZKGT   Frozen subcells per step
C   DFRZMS Residual liquid subcells in tip cell
C   PHI       Geometric advance factor (<= 1)
C   RHOTIP    RTIP/DX
C Geometry:
C   rho = RTIP/DX
C   h/DX = rho - sqrt(rho^2 - 1/4)
C   PHI  = 1 - h/DX = 1 - rho + sqrt(rho^2 - 1/4)
C   rho < 0.5  ->  PHI = 1, DFRZMS = 0
C=====================================================================
      SUBROUTINE XKGTFRZCAP(VELMS, DT, DX, RTIP, 
     &    DFRZKGT, DFRZMS, PHI, RHOTIP)
      IMPLICIT NONE
      REAL VELMS, DT, DX, RTIP, ZUNITS
      REAL DFRZKGT, DFRZMS, PHI, RHOTIP
      REAL ARG, HFRAC
      PARAMETER (ZUNITS=1000.0)
      DFRZKGT = 0.0
      DFRZMS = 0.0
      PHI = 1.0
      RHOTIP = 0.0
      IF (DX .LE. 0.0) RETURN
      IF (ZUNITS .LE. 0.0) RETURN
      RHOTIP = RTIP/DX
      ARG = RHOTIP*RHOTIP - 0.25
      IF (ARG .LT. 0.0) ARG = 0.0
      HFRAC = RHOTIP - SQRT(ARG)
      PHI = 1.0 - HFRAC
      IF (PHI .LT. 0.0) PHI = 0.0
      IF (PHI .GT. 1.0) PHI = 1.0
      DFRZKGT   = ZUNITS*(VELMS*DT/DX)*PHI
      DFRZMS = ZUNITS*HFRAC
      IF (DFRZMS .LT. 0.0) DFRZMS = 0.0
      IF (DFRZMS .GT. ZUNITS) DFRZMS = ZUNITS
      RETURN
      END
C
C ==============================================================
C XDELT0 -- Calculate Solutal Undercooling 
C ==============================================================
      REAL FUNCTION XDELT0(MLIQ, C0, KALPHA)
      IMPLICIT NONE
      REAL MLIQ, C0, KALPHA
      REAL CL0, TL0, TS0
      REAL TLIQBIN
      EXTERNAL TLIQBIN
      XDELT0 = 0.0
      IF (KALPHA .GT. 1.0E-8) THEN
         TL0 = TLIQBIN(C0,0.0,MLIQ,0.0)
         CL0 = C0/KALPHA
         TS0 = TLIQBIN(CL0,0.0,MLIQ,0.0)
         XDELT0 = TL0 - TS0
      ENDIF
      IF (XDELT0 .LT. 0.0) XDELT0 = 0.0
      RETURN
      END
C
C ============================================================
C XRIP -- Calculate Tip Radius with Thermal Correction 
C ============================================================
      REAL FUNCTION XRTIP(V, G, GIBT, DCOEFF, DELT0)
      IMPLICIT NONE
      REAL V, G, GIBT, DCOEFF, DELT0
      REAL SIGSTAR, LCHAR, RTIP0, BB
        PARAMETER (SIGSTAR = 0.0253156)
      LCHAR = GIBT/(SIGSTAR*DELT0)
      RTIP0 = SQRT(DCOEFF*LCHAR/V)    
C Thermal gradient correction factor BB
      BB = 1.0 + (G*GIBT/(V*DELT0*DELT0))
      IF (BB .LT. 1.0) BB = 1.0
      XRTIP = RTIP0/BB
      RETURN
      END
C
C=====================================================================
C  IVANTSOV(P) -- Ivantsov function
C     Iv(P) = P*exp(P)*E1(P)
C * P <  0.67 : E1 series (small x)
C * P >= 0.67 : continued fraction where exp(P)*E1(P) = 1/CF(P)
C               (avoids exp overflow; faster for larger P)
C Ref: see Ivantsov.4th example in J.L. Bezemer's 4tH compiler:
C      https://sourceforge.net/projects/forth-4th/
C=====================================================================
      REAL FUNCTION IVANTSOV(P)
      IMPLICIT NONE
      REAL P, X, EUL, VAL, CF, E1LOW, E1CF
      INTEGER NLOW, NHIGH
      PARAMETER (EUL = 0.57721566)
      PARAMETER (NLOW = 32, NHIGH = 40)
      VAL = 0.0
      X = P
      IF (X .GT. 0.0) THEN
         IF (X .LT. 0.67) THEN
C Low-x: E1(x) = -gamma - ln(x) - sum{n=1..N} [(-x)^n/(n*n!)]
            VAL = X*EXP(X)*E1LOW(X, NLOW, EUL)
         ELSE
C High-x: exp(x)*E1(x) = 1/CF(x) -> Iv(P) = x/CF(x)
            CF = E1CF(X, NHIGH)
            VAL = X/CF
         ENDIF
      ENDIF
      IF (VAL .LT. 0.0) VAL = 0.0
      IVANTSOV = VAL
      END
C
C=====================================================================
C E1LOW(x,N,gamma) -> E1(x) using the (-x)^n/(n*n!) power series
C=====================================================================
      REAL FUNCTION E1LOW(X, N, EUL)
      IMPLICIT NONE
      REAL X, EUL, SUM, FACT, TERM, POWX
      INTEGER N, I
      IF (X .LE. 0.0) THEN
         E1LOW = 1.0E38
      ELSE
         SUM = 0.0
         FACT = 1.0
         POWX = 1.0
         DO I = 1, N
            FACT = FACT*REAL(I)
            POWX = POWX*(-X)
            TERM = POWX/(REAL(I)*FACT)
            SUM = SUM+TERM
         ENDDO
         E1LOW = -EUL - ALOG(X) - SUM
      ENDIF
      END
C
C=====================================================================
C E1CF(x,N): Continued fraction CF(x) such that exp(x)*E1(x) = 1/CF(x)
C=====================================================================
      REAL FUNCTION E1CF(X, N)
      IMPLICIT NONE
      REAL X, CUR
      INTEGER N, I
      CUR = X + REAL(N + 1)
      DO I = N, 1, -1
         CUR = X + REAL(I)/(1.0+ REAL(I)/CUR)
      ENDDO
      E1CF = CUR
      END
C
C ===================================================================
C EMASCL --exponential moving average (EMA) 
C   XAVG(new) = XAVG(old) + alpha(XNEW - XAVG(old))
C     where alpha = DT/(TAU + DT)
C     TAU  - Time const (Larger, more smoothing & slower response)
C     DT   - Time step
C ===================================================================
      SUBROUTINE EMASCL(XNEW, TAU, DT, XAVG, XINIT)
      IMPLICIT NONE
      REAL XNEW, TAU, DT, XAVG, ALPHA
      LOGICAL XINIT
      IF (.NOT. XINIT) THEN
         XAVG = XNEW
         XINIT = .TRUE.
      ELSEIF (DT .GT. 0.0) THEN
         ALPHA = DT/(TAU + DT)
         XAVG = XAVG + ALPHA*(XNEW - XAVG)
      ENDIF
      RETURN
      END  
C
C ===================================================================
C XFINDBAND-- Interface-box adaptive scan window.
C This is safer for center/equiaxed growth than a single tip-centered
C band.  If no partial interface cells exist yet, scan the full domain
C to avoid startup deadlock.  Otherwise scan the bounding box of all
C active interface cells plus a generous halo.
C ====================================================================
      SUBROUTINE XFINDBAND(GRID,XLEN,YLEN,SCUR,
     &     BOX2,IMIN,IMAX,JMIN,JMAX)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX, LCELLS
      INTEGER XLEN, YLEN, SCUR, BOX2
      INTEGER I, J,  IMIN, IMAX, JMIN, JMAX, JMXIF, JMNIF, NIF
      INTEGER IMNIF, IMXIF, HALO
      REAL S1, ZUNITS
      PARAMETER (ZUNITS = 1000.0)
      PARAMETER (NXMAX=1000, NYMAX=1000, LCELLS=1)
      REAL GRID(NXMAX,NYMAX,6,2)
      IMNIF = XLEN
      IMXIF = 1
      JMNIF = YLEN
      JMXIF = 1
      NIF = 0
      DO J = BOX2+1, YLEN-BOX2
      DO I = 1, XLEN
        S1 = GRID(I,J,LCELLS,SCUR)
        IF (S1 .GT. 0.0 .AND. S1 .LT. ZUNITS) THEN
          NIF = NIF + 1
          IF (I .LT. IMNIF) IMNIF = I
          IF (I .GT. IMXIF) IMXIF = I
          IF (J .LT. JMNIF) JMNIF = J
          IF (J .GT. JMXIF) JMXIF = J
        ENDIF
      ENDDO
      ENDDO
      IF (NIF .LE. 0) THEN
        IMIN = 1
        IMAX = XLEN
        JMIN = BOX2 + 1
        JMAX = YLEN - BOX2
      ELSE
        HALO = 80
        IMIN = IMNIF - HALO
        IMAX = IMXIF + HALO
        JMIN = JMNIF - HALO
        JMAX = JMXIF + HALO
      ENDIF
C Full-X fallback if the interface-box X band crosses a side boundary.
      IF (IMIN .LT. 1 .OR. IMAX .GT. XLEN) THEN
        IMIN = 1
        IMAX = XLEN
      ELSE
        IMIN = MAX(1, IMIN)
        IMAX = MIN(XLEN, IMAX)
      ENDIF
      JMIN = MAX(BOX2+1, JMIN)
      JMAX = MIN(YLEN-BOX2, JMAX)
C Safety guards
      IF (JMIN .LT. BOX2+1) JMIN = BOX2 + 1
      IF (JMAX .GT. YLEN-BOX2) JMAX = YLEN - BOX2
      IF (JMIN .GT. JMAX) THEN
        JMIN = BOX2 + 1
        JMAX = YLEN - BOX2
      ENDIF
      IF (IMIN .GT. IMAX) THEN
        IMIN = 1
        IMAX = XLEN
      ENDIF
      RETURN
      END   
C
C=====================================================================
C XADDUPD -- Add one changed cell to the update list.
C            If the list is full, set overflow flag.
C=====================================================================
      SUBROUTINE XADDUPD(I,J,IUPD,JUPD,MAXUPD,NUPD,UPDOVER)
      IMPLICIT NONE
      INTEGER I,J,MAXUPD,NUPD
      INTEGER IUPD(MAXUPD), JUPD(MAXUPD)
      LOGICAL UPDOVER
      IF (NUPD .LT. MAXUPD) THEN
         NUPD = NUPD + 1
         IUPD(NUPD) = I
         JUPD(NUPD) = J
      ELSE
         UPDOVER = .TRUE.
      ENDIF
      RETURN
      END
C
C=====================================================================
C XINITSTAGE -- Initialize per-step staging/tip/update counters.
C=====================================================================
      SUBROUTINE XINITSTAGE(UPDOVER,NUPD,NGROW,
     &     DTBEST,TLBEST,TIPX,TIPY,HAVETIP)
      IMPLICIT NONE
      LOGICAL UPDOVER, HAVETIP
      INTEGER NUPD, NGROW, TIPX, TIPY
      REAL DTBEST, TLBEST
       UPDOVER = .FALSE.
       NUPD = 0
       NGROW = 0
       DTBEST = 1.0E30
       TLBEST = -1.0E30
       TIPX = 1
       TIPY = 1
       HAVETIP = .FALSE.
      RETURN
      END
C
C=====================================================================
C XLOCALKGT -- Compute local projected KGT curvature & freezing scale.
C  Uses the leading-tip thermal-gradient direction to project the
C  local growth orientation.  The projected velocity is then passed
C  to KGT.  If the local KGT calculation is not available, the routine
C  falls back to the input global CURVKGT and DFRZ values.
C
C   Inputs--
C   PTHETA   : local/preferred growth direction [rad]
C   GTXTIP   : x-component of tip thermal-gradient direction/field
C   GTYTIP   : y-component of tip thermal-gradient direction/field
C   GGTIP    : tip thermal-gradient magnitude [K/m]
C   VEL      : imposed velocity [cells/s]
C   GRAD     : fallback gradient [K/cell]
C   DX       : cell size [m]
C   DT       : timestep [s]
C   DFRZ     : fallback freezing increment [subcells]
C   CURVKGT  : fallback KGT curvature [1/m]
C
C   Outputs--
C   CURVKLOC : local KGT curvature [1/m]
C   DFRZLOC  : local freezing increment [subcells]
C=====================================================================
      SUBROUTINE XLOCALKGT(PTHETA,GTXTIP,GTYTIP,GGTIP,
     &     VEL,GRAD,DX,DT,
     &     GIBT,DCOEFF,MLIQ,KALPHA,C0,
     &     CURVKGT,DFRZ,
     &     CURVKLOC,DFRZLOC)
      IMPLICIT NONE
      REAL PTHETA, GTXTIP, GTYTIP, GGTIP
      REAL VEL, GRAD, DX, DT
      REAL GIBT, DCOEFF, MLIQ, KALPHA, C0
      REAL CURVKGT, DFRZ, ZUNITS
      REAL CURVKLOC, DFRZLOC
      REAL TX, TY, GM, PROJAX, VKGTLOC, GGLOC
      REAL DTKGTLOC, RTIPLOC, PELOC
      REAL IVRLOC, DELT0LOC
      REAL DFRZKGTLOC, DFRZMSLOC, PHILOC, RHOLOC
      PARAMETER (ZUNITS=1000.0)
C Safe default
      CURVKLOC = CURVKGT
      DFRZLOC  = DFRZ
C Normalize the local thermal-gradient direction.
      TX = GTXTIP
      TY = GTYTIP
      GM = SQRT(TX*TX + TY*TY)
      IF (GM .GT. 1.0E-12) THEN
        TX = TX/GM
        TY = TY/GM
      ELSE
        TX = 0.0
        TY = 1.0
      ENDIF
C Project local growth direction onto the tip-gradient direction.
      PROJAX = ABS(COS(PTHETA)*TX + SIN(PTHETA)*TY)
      IF (PROJAX .GT. 1.0) PROJAX = 1.0
      VKGTLOC = VEL*DX*PROJAX
C Use measured tip gradient if available; otherwise parameter fallback.
      GGLOC = GGTIP
      IF (GGLOC .LE. 1.0E-12) GGLOC = ABS(GRAD)/DX
C Local KGT state and spherical-cap freezing scale.
      IF (VKGTLOC .GT. 1.0E-12 .AND.
     &    GGLOC .GT. 1.0E-12) THEN
         CALL KGT(VKGTLOC,GGLOC,GIBT,DCOEFF,MLIQ,KALPHA,C0,
     &        DTKGTLOC,RTIPLOC,CURVKLOC,PELOC,IVRLOC,DELT0LOC)
         IF (RTIPLOC .GT. 0.0) THEN
            CALL XKGTFRZCAP(VKGTLOC,DT,DX,RTIPLOC,
     &           DFRZKGTLOC,DFRZMSLOC,PHILOC,RHOLOC)
            DFRZLOC = DFRZMSLOC
            IF (DFRZLOC .LT. 1.0) DFRZLOC = 1.0
            IF (DFRZLOC .GT. 0.8*ZUNITS) DFRZLOC = 0.8*ZUNITS
         ENDIF
      ENDIF
      RETURN
      END
C=====================================================================
C XSTAGE -- Scan active interface band and stage cell-growth updates.
C   Performs interface scan for primary growth and eutectic freezing.
C   Scans only the band selected by XFINDBAND() and then tests each
C   candidate cell for eutectic  or primary solidification.
C For primary growth:
C  - finds neighboring solid/interface support with NNN(),
C  - selects the local growth orientation PTHETA,
C  - obtains projected local KGT curvature scale from XLOCALKGT(),
C  - evaluates local liquidus, curvature undercooling, and CLSTAR,
C  - calls GROW() to update liquid/solid state and mass reservoirs,
C  - records changed cells through XADDUPD().
C Routines called
C   XFINDBAND(): adaptive scan-window selection.
C   XINITSTAGE(): counter/tip-state initialization per time step.
C   XLOCALKGT(): projected local KGT / DFRZ estimate.
C   XADDUPD(): update-list overflow bookkeeping.
C   GROW(): mass-sensitive update path.
C Output:
C ------------
C   NUPD, IUPD, JUPD : list of cells changed this step
C   NGROW            : number of primary-growth updates
C   GINCFRZ          : total frozen-subcell increment
C   GINCPRI          : primary-growth event count
C   GINCEUT          : eutectic-freezing event count
C   DTBEST, TLBEST   : leading local undercooling/tip indicators
C   TIPX, TIPY       : current best tip cell
C   HAVETIP          : true if a growing tip candidate was found
C   UPDOVER          : true if update-list capacity was exceeded
C=====================================================================
      SUBROUTINE XSTAGE(GRID,XADD,XRES,XRND,QLAT,XLEN,YLEN,
     & SCUR,SNEW,
     & BOX2,TOL,
     & TEUT,TL0,MLIQ,KPART,C0,
     & GIBT,DD,DFRZ,CURVKGT,LHEAT,CPHEAT,
     & VEL,GRAD,DX,DCOEFF,DT,
     & GINCFRZ,GINCPRI,GINCEUT,
     & DTBEST,TLBEST,TIPX,TIPY,TTIPREF,HAVETIP,
     & TLX, TLY, TLGX, TLGY,
     & GTXTIP, GTYTIP, GGTIP,
     & TLCLTIP, TLLDIF, TLWID, TLLEN,
     & ITIPLOAD,
     & NSTAGEBOX, ISTAGEBOX, JSTAGEBOX, IFRZSTAGE,
     & IUPD,JUPD,MAXUPD,NUPD,NGROW,UPDOVER, LENBOX, IBCTYPE)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN, SCUR, SNEW, BOX2
      INTEGER LCELLS, LIQUID, TCELL, CELLDIR
      PARAMETER (LCELLS=1, LIQUID=2, TCELL=5, CELLDIR=6)
      REAL GRID(NXMAX,NYMAX,6,2)
      REAL XADD(NXMAX,NYMAX), XRES(NXMAX,NYMAX)
      REAL XRND(NXMAX,NYMAX)
      REAL QLAT(NXMAX,NYMAX)
      REAL LHEAT, CPHEAT
      REAL ZUNITS, TOL
      REAL TEUT, TL0, MLIQ, KPART, C0
      REAL GIBT, DD, DFRZ, CURVKGT
      REAL VEL, GRAD, DX, DCOEFF, DT
      REAL TTIPREF
      INTEGER IBCTYPE
      INTEGER GINCFRZ, GINCPRI, GINCEUT
      REAL DTBEST, TLBEST
      INTEGER TIPX, TIPY
      LOGICAL HAVETIP
      INTEGER MAXUPD, NUPD, NGROW, NSTAGEBOX
      INTEGER ISTAGEBOX(64), JSTAGEBOX(64), IFRZSTAGE(64)
      INTEGER IUPD(MAXUPD), JUPD(MAXUPD)
      LOGICAL UPDOVER
      INTEGER I, J, NC, IP, JP
      INTEGER IMIN, IMAX, JMIN, JMAX
      REAL S1, TL, COMP, LNN
      REAL THETA, CURVUSE, DTCURV
      REAL TEQ, TEQ0, DELT, CLSTAR, CSCELL, AMOUNT
      REAL PTHETA, THUSE, STRENGTH
      REAL ANGDIFF, CANG
      REAL WANG, WTHERM, DTBACK, FRAC
      REAL TLX, TLY, TLGX, TLGY
      REAL GTXTIP, GTYTIP, GGTIP
      REAL TLCLTIP, TLLDIF, TLWID, TLLEN
      INTEGER ITIPLOAD
      REAL PI, HALFPI, FMIN, DT0, DT1
      REAL CURVKLOC, DFRZLOC
      INTEGER LENBOX
      PARAMETER (ZUNITS=1000.0, PI=3.1415927)
      PARAMETER (HALFPI=1.5707963)
      PARAMETER (FMIN=0.15)
      PARAMETER (DT0=0.04)
      PARAMETER (DT1=0.16)
      REAL TLIQBIN, CLIQBIN, CSOLBIN, ALONG, TCURV, THETA4
      EXTERNAL TLIQBIN, CLIQBIN, CSOLBIN, ALONG, TCURV, THETA4
C Select active scan band and initialize per-step staging state.
      CALL XFINDBAND(GRID,XLEN,YLEN,SCUR,BOX2,IMIN,IMAX,JMIN,JMAX)
      CALL XINITSTAGE(UPDOVER,NUPD,NGROW,
     &                            DTBEST,TLBEST,TIPX,TIPY,HAVETIP)
C Scan selected band/window.
      DO J = JMIN, JMAX
      DO I = IMIN, IMAX
         S1 = GRID(I,J,LCELLS,SCUR)
         TL = GRID(I,J,TCELL ,SCUR)
         IF (TL .LE. TEUT .AND. S1 .GT. 0.5) THEN
          CALL EUTFREEZE(GRID,I,J,SCUR,SNEW,TEUT,GINCEUT)
          CALL XADDUPD(I,J,IUPD,JUPD,MAXUPD,NUPD,UPDOVER)
         ELSE IF (S1 .GT. TOL) THEN
          CALL NNN(GRID,I,J,XLEN,YLEN,SCUR,NC,LNN,IP,JP, IBCTYPE)
          IF (NC .GT. 0) THEN
             PTHETA = GRID(I,J,CELLDIR,SCUR)
             IF (S1 .GT. ZUNITS-0.5 .AND. IP .NE. 0) THEN
                PTHETA = GRID(IP,JP,CELLDIR,SCUR)
             ENDIF
C Local projected KGT state for this candidate direction.
             CALL XLOCALKGT(PTHETA,GTXTIP,GTYTIP,GGTIP,
     &            VEL,GRAD,DX,DT, GIBT,DCOEFF,MLIQ,KPART,C0,
     &            CURVKGT,DFRZ, CURVKLOC,DFRZLOC)
C Local liquidus precheck before curvature/interface work.
             COMP = GRID(I,J,LIQUID,SCUR)
             TEQ0 = TLIQBIN(COMP,TL0,MLIQ,C0)
             IF (TEQ0 .GT. TL) THEN
C Interface orientation and curvature weighting.
               CALL BOX(GRID,XLEN,YLEN,SCUR,I,J,LENBOX,THETA,IBCTYPE)
               ANGDIFF = ABS(THETA - PTHETA)
               IF (ANGDIFF .GT. HALFPI) ANGDIFF = PI - ANGDIFF
               CANG = COS(ANGDIFF)
               WANG = CANG*CANG
               DTBACK = TTIPREF - TL
               IF (DTBACK .LE. DT0) THEN
                  WTHERM = 1.0
               ELSE IF (DTBACK .GE. DT1) THEN
                  WTHERM = 0.0
               ELSE
                  FRAC = (DTBACK - DT0)/(DT1 - DT0)
                  WTHERM = 1.0 - FRAC
               ENDIF
               CURVUSE = CURVKLOC*(FMIN + (1.0-FMIN)*WTHERM*WANG)
               DTCURV = TCURV(GIBT, THETA, PTHETA, CURVUSE, DD)
               TEQ = TEQ0 - DTCURV
               DELT = TEQ - TL
               CLSTAR = CLIQBIN(TL + DTCURV,TL0,MLIQ,C0)
               IF (DELT .GT. 0.0) THEN
C Primary growth update.  GROW owns mass-sensitive state changes.
                 CALL TIPUPDATE(I,J,DELT,TL,DTBEST,TLBEST,TIPX,TIPY)
                 HAVETIP = .TRUE.
                 CSCELL = CSOLBIN(COMP,TL,TL0,MLIQ,KPART,C0)
                 AMOUNT = 0.0
                 IF (S1 .GT. 1.0) AMOUNT = (COMP - CSCELL)/S1
                 IF (AMOUNT .LT. 0.0) AMOUNT = 0.0
                 THUSE = THETA4(THETA, PTHETA)
                 STRENGTH = ALONG(THUSE, PTHETA, DD)
                 CALL GROW(GRID,I,J,XLEN,YLEN,SCUR,SNEW,KPART,
     &               CLSTAR,COMP,S1,NC,AMOUNT,STRENGTH,DFRZLOC,
     &               CSCELL,GINCFRZ,GINCPRI,PTHETA,XADD,XRES,
     &               XRND,QLAT,MLIQ,TL0,C0,LHEAT,CPHEAT,
     &               TLX, TLY, TLGX, TLGY,
     &               TLCLTIP, TLLDIF, TLWID, TLLEN,
     &               ITIPLOAD,NSTAGEBOX,
     &               ISTAGEBOX, JSTAGEBOX, IFRZSTAGE, IBCTYPE)
                 NGROW = NGROW + 1
                 CALL XADDUPD(I,J,IUPD,JUPD,MAXUPD,NUPD,UPDOVER)
               ENDIF
             ENDIF
          ENDIF
         ENDIF
      ENDDO
      ENDDO
      RETURN
      END
C
C=====================================================================
C DTINFC -- Interface density count
C=====================================================================
      SUBROUTINE DTINFC(GRID,XLEN,YLEN,SCRN)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      REAL GRID(NXMAX,NYMAX,6,2)
      INTEGER XLEN, YLEN, SCRN, LCELLS
      INTEGER NIFACE
      REAL ZUNITS, TOL
      PARAMETER (LCELLS = 1, ZUNITS=1000.0, TOL=1.0E-6)
      INTEGER I, J
      REAL S1
      NIFACE = 0
      DO J = 1, YLEN
      DO I = 1, XLEN
         S1 = GRID(I,J,LCELLS,SCRN)
         IF (S1 .GT. TOL .AND. S1 .LT. ZUNITS-TOL) THEN
            NIFACE = NIFACE + 1
         ENDIF
      ENDDO
      ENDDO
      RETURN
      END
C
C ====================================================================
C NNN() Find number of solid cardinal & diagonal neighbors for a cell
C  NN   = count of fully solid cardinal neighbors, or if none exist,
C          then 100 times the count of solid diagonal neighbors.
C  LNN  = number of fully liquid cardinal neighbors.
C IP,JP = coordinates of the FIRST solid neighbor encountered.
C ====================================================================
      SUBROUTINE NNN(GRID,I,J,XLEN,YLEN,SCRN,NN,LNN,IP,JP, IBCTYPE)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER I, J, NN, SCRN, XLEN, YLEN
      INTEGER LCELLS, IP, JP, IBCTYPE
      REAL GRID(NXMAX, NYMAX, 6, 2)
      REAL LNN, LIQALL, SLDALL
      REAL S
      INTEGER I2, J2
      PARAMETER (LCELLS=1)
      PARAMETER (LIQALL=999.5, SLDALL=0.5)
        NN = 0
        LNN = 0.0
        IP = 0
        JP = 0
C Cardinals: (I,J+1), (I,J-1), (I+1,J), (I-1,J)
      CALL GETBCIJ(I,J+1,XLEN,YLEN,IBCTYPE,I2,J2)
      S = GRID(I2,J2,LCELLS,SCRN)
      IF (S .GT. LIQALL) LNN = LNN + 1.0
      IF (S .LT. SLDALL) THEN
         NN = NN + 1
         IF (IP .EQ. 0) THEN
            IP = I2
            JP = J2
         ENDIF
      ENDIF
      CALL GETBCIJ(I,J-1,XLEN,YLEN,IBCTYPE,I2,J2)
      S = GRID(I2,J2,LCELLS,SCRN)
      IF (S .GT. LIQALL) LNN = LNN + 1.0
      IF (S .LT. SLDALL) THEN
         NN = NN + 1
         IF (IP .EQ. 0) THEN
            IP = I2
            JP = J2
         ENDIF
      ENDIF
      CALL GETBCIJ(I+1,J,XLEN,YLEN,IBCTYPE,I2,J2)
      S = GRID(I2,J2,LCELLS,SCRN)
      IF (S .GT. LIQALL) LNN = LNN + 1.0
      IF (S .LT. SLDALL) THEN
         NN = NN + 1
         IF (IP .EQ. 0) THEN
            IP = I2
            JP = J2
         ENDIF
      ENDIF
      CALL GETBCIJ(I-1,J,XLEN,YLEN,IBCTYPE,I2,J2)
      S = GRID(I2,J2,LCELLS,SCRN)
      IF (S .GT. LIQALL) LNN = LNN + 1.0
      IF (S .LT. SLDALL) THEN
         NN = NN + 1
         IF (IP .EQ. 0) THEN
            IP = I2
            JP = J2
         ENDIF
      ENDIF
C Diagonals only if no cardinal solids
      IF (NN .EQ. 0) THEN
         CALL GETBCIJ(I+1,J+1,XLEN,YLEN,IBCTYPE,I2,J2)
         S = GRID(I2,J2,LCELLS,SCRN)
         IF (S .LT. SLDALL) THEN
            NN = NN + 100
            IF (IP .EQ. 0) THEN
               IP = I2
               JP = J2
            ENDIF
         ENDIF
         CALL GETBCIJ(I+1,J-1,XLEN,YLEN,IBCTYPE,I2,J2)
         S = GRID(I2,J2,LCELLS,SCRN)
         IF (S .LT. SLDALL) THEN
            NN = NN + 100
            IF (IP .EQ. 0) THEN
               IP = I2
               JP = J2
            ENDIF
         ENDIF
         CALL GETBCIJ(I-1,J+1,XLEN,YLEN,IBCTYPE,I2,J2)
         S = GRID(I2,J2,LCELLS,SCRN)
         IF (S .LT. SLDALL) THEN
            NN = NN + 100
            IF (IP .EQ. 0) THEN
               IP = I2
               JP = J2
            ENDIF
         ENDIF
         CALL GETBCIJ(I-1,J-1,XLEN,YLEN,IBCTYPE,I2,J2)
         S = GRID(I2,J2,LCELLS,SCRN)
         IF (S .LT. SLDALL) THEN
            NN = NN + 100
            IF (IP .EQ. 0) THEN
               IP = I2
               JP = J2
            ENDIF
         ENDIF
      ENDIF
      RETURN
      END
C
C=====================================================================
C TIPUPDATE -- Select tip by: minimum undercooling (DELT) and
C   maximum temperature (TL). Called only for active growth 
C   sites (DELT > 0)
C=====================================================================
      SUBROUTINE TIPUPDATE(I,J,DELT,TL,DTBEST,TLBEST,TIPX,TIPY)
      IMPLICIT NONE
      INTEGER I,J,TIPX,TIPY
      REAL DELT,TL
      REAL DTBEST,TLBEST
      REAL DTTOL
      PARAMETER (DTTOL = 1.0E-6)
C Primary: minimum undercooling
      IF (DELT .LT. DTBEST - DTTOL) THEN
         DTBEST = DELT
         TLBEST = TL
         TIPX = I
         TIPY = J
C Also check maximum temperature
      ELSE IF (ABS(DELT - DTBEST) .LE. DTTOL) THEN
         IF (TL .GT. TLBEST) THEN
            TLBEST = TL
            TIPX = I
            TIPY = J
         ENDIF
      ENDIF
      RETURN
      END      
C
C=====================================================================
C THETA4 -- fold THETA to nearest 4-fold equivalent about PTHETA
C Ensures DELTA = THETA-PTHETA lies in [-PI/4, +PI/4).
C This removes branch/snapping artifacts when THETA comes from BOX().
C=====================================================================
      REAL FUNCTION THETA4(THETA, PTHETA)
      IMPLICIT NONE
      REAL THETA, PTHETA
      REAL PI2, PI4, D
      PARAMETER (PI2 = 2.0*ATAN(1.0))     
      PARAMETER (PI4 = ATAN(1.0))         
      D = THETA - PTHETA
C Fold D into [-PI/4, PI/4) with period PI/2
      D = MOD(D + PI4, PI2)
      IF (D .LT. 0.0) D = D + PI2
      D = D - PI4
      THETA4 = PTHETA + D
      RETURN
      END
C         
C ==============================================================
C ALONG -- 4-fold anisotropy weighting 
C ==============================================================
      REAL FUNCTION ALONG(THETA, PTHETA, DDALONG)
      IMPLICIT NONE
      REAL THETA, PTHETA, DDALONG
      REAL TH, PH
      REAL NORMANG, ANGDIF
      TH = NORMANG(THETA)
      PH = NORMANG(PTHETA)
      ALONG = 1.0 + 15.0*DDALONG*COS(4.0*ANGDIF(TH,PH))
      RETURN
      END
C
C=====================================================================
C GROW -- Freeze one interface cell using local partitioning,
C         integer inventory bookkeeping, and staged solute rejection.
C   - Reads cell state from SCUR; neighbor queries also use SCUR
C   - Writes updated cell state to SNEW
C   - Stages rejected solute outward through EJECT9/XADD
C   - Stores any local ejection remainder in XRES(I,J)
C=====================================================================
      SUBROUTINE GROW(GRID, I, J, XLEN, YLEN, SCUR, SNEW, KPART,
     &      CLSTAR, COMP, S1, NC, AMOUNT, STRENGTH, DFRZ, CSCELL,
     &      GINCFRZ, GINCPRI, PTHETA, XADD, XRES, XRND, QLAT,
     &      MLIQ, TL0,
     &      C0, LHEAT, CPHEAT,
     &      TLX, TLY, TLGX, TLGY,
     &      TLCLTIP, TLLDIF, TLWID, TLLEN,
     &      ITIPLOAD,
     &      NSTAGEBOX, ISTAGEBOX, JSTAGEBOX, IFRZSTAGE, IBCTYPE)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN, SCUR, SNEW, I, J, NC
      INTEGER GINCFRZ, GINCPRI
      INTEGER FRZCELL
      REAL GRID(NXMAX, NYMAX, 6, 2)
      REAL QLAT(NXMAX, NYMAX)
      REAL LHEAT, CPHEAT
      REAL DTLAT      
      REAL KPART, CLSTAR, AMOUNT, STRENGTH
      REAL DFRZ, CSCELL, PTHETA
      REAL COMP, S1
      REAL XADD(NXMAX, NYMAX), XRES(NXMAX, NYMAX)
      REAL XRND(NXMAX, NYMAX)
      REAL MLIQ, TL0, C0
      REAL CSOLBIN
      EXTERNAL CSOLBIN
      REAL ZUNITS, TOL, EPS1
      PARAMETER (ZUNITS=1000.0, TOL=1.0E-6, EPS1=1.0)
      INTEGER LCELLS, LIQUID, SOLID, PHASE, CELLDIR
      PARAMETER (LCELLS=1, LIQUID=2, SOLID=3, PHASE=4, CELLDIR=6)
      REAL PRIMARY, INTERFACE
      PARAMETER (PRIMARY=1.0, INTERFACE=4.0)
      REAL A2, CSINIT
      REAL STR6, H
C Integer inventory carried in subcell*wt% quanta
      INTEGER S1I, SOLI
      INTEGER MTOT, ML, MSI, MS
      INTEGER ICS, ICL, IREJ
      REAL LEFT
      DOUBLE PRECISION MREAL, MERR
      REAL TLX, TLY, TLGX, TLGY
      REAL TLCLTIP, TLLDIF, TLWID, TLLEN
      INTEGER ITIPLOAD
      INTEGER NSTAGEBOX, ISTAGEBOX(64), JSTAGEBOX(64), IFRZSTAGE(64)
      INTEGER IFORCE, KSTG, IBCTYPE, KDO
      LOGICAL FOUND
C --
C (1) Build integer cell inventory from the current SCUR state
      FRZCELL = 0
      CSINIT = GRID(I,J,SOLID,SCUR)
C Snap liquid subcells to an integer count in [0,ZUNITS]
      S1I = NINT(S1)
      IF (S1I .LT. 0) S1I = 0
      IF (S1I .GT. NINT(ZUNITS)) S1I = NINT(ZUNITS)
      S1 = FLOAT(S1I)
C Total solute inventory in this cell (liquid + solid).
C MTOT is the integer inventory used by GROW.  The conversion from
C real-valued grid state to integer inventory creates a signed residual.
C Store that residual in XRND, not XRES.  XRND is included in MASSBAL
C but is never processed by CATRANSFER.
      MREAL = DBLE(FLOAT(S1I))*DBLE(COMP)
     &      + DBLE(FLOAT(NINT(ZUNITS)-S1I))*DBLE(CSINIT)
      MTOT = NINT(SNGL(MREAL))
      IF (MTOT .LT. 0) MTOT = 0
      MERR = DBLE(MTOT) - MREAL
      XRND(I,J) = XRND(I,J) - SNGL(MERR)
C Liquid solute inventory
      ML = NINT( FLOAT(S1I)*COMP )
      IF (ML .LT. 0) ML = 0
      IF (ML .GT. MTOT) ML = MTOT
C Solid solute inventory carried locally in the cell
      MSI = MTOT - ML
      IF (MSI .LT. 0) MSI = 0
C Apply any forced staged freezing assigned to this cell
      IFORCE = 0
      FOUND = .FALSE.
      KSTG = 1
      DO WHILE (KSTG .LE. NSTAGEBOX .AND. .NOT. FOUND)
        IF (ISTAGEBOX(KSTG) .EQ. I .AND.
     &     JSTAGEBOX(KSTG) .EQ. J) THEN
           IFORCE = IFRZSTAGE(KSTG)
           FOUND = .TRUE.
        ENDIF
        KSTG = KSTG + 1
      ENDDO
C Forced freezing step:
C - above DFRZ, partition remains within the cell
C - below DFRZ, rejected solute is ejected outward
      KDO = 0
      DO WHILE (KDO .LT. IFORCE .AND. S1I .GT. 0)
         KDO = KDO + 1
         COMP = ML/FLOAT(S1I)
         CSCELL = CSOLBIN(COMP,0.0,TL0,MLIQ,KPART,C0)
         IF (S1 .GT. (DFRZ + TOL)) THEN
            ICS = NINT(CSCELL)
            IF (ICS .LT. 0) ICS = 0
            IF (ICS .GT. ML) ICS = ML
            ML = ML - ICS
            IF (ML .LT. 0) ML = 0
            MSI = MSI + ICS
         ELSE
            ICS = NINT(CSCELL)
            IF (ICS .LT. 0) ICS = 0
            ICL = (ML + S1I/2)/S1I
            IF (ICL .LT. 0) ICL = 0
            IF (ICL .GT. ML) ICL = ML
            IF (ICS .GT. ICL) ICS = ICL
            IREJ = ICL - ICS
            IF (IREJ .LT. 0) IREJ = 0
            ML = ML - ICL
            IF (ML .LT. 0) ML = 0
            MSI = MSI + ICS
            IF (IREJ .GT. 0) THEN
               CALL EJECT9(I, J, GRID, XLEN, YLEN, SCUR,
     &              FLOAT(IREJ),
     &              XADD, MLIQ, TL0, C0,
     &              TLX, TLY, TLGX, TLGY,
     &              TLCLTIP, TLLDIF, TLWID, TLLEN,
     &              ITIPLOAD, LEFT, IBCTYPE)
               IF (LEFT .GT. 0.0) XRES(I,J) = XRES(I,J) + LEFT
            ENDIF
         ENDIF
         S1I = S1I - 1
         IF (S1I .LT. 0) S1I = 0
         S1 = FLOAT(S1I)
         FRZCELL = FRZCELL + 1
      ENDDO
C Keep liquid composition consistent with the remaining liquid fraction
      IF (S1I .GT. 0) THEN
         COMP = ML/FLOAT(S1I)
      ELSE
         COMP = 0.0
      ENDIF
C --
C (2) Local partitioning while liquid fraction remains above DFRZ
C Driving supersaturation for continued freezing
      A2 = CLSTAR - COMP
      IF (A2 .LT. 0.0) A2 = 0.0
C Directional weighting for local interface advance.
C Prefer nearest-neighbor progression, with limited diagonal
C contribution when aligned with the imposed growth direction.
      STR6 = 1.0
      IF (NC .GT. 0) THEN
         H = 0.5*STRENGTH
         STR6 = H*H
         STR6 = STR6*STR6
         STR6 = STR6*(H*H)
         A2 = A2*STR6
      ENDIF
C Local partition hurdle for continued in-cell freezing
      CSCELL = CSOLBIN(COMP,0.0,TL0,MLIQ,KPART,C0)
      IF (S1 .GT. 1.0) THEN
        AMOUNT = (COMP - CSCELL)/S1
      ELSE
        AMOUNT = 0.0
      ENDIF
      IF (AMOUNT .LT. 0.0) AMOUNT = 0.0
C Main in-cell freezing loop while liquid remains above DFRZ
      DO WHILE (NC.GT.0 .AND. S1I.GT.1 .AND. A2*STRENGTH.GT.AMOUNT)
         COMP = ML/FLOAT(S1I)
         CSCELL = CSOLBIN(COMP,0.0,TL0,MLIQ,KPART,C0)
         ICS = NINT(CSCELL)
         IF (ICS .LT. 0) ICS = 0
         IF (ICS .GT. ML) ICS = ML
C Partition within the cell while S1 > DFRZ:
C the solid takes ICS and the remaining solute stays local
         ML = ML - ICS
         IF (ML .LT. 0) ML = 0
         MSI = MSI + ICS
         S1I = S1I - 1
         IF (S1I .LT. 0) S1I = 0
         S1 = FLOAT(S1I)
         FRZCELL = FRZCELL + 1
         IF (S1I .GT. 0) THEN
          COMP = ML/FLOAT(S1I)
         ELSE
           COMP = 0.0
         ENDIF
         CSCELL = CSOLBIN(COMP,0.0,TL0,MLIQ,KPART,C0)
C Update local partition hurdle and driving supersaturation
         IF (S1 .GT. 1.0) THEN
           AMOUNT = (COMP - CSCELL)/S1
         ELSE
           AMOUNT = 0.0
         ENDIF
         IF (AMOUNT .LT. 0.0) AMOUNT = 0.0
         A2 = CLSTAR - COMP
         IF (A2 .LT. 0.0) A2 = 0.0
         A2 = A2*STR6
C --
C (3) Final freeze-out with staged solute ejection below DFRZ
C Once S1 falls below DFRZ, remaining liquid subcells are frozen
C out one by one. Each subcell partitions locally, and the rejected
C solute is staged outward through EJECT9 into XADD.
C Any local capacity shortfall is stored in XRES for later drainage.
         IF (S1 .LE. (DFRZ + TOL)) THEN
            DO WHILE (S1I .GT. 0)
               COMP = ML/FLOAT(S1I)
               CSCELL = CSOLBIN(COMP,0.0,TL0,MLIQ,KPART,C0)
               ICS = NINT(CSCELL)
               IF (ICS .LT. 0) ICS = 0
C Integer solute content of one remaining liquid subcell
               ICL = (ML + S1I/2)/S1I
               IF (ICL .LT. 0) ICL = 0
               IF (ICL .GT. ML) ICL = ML
               IF (ICS .GT. ICL) ICS = ICL
               IREJ = ICL - ICS
               IF (IREJ .LT. 0) IREJ = 0
C Remove one liquid subcell, add partitioned solute to solid
               ML = ML - ICL
               IF (ML .LT. 0) ML = 0
               MSI = MSI + ICS
               S1I = S1I - 1
               IF (S1I .LT. 0) S1I = 0
               S1 = FLOAT(S1I)
               FRZCELL = FRZCELL + 1
C Stage rejected solute into neighboring liquid cells
               IF (IREJ .GT. 0) THEN
                 CALL EJECT9(I, J, GRID, XLEN, YLEN, SCUR,
     &                FLOAT(IREJ),
     &                XADD, MLIQ, TL0, C0,
     &                TLX, TLY, TLGX, TLGY,
     &                TLCLTIP, TLLDIF, TLWID, TLLEN,
     &                ITIPLOAD, LEFT, IBCTYPE)
C Capacity shortfall remains in the local reservoir
                 IF (LEFT .GT. 0.0) XRES(I,J) = XRES(I,J) + LEFT
               ENDIF
            END DO
C Cell fully captured; outer loop terminates through NC = 0            
            S1I = 0
            S1 = 0.0
            COMP = 0.0
            GINCPRI = GINCPRI + 1
            NC = 0
         ELSE
            NC = NC - 1
         ENDIF
      END DO
C --
C Accumulate latent heat as a temperature increment 
      IF (FRZCELL .GT. 0 .AND. CPHEAT .GT. 0.0) THEN
         DTLAT = (LHEAT/CPHEAT)*FLOAT(FRZCELL)/ZUNITS
         QLAT(I,J) = QLAT(I,J) + DTLAT
      ENDIF
C --      
C (4) Write updated cell state to SNEW
      GINCFRZ = GINCFRZ + FRZCELL
      MS = MSI
      IF (MS .LT. 0) MS = 0
      MTOT = ML + MS
      IF (MTOT .LT. 0) MTOT = 0
      SOLI = NINT(ZUNITS) - S1I
      IF (SOLI .GT. 0) THEN
          GRID(I,J,SOLID,SNEW) = FLOAT(MS)/FLOAT(SOLI)
      ELSE
          GRID(I,J,SOLID,SNEW) = FLOAT(MTOT)/ZUNITS
      ENDIF
      GRID(I,J,LCELLS,SNEW) = FLOAT(S1I)
      GRID(I,J,LIQUID,SNEW) = COMP
C Interface cells retain liquid content and direction;
C fully captured cells are written as primary solid
      IF (S1I .GE. NINT(EPS1)) THEN
          GRID(I,J,PHASE,SNEW) = INTERFACE
          GRID(I,J,CELLDIR,SNEW) = PTHETA
      ELSE
          GRID(I,J,PHASE,SNEW) = PRIMARY
          GRID(I,J,CELLDIR,SNEW) = PTHETA
          GRID(I,J,LCELLS,SNEW) = 0.0
          GRID(I,J,LIQUID,SNEW) = 0.0
      ENDIF
      RETURN
      END
C
C=====================================================================
C EUTFREEZE -- Simple eutectic freeze via a thermal trigger.
C   If TCELL <= TEUT and cell contains liquid, convert full cell
C   (S+L) to solid at cell-average composition.
C=====================================================================
      SUBROUTINE EUTFREEZE(GRID,I,J,SCUR,SNEW,TEUT,GINCEUT)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER I,J,SCUR,SNEW
      REAL GRID(NXMAX, NYMAX, 6, 2)
      REAL TEUT, ZUNITS
      INTEGER GINCEUT
      INTEGER LCELLS, LIQUID, SOLID, PHASE, TCELL, CELLDIR
      PARAMETER (LCELLS=1, LIQUID=2, SOLID=3, PHASE=4, TCELL=5,
     &           CELLDIR=6)
      REAL SECONDARY
      PARAMETER (ZUNITS=1000.0, SECONDARY=2.0)
      REAL TL, CL, CS
      INTEGER S1I, SOLI, MTOT
C Read from SCUR
      TL = GRID(I,J,TCELL,SCUR)
      IF (TL .GT. TEUT) RETURN
      S1I = NINT(GRID(I,J,LCELLS,SCUR))
      IF (S1I .LE. 0) RETURN
C Check S1I to [0..ZUNITS]
      IF (S1I .LT. 0) S1I = 0
      IF (S1I .GT. NINT(ZUNITS)) S1I = NINT(ZUNITS)
      SOLI = NINT(ZUNITS) - S1I
      IF (SOLI .LT. 0) SOLI = 0
      CL = GRID(I,J,LIQUID,SCUR)
      CS = GRID(I,J,SOLID ,SCUR)
C Total solute inventory (integer quanta)
      MTOT = NINT( FLOAT(S1I)*CL + FLOAT(SOLI)*CS )
      IF (MTOT .LT. 0) MTOT = 0
C Write eutectic-closed state to SNEW
      GRID(I,J,LCELLS,SNEW) = 0.0
      GRID(I,J,LIQUID,SNEW) = 0.0
      GRID(I,J,SOLID ,SNEW) = FLOAT(MTOT)/ZUNITS
      GRID(I,J,PHASE ,SNEW) = SECONDARY
      GRID(I,J,CELLDIR,SNEW) = GRID(I,J,CELLDIR,SCUR)
      GINCEUT = GINCEUT + 1
      RETURN
      END
C
C =========================================================
C     Keep angles between 0 and 2PI (CCW) 
C =========================================================
      REAL FUNCTION NORMANG(ANG)
C     Normalize ANG to [0, 2PI) for positive CCW representation.
      IMPLICIT NONE
      REAL ANG, A
      REAL TWOPI
      PARAMETER (TWOPI = 6.28318531)
      A = ANG
      DO WHILE (A .LT. 0.0)
         A = A + TWOPI
      END DO
      DO WHILE (A .GE. TWOPI)
         A = A - TWOPI
      END DO
      NORMANG = A
      RETURN
      END
C    
C =========================================================
C     Min signed difference 
C =========================================================
      REAL FUNCTION ANGDIF(A, B)
C     Return minimal signed difference A - B in [-PI, PI].
      IMPLICIT NONE
      REAL A, B, DA
      REAL PI, TWOPI
      PARAMETER (PI = 3.14159265)
      PARAMETER (TWOPI = 6.28318531)
      DA = A - B
      DO WHILE (DA .GT. PI)
         DA = DA - TWOPI
      END DO
      DO WHILE (DA .LT. -PI)
         DA = DA + TWOPI
      END DO
      ANGDIF = DA
      RETURN
      END
C
C ==================================================================
C KVEFF: Aziz Effective Partition Coefficient Function      
C  Ref: M.J. Aziz, "Influence of the thermodynamic properties of the 
C  interface on the partition coefficient during rapid solidification,"
C  J. Applied Physics, 53(1982) p. 8899; doi: 10.1063/1.332629
C ====================================================================
      REAL FUNCTION KVEFF(V, K0, VD)
      IMPLICIT NONE
      REAL V, K0, VD, VRAT
      IF (VD .LE. 0.0) THEN
        KVEFF = K0
      ELSE
        VRAT = V/VD
        KVEFF = (K0 + VRAT)/(1.0 + VRAT)
      ENDIF
      RETURN
      END
C
C=====================================================================
C EJECT9 -- stage rejected solute into neighboring cells
C   REJMASS = rejected solute mass (wt%*subcells)
C   XADD    = staged solute-deposition buffer
C   LEFT    = remainder not accepted locally
C=====================================================================
      SUBROUTINE EJECT9(I, J, GRID, XLEN, YLEN, SCRN, REJMASS,
     &            XADD, MLIQ, TL0, C0,
     &            TLX, TLY, TLGX, TLGY,
     &            TLCLTIP, TLLDIF, TLWID, TLLEN,
     &            ITIPLOAD, LEFT, IBCTYPE)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER I, J, XLEN, YLEN, SCRN
      INTEGER II, JJ, I2, J2
      INTEGER ITIPLOAD, IBCTYPE
      REAL GRID(NXMAX, NYMAX, 6, 2), XADD(NXMAX, NYMAX)
      REAL REJMASS, MLIQ, TL0, C0, LEFT
      REAL TLX, TLY, TLGX, TLGY
      REAL TLCLTIP, TLLDIF, TLWID, TLLEN
      REAL CAP, CAPSUM, TAKE, TOL
      REAL GETCAP2
      PARAMETER (TOL = 1.0E-8)
C Default: all rejected mass remains unresolved
      LEFT = REJMASS
      IF (REJMASS .LE. TOL) RETURN
C Pass 1: sum admissible receiver capacity in the 3x3 neighborhood
      CAPSUM = 0.0
      DO JJ = -1, 1
      DO II = -1, 1
         IF (II .NE. 0 .OR. JJ .NE. 0) THEN
            CALL GETBCIJ(I+II,J+JJ,XLEN,YLEN,IBCTYPE,I2,J2)
            CAPSUM = CAPSUM + GETCAP2(I2, J2, GRID, SCRN,
     &         MLIQ, TL0, C0, TLX, TLY, TLGX, TLGY,
     &         TLCLTIP, TLLDIF, TLWID, TLLEN, ITIPLOAD)
         ENDIF
      ENDDO
      ENDDO
C No admissible local receiver capacity
      IF (CAPSUM .LE. TOL) RETURN
C Accepted portion of the rejected mass
      TAKE = AMIN1(REJMASS, CAPSUM)
C Pass 2: distribute accepted mass in proportion to local capacity
      DO JJ = -1, 1
      DO II = -1, 1
         IF (II .NE. 0 .OR. JJ .NE. 0) THEN
            CALL GETBCIJ(I+II,J+JJ,XLEN,YLEN,IBCTYPE,I2,J2)
            CAP = GETCAP2(I2, J2, GRID, SCRN,
     &           MLIQ, TL0, C0,
     &           TLX, TLY, TLGX, TLGY,
     &           TLCLTIP, TLLDIF, TLWID, TLLEN, ITIPLOAD)
            IF (CAP .GT. TOL) THEN
               XADD(I2, J2) = XADD(I2, J2) + TAKE*(CAP/CAPSUM)
            ENDIF
         ENDIF
      ENDDO
      ENDDO
C Any remainder stays unresolved and is returned in LEFT
      LEFT = REJMASS - TAKE
      IF (LEFT .LT. TOL) LEFT = 0.0
      RETURN
      END
C
C=====================================================================
C GETBCIJ -- Apply selected 2-D boundary condition to one index pair.
C=====================================================================
      SUBROUTINE GETBCIJ(ITEST,JTEST,XLEN,YLEN,IBCTYPE,I2,J2)
      IMPLICIT NONE
      INTEGER ITEST,JTEST,XLEN,YLEN,IBCTYPE,I2,J2
      LOGICAL XPER, YPER
      XPER = (IBCTYPE .EQ. 1 .OR. IBCTYPE .EQ. 3)
      YPER = (IBCTYPE .EQ. 2 .OR. IBCTYPE .EQ. 3)
      I2 = ITEST
      J2 = JTEST
      IF (XPER) THEN
         DO WHILE (I2 .LT. 1)
            I2 = I2 + XLEN
         END DO
         DO WHILE (I2 .GT. XLEN)
            I2 = I2 - XLEN
         END DO
      ELSE
         IF (I2 .LT. 1)    I2 = 1
         IF (I2 .GT. XLEN) I2 = XLEN
      ENDIF
      IF (YPER) THEN
         DO WHILE (J2 .LT. 1)
            J2 = J2 + YLEN
         END DO
         DO WHILE (J2 .GT. YLEN)
            J2 = J2 - YLEN
         END DO
      ELSE
         IF (J2 .LT. 1)    J2 = 1
         IF (J2 .GT. YLEN) J2 = YLEN
      ENDIF
      RETURN
      END
C
C=====================================================================
C GETCAP2  Local capture driving force (solute uptake)
C  Computes local solute capture term at an interface cell based on
C  the difference between the current liquid composition (CL) and the
C  local equilibrium composition (CLSTAR).
C      GETCAP2 = S1 * max(0, CLSTARUSE - CL)
C  where S1 is the liquid fraction (subcells) in the control volume.
C=====================================================================
      REAL FUNCTION GETCAP2(I2, J2, GRID, SCRN, MLIQ, TL0, C0,
     & TLX, TLY, TLGX, TLGY, TLCLTIP, TLLDIF, TLWID, TLLEN, ITIPLOAD)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER I2, J2, SCRN, ITIPLOAD
      REAL GRID(NXMAX, NYMAX, 6, 2), MLIQ, TL0, C0
      REAL TLX, TLY, TLGX, TLGY
      REAL TLCLTIP, TLLDIF, TLWID, TLLEN
      REAL S1, CL, TL, CLSTARLOC, CLSTARTGT, CLSTARUSE, DELCL
      REAL CLIQBIN
      EXTERNAL CLIQBIN
      REAL DXP, DYP, SPAR, SPER, ARG
      INTEGER LCELLS, LIQUID, TCELL
        PARAMETER (LCELLS=1, LIQUID=2, TCELL=5)
      GETCAP2 = 0.0
      S1 = GRID(I2, J2, LCELLS, SCRN)
      IF (S1 .LE. 0.0) RETURN
      CL = GRID(I2, J2, LIQUID, SCRN)
      TL = GRID(I2, J2, TCELL , SCRN)
      CLSTARLOC = C0
      IF (ABS(MLIQ) .GT. 1.0E-12) THEN
         CLSTARLOC = CLIQBIN(TL,TL0,MLIQ,C0)
      ENDIF
      CLSTARUSE = CLSTARLOC
      IF (ITIPLOAD .EQ. 1 .AND. TLLDIF .GT. 1.0E-6) THEN
         DXP = FLOAT(I2) - TLX
         DYP = FLOAT(J2) - TLY
         SPAR = DXP*TLGX + DYP*TLGY
         SPER = ABS(-DXP*TLGY + DYP*TLGX)
         IF (SPAR .GE. 0.0 .AND. SPAR .LE. TLLEN
     & .AND. SPER .LE. TLWID) THEN
            ARG = -SPAR/TLLDIF
            IF (ARG .LT. -50.0) ARG = -50.0
            CLSTARTGT = C0 + (TLCLTIP - C0)*EXP(ARG)
            CLSTARUSE = AMIN1(CLSTARLOC, CLSTARTGT)
         ENDIF
      ENDIF
      DELCL = CLSTARUSE - CL
      IF (DELCL .GT. 0.0) THEN
         GETCAP2 = S1*DELCL
      ENDIF
      RETURN
      END
C
C ====================================================================
C DIFFUSE -- Conservative explicit liquid-solute diffusion update.
C Transport scheme based on:
C   Beltran-Sanchez, L., Stefanescu, D.M.,
C   Metall. Mater. Trans. A 34 (2003) 367-382.
C   https://doi.org/10.1007/s11661-003-0338-z
C Uses face fluxes weighted by liquid fraction:
C   Q = W*C
C   F = alpha*min(Wi,Wj)*(Cj-Ci)
C   Cnew = Qnew/W
C Faces touching partially-liquid cells are scaled by DMASK.
C Boundary conditions are selected by IBCTYPE through DBCFLAGS.
C alpha = D*dt/dx**2.
C ====================================================================
      SUBROUTINE DIFFUSE(GRID,XLEN,YLEN,SCUR,ALPHA,DMASK,IBCTYPE)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN, SCUR, IBCTYPE
      REAL GRID(NXMAX, NYMAX, 6, 2)
      REAL ALPHA, DMASK
      LOGICAL XPER, YPER
      REAL DMLOCAL
      REAL FX(XLEN, YLEN)
      REAL FY(XLEN, YLEN+1)
      REAL CNEW(XLEN, YLEN)
      CALL DBCFLAGS(IBCTYPE,XPER,YPER)
       DMLOCAL = AMAX1(DMASK,0.05)
      CALL DFLUXX(GRID,XLEN,YLEN,SCUR,ALPHA,DMLOCAL,XPER,FX)
      CALL DFLUXY(GRID,XLEN,YLEN,SCUR,ALPHA,DMLOCAL,YPER,FY)
      CALL DUPDATEC(GRID,XLEN,YLEN,SCUR,XPER,FX,FY,CNEW)
      CALL DCOMMITC(GRID,XLEN,YLEN,SCUR,CNEW)
      RETURN
      END 
C 
C ====================================================================
C DBCFLAGS -- Convert integer boundary-condition code to logical flags.
C ====================================================================
      SUBROUTINE DBCFLAGS(IBCTYPE,XPER,YPER)
      IMPLICIT NONE
      INTEGER IBCTYPE
      LOGICAL XPER, YPER
      XPER = (IBCTYPE .EQ. 1 .OR. IBCTYPE .EQ. 3)
      YPER = (IBCTYPE .EQ. 2 .OR. IBCTYPE .EQ. 3)
      RETURN
      END
C
C ====================================================================
C DFLUXX -- Build conservative X-face liquid-solute fluxes.
C ====================================================================
      SUBROUTINE DFLUXX(GRID,XLEN,YLEN,SCUR,ALPHA,DMLOCAL,XPER,FX)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN, SCUR, LCELLS, LIQUID
      INTEGER I, J, IP
      REAL GRID(NXMAX, NYMAX, 6, 2)
      REAL ALPHA, DMLOCAL
      REAL FX(XLEN,YLEN)
      REAL CC, CN, WC, WN, AFACE
      LOGICAL XPER
      REAL TOLW, Z1, INVN
      PARAMETER (INVN = 1.0/1000.0)
      PARAMETER (TOLW=1.0E-6, Z1=0.9999)
      PARAMETER (LCELLS=1, LIQUID=2)
      DO J = 1, YLEN
      DO I = 1, XLEN
         IF (I .EQ. XLEN .AND. .NOT. XPER) THEN
            FX(I,J) = 0.0
         ELSE
            IP = I + 1
            IF (IP .GT. XLEN) IP = 1
            WC = GRID(I ,J,LCELLS,SCUR)*INVN
            WN = GRID(IP,J,LCELLS,SCUR)*INVN
            AFACE = AMIN1(WC, WN)
            IF (AFACE .GT. TOLW) THEN
               IF (WC .LT. Z1 .OR. WN .LT. Z1) AFACE = AFACE*DMLOCAL
               CC = GRID(I ,J,LIQUID,SCUR)
               CN = GRID(IP,J,LIQUID,SCUR)
               FX(I,J) = ALPHA*AFACE*(CN - CC)
            ELSE
               FX(I,J) = 0.0
            ENDIF
         ENDIF
      END DO
      END DO
      RETURN
      END
C
C ====================================================================
C DFLUXY -- Build conservative Y-face liquid-solute fluxes.
C ====================================================================
      SUBROUTINE DFLUXY(GRID,XLEN,YLEN,SCUR,ALPHA,DMLOCAL,
     &     YPER,FY)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN, SCUR, LCELLS, LIQUID
      INTEGER I, J
      REAL GRID(NXMAX, NYMAX, 6, 2)
      REAL ALPHA, DMLOCAL
      REAL FY(XLEN,YLEN+1)
      REAL CC, CN, WC, WN, AFACE
      LOGICAL YPER
      REAL TOLW, Z1, INVN
      PARAMETER (INVN = 1.0/1000.0)
      PARAMETER (TOLW=1.0E-6, Z1=0.9999)
      PARAMETER (LCELLS=1, LIQUID=2)
C Closed/default boundary initialization.
      DO I = 1, XLEN
         FY(I,1) = 0.0
         FY(I,YLEN+1) = 0.0
      END DO
C Internal Y faces.
      DO J = 1, YLEN-1
      DO I = 1, XLEN
         WC = GRID(I,J  ,LCELLS,SCUR)*INVN
         WN = GRID(I,J+1,LCELLS,SCUR)*INVN
         AFACE = AMIN1(WC, WN)
         IF (AFACE .GT. TOLW) THEN
            IF (WC .LT. Z1 .OR. WN .LT. Z1) AFACE = AFACE*DMLOCAL
            CC = GRID(I,J  ,LIQUID,SCUR)
            CN = GRID(I,J+1,LIQUID,SCUR)
            FY(I,J+1) = ALPHA*AFACE*(CN - CC)
         ELSE
            FY(I,J+1) = 0.0
         ENDIF
      END DO
      END DO
C Periodic Y face between YLEN and 1.
      IF (YPER) THEN
         DO I = 1, XLEN
            WC = GRID(I,YLEN,LCELLS,SCUR)*INVN
            WN = GRID(I,1   ,LCELLS,SCUR)*INVN
            AFACE = AMIN1(WC, WN)
            IF (AFACE .GT. TOLW) THEN
               IF (WC .LT. Z1 .OR. WN .LT. Z1) AFACE = AFACE*DMLOCAL
               CC = GRID(I,YLEN,LIQUID,SCUR)
               CN = GRID(I,1   ,LIQUID,SCUR)
               FY(I,YLEN+1) = ALPHA*AFACE*(CN - CC)
               FY(I,1)      = FY(I,YLEN+1)
            ELSE
               FY(I,YLEN+1) = 0.0
               FY(I,1)      = 0.0
            ENDIF
         END DO
      ENDIF
      RETURN
      END
C
C ====================================================================
C DUPDATEC -- Apply face-flux divergence and compute CNEW.
C ====================================================================
      SUBROUTINE DUPDATEC(GRID,XLEN,YLEN,SCUR,XPER,FX,FY,CNEW)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN, SCUR, LCELLS, LIQUID
      INTEGER I, J, IM
      REAL GRID(NXMAX, NYMAX, 6, 2)
      REAL FX(XLEN,YLEN)
      REAL FY(XLEN,YLEN+1)
      REAL CNEW(XLEN,YLEN)
      LOGICAL XPER
      REAL WC, COLD, Q, QNEW
      REAL TOLW, INVN
      PARAMETER (INVN = 1.0/1000.0)
      PARAMETER (TOLW=1.0E-6)
      PARAMETER (LCELLS=1, LIQUID=2)
      DO J = 1, YLEN
      DO I = 1, XLEN
         WC = GRID(I,J,LCELLS,SCUR)*INVN
         COLD = GRID(I,J,LIQUID,SCUR)
         CNEW(I,J) = COLD
         IF (WC .GT. TOLW) THEN
           IM = I - 1
           IF (IM .LT. 1) THEN
              IF (XPER) THEN
                 IM = XLEN
              ELSE
                 IM = 0
              ENDIF
           ENDIF
           Q = WC*COLD
           QNEW = Q + FX(I,J) + FY(I,J+1) - FY(I,J)
           IF (IM .GT. 0) THEN
              QNEW = QNEW - FX(IM,J)
           ENDIF
           CNEW(I,J) = QNEW/WC
           IF (CNEW(I,J) .LT. 0.0) CNEW(I,J) = 0.0
         ENDIF
      END DO
      END DO
      RETURN
      END
C
C ====================================================================
C DCOMMITC -- Commit updated LIQUID composition back to GRID.
C ====================================================================
      SUBROUTINE DCOMMITC(GRID,XLEN,YLEN,SCUR,CNEW)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN, SCUR, LIQUID
      INTEGER I, J
      REAL GRID(NXMAX, NYMAX, 6, 2)
      REAL CNEW(XLEN,YLEN)
      PARAMETER (LIQUID=2)
      DO J = 1, YLEN
      DO I = 1, XLEN
         GRID(I,J,LIQUID,SCUR) = CNEW(I,J)
      END DO
      END DO
      RETURN
      END
C
C=====================================================================
C COMMA2BLANK -- prepare simple CSV line for list-directed read
C=====================================================================
      SUBROUTINE COMMA2BLANK(LINE)
      IMPLICIT NONE
      CHARACTER*(*) LINE
      INTEGER I, N, IC
        N = LEN(LINE)
      DO I = 1, N
         IC = ICHAR(LINE(I:I))
         IF (LINE(I:I) .EQ. ',') LINE(I:I) = ' '
         IF (IC .EQ. 34) LINE(I:I) = ' '
      ENDDO
      RETURN
      END
C
C=====================================================================
C LOADTHERMCSV -- read imposed macro thermal-history table
C Expected columns:
C   time_s, Tref_C, G_K_per_m, cooling_rate_K_per_s,
C   Vthermal_m_per_s
C=====================================================================
      SUBROUTINE LOADTHERMCSV(FNAME,TTAB,TREFTAB,GTAB,RTAB,VTAB,
     &     MAXTHERM,NTHERM)
      IMPLICIT NONE
      CHARACTER*(*) FNAME
      INTEGER MAXTHERM, NTHERM
      REAL TTAB(MAXTHERM), TREFTAB(MAXTHERM), GTAB(MAXTHERM)
      REAL RTAB(MAXTHERM), VTAB(MAXTHERM)
      INTEGER UNIT, IOS, JSTAT
      CHARACTER*160 LINE
      REAL T, TR, GX, GY, GREAD, RR, VV
      LOGICAL OPENOK
      UNIT = 24
      NTHERM = 0
      OPENOK = .TRUE.
      OPEN(UNIT=UNIT,FILE=FNAME,STATUS='OLD',IOSTAT=IOS)
      IF (IOS .NE. 0) THEN
         WRITE(*,*) 'LOADTHERMCSV: could not open ', FNAME
         OPENOK = .FALSE.
      ENDIF
      IF (OPENOK) THEN
       READ(UNIT,'(A)',IOSTAT=IOS) LINE
       DO WHILE (IOS .EQ. 0)
        IF (LINE .NE. ' ') THEN
           CALL COMMA2BLANK(LINE)
           READ(LINE,*,IOSTAT=JSTAT) T, TR, GX, GY, RR, VV
           IF (JSTAT .EQ. 0) THEN
             GREAD = GY
           ELSE
             READ(LINE,*,IOSTAT=JSTAT) T, TR, GREAD, RR, VV
           ENDIF
           IF (JSTAT .EQ. 0) THEN
            IF (NTHERM .LT. MAXTHERM) THEN
              NTHERM = NTHERM + 1
              TTAB(NTHERM) = T
              TREFTAB(NTHERM) = TR
              GTAB(NTHERM) = GREAD
              RTAB(NTHERM)= RR
              VTAB(NTHERM) = VV
             ENDIF
           ENDIF
        ENDIF
        READ(UNIT,'(A)',IOSTAT=IOS) LINE
       ENDDO
       CLOSE(UNIT)
       WRITE(*,*) 'LOADTHERMCSV loaded ', NTHERM,' rows from ', FNAME
      ENDIF
      RETURN
      END
C
C=====================================================================
C THERMLOOK -- linear interpolation in imposed thermal table
C=====================================================================
      SUBROUTINE THERMLOOK(TIME,TTAB,TREFTAB,GTAB,RTAB,VTAB,NTHERM,
     &     TREF, GKM, RATE, VMS)
      IMPLICIT NONE
      INTEGER NTHERM, K
      REAL TIME, TTAB(*), TREFTAB(*), GTAB(*), RTAB(*), VTAB(*)
      REAL TREF, GKM, RATE, VMS, A, DTAB
      LOGICAL FOUND
      FOUND = .FALSE.
      IF (NTHERM .LE. 0) THEN
       TREF = 0.0
       GKM  = 0.0
       RATE = 0.0
       VMS  = 0.0
       FOUND = .TRUE.
      ENDIF
      IF (.NOT. FOUND) THEN
       IF (TIME .LE. TTAB(1) .OR. NTHERM .EQ. 1) THEN
         TREF = TREFTAB(1)
         GKM  = GTAB(1)
         RATE = RTAB(1)
         VMS  = VTAB(1)
         FOUND = .TRUE.
       ENDIF
      ENDIF
      IF (.NOT. FOUND) THEN
       IF (TIME .GE. TTAB(NTHERM)) THEN
         TREF = TREFTAB(NTHERM)
         GKM  = GTAB(NTHERM)
         RATE = RTAB(NTHERM)
         VMS  = VTAB(NTHERM)
         FOUND = .TRUE.
        ENDIF
      ENDIF
      K = 1
      DO WHILE (K .LE. NTHERM-1 .AND. .NOT. FOUND)
       IF (TIME .GE. TTAB(K) .AND. TIME .LE. TTAB(K+1)) THEN
          DTAB = TTAB(K+1) - TTAB(K)
          IF (DTAB .LE. 0.0) THEN
             A = 0.0
          ELSE
             A = (TIME - TTAB(K))/DTAB
          ENDIF
          TREF = TREFTAB(K) + A*(TREFTAB(K+1)-TREFTAB(K))
          GKM  = GTAB(K) + A*(GTAB(K+1)-GTAB(K))
          RATE = RTAB(K) + A*(RTAB(K+1)-RTAB(K))
          VMS  = VTAB(K) + A*(VTAB(K+1)-VTAB(K))
          FOUND = .TRUE.
       ENDIF
       K = K + 1
      ENDDO
      IF (.NOT. FOUND) THEN
       TREF = TREFTAB(NTHERM)
       GKM  = GTAB(NTHERM)
       RATE = RTAB(NTHERM)
       VMS  = VTAB(NTHERM)
      ENDIF
      RETURN
      END
C=====================================================================
C TMACRO1D -- imposed local macro field
C The CSV Tref is interpreted at reference row JREFTH after an
C optional one-time TSHIFT:
C   Tmacro(J,t) = Tref(t) + TSHIFT + G(t)*(J-JREFTH)*DX
C Embedded absolute mode: JREFTH = center of the CA window, TSHIFT=0.
C Embedded seed-align:    JREFTH = center, but TSHIFT anchors the
C                         initial seed/front at TL0-DTSEED.
C External mode:          JREFTH = initial seed/front row, TSHIFT
C                         anchors the initial front at TL0-DTSEED.
C=====================================================================
      REAL FUNCTION TMACRO1D(J,DX,TREF,GKM,JREFTH,TSHIFT)
      IMPLICIT NONE
      INTEGER J, JREFTH
      REAL DX, TREF, GKM, TSHIFT
      TMACRO1D = TREF + TSHIFT + GKM*(FLOAT(J)-FLOAT(JREFTH))*DX
      RETURN
      END
C
C=====================================================================
C INITTHERMCSV -- initialize macro temperature field
C One thermal-history table is used for both modes.  The difference is
C how the table temperature Tref is anchored to the CA domain:
C   IEMBED = 1 : direct macro-cell anchoring
C                Tref is the imposed macro temperature at the center
C                of the CA window.  No temperature shift is applied.
C                This is the intended absolute embedded mode.
C   IEMBED = 2 : embedded seed-aligned anchoring
C                Tref is still interpreted at the center of the CA
C                window, but the whole macro field is shifted so the
C                initial seed/front row starts at TL0-DTSEED.
C   IEMBED = 0 : seed-aligned standalone anchoring
C                Tref is shifted so the initial seed/front row starts
C                at TL0-DTSEED.  This is used for Bridgman-style
C                standalone tests and comparison to fixed-gradient runs.
C After this call, TCELL is reset to the imposed macro field, so the
C initial perturbation temperature is zero.
C=====================================================================
      SUBROUTINE INITTHERMCSV(GRID,XLEN,YLEN,SCRN,SNEW,DX,
     &     GRAD,VEL,TTAB,TREFTAB,GTAB,RTAB,VTAB,NTHERM,
     &     IEMBED,TL0,DTSEED,JREFTH,TSHIFT)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN, SCRN, SNEW, NTHERM
      INTEGER IEMBED, JREFTH
      INTEGER I, J, JFRONT, TCELL
      PARAMETER (TCELL=5)
      REAL GRID(NXMAX,NYMAX,6,2)
      REAL DX, GRAD, VEL, TL0, DTSEED, TSHIFT
      REAL TTAB(*), TREFTAB(*), GTAB(*), RTAB(*), VTAB(*)
      REAL TREF, GKM, RATE, VMS, TFRONT, TCSVFRONT
      REAL TMACRO1D
      EXTERNAL TMACRO1D
      CALL THERMLOOK(0.0,TTAB,TREFTAB,GTAB,RTAB,VTAB,NTHERM,
     &     TREF,GKM,RATE,VMS)
      IF (IEMBED .EQ. 1) THEN
C Absolute embedded anchoring: trust the CSV directly.
         JREFTH = MAX(1,MIN(YLEN,YLEN/2))
         TSHIFT = 0.0
      ELSE IF (IEMBED .EQ. 2) THEN
C Embedded seed-aligned anchoring: keep the CSV gradient/history
C referenced to the center row, but shift the field so the
C initial seed/front row starts at TL0-DTSEED.
         JREFTH = MAX(1,MIN(YLEN,YLEN/2))
         CALL FNDLOWR(GRID,XLEN,YLEN,SCRN,JFRONT)
         TFRONT = TL0 - DTSEED
         TCSVFRONT = TMACRO1D(JFRONT,DX,TREF,GKM,JREFTH,0.0)
         TSHIFT = TFRONT - TCSVFRONT
      ELSE
C Standalone seed-aligned anchoring: Tref itself is placed on
C the initial seed/front row.
         CALL FINDFRONTROW(GRID,XLEN,YLEN,SCRN,JREFTH)
         TFRONT = TL0 - DTSEED
         TCSVFRONT = TREF
         TSHIFT = TFRONT - TCSVFRONT
      ENDIF
      GRAD = GKM*DX
C For normal directional cases, VMS comes from the CSV as the
C thermal/front velocity.  For uniform cooling, GKM may be zero,
C but VMS is still used as the KGT closure velocity scale.
      IF (DX .GT. 0.0) THEN
         VEL = VMS/DX
      ELSE
         VEL = 0.0
      ENDIF
      GRAD = GKM*DX
      DO J = 1, YLEN
      DO I = 1, XLEN
         GRID(I,J,TCELL,SCRN) = TMACRO1D(J,DX,TREF,GKM,JREFTH,TSHIFT)
         GRID(I,J,TCELL,SNEW) = GRID(I,J,TCELL,SCRN)
      ENDDO
      ENDDO
      WRITE(*,*) 'INITTHERMCSV: TCELL reset to imposed macro field'
      WRITE(*,*) '              IEMBED=', IEMBED,
     &           '  JREF=', JREFTH, '  TSHIFT=', TSHIFT
      WRITE(*,*) '              Tref(csv)=', TREF, '  G[K/m]=', GKM
      WRITE(*,*) '              Grad[K/cell]=', GRAD,
     &           '  Vel[cell/s]=', VEL
      WRITE(*,*) 'INITTHERMCSV debug:'
      WRITE(*,*) '  IEMBED  = ', IEMBED
      WRITE(*,*) '  TREF    = ', TREF
      WRITE(*,*) '  TL0     = ', TL0
      WRITE(*,*) '  DTSEED  = ', DTSEED
      WRITE(*,*) '  TL0-DTSEED = ', TL0-DTSEED
      WRITE(*,*) '  JREFTH  = ', JREFTH
      WRITE(*,*) '  TSHIFT  = ', TSHIFT
      RETURN
      END
C
C=====================================================================
C FINDFRONTROW -- find highest initial primary-solid row
C=====================================================================
      SUBROUTINE FINDFRONTROW(GRID,XLEN,YLEN,SCRN,JREFTH)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN, SCRN, JREFTH
      INTEGER I, J, LCELLS
      PARAMETER (LCELLS=1)
      REAL GRID(NXMAX,NYMAX,6,2)
      JREFTH = YLEN/2
      DO J = 1, YLEN
      DO I = 1, XLEN
         IF (GRID(I,J,LCELLS,SCRN) .LT. 0.5) THEN
            IF (J .GT. JREFTH) JREFTH = J
         ENDIF
      ENDDO
      ENDDO
      IF (JREFTH .LT. 1) JREFTH = 1
      IF (JREFTH .GT. YLEN) JREFTH = YLEN
      RETURN
      END
C
C=====================================================================
C HEATPERT2D -- embedded thermal-history update
C The casting-scale thermal path is imposed by the CSV table.  The CA
C only evolves the local perturbation from latent heat:
C   GRID(:,:,TCELL) = Tmacro + Tpert
C Each CA step:
C   1. recover Tpert from old macro field
C   2. rebuild total T with the new macro field
C   3. diffuse/source Tpert only with Tpert=0 at top/bottom
C   4. rebuild total T = Tmacro + Tpert
C X is periodic.  The Y boundaries are embedded in the macro field,
C so the perturbation is clamped to zero there.
C=====================================================================
      SUBROUTINE HEATPERT2D(GRID,QLAT,XLEN,YLEN,SCRN,DT,DX,ATHERM,
     &     GTIME,VEL,GRAD,TTAB,TREFTAB,GTAB,RTAB,VTAB,NTHERM,
     &     JREFTH,TSHIFT)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN, SCRN, NTHERM
      INTEGER JREFTH
      REAL GRID(NXMAX,NYMAX,6,2), QLAT(NXMAX,NYMAX)
      REAL DT, DX, ATHERM, GTIME, VEL, GRAD, TSHIFT
      REAL TTAB(*), TREFTAB(*), GTAB(*), RTAB(*), VTAB(*)
      INTEGER I, J, IP, IM, SNEW, ISUB, NSUB
      INTEGER TCELL
      PARAMETER (TCELL=5)
      REAL TREF0, GKM0, RATE0, VMS0
      REAL TREF1, GKM1, RATE1, VMS1
      REAL TIME0, TALPHA0, TALPHA, QSCALE
      REAL TEMPOLD
      REAL MC, ME, MW, MN, MS, PC, PE, PW, PN, PS, PNEW
      REAL TMACRO1D
      EXTERNAL TMACRO1D
      SNEW = 3 - SCRN
      TIME0 = GTIME - DT
      IF (TIME0 .LT. 0.0) TIME0 = 0.0
      CALL THERMLOOK(TIME0,TTAB,TREFTAB,GTAB,RTAB,VTAB,NTHERM,
     &     TREF0,GKM0,RATE0,VMS0)
      CALL THERMLOOK(GTIME,TTAB,TREFTAB,GTAB,RTAB,VTAB,NTHERM,
     &     TREF1,GKM1,RATE1,VMS1)
C Update runtime diagnostic/closure values from the imposed table.
      GRAD = GKM1*DX
      IF (DX .GT. 0.0) THEN
         VEL = VMS1/DX
      ENDIF
C Shift from old macro field to new macro field while preserving Tpert.
C Save old total T before overwriting TCELL.
      DO J = 1, YLEN
      DO I = 1, XLEN
         TEMPOLD = GRID(I,J,TCELL,SCRN)
         GRID(I,J,TCELL,SCRN)= TMACRO1D(J,DX,TREF1,GKM1,JREFTH,TSHIFT)
     &      + TEMPOLD - TMACRO1D(J,DX,TREF0,GKM0,JREFTH,TSHIFT)
      ENDDO
      ENDDO
C Select explicit thermal subcycling.
      TALPHA0 = ATHERM*DT/(DX*DX)
      NSUB = 1
      IF (TALPHA0 .GT. 0.20) NSUB = INT(TALPHA0/0.20) + 1
      TALPHA = TALPHA0/FLOAT(NSUB)
      QSCALE = 1.0/FLOAT(NSUB)
      IF (TALPHA .GT. 0.25) THEN
         WRITE(*,*) 'HEATPERT2D warning: TALPHA = ', TALPHA
      ENDIF
      DO ISUB = 1, NSUB
         DO J = 2, YLEN-1
         DO I = 1, XLEN
            IP = I + 1
            IM = I - 1
            IF (IP .GT. XLEN) IP = 1
            IF (IM .LT. 1)    IM = XLEN
            MC = TMACRO1D(J,DX,TREF1,GKM1,JREFTH,TSHIFT)
            ME = MC
            MW = MC
            MN = TMACRO1D(J+1,DX,TREF1,GKM1,JREFTH,TSHIFT)
            MS = TMACRO1D(J-1,DX,TREF1,GKM1,JREFTH,TSHIFT)
            PC = GRID(I ,J,TCELL,SCRN) - MC
            PE = GRID(IP,J,TCELL,SCRN) - ME
            PW = GRID(IM,J,TCELL,SCRN) - MW
            PN = GRID(I,J+1,TCELL,SCRN) - MN
            PS = GRID(I,J-1,TCELL,SCRN) - MS
            PNEW = PC + TALPHA*(PE + PW + PN + PS - 4.0*PC)
     &           + QSCALE*QLAT(I,J)
            GRID(I,J,TCELL,SNEW) = MC + PNEW
         ENDDO
         ENDDO
C Embedded thermal bath: no local perturbation at Y boundaries
         DO I = 1, XLEN
          GRID(I,1,TCELL,SNEW)=TMACRO1D(1,DX,TREF1,GKM1,JREFTH,TSHIFT)
          GRID(I,YLEN,TCELL,SNEW) =
     &         TMACRO1D(YLEN,DX,TREF1,GKM1,JREFTH,TSHIFT)
         ENDDO
         DO J = 1, YLEN
         DO I = 1, XLEN
            GRID(I,J,TCELL,SCRN) = GRID(I,J,TCELL,SNEW)
         ENDDO
         ENDDO
      ENDDO
      DO J = 1, YLEN
      DO I = 1, XLEN
         QLAT(I,J) = 0.0
      ENDDO
      ENDDO
      RETURN
      END
C
C=====================================================================
C XRSUMAVG -- Compute domain-average residual solute in XRES.
C   Returns the average solute stored in the residual reservoir XRES,
C   expressed in wt% per cell over the full CA domain.
C    - XRES stores solute as (wt% * subcells) for each cell.
C    - Division by ZUNITS converts to wt% per cell.
C    - This quantity represents solute temporarily held in the staged
C      reservoir (not yet released to the liquid grid).
C    - Useful for monitoring mass balance and transient staging.
C=====================================================================
      REAL FUNCTION XRSUMAVG(XRES,XLEN,YLEN)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN, I, J
      REAL XRES(NXMAX,NYMAX)
      REAL ZUNITS
      PARAMETER (ZUNITS=1000.0)
      DOUBLE PRECISION SUMXR
       SUMXR = 0.0D0
      DO J=1,YLEN
      DO I=1,XLEN
         SUMXR = SUMXR + DBLE(XRES(I,J))/DBLE(ZUNITS)
      ENDDO
      ENDDO
       XRSUMAVG = REAL(SUMXR/(DBLE(XLEN)*DBLE(YLEN)))
      RETURN
      END
C
C=====================================================================
C FNDLOWR -- find lowest initial primary-solid/interface row
C  Used for IEMBED=2 startup alignment.
C  This finds the first occupied row near the bottom of the seed/front,
C  so the embedded thermal field can be shifted such that this row
C  starts at TL0-DTSEED.
C=====================================================================
      SUBROUTINE FNDLOWR(GRID,XLEN,YLEN,SCRN,JSEED)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN, SCRN, JSEED
      INTEGER I, J, LCELLS
      PARAMETER (LCELLS=1)
      REAL GRID(NXMAX,NYMAX,6,2)
      JSEED = YLEN
      DO J = 1, YLEN
      DO I = 1, XLEN
         IF (GRID(I,J,LCELLS,SCRN) .LT. 0.5) THEN
            IF (J .LT. JSEED) JSEED = J
         ENDIF
      ENDDO
      ENDDO
      IF (JSEED .LT. 1) JSEED = 1
      IF (JSEED .GT. YLEN) JSEED = YLEN
      RETURN
      END      
C
C=====================================================================
C MOVPPM-- write next numbered movie frame: out1.ppm, out2.ppm, ...
C  Each call writes a side-by-side 2X periodic-style PPM image
C=====================================================================
      SUBROUTINE MOVPPM(GRID, XLEN, YLEN, SCRN, C0, CEUT)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN, SCRN
      REAL GRID(NXMAX, NYMAX, 6, 2)
      REAL C0, CEUT, CMIN, CMAX
      INTEGER NFRAME
      CHARACTER*128 FNAME
      SAVE NFRAME
      DATA NFRAME /0/
      NFRAME = NFRAME + 1
      CALL MOVNAME(NFRAME, FNAME)
      CALL MOVCLRANGE(GRID, XLEN, YLEN, SCRN, CMIN, CMAX, C0, CEUT)
      CALL MOVSCALEPPM(GRID, XLEN, YLEN, SCRN, FNAME, CMIN, CMAX)
      RETURN
      END
C
C=====================================================================
C  MOVNAME -- make filename outN.ppm using only F77-safe operations
C=====================================================================
      SUBROUTINE MOVNAME(N, FNAME)
      IMPLICIT NONE
      INTEGER N
      CHARACTER*(*) FNAME
      CHARACTER*12 NUM
      INTEGER I, K, L
      FNAME = ' '
      FNAME(1:3) = 'out'
      WRITE(NUM,'(I12)') N
      K = 4
      DO I = 1, 12
       IF (NUM(I:I) .NE. ' ') THEN
        IF (K .LE. LEN(FNAME)) THEN
          FNAME(K:K) = NUM(I:I)
          K = K + 1
        ENDIF
       ENDIF
      ENDDO
      L = LEN(FNAME)
      IF (K+3 .LE. L) THEN
         FNAME(K:K+3) = '.ppm'
      ENDIF
      RETURN
      END
C
C=====================================================================
C  MOVCLRANGE -- local composition range for liquid cells
C=====================================================================
      SUBROUTINE MOVCLRANGE(GRID, XLEN,YLEN,SCRN, CMIN, CMAX,C0,CEUT)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN, SCRN
      REAL GRID(NXMAX, NYMAX, 6, 2)
      REAL CMIN, CMAX, C0, CEUT
      INTEGER LCELLS, LIQUID
      PARAMETER (LCELLS=1, LIQUID=2)
      INTEGER I, J
      REAL COMP, S1
      CMIN = 999.0
      CMAX = -999.0
      DO J = 1, YLEN
      DO I = 1, XLEN
       S1 = GRID(I, J, LCELLS, SCRN)
       IF (S1 .GT. 0.0) THEN
         COMP = GRID(I, J, LIQUID, SCRN)
         IF (COMP .LT. CMIN) CMIN = COMP
         IF (COMP .GT. CMAX) CMAX = COMP
       ENDIF
      ENDDO
      ENDDO
      IF (CMIN .GT. CMAX) THEN
         CMIN = C0
         CMAX = CEUT
      ENDIF
      IF (CMIN .LT. C0) CMIN = C0
      IF (CMAX .GT. CEUT) CMAX = CEUT
      IF (CMAX .LE. CMIN) CMAX = CMIN + 1.0E-6
      RETURN
      END
C
C=====================================================================
C  MOVSCALEPPM -- write ASCII P3 PPM, side-by-side 2X image
C=====================================================================
      SUBROUTINE MOVSCALEPPM(GRID, XLEN,YLEN, SCRN, FNAME, CMIN, CMAX)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN, SCRN
      REAL GRID(NXMAX, NYMAX, 6, 2)
      CHARACTER*(*) FNAME
      REAL CMIN, CMAX
      INTEGER W, H, LUN, IOUT, I, J, R, G, B
      W = 2*XLEN
      H = YLEN
      LUN = 33
      OPEN(UNIT=LUN, FILE=FNAME, STATUS='UNKNOWN', FORM='FORMATTED')
      WRITE(LUN,'(A)') 'P3'
      WRITE(LUN,'(I10,1X,I10)') W, H
      WRITE(LUN,'(I10)') 255
C Write top image row first.  IOUT runs across the 2X tiled image.
      DO J = YLEN, 1, -1
      DO IOUT = 1, W
        I = MOD(IOUT-1, XLEN) + 1
        CALL MOVCOLOR(GRID, XLEN, YLEN, SCRN, I, J, CMIN,CMAX, R,G,B)
        WRITE(LUN,'(3(I4,1X))') R, G, B
      ENDDO
      ENDDO
      CLOSE(LUN)
      WRITE(*,*) 'PPM movie frame written: ', FNAME
      RETURN
      END
C
C=====================================================================
C  MOVCOLOR -- color rule for one CA cell
C=====================================================================
      SUBROUTINE MOVCOLOR(GRID,XLEN,YLEN, SCRN, I,J, CMIN,CMAX, R,G,B)
      IMPLICIT NONE
      INTEGER NXMAX, NYMAX
      PARAMETER (NXMAX=1000, NYMAX=1000)
      INTEGER XLEN, YLEN, SCRN, I, J, R, G, B
      REAL GRID(NXMAX, NYMAX, 6, 2)
      REAL CMIN, CMAX
      INTEGER LCELLS, LIQUID, PHASE
      PARAMETER (LCELLS=1, LIQUID=2, PHASE=4)
      REAL SECONDARY, DEFECT
      PARAMETER (SECONDARY=2.0, DEFECT=3.0)
      INTEGER II, JJ, I1, I2, J1, J2, RTHICK
      INTEGER GCHAN
      REAL COMP, NORM, GRAY, F1, F2, S1, PHS
      LOGICAL THICKSOL
      PARAMETER (F1=1.45, F2=0.60, RTHICK=1)
C Liquid composition color.
      COMP = GRID(I,J,LIQUID,SCRN)
      IF (COMP .LT. CMIN) COMP = CMIN
      IF (COMP .GT. CMAX) COMP = CMAX
      NORM = (COMP - CMIN)/(CMAX - CMIN)
      GRAY = 255.0*(1.0 - NORM)
      G = MIN(255, NINT(GRAY))
      R = G
      GCHAN = MIN(255, NINT(G*F1))
      B = MIN(255, NINT(G*F2))
C Small cosmetic R=1 solid thickening.
      THICKSOL = .FALSE.
      S1 = GRID(I,J,LCELLS,SCRN)
      IF (S1 .LT. 0.5) THEN
         THICKSOL = .TRUE.
      ELSE
         I1 = MAX(1, I-RTHICK)
         I2 = MIN(XLEN, I+RTHICK)
         J1 = MAX(1, J-RTHICK)
         J2 = MIN(YLEN, J+RTHICK)
         DO JJ = J1, J2
         DO II = I1, I2
          IF (.NOT. THICKSOL) THEN
            IF (GRID(II,JJ,LCELLS,SCRN) .LT. 0.5) THEN
              THICKSOL = .TRUE.
            ENDIF
          ENDIF
         ENDDO
         ENDDO
      ENDIF
      IF (THICKSOL) THEN
         R = 0
         GCHAN = 220
         B = NINT(220*F2)
      ENDIF
C Phase overrides.
      PHS = GRID(I,J,PHASE,SCRN)
      IF (ABS(PHS-SECONDARY) .LE. 0.5) THEN
         R = 120
         GCHAN = 40
         B = 40
      ELSE IF (ABS(PHS-DEFECT) .LE. 0.5) THEN
         R = 200
         GCHAN = 40
         B = 0
      ENDIF
      G = GCHAN
      CALL MOVCLAMP(R, G, B)
      RETURN
      END
C
C=====================================================================
C  MOVCLAMP -- keep RGB values in [0,255]
C=====================================================================
      SUBROUTINE MOVCLAMP(R, G, B)
      IMPLICIT NONE
      INTEGER R, G, B
      IF (R .LT. 0) R = 0
      IF (G .LT. 0) G = 0
      IF (B .LT. 0) B = 0
      IF (R .GT. 255) R = 255
      IF (G .GT. 255) G = 255
      IF (B .GT. 255) B = 255
      RETURN
      END
C
C==============last line========================================