subroutine ocean_redi_init(this, grid, nz_ml)
!! Allocate the Phase-A coefficient arrays. Always allocates
!! (configure runs after init); off-state footprint is the six
!! face-shaped (nsurf) coefficient arrays. Plain host allocation
!! (no `do concurrent` before enter_data).
class(ocean_redi_t), intent(inout) :: this
type(hgrid_t), intent(in) :: grid
integer, intent(in), optional :: nz_ml
integer :: nx, ny, nz, ns
nx = grid%nx_total
ny = grid%ny_total
nz = 1
if (present(nz_ml)) nz = nz_ml
if (nz < 1) nz = 1
! Fail loud: the sweep uses NZ_STACK_MAX-sized column-pair locals.
!
! This used to demand `2*nz + 2 <= NZ_STACK_MAX`, which is ~2x
! stricter than anything here needs. The neutral-surface locals
! (`PoLc`/`PoRc`/`KoLc`/`KoRc` at 2*NZ_STACK_MAX+2, `hEc` at
! 2*NZ_STACK_MAX+1) are declared as MULTIPLES of the constant, so
! they already scale with it: nsurf = 2*nz+2 <= 2*NZ_STACK_MAX+2
! whenever nz <= NZ_STACK_MAX. The binding locals are the
! plain-`NZ_STACK_MAX` column arrays (`hL`/`tcL`/`scL`/...), which
! need `nz`, and the `NZ_STACK_MAX+1` interface arrays, which need
! `nz+1`. So the requirement is the house-wide `nz_stack_required`.
if (.not. nz_stack_is_sufficient(nz)) then
error stop "ocean_redi_init: nz_ml+1 exceeds NZ_STACK_MAX; "// &
"raise -DRDB_NZ_STACK_MAX or disable Redi"
end if
this%nx_total = nx
this%ny_total = ny
this%nz_ml = nz
ns = 2*nz + 2
this%nsurf = ns
allocate (this%uPoL(nx + 1, ny, ns), source=0.0_wp)
allocate (this%uPoR(nx + 1, ny, ns), source=0.0_wp)
allocate (this%uKoL(nx + 1, ny, ns), source=1)
allocate (this%uKoR(nx + 1, ny, ns), source=1)
allocate (this%uhEff(nx + 1, ny, ns - 1), source=0.0_wp)
allocate (this%vPoL(nx, ny + 1, ns), source=0.0_wp)
allocate (this%vPoR(nx, ny + 1, ns), source=0.0_wp)
allocate (this%vKoL(nx, ny + 1, ns), source=1)
allocate (this%vKoR(nx, ny + 1, ns), source=1)
allocate (this%vhEff(nx, ny + 1, ns - 1), source=0.0_wp)
allocate (this%uKb(nx + 1, ny), source=1)
allocate (this%uKt(nx + 1, ny), source=nz)
allocate (this%vKb(nx, ny + 1), source=1)
allocate (this%vKt(nx, ny + 1), source=nz)
allocate (this%khtr_u(nx + 1, ny), source=0.0_wp)
allocate (this%khtr_v(nx, ny + 1), source=0.0_wp)
allocate (this%tr_snap(nx, ny, this%nz_ml), source=0.0_wp)
this%is_init = .true.
end subroutine ocean_redi_init