roquet_spv_point Subroutine

   pure subroutine roquet_spv_point(T_pt, S_sp, p, sv, dsv_dt_model, dsv_ds_model)
      !! Fused Roquet et al. (2015) SpV evaluation at a point, in MODEL
      !! variables (potential temperature `T_pt` degC, practical salinity
      !! `S_sp` PSU, pressure `p` Pa).  Returns specific volume `sv`
      !! (m^3/kg) and the analytic sensitivities w.r.t. the MODEL
      !! variables (`dsv_dt_model` = dSV/dPT, `dsv_ds_model` = dSV/dSP)
      !! so the variant-agnostic consumers (which work in PT, SP) get the
      !! correct chain-ruled derivatives.
      !!
      !! Convention (pinned 2026-06-16 — see EOS_VARIANT_ROQUET_SPV):
      !!   SR = S_sp * (35.16504/35)         [Reference Salinity, g/kg]
      !!   CT = ct_from_pt(SR, T_pt)         [Conservative Temperature]
      !! The polynomial is fit in (CT, SA, p); we feed SR for SA (bounded
      !! anomaly deviation) and CT via the local conversion poly.
      !!
      !! Chain rule back to the model variables.  CT = ct_from_pt(SR, PT)
      !! depends on BOTH PT and SR (= SP·factor), so the SP derivative
      !! carries TWO routes into SV — the direct SA route AND the
      !! CT-via-SR route:
      !!   dSV/dPT = (dSV/dCT) · (dCT/dPT)
      !!   dSV/dSP = [ (dSV/dSA) + (dSV/dCT)·(dCT/dSR) ] · (35.16504/35)
      !! The CT-via-SR term is ~0.5 % of dSV/dSP (dCT/dSR ≈ −0.02 degC per
      !! g/kg); dropping it (the simplified spec formula) leaves a real
      !! ~5e-3 relative error against the true total derivative, so we
      !! keep the full chain — this is what the FD regression locks.
      !!
      !! dCT/dPT and dCT/dSR are analytic derivatives of the ct_from_pt
      !! poly: with yy = 0.025·PT, x2 = sfac·SR, xx = sqrt(x2), each
      !! yy-coefficient c_k(xx) = A + B·xx² + C·xx³ + D·xx⁴ + E·xx⁵, so
      !!   dCT/dPT = 0.025·(dhh/dyy)/cp0
      !!   dCT/dSR = (dhh/dxx)·(sfac/(2·xx))/cp0.
      !!
      !! Transcribed from the verified prototype roquet_spv_eos.py.  One
      !! sqrt for zs (shared by SV and both derivatives) + one sqrt for
      !! the ct_from_pt poly normalisation.
      !$acc routine seq
      real(wp), intent(in) :: T_pt, S_sp, p
      real(wp), intent(out) :: sv, dsv_dt_model, dsv_ds_model

      real(wp) :: SR, CT, dct_dpt, dct_dsr
      real(wp) :: zt, zs, zp
      real(wp) :: sv_ts0, sv_ts1, sv_ts2, sv_ts3, sv_0s0, sv_00p
      real(wp) :: dvdzt0, dvdzt1, dvdzt2, dvdzt3, dsv_dct
      real(wp) :: dvdzs0, dvdzs1, dvdzs2, dvdzs3, dsv_dsa
      real(wp) :: x2, xx, yy, hh, dh_dy, dh_dx
      real(wp) :: c0, c1, c2, c3, c4, c5, c6, c7
      real(wp) :: d0, d1, d2, d3, d4, d5, d6, d7

      SR = S_sp*ROQ_SR_FACTOR

      ! CT = ct_from_pt(SR, PT): the 7-term gsw surface poly is degree 7 in
      ! yy = 0.025*PT.  Each yy-coefficient c_k is a polynomial in xx =
      ! sqrt(sfac*SR): c_k = A + B*xx^2 + C*xx^3 + D*xx^4 + E*xx^5 (x2 = xx^2).
      ! Building c_k and dc_k/dxx (= d_k) gives exact analytic dCT/dPT and
      ! dCT/dSR by Horner — no risk of mis-differentiating the published
      ! nesting.
      ! Floor x2 to a tiny positive so the `dct_dsr = dh_dx/(2*xx)` below
      ! never hits 0/0 at SP=0 (fresh water): dh_dx is itself proportional
      ! to xx (lowest term 2B*xx), so dh_dx/xx has a finite limit — the
      ! floor reproduces it (xx tiny-but-nonzero) instead of NaN.  The
      ! floor only bites at SP < ~1e-9 PSU, so it is bit-identical for any
      ! real-ocean salinity.
      x2 = max(ROQ_CT_SFAC*SR, 1.0e-20_wp)
      xx = sqrt(x2)
      yy = T_pt*0.025_wp
      c0 = 61.01362420681071_wp &
           + x2*(268.5520265845071_wp &
                 + xx*(937.2099110620707_wp &
                       + xx*(-1687.914374187449_wp + xx*246.9598888781377_wp)))
      c1 = 168776.46138048015_wp &
           + x2*(-12019.028203559312_wp &
                 + xx*(588.1802812170108_wp &
                       + xx*(936.3206544460336_wp + xx*123.59576582457964_wp)))
      c2 = -2735.2785605119625_wp &
           + x2*(3734.858026725145_wp &
                 + xx*(248.39476522971285_wp &
                       + xx*(-942.7827304544439_wp + xx*(-48.5891069025409_wp))))
      c3 = 2574.2164453821433_wp &
           + x2*(-2046.7671145057618_wp &
                 + xx*(-3.871557904936333_wp + xx*369.4389437509002_wp))
      c4 = -1536.6644434977543_wp &
           + x2*(465.28655623126450_wp + xx*(-2.6268019854268356_wp + xx*(-33.83664947895248_wp)))
      c5 = 545.7340497931629_wp &
           + x2*(-0.6370820302831379_wp + xx*(-9.987880382780322_wp))
      c6 = -50.91091728474331_wp + x2*(-10.650848542359153_wp)
      c7 = -18.30489878927802_wp
      ! dc_k/dxx = 2B*xx + 3C*xx^2 + 4D*xx^3 + 5E*xx^4 (A and the xx^0/xx^1
      ! terms vanish; B,C,D,E read off the c_k expansions above).
      d0 = xx*(2.0_wp*268.5520265845071_wp &
               + xx*(3.0_wp*937.2099110620707_wp &
                     + xx*(4.0_wp*(-1687.914374187449_wp) + xx*5.0_wp*246.9598888781377_wp)))
      d1 = xx*(2.0_wp*(-12019.028203559312_wp) &
               + xx*(3.0_wp*588.1802812170108_wp &
                     + xx*(4.0_wp*936.3206544460336_wp + xx*5.0_wp*123.59576582457964_wp)))
      d2 = xx*(2.0_wp*3734.858026725145_wp &
               + xx*(3.0_wp*248.39476522971285_wp &
                     + xx*(4.0_wp*(-942.7827304544439_wp) + xx*5.0_wp*(-48.5891069025409_wp))))
      d3 = xx*(2.0_wp*(-2046.7671145057618_wp) &
               + xx*(3.0_wp*(-3.871557904936333_wp) + xx*4.0_wp*369.4389437509002_wp))
      d4 = xx*(2.0_wp*465.28655623126450_wp &
               + xx*(3.0_wp*(-2.6268019854268356_wp) + xx*4.0_wp*(-33.83664947895248_wp)))
      d5 = xx*(2.0_wp*(-0.6370820302831379_wp) + xx*3.0_wp*(-9.987880382780322_wp))
      d6 = xx*2.0_wp*(-10.650848542359153_wp)
      d7 = 0.0_wp
      hh = c0 + yy*(c1 + yy*(c2 + yy*(c3 + yy*(c4 + yy*(c5 + yy*(c6 + yy*c7))))))
      dh_dy = c1 + yy*(2.0_wp*c2 + yy*(3.0_wp*c3 + yy*(4.0_wp*c4 &
                                                       + yy*(5.0_wp*c5 + yy*(6.0_wp*c6 + yy*7.0_wp*c7)))))
      dh_dx = d0 + yy*(d1 + yy*(d2 + yy*(d3 + yy*(d4 + yy*(d5 + yy*(d6 + yy*d7))))))
      CT = hh/ROQ_CP0
      dct_dpt = 0.025_wp*dh_dy/ROQ_CP0
      dct_dsr = dh_dx*(ROQ_CT_SFAC/(2.0_wp*xx))/ROQ_CP0

      zt = CT
      zs = sqrt(abs(S_sp*ROQ_SR_FACTOR + ROQ_RDELTAS)*ROQ_R1_S0)
      zp = p

      ! --- specific volume SV(zs, zt, zp) ---
      sv_ts3 = SPV003 + (zs*SPV103 + zt*SPV013)
      sv_ts2 = SPV002 + (zs*(SPV102 + zs*SPV202) &
                         + zt*(SPV012 + (zs*SPV112 + zt*SPV022)))
      sv_ts1 = SPV001 + (zs*(SPV101 + zs*(SPV201 + zs*(SPV301 + zs*SPV401))) &
                         + zt*(SPV011 + (zs*(SPV111 + zs*(SPV211 + zs*SPV311)) &
                                         + zt*(SPV021 + (zs*(SPV121 + zs*SPV221) &
                                                         + zt*(SPV031 + (zs*SPV131 + zt*SPV041)))))))
      sv_ts0 = zt*(SPV010 &
                   + (zs*(SPV110 + zs*(SPV210 + zs*(SPV310 + zs*(SPV410 + zs*SPV510)))) &
                      + zt*(SPV020 + (zs*(SPV120 + zs*(SPV220 + zs*(SPV320 + zs*SPV420))) &
                                      + zt*(SPV030 + (zs*(SPV130 + zs*(SPV230 + zs*SPV330)) &
                                                      + zt*(SPV040 + (zs*(SPV140 + zs*SPV240) &
                                                                      + zt*(SPV050 + (zs*SPV150 + zt*SPV060))))))))))
      sv_0s0 = SPV000 + zs*(SPV100 + zs*(SPV200 + zs*(SPV300 + zs*(SPV400 &
                                                                   + zs*(SPV500 + zs*SPV600)))))
      sv_00p = zp*(ROQ_V00 + zp*(ROQ_V01 + zp*(ROQ_V02 + zp*(ROQ_V03 &
                                                             + zp*(ROQ_V04 + zp*ROQ_V05)))))
      sv = ((sv_ts0 + sv_0s0) + zp*(sv_ts1 + zp*(sv_ts2 + zp*sv_ts3))) + sv_00p

      ! --- dSV/dCT ---
      dvdzt3 = ALP003
      dvdzt2 = ALP002 + (zs*ALP102 + zt*ALP012)
      dvdzt1 = ALP001 + (zs*(ALP101 + zs*(ALP201 + zs*ALP301)) &
                         + zt*(ALP011 + (zs*(ALP111 + zs*ALP211) &
                                         + zt*(ALP021 + (zs*ALP121 + zt*ALP031)))))
      dvdzt0 = ALP000 + (zs*(ALP100 + zs*(ALP200 + zs*(ALP300 + zs*(ALP400 + zs*ALP500)))) &
                         + zt*(ALP010 + (zs*(ALP110 + zs*(ALP210 + zs*(ALP310 + zs*ALP410))) &
                                         + zt*(ALP020 + (zs*(ALP120 + zs*(ALP220 + zs*ALP320)) &
                                                         + zt*(ALP030 + (zt*(ALP040 + (zs*ALP140 + zt*ALP050)) &
                                                                         + zs*(ALP130 + zs*ALP230))))))))
      dsv_dct = dvdzt0 + zp*(dvdzt1 + zp*(dvdzt2 + zp*dvdzt3))

      ! --- dSV/dSA (per-coef 0.5*r1_S0 folded into BET; residual /zs here) ---
      dvdzs3 = BET003
      dvdzs2 = BET002 + (zs*BET102 + zt*BET012)
      dvdzs1 = BET001 + (zs*(BET101 + zs*(BET201 + zs*BET301)) &
                         + zt*(BET011 + (zs*(BET111 + zs*BET211) &
                                         + zt*(BET021 + (zs*BET121 + zt*BET031)))))
      dvdzs0 = BET000 + (zs*(BET100 + zs*(BET200 + zs*(BET300 + zs*(BET400 + zs*BET500)))) &
                         + zt*(BET010 + (zs*(BET110 + zs*(BET210 + zs*(BET310 + zs*BET410))) &
                                         + zt*(BET020 + (zs*(BET120 + zs*(BET220 + zs*BET320)) &
                                                         + zt*(BET030 + (zt*(BET040 + (zs*BET140 + zt*BET050)) &
                                                                         + zs*(BET130 + zs*BET230))))))))
      dsv_dsa = (dvdzs0 + zp*(dvdzs1 + zp*(dvdzs2 + zp*dvdzs3)))/zs

      ! Chain rule to the model variables (PT, SP).  The SP route includes
      ! the CT-via-SR coupling (CT depends on SR = SP·factor).
      dsv_dt_model = dsv_dct*dct_dpt
      dsv_ds_model = (dsv_dsa + dsv_dct*dct_dsr)*ROQ_SR_FACTOR
   end subroutine roquet_spv_point