ocean_fold_unpack_3d Subroutine

   subroutine ocean_fold_unpack_3d(fld, nxa, nya, nz, stagger, negate, device_resident)
      !! Unpack the next field of the exchanged group: write every north
      !! ghost row (all storage columns) and, for v / corner, the fold-line
      !! row's west half (self-conjugate column → 0 for a vector).
      integer, intent(in) :: nxa, nya, nz
         !! Storage extents of `fld`.
      real(wp), intent(inout) :: fld(nxa, nya, nz)
         !! Field; same stagger and order as at pack.
      integer, intent(in) :: stagger
         !! `FOLD_STAG_*`.
      logical, intent(in) :: negate
         !! `.true.` for a true-vector component.
      logical, intent(in), optional :: device_resident
         !! `.false.` ⇒ host arrays; default device-resident.

      integer :: fam, nrow, dbase, g, g0, ng_e, L, r, p, cap, oT, oU, ib, cls
      logical :: on_device, has_row
      real(wp) :: sgn, val

      if (.not. fx_active) return
      if (.not. fx_open .or. .not. fx_sent) then
         error stop "rdb_ocean_fold_exchange: ocean_fold_unpack before ocean_fold_exchange"
      end if
      on_device = .true.
      if (present(device_resident)) on_device = device_resident

      fam = fold_stagger_family(stagger)
      nrow = fold_stagger_nrows(stagger, fx_ng)
      has_row = (nrow == fx_ng + 1)
      ! Destination row of message row r: ng+nyl+r for every stagger
      ! (`fold_row_map`: T/u ng+nyl+d, v/corner ng+nyl+1+d with d = r-1).
      dbase = fx_ng + fx_nyl
      g0 = fx_r0(fam)
      ng_e = fx_nr(fam)
      cap = fx_cap
      oT = fx_uoff(FOLD_FAM_T)
      oU = fx_uoff(FOLD_FAM_U)
      sgn = 1.0_wp
      if (negate) sgn = -1.0_wp

      if (ng_e > 0) then
         if (on_device) then
            !$acc parallel loop collapse(3) private(p, ib, cls, val) &
            !$acc& present(fx_rbuf, fx_rcol, fx_rpeer, fx_re, fx_rcls, fx_rn, fld)
            do L = 1, nz
               do r = 1, nrow
                  do g = g0 + 1, g0 + ng_e
                     p = fx_rpeer(g)
                     ib = (p - 1)*cap + oT*fx_rn(p, FOLD_FAM_T) + oU*fx_rn(p, FOLD_FAM_U) &
                          + ((L - 1)*nrow + (r - 1))*fx_rn(p, fam) + fx_re(g)
                     val = sgn*fx_rbuf(ib)
                     if (has_row .and. r == 1) then
                        cls = fx_rcls(g)
                        if (cls == FOLD_ROW_WEST) then
                           fld(fx_rcol(g), dbase + 1, L) = val
                        else if (cls == FOLD_ROW_SELF .and. negate) then
                           fld(fx_rcol(g), dbase + 1, L) = 0.0_wp
                        end if
                     else
                        fld(fx_rcol(g), dbase + r, L) = val
                     end if
                  end do
               end do
            end do
         else
            do L = 1, nz
               do r = 1, nrow
                  do g = g0 + 1, g0 + ng_e
                     p = fx_rpeer(g)
                     ib = (p - 1)*cap + oT*fx_rn(p, FOLD_FAM_T) + oU*fx_rn(p, FOLD_FAM_U) &
                          + ((L - 1)*nrow + (r - 1))*fx_rn(p, fam) + fx_re(g)
                     val = sgn*fx_rbuf(ib)
                     if (has_row .and. r == 1) then
                        cls = fx_rcls(g)
                        if (cls == FOLD_ROW_WEST) then
                           fld(fx_rcol(g), dbase + 1, L) = val
                        else if (cls == FOLD_ROW_SELF .and. negate) then
                           fld(fx_rcol(g), dbase + 1, L) = 0.0_wp
                        end if
                     else
                        fld(fx_rcol(g), dbase + r, L) = val
                     end if
                  end do
               end do
            end do
         end if
      end if
      fx_uoff(fam) = fx_uoff(fam) + nrow*nz
   end subroutine ocean_fold_unpack_3d