ice_ride_update_x_impl Subroutine

   pure subroutine ice_ride_update_x_impl(iareaT, uh, tr_flux_x_work, mca, val, dt_adv, &
                                          ncat, nx, ny)
      !! Cell-update pass (race-free): reads `tr_flux_x_work` (the
      !! GATHERED donor VALUE at each face, from `ice_gather_flux_x_impl`)
      !! + its OWN cell's `mca`/`val` (never a neighbour's `val`), applies
      !! SIS2's `advect_scalar_x` flux-form update + the `H_NEGLECT`
      !! conditioning guard (`SIS_tracer_advect.F90:735-760`), writes
      !! `val(i,j,c)`.  `hnew` is the SAME expression the mass-update
      !! kernel (d) uses, so the implied masses agree bitwise.
      !! Bitwise no-op when both faces carry zero flux (`F_W==F_E==0`).
      !!
      !! Per-area transcription of SIS2's cell-integrated algebra
      !! (SIS2 works in `hprev = mca_old*areaT`, `uhh = flux*dt`; dividing
      !! every SIS2 quantity by `areaT` gives the per-area form below —
      !! `dtI = dt_adv*iareaT` plays the role of SIS2's `dt/areaT`):
      !!   `hlst = max(mca_old, 0)`; `hnew = mca_old - dtI*(fe-fw)`.
      !!   `hnew <= 0`            : hold (mass gone).
      !!   `0 < hnew < H_NEGLECT`  : `h_add = H_NEGLECT - hnew`;
      !!     `I_htot = 1/(hlst + dtI*(|fe|+|fw|))` (0 if the denominator
      !!     is 0); `hlst_adj = hlst + h_add*hlst*I_htot`;
      !!     `haddE = h_add*dtI*|fe|*I_htot`, `haddW = h_add*dtI*|fw|*I_htot`
      !!     (SIS2's sign convention: `haddE` is SUBTRACTED from the east
      !!     outflow term, `haddW` is ADDED to the west inflow term — both
      !!     push the effective in/outflow toward "more mass stays");
      !!     `val_new = (val*hlst_adj - ((fe*dtI-haddE)*val_e -
      !!     (fw*dtI+haddW)*val_w)) / H_NEGLECT`.
      !!   `hnew >= H_NEGLECT`     : plain flux form,
      !!     `val_new = (val*mca_old - dtI*(fe*val_e-fw*val_w))/hnew`.
      integer, intent(in) :: ncat, nx, ny
      real(wp), intent(in) :: iareaT(nx, ny)
      real(wp), intent(in) :: uh(nx + 1, ny, ncat)
      real(wp), intent(in) :: tr_flux_x_work(nx + 1, ny, ncat)
      real(wp), intent(in) :: mca(nx, ny, ncat)
      real(wp), intent(inout) :: val(nx, ny, ncat)
      real(wp), intent(in) :: dt_adv
      integer :: i, j, c
      real(wp) :: fw, fe, val_w, val_e, mca_old, hlst, dti, hnew, h_add, denom, i_htot
      real(wp) :: hlst_adj, haddw, hadde, fw_term, fe_term

      do concurrent(j=1:ny, i=1:nx, c=1:ncat) &
         local(fw, fe, val_w, val_e, mca_old, hlst, dti, hnew, h_add, denom, i_htot, &
               hlst_adj, haddw, hadde, fw_term, fe_term)
         fw = uh(i, j, c)
         fe = uh(i + 1, j, c)
         if (fw /= 0.0_wp .or. fe /= 0.0_wp) then
            val_w = tr_flux_x_work(i, j, c)
            val_e = tr_flux_x_work(i + 1, j, c)
            mca_old = mca(i, j, c)
            hlst = max(mca_old, 0.0_wp)
            dti = dt_adv*iareaT(i, j)
            hnew = mca_old - dti*(fe - fw)
            if (hnew <= 0.0_wp) then
               ! mass gone; hold val (inert — PR 4a massless convention)
               continue
            else if (hnew < H_NEGLECT_ICE_TRANSPORT) then
               h_add = H_NEGLECT_ICE_TRANSPORT - hnew
               denom = hlst + dti*(abs(fe) + abs(fw))
               if (denom > 0.0_wp) then
                  i_htot = 1.0_wp/denom
               else
                  i_htot = 0.0_wp
               end if
               hlst_adj = hlst + h_add*hlst*i_htot
               haddw = h_add*dti*abs(fw)*i_htot
               hadde = h_add*dti*abs(fe)*i_htot
               fe_term = fe*dti - hadde
               fw_term = fw*dti + haddw
               val(i, j, c) = (val(i, j, c)*hlst_adj - (fe_term*val_e - fw_term*val_w)) &
                              /H_NEGLECT_ICE_TRANSPORT
            else
               val(i, j, c) = (val(i, j, c)*mca_old - dti*(fe*val_e - fw*val_w))/hnew
            end if
         end if
      end do
   end subroutine ice_ride_update_x_impl