ApplyExtrudedTransferPlan_DGModel3D_t Subroutine

  subroutine ApplyExtrudedTransferPlan_DGModel3D_t(this,plan,nz,interp,cFirst,cLast)
    !! Transfer the staged (or migrated) pre-regrid solution of a model on an EXTRUDED mesh onto
    !! its regridded extruded mesh, driven by the 2-D transfer plan of the surface forest. Fills
    !! the new solution's local columns cFirst..cLast (global surface element indices), which
    !! must be this rank's whole new column range: (cLast-cFirst+1)*nz = solution%nElem.
    !!
    !! An extruded mesh orders its elements column by column, e3 = (e2-1)*nz + k, and adaptation
    !! acts on the surface only, so every layer of a column undergoes the same surface operation:
    !! copy, exact prolongation down the quadtree path, or conservative L2 restriction of a
    !! four-child family, tensored with the identity in zeta. This applies the 2-D plan through
    !! the tested 2-D windowed apply, ApplyTransferPlanWindow of SELF_TransferPlan_2D, on the
    !! field repacked as u(i,j,column, (zeta node, layer, variable)) - (N+1)*nz*nvar "variables"
    !! per surface element - and unpacks the result. Prolongation is therefore exact for data of
    !! degree N in (xi,eta), and restriction conserves the integral over each (xi,eta) cell at
    !! every zeta node; the Jacobian-weighted volume integral is conserved where the children's
    !! columns are nested in the parent's (flat bottom, or a bottom that is linear across the
    !! refined element). Terrain-following levels over a general bathymetry change the geometry
    !! itself on refinement, since the new nodes sample the bathymetry, and then the transfer is
    !! accurate but not exactly conservative.
    !!
    !! The old field is read from the window MigrateOldWindow_DGModel3D_t left in winStage (its
    !! element range must cover whole columns) or, on one rank, from the transferStage of
    !! StageSolutionForTransfer_DGModel3D_t. These are the portable host paths; the extruded AMR
    !! controller calls them directly, by name, on every backend. Runs once per adaptation epoch.
    implicit none
    class(DGModel3D_t),intent(inout),target :: this
    type(TransferPlan2D),intent(in) :: plan
    integer,intent(in) :: nz
    type(Lagrange),intent(in) :: interp
    integer,intent(in) :: cFirst
    integer,intent(in) :: cLast
    ! Local
    integer :: Np,nv,nvc,perElem,o1,oN,c1,c2,c,cl,k,kz,iv
    real(prec),pointer :: src(:,:,:,:,:)
    real(prec),allocatable :: uOld2(:,:,:,:),uNew2(:,:,:,:)

    Np = interp%N+1
    nv = this%nvar
    nvc = Np*nz*nv
    perElem = Np*Np*Np

    if((cLast-cFirst+1)*nz /= this%solution%nElem) then
      print*,__FILE__,':',__LINE__, &
        ' : Error : ApplyExtrudedTransferPlan needs this rank''s whole new column range.'
      stop 1
    endif

    if(this%winStageN > 0) then
      if(allocated(this%transferStage)) then
        print*,__FILE__,':',__LINE__, &
          ' : Error : ApplyExtrudedTransferPlan has both a staged local field and a migrated window.'
        stop 1
      endif
      o1 = this%winStageFirst
      oN = this%winStageN
      src(1:Np,1:Np,1:Np,o1:o1+oN-1,1:nv) => this%winStage(1:perElem*oN*nv)
    elseif(allocated(this%transferStage)) then
      o1 = 1
      oN = size(this%transferStage,4)
      src => this%transferStage
    else
      print*,__FILE__,':',__LINE__, &
        ' : Error : ApplyExtrudedTransferPlan called without a staged solution or a migrated window.'
      stop 1
    endif
    if(mod(o1-1,nz) /= 0 .or. mod(oN,nz) /= 0) then
      print*,__FILE__,':',__LINE__, &
        ' : Error : ApplyExtrudedTransferPlan given an old-element range that splits a column.'
      stop 1
    endif
    c1 = (o1-1)/nz+1
    c2 = c1+oN/nz-1

    allocate(uOld2(1:Np,1:Np,c1:c2,1:nvc),uNew2(1:Np,1:Np,1:cLast-cFirst+1,1:nvc))
    do iv = 1,nv
      do k = 1,nz
        do kz = 1,Np
          do c = c1,c2
            uOld2(1:Np,1:Np,c,kz+Np*((k-1)+nz*(iv-1))) = src(1:Np,1:Np,kz,(c-1)*nz+k,iv)
          enddo
        enddo
      enddo
    enddo
    src => null()

    call ApplyTransferPlanWindow2D(plan,interp,nvc,uOld2,c1,c2,cFirst,cLast,uNew2)

    do iv = 1,nv
      do k = 1,nz
        do kz = 1,Np
          do cl = 1,cLast-cFirst+1
            this%solution%interior(1:Np,1:Np,kz,(cl-1)*nz+k,iv) = uNew2(1:Np,1:Np,cl,kz+Np*((k-1)+nz*(iv-1)))
          enddo
        enddo
      enddo
    enddo
    deallocate(uOld2,uNew2)

    call this%solution%UpdateDevice()

    if(allocated(this%transferStage)) deallocate(this%transferStage)
    this%winStageN = 0

  endsubroutine ApplyExtrudedTransferPlan_DGModel3D_t