ice_evp_dynamics_impl Subroutine

   subroutine ice_evp_dynamics_impl(grid, metrics, f_corner, mis, mice, ci, uo, vo, &
                                    tau_ax, tau_ay, ui, vi, str_d, str_t, str_s, &
                                    fxoc, fyoc, dt_slow, par, periodic_x, periodic_y, nx, ny, &
                                    mis_w, mice_w, ci_w, pres_mice_w, del_sh_min_pr_w, &
                                    sh_dd_w, sh_dt_w, zeta_w, del_sh_w, mask_t_w, &
                                    mi_u_w, mask_u_w, u_tmp_w, mi_v_w, mask_v_w, a_u_w, a_v_w, &
                                    sh_ds_w, mi_ratio_a_q_w, q_w, mask_q_w, &
                                    halo_x, halo_y, land_w, land_e, land_s, land_n, &
                                    dt_transport, n_trunc)
      !! Flat-impl core of `ice_evp_dynamics`: the EVP subcycle body with
      !! the persistent scratch passed as EXPLICIT-SHAPE dummies (memory:
      !! never assumed-shape into a `do concurrent` feeder — NVHPC would
      !! emit descriptor-walk memcpys per launch).  The `*_w` scratch names
      !! mirror the retired module workspace 1:1, so the body below is
      !! unchanged from the pre-slot version.
      !!
      !! **Ghost policy is SPLIT — not "all owned here".**  This routine
      !! owns the ghost policy (periodic wrap or zero) for the workspace
      !! copies of `mis`/`mice`/`ci`, for `ui`/`vi`, and for
      !! `str_d`/`str_t`/`str_s` — it wraps those below.  The `intent(in)`
      !! `uo`/`vo`/`tau_ax`/`tau_ay` are consumed DIRECTLY (the momentum
      !! kernel reads their ghost rows, e.g. `uo(i,j)`, `vo(i-1,j+1)`), so
      !! their ghosts are the CALLER's responsibility: a test in a periodic
      !! config must fill those four with their ghosts already wrapped.
      type(hgrid_t), intent(in) :: grid
      type(ocean_metrics_t), intent(in) :: metrics
      real(wp), intent(in) :: f_corner(:, :)
      real(wp), intent(in) :: mis(:, :), mice(:, :), ci(:, :)
      real(wp), intent(in) :: uo(:, :), vo(:, :)
      real(wp), intent(in) :: tau_ax(:, :), tau_ay(:, :)
      real(wp), intent(inout) :: ui(:, :), vi(:, :)
      real(wp), intent(inout) :: str_d(:, :), str_t(:, :), str_s(:, :)
      real(wp), intent(inout) :: fxoc(:, :), fyoc(:, :)
      real(wp), intent(in) :: dt_slow
      type(ice_evp_params_t), intent(in) :: par
      logical, intent(in) :: periodic_x, periodic_y
      integer, intent(in) :: nx, ny
         !! T-cell extents (declared before the explicit-shape scratch that
         !! uses them — decl-order rule).
      real(wp), intent(inout) :: mis_w(nx, ny), mice_w(nx, ny), ci_w(nx, ny)
      real(wp), intent(inout) :: pres_mice_w(nx, ny), del_sh_min_pr_w(nx, ny)
      real(wp), intent(inout) :: sh_dd_w(nx, ny), sh_dt_w(nx, ny)
      real(wp), intent(inout) :: zeta_w(nx, ny), del_sh_w(nx, ny)
      real(wp), intent(inout) :: mask_t_w(nx, ny)
      real(wp), intent(inout) :: mi_u_w(nx + 1, ny), mask_u_w(nx + 1, ny), u_tmp_w(nx + 1, ny)
      real(wp), intent(inout) :: mi_v_w(nx, ny + 1), mask_v_w(nx, ny + 1)
      real(wp), intent(inout) :: a_u_w(nx + 1, ny), a_v_w(nx, ny + 1)
         !! PR 62: face ice concentration, valid ONLY when `par%a_face_stress`
         !! (uninitialised device memory otherwise — never read off-gate).
      real(wp), intent(inout) :: sh_ds_w(nx + 1, ny + 1), mi_ratio_a_q_w(nx + 1, ny + 1)
      real(wp), intent(inout) :: q_w(nx + 1, ny + 1), mask_q_w(nx + 1, ny + 1)
      logical, intent(in) :: halo_x, halo_y
         !! Axis split across ranks (see `ice_evp_dynamics`).
      logical, intent(in) :: land_w, land_e, land_s, land_n
         !! Ghost band pinned to land (see `ice_evp_dynamics`).
      real(wp), intent(in), optional :: dt_transport
         !! PR 36: the dt TRANSPORT will actually consume
         !! (`ocean_dyn%therm_dt(dt)`), NOT this call's `dt_slow` — EVP
         !! runs every outer step, transport at thermo cadence
         !! (`rdb_driver.F90`). Absent => `dt_slow` (SIS2's own assumption:
         !! `SIS_C_dynamics` and `SIS_transport` share `dt_slow`). Read
         !! only when `par%cfl_trunc > 0`.
      integer, intent(out), optional :: n_trunc
         !! PR 36: count of ice-bearing faces (`mi_u`/`mi_v > m_neglect`)
         !! the FINAL clip touched. `0` when `par%cfl_trunc <= 0`.

      integer :: nx_phys, ny_phys, nghost, n
      real(wp) :: dt, tdamp_eff, dt_2tdamp, i_1pdt_t, ec2, i_ec2
      real(wp) :: cdrho, i_cdrhodt, p0_rho, m_neglect, m_neglect2, m_neglect4
      real(wp) :: dt_tr, dt_cum
      logical :: a_face_on, do_trunc_its, do_trunc_fin
      logical :: wrap_x, wrap_y
      integer :: n_out

      nx_phys = grid%nx_phys
      ny_phys = grid%ny_phys
      nghost = grid%nghost
      ! Local periodic wrap only on an axis the halo does not own (F3):
      ! on a split axis the wrap would copy THIS tile's own interior into
      ! a ghost band the neighbour rank owns.
      wrap_x = periodic_x .and. .not. halo_x
      wrap_y = periodic_y .and. .not. halo_y

      ! ---- Scalar precompute (hoisted before the substep loop) ----
      a_face_on = par%a_face_stress
      dt = dt_slow/real(par%evp_sub_steps, wp)
      if (par%tdamp > 0.0_wp) then
         tdamp_eff = par%tdamp
      else if (par%tdamp == 0.0_wp) then
         tdamp_eff = max(0.2_wp*dt_slow, 3.0_wp*dt)
      else
         tdamp_eff = max(-par%tdamp*dt_slow, 3.0_wp*dt)
      end if
      dt_2tdamp = dt/(2.0_wp*tdamp_eff)
      ec2 = par%ec*par%ec
      i_ec2 = 0.0_wp
      if (ec2 > 0.0_wp) i_ec2 = 1.0_wp/ec2
      i_1pdt_t = 1.0_wp/(1.0_wp + dt_2tdamp)
      cdrho = par%cdw*par%rho_ocean
      i_cdrhodt = 1.0_wp/(par%cdw*par%rho_ocean*dt)
      p0_rho = par%p0/ICE_RHO_ICE
      m_neglect = ICE_RHO_ICE*M_NEGLECT_FACTOR
      m_neglect2 = m_neglect*m_neglect
      m_neglect4 = m_neglect2*m_neglect2

      ! ---- PR 36: CFL-truncation gates + the dt the bound must use.
      ! The bound is against the dt TRANSPORT will consume, NOT this call's
      ! `dt_slow` -- Roundabout decouples EVP (every outer step) from transport
      ! (thermo cadence); SIS2's structure assumes they are the same dt.
      ! Absent `dt_transport` => `dt_slow`, matching SIS2's own assumption. ----
      dt_tr = dt_slow
      if (present(dt_transport)) dt_tr = dt_transport
      do_trunc_its = par%cfl_trunc_dyn_its .and. (par%cfl_trunc > 0.0_wp) .and. (dt_tr > 0.0_wp)
      do_trunc_fin = (par%cfl_trunc > 0.0_wp) .and. (dt_tr > 0.0_wp)

      ! ---- Effective masks (SIS2 mask2dT/Cu/Cv/Bu), built once ----
      call evp_build_masks_impl(metrics%wet_T, mask_t_w, mask_u_w, mask_v_w, mask_q_w, &
                                nx_phys, ny_phys, nghost, wrap_x, wrap_y, &
                                land_w, land_e, land_s, land_n, nx, ny)

      ! ---- Category fields into the workspace, ghost-wrapped/zeroed ----
      call evp_fill_cell_fields_impl(mask_t_w, mis, mice, ci, mis_w, mice_w, ci_w, nx, ny)
      call ocean_periodic_wrap_centre_2d(mis_w, nx, ny, nx_phys, ny_phys, nghost, &
                                         wrap_x, wrap_y)
      call ocean_periodic_wrap_centre_2d(mice_w, nx, ny, nx_phys, ny_phys, nghost, &
                                         wrap_x, wrap_y)
      call ocean_periodic_wrap_centre_2d(ci_w, nx, ny, nx_phys, ny_phys, nghost, &
                                         wrap_x, wrap_y)

      ! ---- Zero ice velocities with no mass (SIS2 :899-907) ----
      call evp_zero_massless_velocity_impl(mask_u_w, mask_v_w, mis_w, ui, vi, nx, ny)

      ! ---- pres_mice + del_sh_min_pr precompute (:877-890) ----
      call evp_pres_mice_impl(metrics%dxT, metrics%dyT, ci_w, p0_rho, par%c0, &
                              par%del_sh_min_scale, tdamp_eff, dt, pres_mice_w, &
                              del_sh_min_pr_w, nx, ny)

      ! ---- mi_u / mi_v (:967-974) ----
      call evp_mi_face_impl(mis_w, mi_u_w, mi_v_w, nx, ny)

      ! ---- PR 62: face ice concentration a_u/a_v, ONLY when a_face_stress.
      ! Reuses evp_mi_face_impl's 0.5-face-average — same expression/edge
      ! convention `ice_ocean_stress_flux_impl`'s `a_u` already uses, so the
      ! momentum budget it closes matches the coupler bit-for-bit (§5.4). ----
      if (a_face_on) call evp_mi_face_impl(ci_w, a_u_w, a_v_w, nx, ny)

      ! ---- q + mi_ratio_A_q (:926-982) ----
      call evp_q_and_mi_ratio_impl(metrics%areaT, f_corner, mask_t_w, mask_u_w, mask_v_w, &
                                   mask_q_w, mis_w, m_neglect, m_neglect2, m_neglect4, &
                                   q_w, mi_ratio_a_q_w, nx, ny)

      ! ---- limit_stresses ONCE before the substep loop (:896) — req (2) ----
      call ice_limit_stresses(metrics%areaT, mask_t_w, pres_mice_w, mice_w, &
                              str_d, str_t, str_s, par%ec, nx, ny)
      call ocean_periodic_wrap_centre_2d(str_d, nx, ny, nx_phys, ny_phys, nghost, &
                                         wrap_x, wrap_y)
      call ocean_periodic_wrap_centre_2d(str_t, nx, ny, nx_phys, ny_phys, nghost, &
                                         wrap_x, wrap_y)
      call evp_wrap_corner_impl(str_s, nx + 1, ny + 1, nx_phys, ny_phys, nghost, &
                                wrap_x, wrap_y)

      ! ---- Zero the subcycle-averaged ice->ocean stress (F1: an explicit
      ! device kernel, NOT a bare host whole-array assignment — fxoc/fyoc
      ! are copyin-mapped, so a host `= 0.0_wp` would leave the DEVICE copy
      ! stale and each outer call would accumulate onto the prior call's
      ! average). ----
      call evp_zero_stress_impl(fxoc, fyoc, nx, ny)

      ! ---- The EVP subcycle loop — device-resident, no per-substep H<->D ----
      dt_cum = 0.0_wp
      do n = 1, par%evp_sub_steps
         ! PR 36: cumulative elapsed time within THIS call, at the TOP of the
         ! loop (SIS2 :1023,1028) => subcycle n sees t_cum = n*dt. Host
         ! scalar by value -- the substep loop is host-driven, so this costs
         ! nothing and breaks no device residency.
         dt_cum = dt_cum + dt

         ! X2: seam ghosts of the ice velocity from the neighbour ranks
         ! (two-pass, so corner ghosts are valid too), then the local wrap
         ! on any periodic axis the halo does not own.
         if (halo_x .or. halo_y) then
            call ocean_halo_face_x(ui)
            call ocean_halo_face_y(vi)
         end if
         call ocean_periodic_wrap_face_x_2d(ui, nx + 1, ny, nx_phys, ny_phys, nghost, &
                                            wrap_x, wrap_y)
         call ocean_periodic_wrap_face_y_2d(vi, nx, ny + 1, nx_phys, ny_phys, nghost, &
                                            wrap_x, wrap_y)

         call evp_sh_ds_impl(metrics%dx_dyBu, metrics%dy_dxBu, metrics%idxCu, metrics%idyCv, &
                             mask_q_w, ui, vi, sh_ds_w, nx, ny)
         call evp_sh_dd_dt_impl(metrics%dy_dxT, metrics%dx_dyT, metrics%iareaT, &
                                metrics%idyCu, metrics%idxCv, metrics%dyCu, metrics%dxCv, &
                                ui, vi, sh_dd_w, sh_dt_w, nx, ny)

         ! ---- PR 36: PROJECT_ICE_CONCENTRATION -- SIS2's position exactly
         ! (:1077 -> :1082): after sh_Dd, before zeta. ci_w is the CALL's
         ! initial (entry-gathered) concentration, never mutated by this --
         ! keep it that way (the projection is `ci*exp(-t_cum*sh_dd)`, not
         ! an incremental accumulator). del_sh_min_pr_w is NOT recomputed
         ! (no ci dependence). ----
         if (par%project_ci) then
            call evp_project_ci_impl(ci_w, sh_dd_w, dt_cum, p0_rho, par%c0, pres_mice_w, nx, ny)
         end if

         call evp_zeta_impl(sh_dd_w, sh_dt_w, sh_ds_w, i_ec2, pres_mice_w, mice_w, &
                            del_sh_min_pr_w, del_sh_w, zeta_w, nx, ny)
         call evp_stress_relax_impl(zeta_w, sh_dd_w, sh_dt_w, pres_mice_w, mice_w, &
                                    i_1pdt_t, dt_2tdamp, i_ec2, str_d, str_t, nx, ny)
         call evp_str_s_relax_impl(metrics%areaT, zeta_w, sh_ds_w, mi_ratio_a_q_w, &
                                   i_1pdt_t, dt_2tdamp, i_ec2, str_s, nx, ny)

         call evp_copy_u_impl(ui, u_tmp_w, nx, ny)

         call evp_u_momentum_impl(metrics%idxCu, metrics%idyCu, metrics%dy2h, &
                                  metrics%dx2q, metrics%iareaCu, mask_u_w, mi_u_w, mi_v_w, &
                                  q_w, str_d, str_t, str_s, uo, vo, tau_ax, ui, vi, &
                                  fxoc, m_neglect, i_cdrhodt, cdrho, dt, nx_phys, ny_phys, &
                                  nghost, nx, ny, a_u_w, a_face_on)
         call evp_v_momentum_impl(metrics%idyCv, metrics%idxCv, metrics%dx2h, &
                                  metrics%dy2q, metrics%iareaCv, mask_v_w, mi_v_w, mi_u_w, &
                                  q_w, str_d, str_t, str_s, uo, vo, tau_ay, u_tmp_w, vi, &
                                  fyoc, m_neglect, i_cdrhodt, cdrho, dt, nx_phys, ny_phys, &
                                  nghost, nx, ny, a_v_w, a_face_on)

         ! ---- PR 36: in-loop CFL clip (cfl_trunc_dyn_its) -- SIS2's
         ! position (:1338, bottom of the loop, after both momentum
         ! solves). Exact bound, no count. No re-wrap needed: the next
         ! iteration's wrap at the top of the loop does it. ----
         if (do_trunc_its) then
            call evp_truncate_velocity_impl(metrics%areaT, metrics%dy_cu, metrics%dx_cv, &
                                            ui, vi, par%cfl_trunc, dt_tr, 1.0_wp, &
                                            nghost, nx_phys, ny_phys, nx, ny)
         end if
      end do

      ! ---- PR 36: FINAL CFL clip (cfl_trunc) -- always on when cfl_trunc>0
      ! and dt_tr>0. 0.95*bound back-off; counts ice-bearing faces touched. ----
      n_out = 0
      if (do_trunc_fin) then
         call evp_truncate_final_impl(metrics%areaT, metrics%dy_cu, metrics%dx_cv, &
                                      mi_u_w, mi_v_w, ui, vi, par%cfl_trunc, dt_tr, &
                                      m_neglect, nghost, nx_phys, ny_phys, nx, ny, n_out, &
                                      count_w=land_w, count_s=land_s)
      end if
      ! X3: the last subcycle's momentum solve (and the clip) wrote PHYSICAL
      ! faces only, so every ghost face still holds the value exchanged at
      ! the TOP of that subcycle -- one subcycle old.  Transport reads them
      ! as donor velocities (and, in the ghost rows, as the x-pass's own
      ! face velocities), so refresh them on EVERY call, clip or not: seam
      ! ghosts first, then the local wrap, same order as X2.  (A periodic
      ! ghost face that kept its unclipped value would also make transport
      ! abort for a reason no test would name.)
      if (halo_x .or. halo_y) then
         call ocean_halo_face_x(ui)
         call ocean_halo_face_y(vi)
      end if
      call ocean_periodic_wrap_face_x_2d(ui, nx + 1, ny, nx_phys, ny_phys, nghost, &
                                         wrap_x, wrap_y)
      call ocean_periodic_wrap_face_y_2d(vi, nx, ny + 1, nx_phys, ny_phys, nghost, &
                                         wrap_x, wrap_y)
      if (present(n_trunc)) n_trunc = n_out

      ! ---- fxoc/fyoc average + mask (:1415-1441) ----
      call evp_average_stress_impl(mask_u_w, mask_v_w, fxoc, fyoc, par%evp_sub_steps, nx, ny)
   end subroutine ice_evp_dynamics_impl