solverUtils_b.f90

Go to the documentation of this file.

!        generated by tapenade     (inria, ecuador team)
!  tapenade 3.16 (develop) - 22 aug 2023 15:51
!
module solverutils_b
  implicit none
 
contains
!  differentiation of timestep_block in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: rhoinf pinfcorr *rev *dtl *p
!                *sfacei *sfacej *sfacek *w *rlv *vol *si *sj *sk
!                *radi *radj *radk
!   with respect to varying inputs: rhoinf pinfcorr *rev *dtl *p
!                *sfacei *sfacej *sfacek *w *rlv *vol *si *sj *sk
!                *radi *radj *radk
!   rw status of diff variables: rhoinf:incr pinfcorr:incr *rev:incr
!                *dtl:in-out *p:incr *sfacei:incr *sfacej:incr
!                *sfacek:incr *w:incr *rlv:incr *vol:incr *si:incr
!                *sj:incr *sk:incr *radi:in-out *radj:in-out *radk:in-out
!   plus diff mem management of: rev:in dtl:in p:in sfacei:in sfacej:in
!                sfacek:in w:in rlv:in vol:in si:in sj:in sk:in
!                radi:in radj:in radk:in
  subroutine timestep_block_b(onlyradii)
!
!       timestep computes the time step, or more precisely the time
!       step divided by the volume per unit cfl, in the owned cells.
!       however, for the artificial dissipation schemes, the spectral
!       radii in the halo's are needed. therefore the loop is taken
!       over the the first level of halo cells. the spectral radii are
!       stored and possibly modified for high aspect ratio cells.
!
    use constants
    use blockpointers, only : ie, je, ke, il, jl, kl, w, wd, p, pd, &
&   rlv, rlvd, rev, revd, radi, radid, radj, radjd, radk, radkd, si, sid&
&   , sj, sjd, sk, skd, sfacei, sfaceid, sfacej, sfacejd, sfacek, &
&   sfacekd, dtl, dtld, gamma, vol, vold, addgridvelocities, sectionid
    use flowvarrefstate, only : timeref, timerefd, eddymodel, gammainf&
&   , pinfcorr, pinfcorrd, viscous, rhoinf, rhoinfd
    use inputdiscretization, only : adis, dirscaling, radiineededcoarse,&
&   radiineededfine, precond, acousticscalefactor
    use inputphysics, only : equationmode
    use iteration, only : groundlevel, currentlevel
    use section, only : sections
    use inputtimespectral, only : ntimeintervalsspectral
    use utils_b, only : terminate
    implicit none
!
!      subroutine argument.
!
    logical, intent(in) :: onlyradii
!
!      local parameters.
!
    real(kind=realtype), parameter :: b=2.0_realtype
!
!      local variables.
!
    integer(kind=inttype) :: i, j, k, ii
    real(kind=realtype) :: plim, rlim, clim2
    real(kind=realtype) :: plimd, clim2d
    real(kind=realtype) :: uux, uuy, uuz, cc2, qsi, qsj, qsk, sx, sy, sz&
&   , rmu
    real(kind=realtype) :: uuxd, uuyd, uuzd, cc2d, qsid, qsjd, qskd, sxd&
&   , syd, szd, rmud
    real(kind=realtype) :: ri, rj, rk, rij, rjk, rki
    real(kind=realtype) :: rid, rjd, rkd, rijd, rjkd, rkid
    real(kind=realtype) :: vsi, vsj, vsk, rfl, dpi, dpj, dpk
    real(kind=realtype) :: vsid, vsjd, vskd, rfld, dpid, dpjd, dpkd
    real(kind=realtype) :: sface, tmp
    real(kind=realtype) :: sfaced
    logical :: radiineeded, doscaling
    intrinsic mod
    intrinsic max
    intrinsic abs
    intrinsic sqrt
    real(kind=realtype) :: abs0
    real(kind=realtype) :: abs0d
    real(kind=realtype) :: abs1
    real(kind=realtype) :: abs1d
    real(kind=realtype) :: abs2
    real(kind=realtype) :: abs2d
    real(kind=realtype) :: abs3
    real(kind=realtype) :: abs3d
    real(kind=realtype) :: abs4
    real(kind=realtype) :: abs4d
    real(kind=realtype) :: abs5
    real(kind=realtype) :: abs5d
    real(kind=realtype) :: temp
    real(kind=realtype) :: tempd
    real(kind=realtype) :: temp0
    real(kind=realtype) :: temp1
    real(kind=realtype) :: tempd0
    real(kind=realtype) :: tempd1
    integer :: branch
! determine whether or not the spectral radii are needed for the
! flux computation.
    radiineeded = radiineededcoarse
    if (currentlevel .le. groundlevel) radiineeded = radiineededfine
! return immediately if only the spectral radii must be computed
! and these are not needed for the flux computation.
    if (.not.(onlyradii .and. (.not.radiineeded))) then
! set the value of plim. to be fully consistent this must have
! the dimension of a pressure. therefore a fraction of pinfcorr
! is used. idem for rlim; compute clim2 as well.
      plim = 0.001_realtype*pinfcorr
      clim2 = 0.000001_realtype*gammainf*pinfcorr/rhoinf
      doscaling = dirscaling .and. currentlevel .le. groundlevel
! initialize sface to zero. this value will be used if the
! block is not moving.
      sface = zero
!
!           inviscid contribution, depending on the preconditioner.
!           compute the cell centered values of the spectral radii.
!
      select case  (precond) 
      case (noprecond) 
        call pushreal8(sface)
!$fwd-of ii-loop 
! no preconditioner. simply the standard spectral radius.
! loop over the cells, including the first level halo.
        do ii=0,ie*je*ke-1
          i = mod(ii, ie) + 1
          j = mod(ii/ie, je) + 1
          k = ii/(ie*je) + 1
! compute the velocities and speed of sound squared.
          uux = w(i, j, k, ivx)
          uuy = w(i, j, k, ivy)
          uuz = w(i, j, k, ivz)
          cc2 = gamma(i, j, k)*p(i, j, k)/w(i, j, k, irho)
          if (cc2 .lt. clim2) then
            cc2 = clim2
          else
            cc2 = cc2
          end if
! set the dot product of the grid velocity and the
! normal in i-direction for a moving face. to avoid
! a number of multiplications by 0.5 simply the sum
! is taken.
          if (addgridvelocities) sface = sfacei(i-1, j, k) + sfacei(i, j&
&             , k)
! spectral radius in i-direction.
          sx = si(i-1, j, k, 1) + si(i, j, k, 1)
          sy = si(i-1, j, k, 2) + si(i, j, k, 2)
          sz = si(i-1, j, k, 3) + si(i, j, k, 3)
          qsi = uux*sx + uuy*sy + uuz*sz - sface
          if (qsi .ge. 0.) then
            abs0 = qsi
          else
            abs0 = -qsi
          end if
          ri = half*(abs0+acousticscalefactor*sqrt(cc2*(sx**2+sy**2+sz**&
&           2)))
! the grid velocity in j-direction.
          if (addgridvelocities) sface = sfacej(i, j-1, k) + sfacej(i, j&
&             , k)
! spectral radius in j-direction.
          sx = sj(i, j-1, k, 1) + sj(i, j, k, 1)
          sy = sj(i, j-1, k, 2) + sj(i, j, k, 2)
          sz = sj(i, j-1, k, 3) + sj(i, j, k, 3)
          qsj = uux*sx + uuy*sy + uuz*sz - sface
          if (qsj .ge. 0.) then
            abs1 = qsj
          else
            abs1 = -qsj
          end if
          rj = half*(abs1+acousticscalefactor*sqrt(cc2*(sx**2+sy**2+sz**&
&           2)))
! the grid velocity in k-direction.
          if (addgridvelocities) sface = sfacek(i, j, k-1) + sfacek(i, j&
&             , k)
! spectral radius in k-direction.
          sx = sk(i, j, k-1, 1) + sk(i, j, k, 1)
          sy = sk(i, j, k-1, 2) + sk(i, j, k, 2)
          sz = sk(i, j, k-1, 3) + sk(i, j, k, 3)
          qsk = uux*sx + uuy*sy + uuz*sz - sface
          if (qsk .ge. 0.) then
            abs2 = qsk
          else
            abs2 = -qsk
          end if
          rk = half*(abs2+acousticscalefactor*sqrt(cc2*(sx**2+sy**2+sz**&
&           2)))
! compute the inviscid contribution to the time step.
          if (.not.onlyradii) dtl(i, j, k) = ri + rj + rk
!
!           adapt the spectral radii if directional scaling must be
!           applied.
!
          if (doscaling) then
            if (ri .lt. eps) then
              ri = eps
            else
              ri = ri
            end if
            if (rj .lt. eps) then
              rj = eps
            else
              rj = rj
            end if
            if (rk .lt. eps) then
              rk = eps
            else
              rk = rk
            end if
! compute the scaling in the three coordinate
! directions.
            rij = (ri/rj)**adis
            rjk = (rj/rk)**adis
            rki = (rk/ri)**adis
! create the scaled versions of the aspect ratios.
! note that the multiplication is done with radi, radj
! and radk, such that the influence of the clipping
! is negligible.
            radi(i, j, k) = ri*(one+one/rij+rki)
            radj(i, j, k) = rj*(one+one/rjk+rij)
            radk(i, j, k) = rk*(one+one/rki+rjk)
          else
            radi(i, j, k) = ri
            radj(i, j, k) = rj
            radk(i, j, k) = rk
          end if
        end do
        call pushcontrol2b(1)
      case (turkel) 
        call pushcontrol2b(2)
      case (choimerkle) 
        call pushcontrol2b(3)
      case default
        call pushcontrol2b(0)
      end select
! the rest of this file can be skipped if only the spectral
! radii need to be computed.
      if (.not.onlyradii) then
! the viscous contribution, if needed.
        if (viscous) then
! loop over the owned cell centers.
          do k=2,kl
            do j=2,jl
              do i=2,il
! compute the effective viscosity coefficient. the
! factor 0.5 is a combination of two things. in the
! standard central discretization of a second
! derivative there is a factor 2 multiplying the
! central node. however in the code below not the
! average but the sum of the left and the right face
! is taken and squared. this leads to a factor 4.
! combining both effects leads to 0.5. furthermore,
! it is divided by the volume and density to obtain
! the correct dimensions and multiplied by the
! non-dimensional factor factvis.
                call pushreal8(rmu)
                rmu = rlv(i, j, k)
                if (eddymodel) then
                  rmu = rmu + rev(i, j, k)
                  call pushcontrol1b(0)
                else
                  call pushcontrol1b(1)
                end if
                call pushreal8(rmu)
                rmu = half*rmu/(w(i, j, k, irho)*vol(i, j, k))
! add the viscous contribution in i-direction to the
! (inverse) of the time step.
                call pushreal8(sx)
                sx = si(i, j, k, 1) + si(i-1, j, k, 1)
                call pushreal8(sy)
                sy = si(i, j, k, 2) + si(i-1, j, k, 2)
                call pushreal8(sz)
                sz = si(i, j, k, 3) + si(i-1, j, k, 3)
                vsi = rmu*(sx*sx+sy*sy+sz*sz)
                dtl(i, j, k) = dtl(i, j, k) + vsi
! add the viscous contribution in j-direction to the
! (inverse) of the time step.
                sx = sj(i, j, k, 1) + sj(i, j-1, k, 1)
                sy = sj(i, j, k, 2) + sj(i, j-1, k, 2)
                sz = sj(i, j, k, 3) + sj(i, j-1, k, 3)
                vsj = rmu*(sx*sx+sy*sy+sz*sz)
                dtl(i, j, k) = dtl(i, j, k) + vsj
! add the viscous contribution in k-direction to the
! (inverse) of the time step.
                sx = sk(i, j, k, 1) + sk(i, j, k-1, 1)
                sy = sk(i, j, k, 2) + sk(i, j, k-1, 2)
                sz = sk(i, j, k, 3) + sk(i, j, k-1, 3)
                vsk = rmu*(sx*sx+sy*sy+sz*sz)
                dtl(i, j, k) = dtl(i, j, k) + vsk
              end do
            end do
          end do
          call pushcontrol1b(0)
        else
          call pushcontrol1b(1)
        end if
! for the spectral mode an additional term term must be
! taken into account, which corresponds to the contribution
! of the highest frequency.
        if (equationmode .eq. timespectral) then
          tmp = ntimeintervalsspectral*pi*timeref/sections(sectionid)%&
&           timeperiod
! loop over the owned cell centers and add the term.
          do k=2,kl
            do j=2,jl
              do i=2,il
                dtl(i, j, k) = dtl(i, j, k) + tmp*vol(i, j, k)
              end do
            end do
          end do
          call pushcontrol1b(1)
        else
          call pushcontrol1b(0)
        end if
! currently the inverse of dt/vol is stored in dtl. invert
! this value such that the time step per unit cfl number is
! stored and correct in cases of high gradients.
        do k=2,kl
          do j=2,jl
            do i=2,il
              if (p(i+1, j, k) - two*p(i, j, k) + p(i-1, j, k) .ge. 0.) &
&             then
                call pushreal8(abs3)
                abs3 = p(i+1, j, k) - two*p(i, j, k) + p(i-1, j, k)
                call pushcontrol1b(0)
              else
                call pushreal8(abs3)
                abs3 = -(p(i+1, j, k)-two*p(i, j, k)+p(i-1, j, k))
                call pushcontrol1b(1)
              end if
              call pushreal8(dpi)
              dpi = abs3/(p(i+1, j, k)+two*p(i, j, k)+p(i-1, j, k)+plim)
              if (p(i, j+1, k) - two*p(i, j, k) + p(i, j-1, k) .ge. 0.) &
&             then
                call pushreal8(abs4)
                abs4 = p(i, j+1, k) - two*p(i, j, k) + p(i, j-1, k)
                call pushcontrol1b(0)
              else
                call pushreal8(abs4)
                abs4 = -(p(i, j+1, k)-two*p(i, j, k)+p(i, j-1, k))
                call pushcontrol1b(1)
              end if
              call pushreal8(dpj)
              dpj = abs4/(p(i, j+1, k)+two*p(i, j, k)+p(i, j-1, k)+plim)
              if (p(i, j, k+1) - two*p(i, j, k) + p(i, j, k-1) .ge. 0.) &
&             then
                call pushreal8(abs5)
                abs5 = p(i, j, k+1) - two*p(i, j, k) + p(i, j, k-1)
                call pushcontrol1b(0)
              else
                call pushreal8(abs5)
                abs5 = -(p(i, j, k+1)-two*p(i, j, k)+p(i, j, k-1))
                call pushcontrol1b(1)
              end if
              call pushreal8(dpk)
              dpk = abs5/(p(i, j, k+1)+two*p(i, j, k)+p(i, j, k-1)+plim)
              rfl = one/(one+b*(dpi+dpj+dpk))
              call pushreal8(dtl(i, j, k))
              dtl(i, j, k) = rfl/dtl(i, j, k)
            end do
          end do
        end do
        plimd = 0.0_8
        do k=kl,2,-1
          do j=jl,2,-1
            do i=il,2,-1
              rfl = one/(one+b*(dpi+dpj+dpk))
              call popreal8(dtl(i, j, k))
              tempd0 = dtld(i, j, k)/dtl(i, j, k)
              dtld(i, j, k) = -(rfl*tempd0/dtl(i, j, k))
              rfld = tempd0
              temp1 = one + b*(dpi+dpj+dpk)
              tempd1 = -(b*one*rfld/temp1**2)
              dpid = tempd1
              dpjd = tempd1
              dpkd = tempd1
              call popreal8(dpk)
              temp1 = p(i, j, k+1) + two*p(i, j, k) + p(i, j, k-1) + &
&               plim
              abs5d = dpkd/temp1
              tempd0 = -(abs5*dpkd/temp1**2)
              pd(i, j, k+1) = pd(i, j, k+1) + tempd0
              pd(i, j, k) = pd(i, j, k) + two*tempd0
              pd(i, j, k-1) = pd(i, j, k-1) + tempd0
              plimd = plimd + tempd0
              call popcontrol1b(branch)
              if (branch .eq. 0) then
                call popreal8(abs5)
                pd(i, j, k+1) = pd(i, j, k+1) + abs5d
                pd(i, j, k) = pd(i, j, k) - two*abs5d
                pd(i, j, k-1) = pd(i, j, k-1) + abs5d
              else
                call popreal8(abs5)
                pd(i, j, k) = pd(i, j, k) + two*abs5d
                pd(i, j, k+1) = pd(i, j, k+1) - abs5d
                pd(i, j, k-1) = pd(i, j, k-1) - abs5d
              end if
              call popreal8(dpj)
              temp1 = p(i, j+1, k) + two*p(i, j, k) + p(i, j-1, k) + &
&               plim
              abs4d = dpjd/temp1
              tempd0 = -(abs4*dpjd/temp1**2)
              pd(i, j+1, k) = pd(i, j+1, k) + tempd0
              pd(i, j, k) = pd(i, j, k) + two*tempd0
              pd(i, j-1, k) = pd(i, j-1, k) + tempd0
              plimd = plimd + tempd0
              call popcontrol1b(branch)
              if (branch .eq. 0) then
                call popreal8(abs4)
                pd(i, j+1, k) = pd(i, j+1, k) + abs4d
                pd(i, j, k) = pd(i, j, k) - two*abs4d
                pd(i, j-1, k) = pd(i, j-1, k) + abs4d
              else
                call popreal8(abs4)
                pd(i, j, k) = pd(i, j, k) + two*abs4d
                pd(i, j+1, k) = pd(i, j+1, k) - abs4d
                pd(i, j-1, k) = pd(i, j-1, k) - abs4d
              end if
              call popreal8(dpi)
              temp1 = p(i+1, j, k) + two*p(i, j, k) + p(i-1, j, k) + &
&               plim
              abs3d = dpid/temp1
              tempd0 = -(abs3*dpid/temp1**2)
              pd(i+1, j, k) = pd(i+1, j, k) + tempd0
              pd(i, j, k) = pd(i, j, k) + two*tempd0
              pd(i-1, j, k) = pd(i-1, j, k) + tempd0
              plimd = plimd + tempd0
              call popcontrol1b(branch)
              if (branch .eq. 0) then
                call popreal8(abs3)
                pd(i+1, j, k) = pd(i+1, j, k) + abs3d
                pd(i, j, k) = pd(i, j, k) - two*abs3d
                pd(i-1, j, k) = pd(i-1, j, k) + abs3d
              else
                call popreal8(abs3)
                pd(i, j, k) = pd(i, j, k) + two*abs3d
                pd(i+1, j, k) = pd(i+1, j, k) - abs3d
                pd(i-1, j, k) = pd(i-1, j, k) - abs3d
              end if
            end do
          end do
        end do
        call popcontrol1b(branch)
        if (branch .ne. 0) then
          do k=kl,2,-1
            do j=jl,2,-1
              do i=il,2,-1
                vold(i, j, k) = vold(i, j, k) + tmp*dtld(i, j, k)
              end do
            end do
          end do
        end if
        call popcontrol1b(branch)
        if (branch .eq. 0) then
          do k=kl,2,-1
            do j=jl,2,-1
              do i=il,2,-1
                vskd = dtld(i, j, k)
                rmud = (sx**2+sy**2+sz**2)*vskd
                tempd0 = rmu*vskd
                sxd = 2*sx*tempd0
                syd = 2*sy*tempd0
                szd = 2*sz*tempd0
                skd(i, j, k, 3) = skd(i, j, k, 3) + szd
                skd(i, j, k-1, 3) = skd(i, j, k-1, 3) + szd
                skd(i, j, k, 2) = skd(i, j, k, 2) + syd
                skd(i, j, k-1, 2) = skd(i, j, k-1, 2) + syd
                skd(i, j, k, 1) = skd(i, j, k, 1) + sxd
                skd(i, j, k-1, 1) = skd(i, j, k-1, 1) + sxd
                sy = sj(i, j, k, 2) + sj(i, j-1, k, 2)
                sz = sj(i, j, k, 3) + sj(i, j-1, k, 3)
                sx = sj(i, j, k, 1) + sj(i, j-1, k, 1)
                vsjd = dtld(i, j, k)
                rmud = rmud + (sx**2+sy**2+sz**2)*vsjd
                tempd0 = rmu*vsjd
                sxd = 2*sx*tempd0
                syd = 2*sy*tempd0
                szd = 2*sz*tempd0
                sjd(i, j, k, 3) = sjd(i, j, k, 3) + szd
                sjd(i, j-1, k, 3) = sjd(i, j-1, k, 3) + szd
                sjd(i, j, k, 2) = sjd(i, j, k, 2) + syd
                sjd(i, j-1, k, 2) = sjd(i, j-1, k, 2) + syd
                sjd(i, j, k, 1) = sjd(i, j, k, 1) + sxd
                sjd(i, j-1, k, 1) = sjd(i, j-1, k, 1) + sxd
                sy = si(i, j, k, 2) + si(i-1, j, k, 2)
                sz = si(i, j, k, 3) + si(i-1, j, k, 3)
                sx = si(i, j, k, 1) + si(i-1, j, k, 1)
                vsid = dtld(i, j, k)
                rmud = rmud + (sx**2+sy**2+sz**2)*vsid
                tempd0 = rmu*vsid
                sxd = 2*sx*tempd0
                syd = 2*sy*tempd0
                szd = 2*sz*tempd0
                call popreal8(sz)
                sid(i, j, k, 3) = sid(i, j, k, 3) + szd
                sid(i-1, j, k, 3) = sid(i-1, j, k, 3) + szd
                call popreal8(sy)
                sid(i, j, k, 2) = sid(i, j, k, 2) + syd
                sid(i-1, j, k, 2) = sid(i-1, j, k, 2) + syd
                call popreal8(sx)
                sid(i, j, k, 1) = sid(i, j, k, 1) + sxd
                sid(i-1, j, k, 1) = sid(i-1, j, k, 1) + sxd
                call popreal8(rmu)
                temp = w(i, j, k, irho)
                temp0 = temp*vol(i, j, k)
                tempd0 = half*rmud/temp0
                rmud = tempd0
                tempd1 = -(rmu*tempd0/temp0)
                wd(i, j, k, irho) = wd(i, j, k, irho) + vol(i, j, k)*&
&                 tempd1
                vold(i, j, k) = vold(i, j, k) + temp*tempd1
                call popcontrol1b(branch)
                if (branch .eq. 0) revd(i, j, k) = revd(i, j, k) + rmud
                call popreal8(rmu)
                rlvd(i, j, k) = rlvd(i, j, k) + rmud
              end do
            end do
          end do
        end if
      else
        plimd = 0.0_8
      end if
      call popcontrol2b(branch)
      if (branch .lt. 2) then
        if (branch .eq. 0) then
          clim2d = 0.0_8
        else
          clim2d = 0.0_8
          sfaced = 0.0_8
          call popreal8(sface)
!$bwd-of ii-loop 
          do ii=0,ie*je*ke-1
            i = mod(ii, ie) + 1
            j = mod(ii/ie, je) + 1
            k = ii/(ie*je) + 1
! compute the velocities and speed of sound squared.
            uux = w(i, j, k, ivx)
            uuy = w(i, j, k, ivy)
            uuz = w(i, j, k, ivz)
            cc2 = gamma(i, j, k)*p(i, j, k)/w(i, j, k, irho)
            if (cc2 .lt. clim2) then
              cc2 = clim2
              call pushcontrol1b(0)
            else
              call pushcontrol1b(1)
              cc2 = cc2
            end if
! set the dot product of the grid velocity and the
! normal in i-direction for a moving face. to avoid
! a number of multiplications by 0.5 simply the sum
! is taken.
            if (addgridvelocities) then
              sface = sfacei(i-1, j, k) + sfacei(i, j, k)
              call pushcontrol1b(1)
            else
              call pushcontrol1b(0)
            end if
! spectral radius in i-direction.
            sx = si(i-1, j, k, 1) + si(i, j, k, 1)
            sy = si(i-1, j, k, 2) + si(i, j, k, 2)
            sz = si(i-1, j, k, 3) + si(i, j, k, 3)
            qsi = uux*sx + uuy*sy + uuz*sz - sface
            if (qsi .ge. 0.) then
              abs0 = qsi
              call pushcontrol1b(0)
            else
              abs0 = -qsi
              call pushcontrol1b(1)
            end if
            ri = half*(abs0+acousticscalefactor*sqrt(cc2*(sx**2+sy**2+sz&
&             **2)))
! the grid velocity in j-direction.
            if (addgridvelocities) then
              sface = sfacej(i, j-1, k) + sfacej(i, j, k)
              call pushcontrol1b(1)
            else
              call pushcontrol1b(0)
            end if
! spectral radius in j-direction.
            sx = sj(i, j-1, k, 1) + sj(i, j, k, 1)
            sy = sj(i, j-1, k, 2) + sj(i, j, k, 2)
            sz = sj(i, j-1, k, 3) + sj(i, j, k, 3)
            qsj = uux*sx + uuy*sy + uuz*sz - sface
            if (qsj .ge. 0.) then
              abs1 = qsj
              call pushcontrol1b(0)
            else
              abs1 = -qsj
              call pushcontrol1b(1)
            end if
            rj = half*(abs1+acousticscalefactor*sqrt(cc2*(sx**2+sy**2+sz&
&             **2)))
! the grid velocity in k-direction.
            if (addgridvelocities) then
              sface = sfacek(i, j, k-1) + sfacek(i, j, k)
              call pushcontrol1b(1)
            else
              call pushcontrol1b(0)
            end if
! spectral radius in k-direction.
            sx = sk(i, j, k-1, 1) + sk(i, j, k, 1)
            sy = sk(i, j, k-1, 2) + sk(i, j, k, 2)
            sz = sk(i, j, k-1, 3) + sk(i, j, k, 3)
            qsk = uux*sx + uuy*sy + uuz*sz - sface
            if (qsk .ge. 0.) then
              abs2 = qsk
              call pushcontrol1b(0)
            else
              abs2 = -qsk
              call pushcontrol1b(1)
            end if
            rk = half*(abs2+acousticscalefactor*sqrt(cc2*(sx**2+sy**2+sz&
&             **2)))
! compute the inviscid contribution to the time step.
            if (.not.onlyradii) then
              call pushcontrol1b(0)
            else
              call pushcontrol1b(1)
            end if
!
!           adapt the spectral radii if directional scaling must be
!           applied.
!
            if (doscaling) then
              if (ri .lt. eps) then
                ri = eps
                call pushcontrol1b(0)
              else
                call pushcontrol1b(1)
                ri = ri
              end if
              if (rj .lt. eps) then
                rj = eps
                call pushcontrol1b(0)
              else
                call pushcontrol1b(1)
                rj = rj
              end if
              if (rk .lt. eps) then
                rk = eps
                call pushcontrol1b(0)
              else
                call pushcontrol1b(1)
                rk = rk
              end if
! compute the scaling in the three coordinate
! directions.
              rij = (ri/rj)**adis
              rjk = (rj/rk)**adis
              rki = (rk/ri)**adis
! create the scaled versions of the aspect ratios.
! note that the multiplication is done with radi, radj
! and radk, such that the influence of the clipping
! is negligible.
              rkid = ri*radid(i, j, k) - one*rk*radkd(i, j, k)/rki**2
              temp1 = rk/ri
              if (temp1 .le. 0.0_8 .and. (adis .eq. 0.0_8 .or. adis .ne.&
&                 int(adis))) then
                tempd0 = 0.0_8
              else
                tempd0 = adis*temp1**(adis-1)*rkid/ri
              end if
              rkd = (one+one/rki+rjk)*radkd(i, j, k) + tempd0
              rjkd = rk*radkd(i, j, k) - one*rj*radjd(i, j, k)/rjk**2
              radkd(i, j, k) = 0.0_8
              rijd = rj*radjd(i, j, k) - one*ri*radid(i, j, k)/rij**2
              rid = (one+one/rij+rki)*radid(i, j, k) - temp1*tempd0
              radid(i, j, k) = 0.0_8
              temp1 = rj/rk
              if (temp1 .le. 0.0_8 .and. (adis .eq. 0.0_8 .or. adis .ne.&
&                 int(adis))) then
                tempd0 = 0.0_8
              else
                tempd0 = adis*temp1**(adis-1)*rjkd/rk
              end if
              rjd = (one+one/rjk+rij)*radjd(i, j, k) + tempd0
              radjd(i, j, k) = 0.0_8
              rkd = rkd - temp1*tempd0
              temp1 = ri/rj
              if (temp1 .le. 0.0_8 .and. (adis .eq. 0.0_8 .or. adis .ne.&
&                 int(adis))) then
                tempd0 = 0.0_8
              else
                tempd0 = adis*temp1**(adis-1)*rijd/rj
              end if
              rid = rid + tempd0
              rjd = rjd - temp1*tempd0
              call popcontrol1b(branch)
              if (branch .eq. 0) rkd = 0.0_8
              call popcontrol1b(branch)
              if (branch .eq. 0) rjd = 0.0_8
              call popcontrol1b(branch)
              if (branch .eq. 0) rid = 0.0_8
            else
              rkd = radkd(i, j, k)
              radkd(i, j, k) = 0.0_8
              rjd = radjd(i, j, k)
              radjd(i, j, k) = 0.0_8
              rid = radid(i, j, k)
              radid(i, j, k) = 0.0_8
            end if
            call popcontrol1b(branch)
            if (branch .eq. 0) then
              rid = rid + dtld(i, j, k)
              rjd = rjd + dtld(i, j, k)
              rkd = rkd + dtld(i, j, k)
              dtld(i, j, k) = 0.0_8
            end if
            temp1 = sx*sx + sy*sy + sz*sz
            abs2d = half*rkd
            if (cc2*temp1 .eq. 0.0_8) then
              tempd1 = 0.0_8
            else
              tempd1 = acousticscalefactor*half*rkd/(2.0*sqrt(cc2*temp1)&
&               )
            end if
            cc2d = temp1*tempd1
            tempd0 = cc2*tempd1
            sxd = 2*sx*tempd0
            syd = 2*sy*tempd0
            szd = 2*sz*tempd0
            call popcontrol1b(branch)
            if (branch .eq. 0) then
              qskd = abs2d
            else
              qskd = -abs2d
            end if
            uuxd = sx*qskd
            sxd = sxd + uux*qskd
            uuyd = sy*qskd
            syd = syd + uuy*qskd
            uuzd = sz*qskd
            szd = szd + uuz*qskd
            sfaced = sfaced - qskd
            skd(i, j, k-1, 3) = skd(i, j, k-1, 3) + szd
            skd(i, j, k, 3) = skd(i, j, k, 3) + szd
            skd(i, j, k-1, 2) = skd(i, j, k-1, 2) + syd
            skd(i, j, k, 2) = skd(i, j, k, 2) + syd
            skd(i, j, k-1, 1) = skd(i, j, k-1, 1) + sxd
            skd(i, j, k, 1) = skd(i, j, k, 1) + sxd
            call popcontrol1b(branch)
            if (branch .ne. 0) then
              sfacekd(i, j, k-1) = sfacekd(i, j, k-1) + sfaced
              sfacekd(i, j, k) = sfacekd(i, j, k) + sfaced
              sfaced = 0.0_8
            end if
            sx = sj(i, j-1, k, 1) + sj(i, j, k, 1)
            sy = sj(i, j-1, k, 2) + sj(i, j, k, 2)
            sz = sj(i, j-1, k, 3) + sj(i, j, k, 3)
            temp1 = sx*sx + sy*sy + sz*sz
            abs1d = half*rjd
            if (cc2*temp1 .eq. 0.0_8) then
              tempd1 = 0.0_8
            else
              tempd1 = acousticscalefactor*half*rjd/(2.0*sqrt(cc2*temp1)&
&               )
            end if
            cc2d = cc2d + temp1*tempd1
            tempd0 = cc2*tempd1
            sxd = 2*sx*tempd0
            syd = 2*sy*tempd0
            szd = 2*sz*tempd0
            call popcontrol1b(branch)
            if (branch .eq. 0) then
              qsjd = abs1d
            else
              qsjd = -abs1d
            end if
            uuxd = uuxd + sx*qsjd
            sxd = sxd + uux*qsjd
            uuyd = uuyd + sy*qsjd
            syd = syd + uuy*qsjd
            uuzd = uuzd + sz*qsjd
            szd = szd + uuz*qsjd
            sfaced = sfaced - qsjd
            sjd(i, j-1, k, 3) = sjd(i, j-1, k, 3) + szd
            sjd(i, j, k, 3) = sjd(i, j, k, 3) + szd
            sjd(i, j-1, k, 2) = sjd(i, j-1, k, 2) + syd
            sjd(i, j, k, 2) = sjd(i, j, k, 2) + syd
            sjd(i, j-1, k, 1) = sjd(i, j-1, k, 1) + sxd
            sjd(i, j, k, 1) = sjd(i, j, k, 1) + sxd
            call popcontrol1b(branch)
            if (branch .ne. 0) then
              sfacejd(i, j-1, k) = sfacejd(i, j-1, k) + sfaced
              sfacejd(i, j, k) = sfacejd(i, j, k) + sfaced
              sfaced = 0.0_8
            end if
            sx = si(i-1, j, k, 1) + si(i, j, k, 1)
            sy = si(i-1, j, k, 2) + si(i, j, k, 2)
            sz = si(i-1, j, k, 3) + si(i, j, k, 3)
            temp1 = sx*sx + sy*sy + sz*sz
            abs0d = half*rid
            if (cc2*temp1 .eq. 0.0_8) then
              tempd1 = 0.0_8
            else
              tempd1 = acousticscalefactor*half*rid/(2.0*sqrt(cc2*temp1)&
&               )
            end if
            cc2d = cc2d + temp1*tempd1
            tempd0 = cc2*tempd1
            sxd = 2*sx*tempd0
            syd = 2*sy*tempd0
            szd = 2*sz*tempd0
            call popcontrol1b(branch)
            if (branch .eq. 0) then
              qsid = abs0d
            else
              qsid = -abs0d
            end if
            uuxd = uuxd + sx*qsid
            sxd = sxd + uux*qsid
            uuyd = uuyd + sy*qsid
            syd = syd + uuy*qsid
            uuzd = uuzd + sz*qsid
            szd = szd + uuz*qsid
            sfaced = sfaced - qsid
            sid(i-1, j, k, 3) = sid(i-1, j, k, 3) + szd
            sid(i, j, k, 3) = sid(i, j, k, 3) + szd
            sid(i-1, j, k, 2) = sid(i-1, j, k, 2) + syd
            sid(i, j, k, 2) = sid(i, j, k, 2) + syd
            sid(i-1, j, k, 1) = sid(i-1, j, k, 1) + sxd
            sid(i, j, k, 1) = sid(i, j, k, 1) + sxd
            call popcontrol1b(branch)
            if (branch .ne. 0) then
              sfaceid(i-1, j, k) = sfaceid(i-1, j, k) + sfaced
              sfaceid(i, j, k) = sfaceid(i, j, k) + sfaced
              sfaced = 0.0_8
            end if
            call popcontrol1b(branch)
            if (branch .eq. 0) then
              clim2d = clim2d + cc2d
              cc2d = 0.0_8
            end if
            temp1 = w(i, j, k, irho)
            tempd1 = gamma(i, j, k)*cc2d/temp1
            pd(i, j, k) = pd(i, j, k) + tempd1
            wd(i, j, k, irho) = wd(i, j, k, irho) - p(i, j, k)*tempd1/&
&             temp1
            wd(i, j, k, ivz) = wd(i, j, k, ivz) + uuzd
            wd(i, j, k, ivy) = wd(i, j, k, ivy) + uuyd
            wd(i, j, k, ivx) = wd(i, j, k, ivx) + uuxd
          end do
        end if
      else if (branch .eq. 2) then
        clim2d = 0.0_8
      else
        clim2d = 0.0_8
      end if
      tempd = gammainf*0.000001_realtype*clim2d/rhoinf
      pinfcorrd = pinfcorrd + tempd + 0.001_realtype*plimd
      rhoinfd = rhoinfd - pinfcorr*tempd/rhoinf
    end if
  end subroutine timestep_block_b
 
  subroutine timestep_block(onlyradii)
!
!       timestep computes the time step, or more precisely the time
!       step divided by the volume per unit cfl, in the owned cells.
!       however, for the artificial dissipation schemes, the spectral
!       radii in the halo's are needed. therefore the loop is taken
!       over the the first level of halo cells. the spectral radii are
!       stored and possibly modified for high aspect ratio cells.
!
    use constants
    use blockpointers, only : ie, je, ke, il, jl, kl, w, p, rlv, rev, &
&   radi, radj, radk, si, sj, sk, sfacei, sfacej, sfacek, dtl, gamma, &
&   vol, addgridvelocities, sectionid
    use flowvarrefstate, only : timeref, eddymodel, gammainf, pinfcorr&
&   , viscous, rhoinf
    use inputdiscretization, only : adis, dirscaling, radiineededcoarse,&
&   radiineededfine, precond, acousticscalefactor
    use inputphysics, only : equationmode
    use iteration, only : groundlevel, currentlevel
    use section, only : sections
    use inputtimespectral, only : ntimeintervalsspectral
    use utils_b, only : terminate
    implicit none
!
!      subroutine argument.
!
    logical, intent(in) :: onlyradii
!
!      local parameters.
!
    real(kind=realtype), parameter :: b=2.0_realtype
!
!      local variables.
!
    integer(kind=inttype) :: i, j, k, ii
    real(kind=realtype) :: plim, rlim, clim2
    real(kind=realtype) :: uux, uuy, uuz, cc2, qsi, qsj, qsk, sx, sy, sz&
&   , rmu
    real(kind=realtype) :: ri, rj, rk, rij, rjk, rki
    real(kind=realtype) :: vsi, vsj, vsk, rfl, dpi, dpj, dpk
    real(kind=realtype) :: sface, tmp
    logical :: radiineeded, doscaling
    intrinsic mod
    intrinsic max
    intrinsic abs
    intrinsic sqrt
    real(kind=realtype) :: abs0
    real(kind=realtype) :: abs1
    real(kind=realtype) :: abs2
    real(kind=realtype) :: abs3
    real(kind=realtype) :: abs4
    real(kind=realtype) :: abs5
! determine whether or not the spectral radii are needed for the
! flux computation.
    radiineeded = radiineededcoarse
    if (currentlevel .le. groundlevel) radiineeded = radiineededfine
! return immediately if only the spectral radii must be computed
! and these are not needed for the flux computation.
    if (onlyradii .and. (.not.radiineeded)) then
      return
    else
! set the value of plim. to be fully consistent this must have
! the dimension of a pressure. therefore a fraction of pinfcorr
! is used. idem for rlim; compute clim2 as well.
      plim = 0.001_realtype*pinfcorr
      rlim = 0.001_realtype*rhoinf
      clim2 = 0.000001_realtype*gammainf*pinfcorr/rhoinf
      doscaling = dirscaling .and. currentlevel .le. groundlevel
! initialize sface to zero. this value will be used if the
! block is not moving.
      sface = zero
!
!           inviscid contribution, depending on the preconditioner.
!           compute the cell centered values of the spectral radii.
!
      select case  (precond) 
      case (noprecond) 
!$ad ii-loop
! no preconditioner. simply the standard spectral radius.
! loop over the cells, including the first level halo.
        do ii=0,ie*je*ke-1
          i = mod(ii, ie) + 1
          j = mod(ii/ie, je) + 1
          k = ii/(ie*je) + 1
! compute the velocities and speed of sound squared.
          uux = w(i, j, k, ivx)
          uuy = w(i, j, k, ivy)
          uuz = w(i, j, k, ivz)
          cc2 = gamma(i, j, k)*p(i, j, k)/w(i, j, k, irho)
          if (cc2 .lt. clim2) then
            cc2 = clim2
          else
            cc2 = cc2
          end if
! set the dot product of the grid velocity and the
! normal in i-direction for a moving face. to avoid
! a number of multiplications by 0.5 simply the sum
! is taken.
          if (addgridvelocities) sface = sfacei(i-1, j, k) + sfacei(i, j&
&             , k)
! spectral radius in i-direction.
          sx = si(i-1, j, k, 1) + si(i, j, k, 1)
          sy = si(i-1, j, k, 2) + si(i, j, k, 2)
          sz = si(i-1, j, k, 3) + si(i, j, k, 3)
          qsi = uux*sx + uuy*sy + uuz*sz - sface
          if (qsi .ge. 0.) then
            abs0 = qsi
          else
            abs0 = -qsi
          end if
          ri = half*(abs0+acousticscalefactor*sqrt(cc2*(sx**2+sy**2+sz**&
&           2)))
! the grid velocity in j-direction.
          if (addgridvelocities) sface = sfacej(i, j-1, k) + sfacej(i, j&
&             , k)
! spectral radius in j-direction.
          sx = sj(i, j-1, k, 1) + sj(i, j, k, 1)
          sy = sj(i, j-1, k, 2) + sj(i, j, k, 2)
          sz = sj(i, j-1, k, 3) + sj(i, j, k, 3)
          qsj = uux*sx + uuy*sy + uuz*sz - sface
          if (qsj .ge. 0.) then
            abs1 = qsj
          else
            abs1 = -qsj
          end if
          rj = half*(abs1+acousticscalefactor*sqrt(cc2*(sx**2+sy**2+sz**&
&           2)))
! the grid velocity in k-direction.
          if (addgridvelocities) sface = sfacek(i, j, k-1) + sfacek(i, j&
&             , k)
! spectral radius in k-direction.
          sx = sk(i, j, k-1, 1) + sk(i, j, k, 1)
          sy = sk(i, j, k-1, 2) + sk(i, j, k, 2)
          sz = sk(i, j, k-1, 3) + sk(i, j, k, 3)
          qsk = uux*sx + uuy*sy + uuz*sz - sface
          if (qsk .ge. 0.) then
            abs2 = qsk
          else
            abs2 = -qsk
          end if
          rk = half*(abs2+acousticscalefactor*sqrt(cc2*(sx**2+sy**2+sz**&
&           2)))
! compute the inviscid contribution to the time step.
          if (.not.onlyradii) dtl(i, j, k) = ri + rj + rk
!
!           adapt the spectral radii if directional scaling must be
!           applied.
!
          if (doscaling) then
            if (ri .lt. eps) then
              ri = eps
            else
              ri = ri
            end if
            if (rj .lt. eps) then
              rj = eps
            else
              rj = rj
            end if
            if (rk .lt. eps) then
              rk = eps
            else
              rk = rk
            end if
! compute the scaling in the three coordinate
! directions.
            rij = (ri/rj)**adis
            rjk = (rj/rk)**adis
            rki = (rk/ri)**adis
! create the scaled versions of the aspect ratios.
! note that the multiplication is done with radi, radj
! and radk, such that the influence of the clipping
! is negligible.
            radi(i, j, k) = ri*(one+one/rij+rki)
            radj(i, j, k) = rj*(one+one/rjk+rij)
            radk(i, j, k) = rk*(one+one/rki+rjk)
          else
            radi(i, j, k) = ri
            radj(i, j, k) = rj
            radk(i, j, k) = rk
          end if
        end do
      case (turkel) 
        call terminate('timestep', &
&                'turkel preconditioner not implemented yet')
      case (choimerkle) 
        call terminate('timestep', &
&                'choi merkle preconditioner not implemented yet')
      end select
! the rest of this file can be skipped if only the spectral
! radii need to be computed.
      if (.not.onlyradii) then
! the viscous contribution, if needed.
        if (viscous) then
! loop over the owned cell centers.
          do k=2,kl
            do j=2,jl
              do i=2,il
! compute the effective viscosity coefficient. the
! factor 0.5 is a combination of two things. in the
! standard central discretization of a second
! derivative there is a factor 2 multiplying the
! central node. however in the code below not the
! average but the sum of the left and the right face
! is taken and squared. this leads to a factor 4.
! combining both effects leads to 0.5. furthermore,
! it is divided by the volume and density to obtain
! the correct dimensions and multiplied by the
! non-dimensional factor factvis.
                rmu = rlv(i, j, k)
                if (eddymodel) rmu = rmu + rev(i, j, k)
                rmu = half*rmu/(w(i, j, k, irho)*vol(i, j, k))
! add the viscous contribution in i-direction to the
! (inverse) of the time step.
                sx = si(i, j, k, 1) + si(i-1, j, k, 1)
                sy = si(i, j, k, 2) + si(i-1, j, k, 2)
                sz = si(i, j, k, 3) + si(i-1, j, k, 3)
                vsi = rmu*(sx*sx+sy*sy+sz*sz)
                dtl(i, j, k) = dtl(i, j, k) + vsi
! add the viscous contribution in j-direction to the
! (inverse) of the time step.
                sx = sj(i, j, k, 1) + sj(i, j-1, k, 1)
                sy = sj(i, j, k, 2) + sj(i, j-1, k, 2)
                sz = sj(i, j, k, 3) + sj(i, j-1, k, 3)
                vsj = rmu*(sx*sx+sy*sy+sz*sz)
                dtl(i, j, k) = dtl(i, j, k) + vsj
! add the viscous contribution in k-direction to the
! (inverse) of the time step.
                sx = sk(i, j, k, 1) + sk(i, j, k-1, 1)
                sy = sk(i, j, k, 2) + sk(i, j, k-1, 2)
                sz = sk(i, j, k, 3) + sk(i, j, k-1, 3)
                vsk = rmu*(sx*sx+sy*sy+sz*sz)
                dtl(i, j, k) = dtl(i, j, k) + vsk
              end do
            end do
          end do
        end if
! for the spectral mode an additional term term must be
! taken into account, which corresponds to the contribution
! of the highest frequency.
        if (equationmode .eq. timespectral) then
          tmp = ntimeintervalsspectral*pi*timeref/sections(sectionid)%&
&           timeperiod
! loop over the owned cell centers and add the term.
          do k=2,kl
            do j=2,jl
              do i=2,il
                dtl(i, j, k) = dtl(i, j, k) + tmp*vol(i, j, k)
              end do
            end do
          end do
        end if
! currently the inverse of dt/vol is stored in dtl. invert
! this value such that the time step per unit cfl number is
! stored and correct in cases of high gradients.
        do k=2,kl
          do j=2,jl
            do i=2,il
              if (p(i+1, j, k) - two*p(i, j, k) + p(i-1, j, k) .ge. 0.) &
&             then
                abs3 = p(i+1, j, k) - two*p(i, j, k) + p(i-1, j, k)
              else
                abs3 = -(p(i+1, j, k)-two*p(i, j, k)+p(i-1, j, k))
              end if
              dpi = abs3/(p(i+1, j, k)+two*p(i, j, k)+p(i-1, j, k)+plim)
              if (p(i, j+1, k) - two*p(i, j, k) + p(i, j-1, k) .ge. 0.) &
&             then
                abs4 = p(i, j+1, k) - two*p(i, j, k) + p(i, j-1, k)
              else
                abs4 = -(p(i, j+1, k)-two*p(i, j, k)+p(i, j-1, k))
              end if
              dpj = abs4/(p(i, j+1, k)+two*p(i, j, k)+p(i, j-1, k)+plim)
              if (p(i, j, k+1) - two*p(i, j, k) + p(i, j, k-1) .ge. 0.) &
&             then
                abs5 = p(i, j, k+1) - two*p(i, j, k) + p(i, j, k-1)
              else
                abs5 = -(p(i, j, k+1)-two*p(i, j, k)+p(i, j, k-1))
              end if
              dpk = abs5/(p(i, j, k+1)+two*p(i, j, k)+p(i, j, k-1)+plim)
              rfl = one/(one+b*(dpi+dpj+dpk))
              dtl(i, j, k) = rfl/dtl(i, j, k)
            end do
          end do
        end do
      end if
    end if
  end subroutine timestep_block
 
!  differentiation of gridvelocitiesfinelevel_block in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: *(flowdoms.x) veldirfreestream
!                *sfacei *sfacej *sfacek *s *si *sj *sk
!   with respect to varying inputs: pinf timeref rhoinf *(flowdoms.x)
!                veldirfreestream machgrid *sfacei *sfacej *sfacek
!                *s *si *sj *sk
!   rw status of diff variables: pinf:out timeref:out rhoinf:out
!                *(flowdoms.x):incr veldirfreestream:incr machgrid:out
!                *sfacei:in-out *sfacej:in-out *sfacek:in-out *s:in-out
!                *si:incr *sj:incr *sk:incr
!   plus diff mem management of: flowdoms.x:in sfacei:in sfacej:in
!                sfacek:in s:in si:in sj:in sk:in
  subroutine gridvelocitiesfinelevel_block_b(useoldcoor, t, sps, nn)
!
!       gridvelocitiesfinelevel computes the grid velocities for
!       the cell centers and the normal grid velocities for the faces
!       of moving blocks for the currently finest grid, i.e.
!       groundlevel. the velocities are computed at time t for
!       spectral mode sps. if useoldcoor is .true. the velocities
!       are determined using the unsteady time integrator in
!       combination with the old coordinates; otherwise the analytic
!       form is used.
!
    use blockpointers
    use cgnsgrid
    use flowvarrefstate
    use inputmotion
    use inputunsteady
    use iteration
    use inputphysics
    use inputtsstabderiv
    use monitor
    use communication
    use flowutils_b, only : derivativerotmatrixrigid, &
&   derivativerotmatrixrigid_b, getdirvector
    use utils_b, only : setcoeftimeintegrator, tsalpha, tsbeta, tsmach, &
&   terminate, rotmatrixrigidbody, getdirangle
    implicit none
!
!      subroutine arguments.
!
    integer(kind=inttype), intent(in) :: sps, nn
    logical, intent(in) :: useoldcoor
    real(kind=realtype), dimension(*), intent(in) :: t
!
!      local variables.
!
    integer(kind=inttype) :: mm
    integer(kind=inttype) :: i, j, k, ii, iie, jje, kke
    real(kind=realtype) :: oneover4dt, oneover8dt
    real(kind=realtype) :: velxgrid, velygrid, velzgrid, ainf
    real(kind=realtype) :: velxgridd, velygridd, velzgridd, ainfd
    real(kind=realtype) :: velxgrid0, velygrid0, velzgrid0
    real(kind=realtype) :: velxgrid0d, velygrid0d, velzgrid0d
    real(kind=realtype), dimension(3) :: sc, xc, xxc
    real(kind=realtype), dimension(3) :: scd, xcd, xxcd
    real(kind=realtype), dimension(3) :: rotcenter, rotrate
    real(kind=realtype), dimension(3) :: rotcenterd, rotrated
    real(kind=realtype), dimension(3) :: rotationpoint
    real(kind=realtype), dimension(3, 3) :: rotationmatrix, &
&   derivrotationmatrix
    real(kind=realtype), dimension(3, 3) :: derivrotationmatrixd
    real(kind=realtype) :: tnew, told
    real(kind=realtype), dimension(:, :), pointer :: sface
    real(kind=realtype), dimension(:, :, :), pointer :: xx, ss
    real(kind=realtype), dimension(:, :, :, :), pointer :: xxold
    real(kind=realtype) :: intervalmach, alphats, alphaincrement, betats&
&   , betaincrement
    real(kind=realtype), dimension(3) :: veldir
    real(kind=realtype), dimension(3) :: refdirection
    intrinsic sqrt
    real(kind=realtype) :: tempd
! compute the mesh velocity from the given mesh mach number.
! vel{x,y,z}grid0 is the actual velocity you want at the
! geometry.
    ainf = sqrt(gammainf*pinf/rhoinf)
    velxgrid0 = ainf*machgrid*(-veldirfreestream(1))
    velygrid0 = ainf*machgrid*(-veldirfreestream(2))
    velzgrid0 = ainf*machgrid*(-veldirfreestream(3))
! compute the derivative of the rotation matrix and the rotation
! point; needed for velocity due to the rigid body rotation of
! the entire grid. it is assumed that the rigid body motion of
! the grid is only specified if there is only 1 section present.
    call derivativerotmatrixrigid(derivrotationmatrix, rotationpoint, t(&
&                           1))
!compute the rotation matrix to update the velocities for the time
!spectral stability derivative case...
    if (blockismoving) then
! determine the situation we are having here.
      if (useoldcoor) then
        timerefd = 0.0_8
        velxgrid0d = 0.0_8
        velzgrid0d = 0.0_8
        derivrotationmatrixd = 0.0_8
        velygrid0d = 0.0_8
      else
!
!             the velocities must be determined analytically.
!
! store the rotation center and determine the
! nondimensional rotation rate of this block. as the
! reference length is 1 timeref == 1/uref and at the end
! the nondimensional velocity is computed.
        j = nbkglobal
        rotcenter = cgnsdoms(j)%rotcenter
        rotrate = timeref*cgnsdoms(j)%rotrate
        velxgrid = velxgrid0
        velygrid = velygrid0
        velzgrid = velzgrid0
!
!             grid velocities of the cell centers, including the
!             1st level halo cells.
!
! loop over the cells, including the 1st level halo's.
        do k=1,ke
          call pushinteger4(j)
          do j=1,je
            do i=1,ie
! determine the coordinates of the cell center,
! which are stored in xc.
              xc(1) = eighth*(flowdoms(nn, groundlevel, sps)%x(i-1, j-1&
&               , k-1, 1)+flowdoms(nn, groundlevel, sps)%x(i, j-1, k-1, &
&               1)+flowdoms(nn, groundlevel, sps)%x(i-1, j, k-1, 1)+&
&               flowdoms(nn, groundlevel, sps)%x(i, j, k-1, 1)+flowdoms(&
&               nn, groundlevel, sps)%x(i-1, j-1, k, 1)+flowdoms(nn, &
&               groundlevel, sps)%x(i, j-1, k, 1)+flowdoms(nn, &
&               groundlevel, sps)%x(i-1, j, k, 1)+flowdoms(nn, &
&               groundlevel, sps)%x(i, j, k, 1))
              xc(2) = eighth*(flowdoms(nn, groundlevel, sps)%x(i-1, j-1&
&               , k-1, 2)+flowdoms(nn, groundlevel, sps)%x(i, j-1, k-1, &
&               2)+flowdoms(nn, groundlevel, sps)%x(i-1, j, k-1, 2)+&
&               flowdoms(nn, groundlevel, sps)%x(i, j, k-1, 2)+flowdoms(&
&               nn, groundlevel, sps)%x(i-1, j-1, k, 2)+flowdoms(nn, &
&               groundlevel, sps)%x(i, j-1, k, 2)+flowdoms(nn, &
&               groundlevel, sps)%x(i-1, j, k, 2)+flowdoms(nn, &
&               groundlevel, sps)%x(i, j, k, 2))
              xc(3) = eighth*(flowdoms(nn, groundlevel, sps)%x(i-1, j-1&
&               , k-1, 3)+flowdoms(nn, groundlevel, sps)%x(i, j-1, k-1, &
&               3)+flowdoms(nn, groundlevel, sps)%x(i-1, j, k-1, 3)+&
&               flowdoms(nn, groundlevel, sps)%x(i, j, k-1, 3)+flowdoms(&
&               nn, groundlevel, sps)%x(i-1, j-1, k, 3)+flowdoms(nn, &
&               groundlevel, sps)%x(i, j-1, k, 3)+flowdoms(nn, &
&               groundlevel, sps)%x(i-1, j, k, 3)+flowdoms(nn, &
&               groundlevel, sps)%x(i, j, k, 3))
! determine the coordinates relative to the
! center of rotation.
              call pushreal8(xxc(1))
              xxc(1) = xc(1) - rotcenter(1)
              call pushreal8(xxc(2))
              xxc(2) = xc(2) - rotcenter(2)
              call pushreal8(xxc(3))
              xxc(3) = xc(3) - rotcenter(3)
! determine the rotation speed of the cell center,
! which is omega*r.
              sc(1) = rotrate(2)*xxc(3) - rotrate(3)*xxc(2)
              sc(2) = rotrate(3)*xxc(1) - rotrate(1)*xxc(3)
              sc(3) = rotrate(1)*xxc(2) - rotrate(2)*xxc(1)
! determine the coordinates relative to the
! rigid body rotation point.
              call pushreal8(xxc(1))
              xxc(1) = xc(1) - rotationpoint(1)
              call pushreal8(xxc(2))
              xxc(2) = xc(2) - rotationpoint(2)
              call pushreal8(xxc(3))
              xxc(3) = xc(3) - rotationpoint(3)
! determine the total velocity of the cell center.
! this is a combination of rotation speed of this
! block and the entire rigid body rotation.
            end do
          end do
        end do
!
!             normal grid velocities of the faces.
!
! loop over the three directions.
! the original code is elegant but the tapenade has a difficult time
! to understand it. thus, we unfold it and make it easier for the
! tapenade.
! i-direction
        do k=1,ke
          call pushinteger4(j)
          do j=1,je
            do i=0,ie
! determine the coordinates of the face center,
! which are stored in xc.
              call pushreal8(xc(1))
              xc(1) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, j-1, k&
&               -1, 1)+flowdoms(nn, groundlevel, sps)%x(i, j, k-1, 1)+&
&               flowdoms(nn, groundlevel, sps)%x(i, j-1, k, 1)+flowdoms(&
&               nn, groundlevel, sps)%x(i, j, k, 1))
              call pushreal8(xc(2))
              xc(2) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, j-1, k&
&               -1, 2)+flowdoms(nn, groundlevel, sps)%x(i, j, k-1, 2)+&
&               flowdoms(nn, groundlevel, sps)%x(i, j-1, k, 2)+flowdoms(&
&               nn, groundlevel, sps)%x(i, j, k, 2))
              call pushreal8(xc(3))
              xc(3) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, j-1, k&
&               -1, 3)+flowdoms(nn, groundlevel, sps)%x(i, j, k-1, 3)+&
&               flowdoms(nn, groundlevel, sps)%x(i, j-1, k, 3)+flowdoms(&
&               nn, groundlevel, sps)%x(i, j, k, 3))
              call pushreal8array(sc, 3)
              call cellfacevelocities(xc, rotcenter, rotrate, velxgrid, &
&                               velygrid, velzgrid, derivrotationmatrix&
&                               , sc)
! store the dot product of grid velocity sc and
! the normal ss in sface.
            end do
          end do
        end do
! j-direction
        do k=1,ke
          call pushinteger4(j)
          do j=0,je
            do i=1,ie
! determine the coordinates of the face center,
! which are stored in xc.
              call pushreal8(xc(1))
              xc(1) = fourth*(flowdoms(nn, groundlevel, sps)%x(i-1, j, k&
&               , 1)+flowdoms(nn, groundlevel, sps)%x(i, j, k-1, 1)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, j, k-1, 1)+&
&               flowdoms(nn, groundlevel, sps)%x(i, j, k, 1))
              call pushreal8(xc(2))
              xc(2) = fourth*(flowdoms(nn, groundlevel, sps)%x(i-1, j, k&
&               , 2)+flowdoms(nn, groundlevel, sps)%x(i, j, k-1, 2)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, j, k-1, 2)+&
&               flowdoms(nn, groundlevel, sps)%x(i, j, k, 2))
              call pushreal8(xc(3))
              xc(3) = fourth*(flowdoms(nn, groundlevel, sps)%x(i-1, j, k&
&               , 3)+flowdoms(nn, groundlevel, sps)%x(i, j, k-1, 3)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, j, k-1, 3)+&
&               flowdoms(nn, groundlevel, sps)%x(i, j, k, 3))
              call pushreal8array(sc, 3)
              call cellfacevelocities(xc, rotcenter, rotrate, velxgrid, &
&                               velygrid, velzgrid, derivrotationmatrix&
&                               , sc)
! store the dot product of grid velocity sc and
! the normal ss in sface.
            end do
          end do
        end do
! k-direction
        do k=0,ke
          call pushinteger4(j)
          do j=1,je
            do i=1,ie
! determine the coordinates of the face center,
! which are stored in xc.
              call pushreal8(xc(1))
              xc(1) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, j-1, k&
&               , 1)+flowdoms(nn, groundlevel, sps)%x(i-1, j, k, 1)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, j-1, k, 1)+&
&               flowdoms(nn, groundlevel, sps)%x(i, j, k, 1))
              call pushreal8(xc(2))
              xc(2) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, j-1, k&
&               , 2)+flowdoms(nn, groundlevel, sps)%x(i-1, j, k, 2)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, j-1, k, 2)+&
&               flowdoms(nn, groundlevel, sps)%x(i, j, k, 2))
              call pushreal8(xc(3))
              xc(3) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, j-1, k&
&               , 3)+flowdoms(nn, groundlevel, sps)%x(i-1, j, k, 3)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, j-1, k, 3)+&
&               flowdoms(nn, groundlevel, sps)%x(i, j, k, 3))
              call pushreal8array(sc, 3)
              call cellfacevelocities(xc, rotcenter, rotrate, velxgrid, &
&                               velygrid, velzgrid, derivrotationmatrix&
&                               , sc)
! store the dot product of grid velocity sc and
! the normal ss in sface.
            end do
          end do
        end do
        rotrated = 0.0_8
        velygridd = 0.0_8
        xcd = 0.0_8
        velzgridd = 0.0_8
        derivrotationmatrixd = 0.0_8
        scd = 0.0_8
        velxgridd = 0.0_8
        do k=ke,0,-1
          do j=je,1,-1
            do i=ie,1,-1
              scd(1) = scd(1) + sk(i, j, k, 1)*sfacekd(i, j, k)
              skd(i, j, k, 1) = skd(i, j, k, 1) + sc(1)*sfacekd(i, j, k)
              scd(2) = scd(2) + sk(i, j, k, 2)*sfacekd(i, j, k)
              skd(i, j, k, 2) = skd(i, j, k, 2) + sc(2)*sfacekd(i, j, k)
              scd(3) = scd(3) + sk(i, j, k, 3)*sfacekd(i, j, k)
              skd(i, j, k, 3) = skd(i, j, k, 3) + sc(3)*sfacekd(i, j, k)
              sfacekd(i, j, k) = 0.0_8
              call popreal8array(sc, 3)
              call cellfacevelocities_b(xc, xcd, rotcenter, rotcenterd, &
&                                 rotrate, rotrated, velxgrid, velxgridd&
&                                 , velygrid, velygridd, velzgrid, &
&                                 velzgridd, derivrotationmatrix, &
&                                 derivrotationmatrixd, sc, scd)
              call popreal8(xc(3))
              tempd = fourth*xcd(3)
              xcd(3) = 0.0_8
              flowdomsd(nn, groundlevel, sps)%x(i, j-1, k, 3) = &
&               flowdomsd(nn, groundlevel, sps)%x(i, j-1, k, 3) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i-1, j, k, 3) = &
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, k, 3) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, k, 3) = &
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, k, 3) + &
&               tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j, k, 3) = flowdomsd(&
&               nn, groundlevel, sps)%x(i, j, k, 3) + tempd
              call popreal8(xc(2))
              tempd = fourth*xcd(2)
              xcd(2) = 0.0_8
              flowdomsd(nn, groundlevel, sps)%x(i, j-1, k, 2) = &
&               flowdomsd(nn, groundlevel, sps)%x(i, j-1, k, 2) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i-1, j, k, 2) = &
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, k, 2) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, k, 2) = &
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, k, 2) + &
&               tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j, k, 2) = flowdomsd(&
&               nn, groundlevel, sps)%x(i, j, k, 2) + tempd
              call popreal8(xc(1))
              tempd = fourth*xcd(1)
              xcd(1) = 0.0_8
              flowdomsd(nn, groundlevel, sps)%x(i, j-1, k, 1) = &
&               flowdomsd(nn, groundlevel, sps)%x(i, j-1, k, 1) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i-1, j, k, 1) = &
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, k, 1) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, k, 1) = &
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, k, 1) + &
&               tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j, k, 1) = flowdomsd(&
&               nn, groundlevel, sps)%x(i, j, k, 1) + tempd
            end do
          end do
          call popinteger4(j)
        end do
        do k=ke,1,-1
          do j=je,0,-1
            do i=ie,1,-1
              scd(1) = scd(1) + sj(i, j, k, 1)*sfacejd(i, j, k)
              sjd(i, j, k, 1) = sjd(i, j, k, 1) + sc(1)*sfacejd(i, j, k)
              scd(2) = scd(2) + sj(i, j, k, 2)*sfacejd(i, j, k)
              sjd(i, j, k, 2) = sjd(i, j, k, 2) + sc(2)*sfacejd(i, j, k)
              scd(3) = scd(3) + sj(i, j, k, 3)*sfacejd(i, j, k)
              sjd(i, j, k, 3) = sjd(i, j, k, 3) + sc(3)*sfacejd(i, j, k)
              sfacejd(i, j, k) = 0.0_8
              call popreal8array(sc, 3)
              call cellfacevelocities_b(xc, xcd, rotcenter, rotcenterd, &
&                                 rotrate, rotrated, velxgrid, velxgridd&
&                                 , velygrid, velygridd, velzgrid, &
&                                 velzgridd, derivrotationmatrix, &
&                                 derivrotationmatrixd, sc, scd)
              call popreal8(xc(3))
              tempd = fourth*xcd(3)
              xcd(3) = 0.0_8
              flowdomsd(nn, groundlevel, sps)%x(i-1, j, k, 3) = &
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, k, 3) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 3) = &
&               flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 3) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i-1, j, k-1, 3) = &
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, k-1, 3) + &
&               tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j, k, 3) = flowdomsd(&
&               nn, groundlevel, sps)%x(i, j, k, 3) + tempd
              call popreal8(xc(2))
              tempd = fourth*xcd(2)
              xcd(2) = 0.0_8
              flowdomsd(nn, groundlevel, sps)%x(i-1, j, k, 2) = &
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, k, 2) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 2) = &
&               flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 2) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i-1, j, k-1, 2) = &
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, k-1, 2) + &
&               tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j, k, 2) = flowdomsd(&
&               nn, groundlevel, sps)%x(i, j, k, 2) + tempd
              call popreal8(xc(1))
              tempd = fourth*xcd(1)
              xcd(1) = 0.0_8
              flowdomsd(nn, groundlevel, sps)%x(i-1, j, k, 1) = &
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, k, 1) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 1) = &
&               flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 1) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i-1, j, k-1, 1) = &
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, k-1, 1) + &
&               tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j, k, 1) = flowdomsd(&
&               nn, groundlevel, sps)%x(i, j, k, 1) + tempd
            end do
          end do
          call popinteger4(j)
        end do
        do k=ke,1,-1
          do j=je,1,-1
            do i=ie,0,-1
              scd(1) = scd(1) + si(i, j, k, 1)*sfaceid(i, j, k)
              sid(i, j, k, 1) = sid(i, j, k, 1) + sc(1)*sfaceid(i, j, k)
              scd(2) = scd(2) + si(i, j, k, 2)*sfaceid(i, j, k)
              sid(i, j, k, 2) = sid(i, j, k, 2) + sc(2)*sfaceid(i, j, k)
              scd(3) = scd(3) + si(i, j, k, 3)*sfaceid(i, j, k)
              sid(i, j, k, 3) = sid(i, j, k, 3) + sc(3)*sfaceid(i, j, k)
              sfaceid(i, j, k) = 0.0_8
              call popreal8array(sc, 3)
              call cellfacevelocities_b(xc, xcd, rotcenter, rotcenterd, &
&                                 rotrate, rotrated, velxgrid, velxgridd&
&                                 , velygrid, velygridd, velzgrid, &
&                                 velzgridd, derivrotationmatrix, &
&                                 derivrotationmatrixd, sc, scd)
              call popreal8(xc(3))
              tempd = fourth*xcd(3)
              xcd(3) = 0.0_8
              flowdomsd(nn, groundlevel, sps)%x(i, j-1, k-1, 3) = &
&               flowdomsd(nn, groundlevel, sps)%x(i, j-1, k-1, 3) + &
&               tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 3) = &
&               flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 3) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j-1, k, 3) = &
&               flowdomsd(nn, groundlevel, sps)%x(i, j-1, k, 3) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j, k, 3) = flowdomsd(&
&               nn, groundlevel, sps)%x(i, j, k, 3) + tempd
              call popreal8(xc(2))
              tempd = fourth*xcd(2)
              xcd(2) = 0.0_8
              flowdomsd(nn, groundlevel, sps)%x(i, j-1, k-1, 2) = &
&               flowdomsd(nn, groundlevel, sps)%x(i, j-1, k-1, 2) + &
&               tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 2) = &
&               flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 2) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j-1, k, 2) = &
&               flowdomsd(nn, groundlevel, sps)%x(i, j-1, k, 2) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j, k, 2) = flowdomsd(&
&               nn, groundlevel, sps)%x(i, j, k, 2) + tempd
              call popreal8(xc(1))
              tempd = fourth*xcd(1)
              xcd(1) = 0.0_8
              flowdomsd(nn, groundlevel, sps)%x(i, j-1, k-1, 1) = &
&               flowdomsd(nn, groundlevel, sps)%x(i, j-1, k-1, 1) + &
&               tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 1) = &
&               flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 1) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j-1, k, 1) = &
&               flowdomsd(nn, groundlevel, sps)%x(i, j-1, k, 1) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j, k, 1) = flowdomsd(&
&               nn, groundlevel, sps)%x(i, j, k, 1) + tempd
            end do
          end do
          call popinteger4(j)
        end do
        xxcd = 0.0_8
        do k=ke,1,-1
          do j=je,1,-1
            do i=ie,1,-1
              scd(3) = scd(3) + sd(i, j, k, 3)
              velzgridd = velzgridd + sd(i, j, k, 3)
              derivrotationmatrixd(3, 1) = derivrotationmatrixd(3, 1) + &
&               xxc(1)*sd(i, j, k, 3)
              xxcd(1) = xxcd(1) + derivrotationmatrix(3, 1)*sd(i, j, k, &
&               3)
              derivrotationmatrixd(3, 2) = derivrotationmatrixd(3, 2) + &
&               xxc(2)*sd(i, j, k, 3)
              xxcd(2) = xxcd(2) + derivrotationmatrix(3, 2)*sd(i, j, k, &
&               3)
              derivrotationmatrixd(3, 3) = derivrotationmatrixd(3, 3) + &
&               xxc(3)*sd(i, j, k, 3)
              xxcd(3) = xxcd(3) + derivrotationmatrix(3, 3)*sd(i, j, k, &
&               3)
              sd(i, j, k, 3) = 0.0_8
              scd(2) = scd(2) + sd(i, j, k, 2)
              velygridd = velygridd + sd(i, j, k, 2)
              derivrotationmatrixd(2, 1) = derivrotationmatrixd(2, 1) + &
&               xxc(1)*sd(i, j, k, 2)
              xxcd(1) = xxcd(1) + derivrotationmatrix(2, 1)*sd(i, j, k, &
&               2)
              derivrotationmatrixd(2, 2) = derivrotationmatrixd(2, 2) + &
&               xxc(2)*sd(i, j, k, 2)
              xxcd(2) = xxcd(2) + derivrotationmatrix(2, 2)*sd(i, j, k, &
&               2)
              derivrotationmatrixd(2, 3) = derivrotationmatrixd(2, 3) + &
&               xxc(3)*sd(i, j, k, 2)
              xxcd(3) = xxcd(3) + derivrotationmatrix(2, 3)*sd(i, j, k, &
&               2)
              sd(i, j, k, 2) = 0.0_8
              scd(1) = scd(1) + sd(i, j, k, 1)
              velxgridd = velxgridd + sd(i, j, k, 1)
              derivrotationmatrixd(1, 1) = derivrotationmatrixd(1, 1) + &
&               xxc(1)*sd(i, j, k, 1)
              xxcd(1) = xxcd(1) + derivrotationmatrix(1, 1)*sd(i, j, k, &
&               1)
              derivrotationmatrixd(1, 2) = derivrotationmatrixd(1, 2) + &
&               xxc(2)*sd(i, j, k, 1)
              xxcd(2) = xxcd(2) + derivrotationmatrix(1, 2)*sd(i, j, k, &
&               1)
              derivrotationmatrixd(1, 3) = derivrotationmatrixd(1, 3) + &
&               xxc(3)*sd(i, j, k, 1)
              xxcd(3) = xxcd(3) + derivrotationmatrix(1, 3)*sd(i, j, k, &
&               1)
              sd(i, j, k, 1) = 0.0_8
              call popreal8(xxc(3))
              xcd(3) = xcd(3) + xxcd(3)
              xxcd(3) = rotrate(2)*scd(1) - rotrate(1)*scd(2)
              call popreal8(xxc(2))
              xcd(2) = xcd(2) + xxcd(2)
              xxcd(2) = rotrate(1)*scd(3) - rotrate(3)*scd(1)
              call popreal8(xxc(1))
              xcd(1) = xcd(1) + xxcd(1)
              xxcd(1) = rotrate(3)*scd(2) - rotrate(2)*scd(3)
              rotrated(1) = rotrated(1) + xxc(2)*scd(3) - xxc(3)*scd(2)
              rotrated(2) = rotrated(2) + xxc(3)*scd(1) - xxc(1)*scd(3)
              scd(3) = 0.0_8
              rotrated(3) = rotrated(3) + xxc(1)*scd(2) - xxc(2)*scd(1)
              scd(2) = 0.0_8
              scd(1) = 0.0_8
              call popreal8(xxc(3))
              xcd(3) = xcd(3) + xxcd(3)
              xxcd(3) = 0.0_8
              call popreal8(xxc(2))
              xcd(2) = xcd(2) + xxcd(2)
              xxcd(2) = 0.0_8
              call popreal8(xxc(1))
              xcd(1) = xcd(1) + xxcd(1)
              xxcd(1) = 0.0_8
              tempd = eighth*xcd(3)
              xcd(3) = 0.0_8
              flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, k-1, 3) = &
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, k-1, 3) + &
&               tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j-1, k-1, 3) = &
&               flowdomsd(nn, groundlevel, sps)%x(i, j-1, k-1, 3) + &
&               tempd
              flowdomsd(nn, groundlevel, sps)%x(i-1, j, k-1, 3) = &
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, k-1, 3) + &
&               tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 3) = &
&               flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 3) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, k, 3) = &
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, k, 3) + &
&               tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j-1, k, 3) = &
&               flowdomsd(nn, groundlevel, sps)%x(i, j-1, k, 3) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i-1, j, k, 3) = &
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, k, 3) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j, k, 3) = flowdomsd(&
&               nn, groundlevel, sps)%x(i, j, k, 3) + tempd
              tempd = eighth*xcd(2)
              xcd(2) = 0.0_8
              flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, k-1, 2) = &
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, k-1, 2) + &
&               tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j-1, k-1, 2) = &
&               flowdomsd(nn, groundlevel, sps)%x(i, j-1, k-1, 2) + &
&               tempd
              flowdomsd(nn, groundlevel, sps)%x(i-1, j, k-1, 2) = &
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, k-1, 2) + &
&               tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 2) = &
&               flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 2) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, k, 2) = &
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, k, 2) + &
&               tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j-1, k, 2) = &
&               flowdomsd(nn, groundlevel, sps)%x(i, j-1, k, 2) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i-1, j, k, 2) = &
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, k, 2) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j, k, 2) = flowdomsd(&
&               nn, groundlevel, sps)%x(i, j, k, 2) + tempd
              tempd = eighth*xcd(1)
              xcd(1) = 0.0_8
              flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, k-1, 1) = &
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, k-1, 1) + &
&               tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j-1, k-1, 1) = &
&               flowdomsd(nn, groundlevel, sps)%x(i, j-1, k-1, 1) + &
&               tempd
              flowdomsd(nn, groundlevel, sps)%x(i-1, j, k-1, 1) = &
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, k-1, 1) + &
&               tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 1) = &
&               flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 1) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, k, 1) = &
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, k, 1) + &
&               tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j-1, k, 1) = &
&               flowdomsd(nn, groundlevel, sps)%x(i, j-1, k, 1) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i-1, j, k, 1) = &
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, k, 1) + tempd
              flowdomsd(nn, groundlevel, sps)%x(i, j, k, 1) = flowdomsd(&
&               nn, groundlevel, sps)%x(i, j, k, 1) + tempd
            end do
          end do
          call popinteger4(j)
        end do
        velzgrid0d = velzgridd
        velygrid0d = velygridd
        velxgrid0d = velxgridd
        timerefd = sum(cgnsdoms(j)%rotrate*rotrated)
      end if
    else
      timerefd = 0.0_8
      velxgrid0d = 0.0_8
      velzgrid0d = 0.0_8
      derivrotationmatrixd = 0.0_8
      velygrid0d = 0.0_8
    end if
    call derivativerotmatrixrigid_b(derivrotationmatrix, &
&                             derivrotationmatrixd, rotationpoint, t(1))
    ainfd = -(machgrid*veldirfreestream(3)*velzgrid0d) - machgrid*&
&     veldirfreestream(2)*velygrid0d - machgrid*veldirfreestream(1)*&
&     velxgrid0d
    machgridd = -(ainf*veldirfreestream(3)*velzgrid0d) - ainf*&
&     veldirfreestream(2)*velygrid0d - ainf*veldirfreestream(1)*&
&     velxgrid0d
    veldirfreestreamd(3) = veldirfreestreamd(3) - ainf*machgrid*&
&     velzgrid0d
    veldirfreestreamd(2) = veldirfreestreamd(2) - ainf*machgrid*&
&     velygrid0d
    veldirfreestreamd(1) = veldirfreestreamd(1) - ainf*machgrid*&
&     velxgrid0d
    if (gammainf*(pinf/rhoinf) .eq. 0.0_8) then
      tempd = 0.0_8
    else
      tempd = gammainf*ainfd/(rhoinf*2.0*sqrt(gammainf*(pinf/rhoinf)))
    end if
    pinfd = tempd
    rhoinfd = -(pinf*tempd/rhoinf)
  end subroutine gridvelocitiesfinelevel_block_b
 
  subroutine gridvelocitiesfinelevel_block(useoldcoor, t, sps, nn)
!
!       gridvelocitiesfinelevel computes the grid velocities for
!       the cell centers and the normal grid velocities for the faces
!       of moving blocks for the currently finest grid, i.e.
!       groundlevel. the velocities are computed at time t for
!       spectral mode sps. if useoldcoor is .true. the velocities
!       are determined using the unsteady time integrator in
!       combination with the old coordinates; otherwise the analytic
!       form is used.
!
    use blockpointers
    use cgnsgrid
    use flowvarrefstate
    use inputmotion
    use inputunsteady
    use iteration
    use inputphysics
    use inputtsstabderiv
    use monitor
    use communication
    use flowutils_b, only : derivativerotmatrixrigid, getdirvector
    use utils_b, only : setcoeftimeintegrator, tsalpha, tsbeta, tsmach, &
&   terminate, rotmatrixrigidbody, getdirangle
    implicit none
!
!      subroutine arguments.
!
    integer(kind=inttype), intent(in) :: sps, nn
    logical, intent(in) :: useoldcoor
    real(kind=realtype), dimension(*), intent(in) :: t
!
!      local variables.
!
    integer(kind=inttype) :: mm
    integer(kind=inttype) :: i, j, k, ii, iie, jje, kke
    real(kind=realtype) :: oneover4dt, oneover8dt
    real(kind=realtype) :: velxgrid, velygrid, velzgrid, ainf
    real(kind=realtype) :: velxgrid0, velygrid0, velzgrid0
    real(kind=realtype), dimension(3) :: sc, xc, xxc
    real(kind=realtype), dimension(3) :: rotcenter, rotrate
    real(kind=realtype), dimension(3) :: rotationpoint
    real(kind=realtype), dimension(3, 3) :: rotationmatrix, &
&   derivrotationmatrix
    real(kind=realtype) :: tnew, told
    real(kind=realtype), dimension(:, :), pointer :: sface
    real(kind=realtype), dimension(:, :, :), pointer :: xx, ss
    real(kind=realtype), dimension(:, :, :, :), pointer :: xxold
    real(kind=realtype) :: intervalmach, alphats, alphaincrement, betats&
&   , betaincrement
    real(kind=realtype), dimension(3) :: veldir
    real(kind=realtype), dimension(3) :: refdirection
    intrinsic sqrt
! compute the mesh velocity from the given mesh mach number.
! vel{x,y,z}grid0 is the actual velocity you want at the
! geometry.
    ainf = sqrt(gammainf*pinf/rhoinf)
    velxgrid0 = ainf*machgrid*(-veldirfreestream(1))
    velygrid0 = ainf*machgrid*(-veldirfreestream(2))
    velzgrid0 = ainf*machgrid*(-veldirfreestream(3))
! compute the derivative of the rotation matrix and the rotation
! point; needed for velocity due to the rigid body rotation of
! the entire grid. it is assumed that the rigid body motion of
! the grid is only specified if there is only 1 section present.
    call derivativerotmatrixrigid(derivrotationmatrix, rotationpoint, t(&
&                           1))
!compute the rotation matrix to update the velocities for the time
!spectral stability derivative case...
    if (blockismoving) then
! determine the situation we are having here.
      if (.not.useoldcoor) then
!
!             the velocities must be determined analytically.
!
! store the rotation center and determine the
! nondimensional rotation rate of this block. as the
! reference length is 1 timeref == 1/uref and at the end
! the nondimensional velocity is computed.
        j = nbkglobal
        rotcenter = cgnsdoms(j)%rotcenter
        rotrate = timeref*cgnsdoms(j)%rotrate
        velxgrid = velxgrid0
        velygrid = velygrid0
        velzgrid = velzgrid0
!
!             grid velocities of the cell centers, including the
!             1st level halo cells.
!
! loop over the cells, including the 1st level halo's.
        do k=1,ke
          do j=1,je
            do i=1,ie
! determine the coordinates of the cell center,
! which are stored in xc.
              xc(1) = eighth*(flowdoms(nn, groundlevel, sps)%x(i-1, j-1&
&               , k-1, 1)+flowdoms(nn, groundlevel, sps)%x(i, j-1, k-1, &
&               1)+flowdoms(nn, groundlevel, sps)%x(i-1, j, k-1, 1)+&
&               flowdoms(nn, groundlevel, sps)%x(i, j, k-1, 1)+flowdoms(&
&               nn, groundlevel, sps)%x(i-1, j-1, k, 1)+flowdoms(nn, &
&               groundlevel, sps)%x(i, j-1, k, 1)+flowdoms(nn, &
&               groundlevel, sps)%x(i-1, j, k, 1)+flowdoms(nn, &
&               groundlevel, sps)%x(i, j, k, 1))
              xc(2) = eighth*(flowdoms(nn, groundlevel, sps)%x(i-1, j-1&
&               , k-1, 2)+flowdoms(nn, groundlevel, sps)%x(i, j-1, k-1, &
&               2)+flowdoms(nn, groundlevel, sps)%x(i-1, j, k-1, 2)+&
&               flowdoms(nn, groundlevel, sps)%x(i, j, k-1, 2)+flowdoms(&
&               nn, groundlevel, sps)%x(i-1, j-1, k, 2)+flowdoms(nn, &
&               groundlevel, sps)%x(i, j-1, k, 2)+flowdoms(nn, &
&               groundlevel, sps)%x(i-1, j, k, 2)+flowdoms(nn, &
&               groundlevel, sps)%x(i, j, k, 2))
              xc(3) = eighth*(flowdoms(nn, groundlevel, sps)%x(i-1, j-1&
&               , k-1, 3)+flowdoms(nn, groundlevel, sps)%x(i, j-1, k-1, &
&               3)+flowdoms(nn, groundlevel, sps)%x(i-1, j, k-1, 3)+&
&               flowdoms(nn, groundlevel, sps)%x(i, j, k-1, 3)+flowdoms(&
&               nn, groundlevel, sps)%x(i-1, j-1, k, 3)+flowdoms(nn, &
&               groundlevel, sps)%x(i, j-1, k, 3)+flowdoms(nn, &
&               groundlevel, sps)%x(i-1, j, k, 3)+flowdoms(nn, &
&               groundlevel, sps)%x(i, j, k, 3))
! determine the coordinates relative to the
! center of rotation.
              xxc(1) = xc(1) - rotcenter(1)
              xxc(2) = xc(2) - rotcenter(2)
              xxc(3) = xc(3) - rotcenter(3)
! determine the rotation speed of the cell center,
! which is omega*r.
              sc(1) = rotrate(2)*xxc(3) - rotrate(3)*xxc(2)
              sc(2) = rotrate(3)*xxc(1) - rotrate(1)*xxc(3)
              sc(3) = rotrate(1)*xxc(2) - rotrate(2)*xxc(1)
! determine the coordinates relative to the
! rigid body rotation point.
              xxc(1) = xc(1) - rotationpoint(1)
              xxc(2) = xc(2) - rotationpoint(2)
              xxc(3) = xc(3) - rotationpoint(3)
! determine the total velocity of the cell center.
! this is a combination of rotation speed of this
! block and the entire rigid body rotation.
              s(i, j, k, 1) = sc(1) + velxgrid + derivrotationmatrix(1, &
&               1)*xxc(1) + derivrotationmatrix(1, 2)*xxc(2) + &
&               derivrotationmatrix(1, 3)*xxc(3)
              s(i, j, k, 2) = sc(2) + velygrid + derivrotationmatrix(2, &
&               1)*xxc(1) + derivrotationmatrix(2, 2)*xxc(2) + &
&               derivrotationmatrix(2, 3)*xxc(3)
              s(i, j, k, 3) = sc(3) + velzgrid + derivrotationmatrix(3, &
&               1)*xxc(1) + derivrotationmatrix(3, 2)*xxc(2) + &
&               derivrotationmatrix(3, 3)*xxc(3)
            end do
          end do
        end do
!
!             normal grid velocities of the faces.
!
! loop over the three directions.
! the original code is elegant but the tapenade has a difficult time
! to understand it. thus, we unfold it and make it easier for the
! tapenade.
! i-direction
        do k=1,ke
          do j=1,je
            do i=0,ie
! determine the coordinates of the face center,
! which are stored in xc.
              xc(1) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, j-1, k&
&               -1, 1)+flowdoms(nn, groundlevel, sps)%x(i, j, k-1, 1)+&
&               flowdoms(nn, groundlevel, sps)%x(i, j-1, k, 1)+flowdoms(&
&               nn, groundlevel, sps)%x(i, j, k, 1))
              xc(2) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, j-1, k&
&               -1, 2)+flowdoms(nn, groundlevel, sps)%x(i, j, k-1, 2)+&
&               flowdoms(nn, groundlevel, sps)%x(i, j-1, k, 2)+flowdoms(&
&               nn, groundlevel, sps)%x(i, j, k, 2))
              xc(3) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, j-1, k&
&               -1, 3)+flowdoms(nn, groundlevel, sps)%x(i, j, k-1, 3)+&
&               flowdoms(nn, groundlevel, sps)%x(i, j-1, k, 3)+flowdoms(&
&               nn, groundlevel, sps)%x(i, j, k, 3))
              call cellfacevelocities(xc, rotcenter, rotrate, velxgrid, &
&                               velygrid, velzgrid, derivrotationmatrix&
&                               , sc)
! store the dot product of grid velocity sc and
! the normal ss in sface.
              sfacei(i, j, k) = sc(1)*si(i, j, k, 1) + sc(2)*si(i, j, k&
&               , 2) + sc(3)*si(i, j, k, 3)
            end do
          end do
        end do
! j-direction
        do k=1,ke
          do j=0,je
            do i=1,ie
! determine the coordinates of the face center,
! which are stored in xc.
              xc(1) = fourth*(flowdoms(nn, groundlevel, sps)%x(i-1, j, k&
&               , 1)+flowdoms(nn, groundlevel, sps)%x(i, j, k-1, 1)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, j, k-1, 1)+&
&               flowdoms(nn, groundlevel, sps)%x(i, j, k, 1))
              xc(2) = fourth*(flowdoms(nn, groundlevel, sps)%x(i-1, j, k&
&               , 2)+flowdoms(nn, groundlevel, sps)%x(i, j, k-1, 2)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, j, k-1, 2)+&
&               flowdoms(nn, groundlevel, sps)%x(i, j, k, 2))
              xc(3) = fourth*(flowdoms(nn, groundlevel, sps)%x(i-1, j, k&
&               , 3)+flowdoms(nn, groundlevel, sps)%x(i, j, k-1, 3)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, j, k-1, 3)+&
&               flowdoms(nn, groundlevel, sps)%x(i, j, k, 3))
              call cellfacevelocities(xc, rotcenter, rotrate, velxgrid, &
&                               velygrid, velzgrid, derivrotationmatrix&
&                               , sc)
! store the dot product of grid velocity sc and
! the normal ss in sface.
              sfacej(i, j, k) = sc(1)*sj(i, j, k, 1) + sc(2)*sj(i, j, k&
&               , 2) + sc(3)*sj(i, j, k, 3)
            end do
          end do
        end do
! k-direction
        do k=0,ke
          do j=1,je
            do i=1,ie
! determine the coordinates of the face center,
! which are stored in xc.
              xc(1) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, j-1, k&
&               , 1)+flowdoms(nn, groundlevel, sps)%x(i-1, j, k, 1)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, j-1, k, 1)+&
&               flowdoms(nn, groundlevel, sps)%x(i, j, k, 1))
              xc(2) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, j-1, k&
&               , 2)+flowdoms(nn, groundlevel, sps)%x(i-1, j, k, 2)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, j-1, k, 2)+&
&               flowdoms(nn, groundlevel, sps)%x(i, j, k, 2))
              xc(3) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, j-1, k&
&               , 3)+flowdoms(nn, groundlevel, sps)%x(i-1, j, k, 3)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, j-1, k, 3)+&
&               flowdoms(nn, groundlevel, sps)%x(i, j, k, 3))
              call cellfacevelocities(xc, rotcenter, rotrate, velxgrid, &
&                               velygrid, velzgrid, derivrotationmatrix&
&                               , sc)
! store the dot product of grid velocity sc and
! the normal ss in sface.
              sfacek(i, j, k) = sc(1)*sk(i, j, k, 1) + sc(2)*sk(i, j, k&
&               , 2) + sc(3)*sk(i, j, k, 3)
            end do
          end do
        end do
      end if
    end if
  end subroutine gridvelocitiesfinelevel_block
 
!  differentiation of cellfacevelocities in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: rotrate velygrid xc velzgrid
!                derivrotationmatrix sc velxgrid
!   with respect to varying inputs: rotrate velygrid xc velzgrid
!                derivrotationmatrix sc rotcenter velxgrid
!   rw status of diff variables: rotrate:incr velygrid:incr xc:incr
!                velzgrid:incr derivrotationmatrix:incr sc:in-out
!                rotcenter:out velxgrid:incr
  subroutine cellfacevelocities_b(xc, xcd, rotcenter, rotcenterd, &
&   rotrate, rotrated, velxgrid, velxgridd, velygrid, velygridd, &
&   velzgrid, velzgridd, derivrotationmatrix, derivrotationmatrixd, sc, &
&   scd)
!
!  returns the cell face velocities for a given face center
!
    use constants
    implicit none
!
!      subroutine arguments.
!
    real(kind=realtype), dimension(3), intent(in) :: xc, rotcenter, &
&   rotrate
    real(kind=realtype), dimension(3) :: xcd, rotcenterd, rotrated
    real(kind=realtype), intent(in) :: velxgrid, velygrid, velzgrid
    real(kind=realtype) :: velxgridd, velygridd, velzgridd
    real(kind=realtype), dimension(3, 3), intent(in) :: &
&   derivrotationmatrix
    real(kind=realtype), dimension(3, 3) :: derivrotationmatrixd
    real(kind=realtype), dimension(3) :: sc
    real(kind=realtype), dimension(3) :: scd
!
!      local variables.
!
    real(kind=realtype), dimension(3) :: rotationpoint, xxc
    real(kind=realtype), dimension(3) :: xxcd
! determine the coordinates relative to the
! center of rotation.
    xxc(1) = xc(1) - rotcenter(1)
    xxc(2) = xc(2) - rotcenter(2)
    xxc(3) = xc(3) - rotcenter(3)
! determine the rotation speed of the face center,
! which is omega*r.
! determine the coordinates relative to the
! rigid body rotation point.
    call pushreal8(xxc(1))
    xxc(1) = xc(1) - rotationpoint(1)
    call pushreal8(xxc(2))
    xxc(2) = xc(2) - rotationpoint(2)
    call pushreal8(xxc(3))
    xxc(3) = xc(3) - rotationpoint(3)
! determine the total velocity of the cell face.
! this is a combination of rotation speed of this
! block and the entire rigid body rotation.
    xxcd = 0.0_8
    velzgridd = velzgridd + scd(3)
    derivrotationmatrixd(3, 1) = derivrotationmatrixd(3, 1) + xxc(1)*scd&
&     (3)
    xxcd(1) = xxcd(1) + derivrotationmatrix(3, 1)*scd(3) + &
&     derivrotationmatrix(2, 1)*scd(2) + derivrotationmatrix(1, 1)*scd(1&
&     )
    derivrotationmatrixd(3, 2) = derivrotationmatrixd(3, 2) + xxc(2)*scd&
&     (3)
    xxcd(2) = xxcd(2) + derivrotationmatrix(3, 2)*scd(3) + &
&     derivrotationmatrix(2, 2)*scd(2) + derivrotationmatrix(1, 2)*scd(1&
&     )
    derivrotationmatrixd(3, 3) = derivrotationmatrixd(3, 3) + xxc(3)*scd&
&     (3)
    xxcd(3) = xxcd(3) + derivrotationmatrix(3, 3)*scd(3) + &
&     derivrotationmatrix(2, 3)*scd(2) + derivrotationmatrix(1, 3)*scd(1&
&     )
    velygridd = velygridd + scd(2)
    derivrotationmatrixd(2, 1) = derivrotationmatrixd(2, 1) + xxc(1)*scd&
&     (2)
    derivrotationmatrixd(2, 2) = derivrotationmatrixd(2, 2) + xxc(2)*scd&
&     (2)
    derivrotationmatrixd(2, 3) = derivrotationmatrixd(2, 3) + xxc(3)*scd&
&     (2)
    velxgridd = velxgridd + scd(1)
    derivrotationmatrixd(1, 1) = derivrotationmatrixd(1, 1) + xxc(1)*scd&
&     (1)
    derivrotationmatrixd(1, 2) = derivrotationmatrixd(1, 2) + xxc(2)*scd&
&     (1)
    derivrotationmatrixd(1, 3) = derivrotationmatrixd(1, 3) + xxc(3)*scd&
&     (1)
    call popreal8(xxc(3))
    xcd(3) = xcd(3) + xxcd(3)
    xxcd(3) = rotrate(2)*scd(1) - rotrate(1)*scd(2)
    call popreal8(xxc(2))
    xcd(2) = xcd(2) + xxcd(2)
    xxcd(2) = rotrate(1)*scd(3) - rotrate(3)*scd(1)
    call popreal8(xxc(1))
    xcd(1) = xcd(1) + xxcd(1)
    xxcd(1) = rotrate(3)*scd(2) - rotrate(2)*scd(3)
    rotrated(1) = rotrated(1) + xxc(2)*scd(3) - xxc(3)*scd(2)
    rotrated(2) = rotrated(2) + xxc(3)*scd(1) - xxc(1)*scd(3)
    scd(3) = 0.0_8
    rotrated(3) = rotrated(3) + xxc(1)*scd(2) - xxc(2)*scd(1)
    scd(2) = 0.0_8
    scd(1) = 0.0_8
    rotcenterd = 0.0_8
    xcd(3) = xcd(3) + xxcd(3)
    rotcenterd(3) = rotcenterd(3) - xxcd(3)
    xxcd(3) = 0.0_8
    xcd(2) = xcd(2) + xxcd(2)
    rotcenterd(2) = rotcenterd(2) - xxcd(2)
    xxcd(2) = 0.0_8
    xcd(1) = xcd(1) + xxcd(1)
    rotcenterd(1) = rotcenterd(1) - xxcd(1)
  end subroutine cellfacevelocities_b
 
  subroutine cellfacevelocities(xc, rotcenter, rotrate, velxgrid, &
&   velygrid, velzgrid, derivrotationmatrix, sc)
!
!  returns the cell face velocities for a given face center
!
    use constants
    implicit none
!
!      subroutine arguments.
!
    real(kind=realtype), dimension(3), intent(in) :: xc, rotcenter, &
&   rotrate
    real(kind=realtype), intent(in) :: velxgrid, velygrid, velzgrid
    real(kind=realtype), dimension(3, 3), intent(in) :: &
&   derivrotationmatrix
    real(kind=realtype), dimension(3), intent(out) :: sc
!
!      local variables.
!
    real(kind=realtype), dimension(3) :: rotationpoint, xxc
! determine the coordinates relative to the
! center of rotation.
    xxc(1) = xc(1) - rotcenter(1)
    xxc(2) = xc(2) - rotcenter(2)
    xxc(3) = xc(3) - rotcenter(3)
! determine the rotation speed of the face center,
! which is omega*r.
    sc(1) = rotrate(2)*xxc(3) - rotrate(3)*xxc(2)
    sc(2) = rotrate(3)*xxc(1) - rotrate(1)*xxc(3)
    sc(3) = rotrate(1)*xxc(2) - rotrate(2)*xxc(1)
! determine the coordinates relative to the
! rigid body rotation point.
    xxc(1) = xc(1) - rotationpoint(1)
    xxc(2) = xc(2) - rotationpoint(2)
    xxc(3) = xc(3) - rotationpoint(3)
! determine the total velocity of the cell face.
! this is a combination of rotation speed of this
! block and the entire rigid body rotation.
    sc(1) = sc(1) + velxgrid + derivrotationmatrix(1, 1)*xxc(1) + &
&     derivrotationmatrix(1, 2)*xxc(2) + derivrotationmatrix(1, 3)*xxc(3&
&     )
    sc(2) = sc(2) + velygrid + derivrotationmatrix(2, 1)*xxc(1) + &
&     derivrotationmatrix(2, 2)*xxc(2) + derivrotationmatrix(2, 3)*xxc(3&
&     )
    sc(3) = sc(3) + velzgrid + derivrotationmatrix(3, 1)*xxc(1) + &
&     derivrotationmatrix(3, 2)*xxc(2) + derivrotationmatrix(3, 3)*xxc(3&
&     )
  end subroutine cellfacevelocities
 
!  differentiation of slipvelocitiesfinelevel_block in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: *(flowdoms.x) veldirfreestream
!                *(*bcdata.uslip)
!   with respect to varying inputs: pinf timeref rhoinf *(flowdoms.x)
!                veldirfreestream machgrid *(*bcdata.uslip)
!   rw status of diff variables: pinf:out timeref:out rhoinf:out
!                *(flowdoms.x):incr veldirfreestream:incr machgrid:out
!                *(*bcdata.uslip):in-out
!   plus diff mem management of: flowdoms.x:in bcdata:in *bcdata.uslip:in
  subroutine slipvelocitiesfinelevel_block_b(useoldcoor, t, sps, nn)
!
!       slipvelocitiesfinelevel computes the slip velocities for
!       viscous subfaces on all viscous boundaries on groundlevel for
!       the given spectral solution. if useoldcoor is .true. the
!       velocities are determined using the unsteady time integrator;
!       otherwise the analytic form is used.
!
    use constants
    use inputtimespectral
    use blockpointers
    use cgnsgrid
    use flowvarrefstate
    use inputmotion
    use inputunsteady
    use iteration
    use inputphysics
    use inputtsstabderiv
    use monitor
    use communication
    use flowutils_b, only : derivativerotmatrixrigid, &
&   derivativerotmatrixrigid_b, getdirvector
    use utils_b, only : tsalpha, tsbeta, tsmach, terminate, &
&   rotmatrixrigidbody, setcoeftimeintegrator, getdirangle
    implicit none
!
!      subroutine arguments.
!
    integer(kind=inttype), intent(in) :: sps, nn
    logical, intent(in) :: useoldcoor
    real(kind=realtype), dimension(*), intent(in) :: t
!
!      local variables.
!
    integer(kind=inttype) :: mm, i, j, level, ii
    real(kind=realtype) :: oneover4dt
    real(kind=realtype) :: velxgrid, velygrid, velzgrid, ainf
    real(kind=realtype) :: velxgridd, velygridd, velzgridd, ainfd
    real(kind=realtype) :: velxgrid0, velygrid0, velzgrid0
    real(kind=realtype) :: velxgrid0d, velygrid0d, velzgrid0d
    real(kind=realtype), dimension(3) :: xc, xxc
    real(kind=realtype), dimension(3) :: xcd, xxcd
    real(kind=realtype), dimension(3) :: rotcenter, rotrate
    real(kind=realtype), dimension(3) :: rotrated
    real(kind=realtype), dimension(3) :: rotationpoint
    real(kind=realtype), dimension(3, 3) :: rotationmatrix, &
&   derivrotationmatrix
    real(kind=realtype), dimension(3, 3) :: derivrotationmatrixd
    real(kind=realtype) :: tnew, told
    real(kind=realtype), dimension(:, :, :), pointer :: uslip
    real(kind=realtype), dimension(:, :, :), pointer :: xface
    real(kind=realtype), dimension(:, :, :, :), pointer :: xfaceold
    real(kind=realtype) :: intervalmach, alphats, alphaincrement, betats&
&   , betaincrement
    real(kind=realtype), dimension(3) :: veldir
    real(kind=realtype), dimension(3) :: refdirection
    intrinsic sqrt
    real(kind=realtype) :: tempd
    integer :: ad_from
    integer :: ad_to
    integer :: ad_from0
    integer :: ad_to0
    integer :: ad_from1
    integer :: ad_to1
    integer :: ad_from2
    integer :: ad_to2
    integer :: ad_from3
    integer :: ad_to3
    integer :: ad_from4
    integer :: ad_to4
    integer :: ad_from5
    integer :: ad_to5
    integer :: ad_from6
    integer :: ad_to6
    integer :: ad_from7
    integer :: ad_to7
    integer :: ad_from8
    integer :: ad_to8
    integer :: ad_from9
    integer :: ad_to9
    integer :: ad_from10
    integer :: ad_to10
    integer :: branch
! determine the situation we are having here.
    if (useoldcoor) then
      pinfd = 0.0_8
      timerefd = 0.0_8
      rhoinfd = 0.0_8
      machgridd = 0.0_8
    else
! the velocities must be determined analytically.
! compute the mesh velocity from the given mesh mach number.
!  ainf = sqrt(gammainf*pinf/rhoinf)
!  velxgrid = ainf*machgrid(1)
!  velygrid = ainf*machgrid(2)
!  velzgrid = ainf*machgrid(3)
      ainf = sqrt(gammainf*pinf/rhoinf)
! compute the derivative of the rotation matrix and the rotation
! point; needed for velocity due to the rigid body rotation of
! the entire grid. it is assumed that the rigid body motion of
! the grid is only specified if there is only 1 section present.
      call derivativerotmatrixrigid(derivrotationmatrix, rotationpoint, &
&                             t(1))
!compute the rotation matrix to update the velocities for the time
!spectral stability derivative case...
! loop over the number of viscous subfaces.
bocoloop2:do mm=1,nviscbocos
! store the rotation center and the rotation rate
! for this subface.
        call pushinteger4(ii)
        ii = cgnssubface(mm)
        rotcenter = cgnsdoms(nbkglobal)%bocoinfo(ii)%rotcenter
        call pushreal8array(rotrate, 3)
        rotrate = timeref*cgnsdoms(nbkglobal)%bocoinfo(ii)%rotrate
! usewindaxis should go back here!
! loop over the quadrilateral faces of the viscous
! subface.
! the new procedure is less elegant as the previous one.
! but the new stands up to tapenade.
        if (bcfaceid(mm) .eq. imin) then
          ad_from0 = bcdata(mm)%jcbeg
          do j=ad_from0,bcdata(mm)%jcend
            ad_from = bcdata(mm)%icbeg
            do i=ad_from,bcdata(mm)%icend
! compute the coordinates of the centroid of the face.
! normally this would be an average of i-1 and i, but
! due to the usage of the pointer xface and the fact
! that x starts at index 0 this is shifted 1 index.
              xc(1) = fourth*(flowdoms(nn, groundlevel, sps)%x(1, i, j, &
&               1)+flowdoms(nn, groundlevel, sps)%x(1, i, j-1, 1)+&
&               flowdoms(nn, groundlevel, sps)%x(1, i-1, j, 1)+flowdoms(&
&               nn, groundlevel, sps)%x(1, i-1, j-1, 1))
              xc(2) = fourth*(flowdoms(nn, groundlevel, sps)%x(1, i, j, &
&               2)+flowdoms(nn, groundlevel, sps)%x(1, i, j-1, 2)+&
&               flowdoms(nn, groundlevel, sps)%x(1, i-1, j, 2)+flowdoms(&
&               nn, groundlevel, sps)%x(1, i-1, j-1, 2))
              xc(3) = fourth*(flowdoms(nn, groundlevel, sps)%x(1, i, j, &
&               3)+flowdoms(nn, groundlevel, sps)%x(1, i, j-1, 3)+&
&               flowdoms(nn, groundlevel, sps)%x(1, i-1, j, 3)+flowdoms(&
&               nn, groundlevel, sps)%x(1, i-1, j-1, 3))
! determine the coordinates relative to the center
! of rotation.
              call pushreal8(xxc(1))
              xxc(1) = xc(1) - rotcenter(1)
              call pushreal8(xxc(2))
              xxc(2) = xc(2) - rotcenter(2)
              call pushreal8(xxc(3))
              xxc(3) = xc(3) - rotcenter(3)
! compute the velocity, which is the cross product
! of rotrate and xc.
! determine the coordinates relative to the
! rigid body rotation point.
              call pushreal8(xxc(1))
              xxc(1) = xc(1) - rotationpoint(1)
              call pushreal8(xxc(2))
              xxc(2) = xc(2) - rotationpoint(2)
              call pushreal8(xxc(3))
              xxc(3) = xc(3) - rotationpoint(3)
! determine the total velocity of the cell center.
! this is a combination of rotation speed of this
! block and the entire rigid body rotation.
            end do
            call pushinteger4(i - 1)
            call pushinteger4(ad_from)
          end do
          call pushinteger4(j - 1)
          call pushinteger4(ad_from0)
          call pushcontrol3b(6)
        else if (bcfaceid(mm) .eq. imax) then
          ad_from2 = bcdata(mm)%jcbeg
          do j=ad_from2,bcdata(mm)%jcend
            ad_from1 = bcdata(mm)%icbeg
            do i=ad_from1,bcdata(mm)%icend
! compute the coordinates of the centroid of the face.
! normally this would be an average of i-1 and i, but
! due to the usage of the pointer xface and the fact
! that x starts at index 0 this is shifted 1 index.
              xc(1) = fourth*(flowdoms(nn, groundlevel, sps)%x(il, i, j&
&               , 1)+flowdoms(nn, groundlevel, sps)%x(il, i, j-1, 1)+&
&               flowdoms(nn, groundlevel, sps)%x(il, i-1, j, 1)+flowdoms&
&               (nn, groundlevel, sps)%x(il, i-1, j-1, 1))
              xc(2) = fourth*(flowdoms(nn, groundlevel, sps)%x(il, i, j&
&               , 2)+flowdoms(nn, groundlevel, sps)%x(il, i, j-1, 2)+&
&               flowdoms(nn, groundlevel, sps)%x(il, i-1, j, 2)+flowdoms&
&               (nn, groundlevel, sps)%x(il, i-1, j-1, 2))
              xc(3) = fourth*(flowdoms(nn, groundlevel, sps)%x(il, i, j&
&               , 3)+flowdoms(nn, groundlevel, sps)%x(il, i, j-1, 3)+&
&               flowdoms(nn, groundlevel, sps)%x(il, i-1, j, 3)+flowdoms&
&               (nn, groundlevel, sps)%x(il, i-1, j-1, 3))
! determine the coordinates relative to the center
! of rotation.
              call pushreal8(xxc(1))
              xxc(1) = xc(1) - rotcenter(1)
              call pushreal8(xxc(2))
              xxc(2) = xc(2) - rotcenter(2)
              call pushreal8(xxc(3))
              xxc(3) = xc(3) - rotcenter(3)
! compute the velocity, which is the cross product
! of rotrate and xc.
! determine the coordinates relative to the
! rigid body rotation point.
              call pushreal8(xxc(1))
              xxc(1) = xc(1) - rotationpoint(1)
              call pushreal8(xxc(2))
              xxc(2) = xc(2) - rotationpoint(2)
              call pushreal8(xxc(3))
              xxc(3) = xc(3) - rotationpoint(3)
! determine the total velocity of the cell center.
! this is a combination of rotation speed of this
! block and the entire rigid body rotation.
            end do
            call pushinteger4(i - 1)
            call pushinteger4(ad_from1)
          end do
          call pushinteger4(j - 1)
          call pushinteger4(ad_from2)
          call pushcontrol3b(5)
        else if (bcfaceid(mm) .eq. jmin) then
          ad_from4 = bcdata(mm)%jcbeg
          do j=ad_from4,bcdata(mm)%jcend
            ad_from3 = bcdata(mm)%icbeg
            do i=ad_from3,bcdata(mm)%icend
! compute the coordinates of the centroid of the face.
! normally this would be an average of i-1 and i, but
! due to the usage of the pointer xface and the fact
! that x starts at index 0 this is shifted 1 index.
              xc(1) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, 1, j, &
&               1)+flowdoms(nn, groundlevel, sps)%x(i, 1, j-1, 1)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, 1, j, 1)+flowdoms(&
&               nn, groundlevel, sps)%x(i-1, 1, j-1, 1))
              xc(2) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, 1, j, &
&               2)+flowdoms(nn, groundlevel, sps)%x(i, 1, j-1, 2)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, 1, j, 2)+flowdoms(&
&               nn, groundlevel, sps)%x(i-1, 1, j-1, 2))
              xc(3) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, 1, j, &
&               3)+flowdoms(nn, groundlevel, sps)%x(i, 1, j-1, 3)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, 1, j, 3)+flowdoms(&
&               nn, groundlevel, sps)%x(i-1, 1, j-1, 3))
! determine the coordinates relative to the center
! of rotation.
              call pushreal8(xxc(1))
              xxc(1) = xc(1) - rotcenter(1)
              call pushreal8(xxc(2))
              xxc(2) = xc(2) - rotcenter(2)
              call pushreal8(xxc(3))
              xxc(3) = xc(3) - rotcenter(3)
! compute the velocity, which is the cross product
! of rotrate and xc.
! determine the coordinates relative to the
! rigid body rotation point.
              call pushreal8(xxc(1))
              xxc(1) = xc(1) - rotationpoint(1)
              call pushreal8(xxc(2))
              xxc(2) = xc(2) - rotationpoint(2)
              call pushreal8(xxc(3))
              xxc(3) = xc(3) - rotationpoint(3)
! determine the total velocity of the cell center.
! this is a combination of rotation speed of this
! block and the entire rigid body rotation.
            end do
            call pushinteger4(i - 1)
            call pushinteger4(ad_from3)
          end do
          call pushinteger4(j - 1)
          call pushinteger4(ad_from4)
          call pushcontrol3b(4)
        else if (bcfaceid(mm) .eq. jmax) then
          ad_from6 = bcdata(mm)%jcbeg
          do j=ad_from6,bcdata(mm)%jcend
            ad_from5 = bcdata(mm)%icbeg
            do i=ad_from5,bcdata(mm)%icend
! compute the coordinates of the centroid of the face.
! normally this would be an average of i-1 and i, but
! due to the usage of the pointer xface and the fact
! that x starts at index 0 this is shifted 1 index.
              xc(1) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, jl, j&
&               , 1)+flowdoms(nn, groundlevel, sps)%x(i, jl, j-1, 1)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, jl, j, 1)+flowdoms&
&               (nn, groundlevel, sps)%x(i-1, jl, j-1, 1))
              xc(2) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, jl, j&
&               , 2)+flowdoms(nn, groundlevel, sps)%x(i, jl, j-1, 2)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, jl, j, 2)+flowdoms&
&               (nn, groundlevel, sps)%x(i-1, jl, j-1, 2))
              xc(3) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, jl, j&
&               , 3)+flowdoms(nn, groundlevel, sps)%x(i, jl, j-1, 3)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, jl, j, 3)+flowdoms&
&               (nn, groundlevel, sps)%x(i-1, jl, j-1, 3))
! determine the coordinates relative to the center
! of rotation.
              call pushreal8(xxc(1))
              xxc(1) = xc(1) - rotcenter(1)
              call pushreal8(xxc(2))
              xxc(2) = xc(2) - rotcenter(2)
              call pushreal8(xxc(3))
              xxc(3) = xc(3) - rotcenter(3)
! compute the velocity, which is the cross product
! of rotrate and xc.
! determine the coordinates relative to the
! rigid body rotation point.
              call pushreal8(xxc(1))
              xxc(1) = xc(1) - rotationpoint(1)
              call pushreal8(xxc(2))
              xxc(2) = xc(2) - rotationpoint(2)
              call pushreal8(xxc(3))
              xxc(3) = xc(3) - rotationpoint(3)
! determine the total velocity of the cell center.
! this is a combination of rotation speed of this
! block and the entire rigid body rotation.
            end do
            call pushinteger4(i - 1)
            call pushinteger4(ad_from5)
          end do
          call pushinteger4(j - 1)
          call pushinteger4(ad_from6)
          call pushcontrol3b(3)
        else if (bcfaceid(mm) .eq. kmin) then
          ad_from8 = bcdata(mm)%jcbeg
          do j=ad_from8,bcdata(mm)%jcend
            ad_from7 = bcdata(mm)%icbeg
            do i=ad_from7,bcdata(mm)%icend
! compute the coordinates of the centroid of the face.
! normally this would be an average of i-1 and i, but
! due to the usage of the pointer xface and the fact
! that x starts at index 0 this is shifted 1 index.
              xc(1) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, j, 1, &
&               1)+flowdoms(nn, groundlevel, sps)%x(i, j-1, 1, 1)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, j, 1, 1)+flowdoms(&
&               nn, groundlevel, sps)%x(i-1, j-1, 1, 1))
              xc(2) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, j, 1, &
&               2)+flowdoms(nn, groundlevel, sps)%x(i, j-1, 1, 2)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, j, 1, 2)+flowdoms(&
&               nn, groundlevel, sps)%x(i-1, j-1, 1, 2))
              xc(3) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, j, 1, &
&               3)+flowdoms(nn, groundlevel, sps)%x(i, j-1, 1, 3)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, j, 1, 3)+flowdoms(&
&               nn, groundlevel, sps)%x(i-1, j-1, 1, 3))
! determine the coordinates relative to the center
! of rotation.
              call pushreal8(xxc(1))
              xxc(1) = xc(1) - rotcenter(1)
              call pushreal8(xxc(2))
              xxc(2) = xc(2) - rotcenter(2)
              call pushreal8(xxc(3))
              xxc(3) = xc(3) - rotcenter(3)
! compute the velocity, which is the cross product
! of rotrate and xc.
! determine the coordinates relative to the
! rigid body rotation point.
              call pushreal8(xxc(1))
              xxc(1) = xc(1) - rotationpoint(1)
              call pushreal8(xxc(2))
              xxc(2) = xc(2) - rotationpoint(2)
              call pushreal8(xxc(3))
              xxc(3) = xc(3) - rotationpoint(3)
! determine the total velocity of the cell center.
! this is a combination of rotation speed of this
! block and the entire rigid body rotation.
            end do
            call pushinteger4(i - 1)
            call pushinteger4(ad_from7)
          end do
          call pushinteger4(j - 1)
          call pushinteger4(ad_from8)
          call pushcontrol3b(2)
        else if (bcfaceid(mm) .eq. kmax) then
          ad_from10 = bcdata(mm)%jcbeg
          do j=ad_from10,bcdata(mm)%jcend
            ad_from9 = bcdata(mm)%icbeg
            do i=ad_from9,bcdata(mm)%icend
! compute the coordinates of the centroid of the face.
! normally this would be an average of i-1 and i, but
! due to the usage of the pointer xface and the fact
! that x starts at index 0 this is shifted 1 index.
              xc(1) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, j, kl&
&               , 1)+flowdoms(nn, groundlevel, sps)%x(i, j-1, kl, 1)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, j, kl, 1)+flowdoms&
&               (nn, groundlevel, sps)%x(i-1, j-1, kl, 1))
              xc(2) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, j, kl&
&               , 2)+flowdoms(nn, groundlevel, sps)%x(i, j-1, kl, 2)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, j, kl, 2)+flowdoms&
&               (nn, groundlevel, sps)%x(i-1, j-1, kl, 2))
              xc(3) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, j, kl&
&               , 3)+flowdoms(nn, groundlevel, sps)%x(i, j-1, kl, 3)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, j, kl, 3)+flowdoms&
&               (nn, groundlevel, sps)%x(i-1, j-1, kl, 3))
! determine the coordinates relative to the center
! of rotation.
              call pushreal8(xxc(1))
              xxc(1) = xc(1) - rotcenter(1)
              call pushreal8(xxc(2))
              xxc(2) = xc(2) - rotcenter(2)
              call pushreal8(xxc(3))
              xxc(3) = xc(3) - rotcenter(3)
! compute the velocity, which is the cross product
! of rotrate and xc.
! determine the coordinates relative to the
! rigid body rotation point.
              call pushreal8(xxc(1))
              xxc(1) = xc(1) - rotationpoint(1)
              call pushreal8(xxc(2))
              xxc(2) = xc(2) - rotationpoint(2)
              call pushreal8(xxc(3))
              xxc(3) = xc(3) - rotationpoint(3)
! determine the total velocity of the cell center.
! this is a combination of rotation speed of this
! block and the entire rigid body rotation.
            end do
            call pushinteger4(i - 1)
            call pushinteger4(ad_from9)
          end do
          call pushinteger4(j - 1)
          call pushinteger4(ad_from10)
          call pushcontrol3b(1)
        else
          call pushcontrol3b(0)
        end if
      end do bocoloop2
      timerefd = 0.0_8
      velxgrid0d = 0.0_8
      velzgrid0d = 0.0_8
      xcd = 0.0_8
      xxcd = 0.0_8
      derivrotationmatrixd = 0.0_8
      velygrid0d = 0.0_8
      do mm=nviscbocos,1,-1
        call popcontrol3b(branch)
        if (branch .lt. 3) then
          if (branch .eq. 0) then
            rotrated = 0.0_8
            velygridd = 0.0_8
            velzgridd = 0.0_8
            velxgridd = 0.0_8
          else if (branch .eq. 1) then
            rotrated = 0.0_8
            velygridd = 0.0_8
            velzgridd = 0.0_8
            velxgridd = 0.0_8
            call popinteger4(ad_from10)
            call popinteger4(ad_to10)
            do j=ad_to10,ad_from10,-1
              call popinteger4(ad_from9)
              call popinteger4(ad_to9)
              do i=ad_to9,ad_from9,-1
                velzgridd = velzgridd + bcdatad(mm)%uslip(i, j, 3)
                derivrotationmatrixd(3, 1) = derivrotationmatrixd(3, 1) &
&                 + xxc(1)*bcdatad(mm)%uslip(i, j, 3)
                xxcd(1) = xxcd(1) + derivrotationmatrix(3, 1)*bcdatad(mm&
&                 )%uslip(i, j, 3) + derivrotationmatrix(2, 1)*bcdatad(&
&                 mm)%uslip(i, j, 2) + derivrotationmatrix(1, 1)*bcdatad&
&                 (mm)%uslip(i, j, 1)
                derivrotationmatrixd(3, 2) = derivrotationmatrixd(3, 2) &
&                 + xxc(2)*bcdatad(mm)%uslip(i, j, 3)
                xxcd(2) = xxcd(2) + derivrotationmatrix(3, 2)*bcdatad(mm&
&                 )%uslip(i, j, 3) + derivrotationmatrix(2, 2)*bcdatad(&
&                 mm)%uslip(i, j, 2) + derivrotationmatrix(1, 2)*bcdatad&
&                 (mm)%uslip(i, j, 1)
                derivrotationmatrixd(3, 3) = derivrotationmatrixd(3, 3) &
&                 + xxc(3)*bcdatad(mm)%uslip(i, j, 3)
                xxcd(3) = xxcd(3) + derivrotationmatrix(3, 3)*bcdatad(mm&
&                 )%uslip(i, j, 3) + derivrotationmatrix(2, 3)*bcdatad(&
&                 mm)%uslip(i, j, 2) + derivrotationmatrix(1, 3)*bcdatad&
&                 (mm)%uslip(i, j, 1)
                velygridd = velygridd + bcdatad(mm)%uslip(i, j, 2)
                derivrotationmatrixd(2, 1) = derivrotationmatrixd(2, 1) &
&                 + xxc(1)*bcdatad(mm)%uslip(i, j, 2)
                derivrotationmatrixd(2, 2) = derivrotationmatrixd(2, 2) &
&                 + xxc(2)*bcdatad(mm)%uslip(i, j, 2)
                derivrotationmatrixd(2, 3) = derivrotationmatrixd(2, 3) &
&                 + xxc(3)*bcdatad(mm)%uslip(i, j, 2)
                velxgridd = velxgridd + bcdatad(mm)%uslip(i, j, 1)
                derivrotationmatrixd(1, 1) = derivrotationmatrixd(1, 1) &
&                 + xxc(1)*bcdatad(mm)%uslip(i, j, 1)
                derivrotationmatrixd(1, 2) = derivrotationmatrixd(1, 2) &
&                 + xxc(2)*bcdatad(mm)%uslip(i, j, 1)
                derivrotationmatrixd(1, 3) = derivrotationmatrixd(1, 3) &
&                 + xxc(3)*bcdatad(mm)%uslip(i, j, 1)
                call popreal8(xxc(3))
                xcd(3) = xcd(3) + xxcd(3)
                call popreal8(xxc(2))
                xcd(2) = xcd(2) + xxcd(2)
                xxcd(2) = rotrate(1)*bcdatad(mm)%uslip(i, j, 3)
                call popreal8(xxc(1))
                xcd(1) = xcd(1) + xxcd(1)
                xxcd(1) = -(rotrate(2)*bcdatad(mm)%uslip(i, j, 3))
                rotrated(1) = rotrated(1) + xxc(2)*bcdatad(mm)%uslip(i, &
&                 j, 3)
                rotrated(2) = rotrated(2) - xxc(1)*bcdatad(mm)%uslip(i, &
&                 j, 3)
                bcdatad(mm)%uslip(i, j, 3) = 0.0_8
                xxcd(3) = -(rotrate(1)*bcdatad(mm)%uslip(i, j, 2))
                rotrated(3) = rotrated(3) + xxc(1)*bcdatad(mm)%uslip(i, &
&                 j, 2)
                xxcd(1) = xxcd(1) + rotrate(3)*bcdatad(mm)%uslip(i, j, 2&
&                 )
                rotrated(1) = rotrated(1) - xxc(3)*bcdatad(mm)%uslip(i, &
&                 j, 2)
                bcdatad(mm)%uslip(i, j, 2) = 0.0_8
                rotrated(2) = rotrated(2) + xxc(3)*bcdatad(mm)%uslip(i, &
&                 j, 1)
                xxcd(3) = xxcd(3) + rotrate(2)*bcdatad(mm)%uslip(i, j, 1&
&                 )
                rotrated(3) = rotrated(3) - xxc(2)*bcdatad(mm)%uslip(i, &
&                 j, 1)
                xxcd(2) = xxcd(2) - rotrate(3)*bcdatad(mm)%uslip(i, j, 1&
&                 )
                bcdatad(mm)%uslip(i, j, 1) = 0.0_8
                call popreal8(xxc(3))
                xcd(3) = xcd(3) + xxcd(3)
                xxcd(3) = 0.0_8
                call popreal8(xxc(2))
                xcd(2) = xcd(2) + xxcd(2)
                xxcd(2) = 0.0_8
                call popreal8(xxc(1))
                xcd(1) = xcd(1) + xxcd(1)
                xxcd(1) = 0.0_8
                tempd = fourth*xcd(3)
                xcd(3) = 0.0_8
                flowdomsd(nn, groundlevel, sps)%x(i, j, kl, 3) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i, j, kl, 3) + tempd
                flowdomsd(nn, groundlevel, sps)%x(i, j-1, kl, 3) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i, j-1, kl, 3) + &
&                 tempd
                flowdomsd(nn, groundlevel, sps)%x(i-1, j, kl, 3) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i-1, j, kl, 3) + &
&                 tempd
                flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, kl, 3) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, kl, 3) + &
&                 tempd
                tempd = fourth*xcd(2)
                xcd(2) = 0.0_8
                flowdomsd(nn, groundlevel, sps)%x(i, j, kl, 2) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i, j, kl, 2) + tempd
                flowdomsd(nn, groundlevel, sps)%x(i, j-1, kl, 2) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i, j-1, kl, 2) + &
&                 tempd
                flowdomsd(nn, groundlevel, sps)%x(i-1, j, kl, 2) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i-1, j, kl, 2) + &
&                 tempd
                flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, kl, 2) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, kl, 2) + &
&                 tempd
                tempd = fourth*xcd(1)
                xcd(1) = 0.0_8
                flowdomsd(nn, groundlevel, sps)%x(i, j, kl, 1) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i, j, kl, 1) + tempd
                flowdomsd(nn, groundlevel, sps)%x(i, j-1, kl, 1) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i, j-1, kl, 1) + &
&                 tempd
                flowdomsd(nn, groundlevel, sps)%x(i-1, j, kl, 1) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i-1, j, kl, 1) + &
&                 tempd
                flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, kl, 1) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, kl, 1) + &
&                 tempd
              end do
            end do
          else
            rotrated = 0.0_8
            velygridd = 0.0_8
            velzgridd = 0.0_8
            velxgridd = 0.0_8
            call popinteger4(ad_from8)
            call popinteger4(ad_to8)
            do j=ad_to8,ad_from8,-1
              call popinteger4(ad_from7)
              call popinteger4(ad_to7)
              do i=ad_to7,ad_from7,-1
                velzgridd = velzgridd + bcdatad(mm)%uslip(i, j, 3)
                derivrotationmatrixd(3, 1) = derivrotationmatrixd(3, 1) &
&                 + xxc(1)*bcdatad(mm)%uslip(i, j, 3)
                xxcd(1) = xxcd(1) + derivrotationmatrix(3, 1)*bcdatad(mm&
&                 )%uslip(i, j, 3) + derivrotationmatrix(2, 1)*bcdatad(&
&                 mm)%uslip(i, j, 2) + derivrotationmatrix(1, 1)*bcdatad&
&                 (mm)%uslip(i, j, 1)
                derivrotationmatrixd(3, 2) = derivrotationmatrixd(3, 2) &
&                 + xxc(2)*bcdatad(mm)%uslip(i, j, 3)
                xxcd(2) = xxcd(2) + derivrotationmatrix(3, 2)*bcdatad(mm&
&                 )%uslip(i, j, 3) + derivrotationmatrix(2, 2)*bcdatad(&
&                 mm)%uslip(i, j, 2) + derivrotationmatrix(1, 2)*bcdatad&
&                 (mm)%uslip(i, j, 1)
                derivrotationmatrixd(3, 3) = derivrotationmatrixd(3, 3) &
&                 + xxc(3)*bcdatad(mm)%uslip(i, j, 3)
                xxcd(3) = xxcd(3) + derivrotationmatrix(3, 3)*bcdatad(mm&
&                 )%uslip(i, j, 3) + derivrotationmatrix(2, 3)*bcdatad(&
&                 mm)%uslip(i, j, 2) + derivrotationmatrix(1, 3)*bcdatad&
&                 (mm)%uslip(i, j, 1)
                velygridd = velygridd + bcdatad(mm)%uslip(i, j, 2)
                derivrotationmatrixd(2, 1) = derivrotationmatrixd(2, 1) &
&                 + xxc(1)*bcdatad(mm)%uslip(i, j, 2)
                derivrotationmatrixd(2, 2) = derivrotationmatrixd(2, 2) &
&                 + xxc(2)*bcdatad(mm)%uslip(i, j, 2)
                derivrotationmatrixd(2, 3) = derivrotationmatrixd(2, 3) &
&                 + xxc(3)*bcdatad(mm)%uslip(i, j, 2)
                velxgridd = velxgridd + bcdatad(mm)%uslip(i, j, 1)
                derivrotationmatrixd(1, 1) = derivrotationmatrixd(1, 1) &
&                 + xxc(1)*bcdatad(mm)%uslip(i, j, 1)
                derivrotationmatrixd(1, 2) = derivrotationmatrixd(1, 2) &
&                 + xxc(2)*bcdatad(mm)%uslip(i, j, 1)
                derivrotationmatrixd(1, 3) = derivrotationmatrixd(1, 3) &
&                 + xxc(3)*bcdatad(mm)%uslip(i, j, 1)
                call popreal8(xxc(3))
                xcd(3) = xcd(3) + xxcd(3)
                call popreal8(xxc(2))
                xcd(2) = xcd(2) + xxcd(2)
                xxcd(2) = rotrate(1)*bcdatad(mm)%uslip(i, j, 3)
                call popreal8(xxc(1))
                xcd(1) = xcd(1) + xxcd(1)
                xxcd(1) = -(rotrate(2)*bcdatad(mm)%uslip(i, j, 3))
                rotrated(1) = rotrated(1) + xxc(2)*bcdatad(mm)%uslip(i, &
&                 j, 3)
                rotrated(2) = rotrated(2) - xxc(1)*bcdatad(mm)%uslip(i, &
&                 j, 3)
                bcdatad(mm)%uslip(i, j, 3) = 0.0_8
                xxcd(3) = -(rotrate(1)*bcdatad(mm)%uslip(i, j, 2))
                rotrated(3) = rotrated(3) + xxc(1)*bcdatad(mm)%uslip(i, &
&                 j, 2)
                xxcd(1) = xxcd(1) + rotrate(3)*bcdatad(mm)%uslip(i, j, 2&
&                 )
                rotrated(1) = rotrated(1) - xxc(3)*bcdatad(mm)%uslip(i, &
&                 j, 2)
                bcdatad(mm)%uslip(i, j, 2) = 0.0_8
                rotrated(2) = rotrated(2) + xxc(3)*bcdatad(mm)%uslip(i, &
&                 j, 1)
                xxcd(3) = xxcd(3) + rotrate(2)*bcdatad(mm)%uslip(i, j, 1&
&                 )
                rotrated(3) = rotrated(3) - xxc(2)*bcdatad(mm)%uslip(i, &
&                 j, 1)
                xxcd(2) = xxcd(2) - rotrate(3)*bcdatad(mm)%uslip(i, j, 1&
&                 )
                bcdatad(mm)%uslip(i, j, 1) = 0.0_8
                call popreal8(xxc(3))
                xcd(3) = xcd(3) + xxcd(3)
                xxcd(3) = 0.0_8
                call popreal8(xxc(2))
                xcd(2) = xcd(2) + xxcd(2)
                xxcd(2) = 0.0_8
                call popreal8(xxc(1))
                xcd(1) = xcd(1) + xxcd(1)
                xxcd(1) = 0.0_8
                tempd = fourth*xcd(3)
                xcd(3) = 0.0_8
                flowdomsd(nn, groundlevel, sps)%x(i, j, 1, 3) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i, j, 1, 3) + tempd
                flowdomsd(nn, groundlevel, sps)%x(i, j-1, 1, 3) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i, j-1, 1, 3) + &
&                 tempd
                flowdomsd(nn, groundlevel, sps)%x(i-1, j, 1, 3) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i-1, j, 1, 3) + &
&                 tempd
                flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, 1, 3) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, 1, 3) + &
&                 tempd
                tempd = fourth*xcd(2)
                xcd(2) = 0.0_8
                flowdomsd(nn, groundlevel, sps)%x(i, j, 1, 2) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i, j, 1, 2) + tempd
                flowdomsd(nn, groundlevel, sps)%x(i, j-1, 1, 2) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i, j-1, 1, 2) + &
&                 tempd
                flowdomsd(nn, groundlevel, sps)%x(i-1, j, 1, 2) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i-1, j, 1, 2) + &
&                 tempd
                flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, 1, 2) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, 1, 2) + &
&                 tempd
                tempd = fourth*xcd(1)
                xcd(1) = 0.0_8
                flowdomsd(nn, groundlevel, sps)%x(i, j, 1, 1) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i, j, 1, 1) + tempd
                flowdomsd(nn, groundlevel, sps)%x(i, j-1, 1, 1) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i, j-1, 1, 1) + &
&                 tempd
                flowdomsd(nn, groundlevel, sps)%x(i-1, j, 1, 1) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i-1, j, 1, 1) + &
&                 tempd
                flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, 1, 1) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, 1, 1) + &
&                 tempd
              end do
            end do
          end if
        else if (branch .lt. 5) then
          if (branch .eq. 3) then
            rotrated = 0.0_8
            velygridd = 0.0_8
            velzgridd = 0.0_8
            velxgridd = 0.0_8
            call popinteger4(ad_from6)
            call popinteger4(ad_to6)
            do j=ad_to6,ad_from6,-1
              call popinteger4(ad_from5)
              call popinteger4(ad_to5)
              do i=ad_to5,ad_from5,-1
                velzgridd = velzgridd + bcdatad(mm)%uslip(i, j, 3)
                derivrotationmatrixd(3, 1) = derivrotationmatrixd(3, 1) &
&                 + xxc(1)*bcdatad(mm)%uslip(i, j, 3)
                xxcd(1) = xxcd(1) + derivrotationmatrix(3, 1)*bcdatad(mm&
&                 )%uslip(i, j, 3) + derivrotationmatrix(2, 1)*bcdatad(&
&                 mm)%uslip(i, j, 2) + derivrotationmatrix(1, 1)*bcdatad&
&                 (mm)%uslip(i, j, 1)
                derivrotationmatrixd(3, 2) = derivrotationmatrixd(3, 2) &
&                 + xxc(2)*bcdatad(mm)%uslip(i, j, 3)
                xxcd(2) = xxcd(2) + derivrotationmatrix(3, 2)*bcdatad(mm&
&                 )%uslip(i, j, 3) + derivrotationmatrix(2, 2)*bcdatad(&
&                 mm)%uslip(i, j, 2) + derivrotationmatrix(1, 2)*bcdatad&
&                 (mm)%uslip(i, j, 1)
                derivrotationmatrixd(3, 3) = derivrotationmatrixd(3, 3) &
&                 + xxc(3)*bcdatad(mm)%uslip(i, j, 3)
                xxcd(3) = xxcd(3) + derivrotationmatrix(3, 3)*bcdatad(mm&
&                 )%uslip(i, j, 3) + derivrotationmatrix(2, 3)*bcdatad(&
&                 mm)%uslip(i, j, 2) + derivrotationmatrix(1, 3)*bcdatad&
&                 (mm)%uslip(i, j, 1)
                velygridd = velygridd + bcdatad(mm)%uslip(i, j, 2)
                derivrotationmatrixd(2, 1) = derivrotationmatrixd(2, 1) &
&                 + xxc(1)*bcdatad(mm)%uslip(i, j, 2)
                derivrotationmatrixd(2, 2) = derivrotationmatrixd(2, 2) &
&                 + xxc(2)*bcdatad(mm)%uslip(i, j, 2)
                derivrotationmatrixd(2, 3) = derivrotationmatrixd(2, 3) &
&                 + xxc(3)*bcdatad(mm)%uslip(i, j, 2)
                velxgridd = velxgridd + bcdatad(mm)%uslip(i, j, 1)
                derivrotationmatrixd(1, 1) = derivrotationmatrixd(1, 1) &
&                 + xxc(1)*bcdatad(mm)%uslip(i, j, 1)
                derivrotationmatrixd(1, 2) = derivrotationmatrixd(1, 2) &
&                 + xxc(2)*bcdatad(mm)%uslip(i, j, 1)
                derivrotationmatrixd(1, 3) = derivrotationmatrixd(1, 3) &
&                 + xxc(3)*bcdatad(mm)%uslip(i, j, 1)
                call popreal8(xxc(3))
                xcd(3) = xcd(3) + xxcd(3)
                call popreal8(xxc(2))
                xcd(2) = xcd(2) + xxcd(2)
                xxcd(2) = rotrate(1)*bcdatad(mm)%uslip(i, j, 3)
                call popreal8(xxc(1))
                xcd(1) = xcd(1) + xxcd(1)
                xxcd(1) = -(rotrate(2)*bcdatad(mm)%uslip(i, j, 3))
                rotrated(1) = rotrated(1) + xxc(2)*bcdatad(mm)%uslip(i, &
&                 j, 3)
                rotrated(2) = rotrated(2) - xxc(1)*bcdatad(mm)%uslip(i, &
&                 j, 3)
                bcdatad(mm)%uslip(i, j, 3) = 0.0_8
                xxcd(3) = -(rotrate(1)*bcdatad(mm)%uslip(i, j, 2))
                rotrated(3) = rotrated(3) + xxc(1)*bcdatad(mm)%uslip(i, &
&                 j, 2)
                xxcd(1) = xxcd(1) + rotrate(3)*bcdatad(mm)%uslip(i, j, 2&
&                 )
                rotrated(1) = rotrated(1) - xxc(3)*bcdatad(mm)%uslip(i, &
&                 j, 2)
                bcdatad(mm)%uslip(i, j, 2) = 0.0_8
                rotrated(2) = rotrated(2) + xxc(3)*bcdatad(mm)%uslip(i, &
&                 j, 1)
                xxcd(3) = xxcd(3) + rotrate(2)*bcdatad(mm)%uslip(i, j, 1&
&                 )
                rotrated(3) = rotrated(3) - xxc(2)*bcdatad(mm)%uslip(i, &
&                 j, 1)
                xxcd(2) = xxcd(2) - rotrate(3)*bcdatad(mm)%uslip(i, j, 1&
&                 )
                bcdatad(mm)%uslip(i, j, 1) = 0.0_8
                call popreal8(xxc(3))
                xcd(3) = xcd(3) + xxcd(3)
                xxcd(3) = 0.0_8
                call popreal8(xxc(2))
                xcd(2) = xcd(2) + xxcd(2)
                xxcd(2) = 0.0_8
                call popreal8(xxc(1))
                xcd(1) = xcd(1) + xxcd(1)
                xxcd(1) = 0.0_8
                tempd = fourth*xcd(3)
                xcd(3) = 0.0_8
                flowdomsd(nn, groundlevel, sps)%x(i, jl, j, 3) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i, jl, j, 3) + tempd
                flowdomsd(nn, groundlevel, sps)%x(i, jl, j-1, 3) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i, jl, j-1, 3) + &
&                 tempd
                flowdomsd(nn, groundlevel, sps)%x(i-1, jl, j, 3) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i-1, jl, j, 3) + &
&                 tempd
                flowdomsd(nn, groundlevel, sps)%x(i-1, jl, j-1, 3) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i-1, jl, j-1, 3) + &
&                 tempd
                tempd = fourth*xcd(2)
                xcd(2) = 0.0_8
                flowdomsd(nn, groundlevel, sps)%x(i, jl, j, 2) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i, jl, j, 2) + tempd
                flowdomsd(nn, groundlevel, sps)%x(i, jl, j-1, 2) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i, jl, j-1, 2) + &
&                 tempd
                flowdomsd(nn, groundlevel, sps)%x(i-1, jl, j, 2) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i-1, jl, j, 2) + &
&                 tempd
                flowdomsd(nn, groundlevel, sps)%x(i-1, jl, j-1, 2) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i-1, jl, j-1, 2) + &
&                 tempd
                tempd = fourth*xcd(1)
                xcd(1) = 0.0_8
                flowdomsd(nn, groundlevel, sps)%x(i, jl, j, 1) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i, jl, j, 1) + tempd
                flowdomsd(nn, groundlevel, sps)%x(i, jl, j-1, 1) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i, jl, j-1, 1) + &
&                 tempd
                flowdomsd(nn, groundlevel, sps)%x(i-1, jl, j, 1) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i-1, jl, j, 1) + &
&                 tempd
                flowdomsd(nn, groundlevel, sps)%x(i-1, jl, j-1, 1) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i-1, jl, j-1, 1) + &
&                 tempd
              end do
            end do
          else
            rotrated = 0.0_8
            velygridd = 0.0_8
            velzgridd = 0.0_8
            velxgridd = 0.0_8
            call popinteger4(ad_from4)
            call popinteger4(ad_to4)
            do j=ad_to4,ad_from4,-1
              call popinteger4(ad_from3)
              call popinteger4(ad_to3)
              do i=ad_to3,ad_from3,-1
                velzgridd = velzgridd + bcdatad(mm)%uslip(i, j, 3)
                derivrotationmatrixd(3, 1) = derivrotationmatrixd(3, 1) &
&                 + xxc(1)*bcdatad(mm)%uslip(i, j, 3)
                xxcd(1) = xxcd(1) + derivrotationmatrix(3, 1)*bcdatad(mm&
&                 )%uslip(i, j, 3) + derivrotationmatrix(2, 1)*bcdatad(&
&                 mm)%uslip(i, j, 2) + derivrotationmatrix(1, 1)*bcdatad&
&                 (mm)%uslip(i, j, 1)
                derivrotationmatrixd(3, 2) = derivrotationmatrixd(3, 2) &
&                 + xxc(2)*bcdatad(mm)%uslip(i, j, 3)
                xxcd(2) = xxcd(2) + derivrotationmatrix(3, 2)*bcdatad(mm&
&                 )%uslip(i, j, 3) + derivrotationmatrix(2, 2)*bcdatad(&
&                 mm)%uslip(i, j, 2) + derivrotationmatrix(1, 2)*bcdatad&
&                 (mm)%uslip(i, j, 1)
                derivrotationmatrixd(3, 3) = derivrotationmatrixd(3, 3) &
&                 + xxc(3)*bcdatad(mm)%uslip(i, j, 3)
                xxcd(3) = xxcd(3) + derivrotationmatrix(3, 3)*bcdatad(mm&
&                 )%uslip(i, j, 3) + derivrotationmatrix(2, 3)*bcdatad(&
&                 mm)%uslip(i, j, 2) + derivrotationmatrix(1, 3)*bcdatad&
&                 (mm)%uslip(i, j, 1)
                velygridd = velygridd + bcdatad(mm)%uslip(i, j, 2)
                derivrotationmatrixd(2, 1) = derivrotationmatrixd(2, 1) &
&                 + xxc(1)*bcdatad(mm)%uslip(i, j, 2)
                derivrotationmatrixd(2, 2) = derivrotationmatrixd(2, 2) &
&                 + xxc(2)*bcdatad(mm)%uslip(i, j, 2)
                derivrotationmatrixd(2, 3) = derivrotationmatrixd(2, 3) &
&                 + xxc(3)*bcdatad(mm)%uslip(i, j, 2)
                velxgridd = velxgridd + bcdatad(mm)%uslip(i, j, 1)
                derivrotationmatrixd(1, 1) = derivrotationmatrixd(1, 1) &
&                 + xxc(1)*bcdatad(mm)%uslip(i, j, 1)
                derivrotationmatrixd(1, 2) = derivrotationmatrixd(1, 2) &
&                 + xxc(2)*bcdatad(mm)%uslip(i, j, 1)
                derivrotationmatrixd(1, 3) = derivrotationmatrixd(1, 3) &
&                 + xxc(3)*bcdatad(mm)%uslip(i, j, 1)
                call popreal8(xxc(3))
                xcd(3) = xcd(3) + xxcd(3)
                call popreal8(xxc(2))
                xcd(2) = xcd(2) + xxcd(2)
                xxcd(2) = rotrate(1)*bcdatad(mm)%uslip(i, j, 3)
                call popreal8(xxc(1))
                xcd(1) = xcd(1) + xxcd(1)
                xxcd(1) = -(rotrate(2)*bcdatad(mm)%uslip(i, j, 3))
                rotrated(1) = rotrated(1) + xxc(2)*bcdatad(mm)%uslip(i, &
&                 j, 3)
                rotrated(2) = rotrated(2) - xxc(1)*bcdatad(mm)%uslip(i, &
&                 j, 3)
                bcdatad(mm)%uslip(i, j, 3) = 0.0_8
                xxcd(3) = -(rotrate(1)*bcdatad(mm)%uslip(i, j, 2))
                rotrated(3) = rotrated(3) + xxc(1)*bcdatad(mm)%uslip(i, &
&                 j, 2)
                xxcd(1) = xxcd(1) + rotrate(3)*bcdatad(mm)%uslip(i, j, 2&
&                 )
                rotrated(1) = rotrated(1) - xxc(3)*bcdatad(mm)%uslip(i, &
&                 j, 2)
                bcdatad(mm)%uslip(i, j, 2) = 0.0_8
                rotrated(2) = rotrated(2) + xxc(3)*bcdatad(mm)%uslip(i, &
&                 j, 1)
                xxcd(3) = xxcd(3) + rotrate(2)*bcdatad(mm)%uslip(i, j, 1&
&                 )
                rotrated(3) = rotrated(3) - xxc(2)*bcdatad(mm)%uslip(i, &
&                 j, 1)
                xxcd(2) = xxcd(2) - rotrate(3)*bcdatad(mm)%uslip(i, j, 1&
&                 )
                bcdatad(mm)%uslip(i, j, 1) = 0.0_8
                call popreal8(xxc(3))
                xcd(3) = xcd(3) + xxcd(3)
                xxcd(3) = 0.0_8
                call popreal8(xxc(2))
                xcd(2) = xcd(2) + xxcd(2)
                xxcd(2) = 0.0_8
                call popreal8(xxc(1))
                xcd(1) = xcd(1) + xxcd(1)
                xxcd(1) = 0.0_8
                tempd = fourth*xcd(3)
                xcd(3) = 0.0_8
                flowdomsd(nn, groundlevel, sps)%x(i, 1, j, 3) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i, 1, j, 3) + tempd
                flowdomsd(nn, groundlevel, sps)%x(i, 1, j-1, 3) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i, 1, j-1, 3) + &
&                 tempd
                flowdomsd(nn, groundlevel, sps)%x(i-1, 1, j, 3) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i-1, 1, j, 3) + &
&                 tempd
                flowdomsd(nn, groundlevel, sps)%x(i-1, 1, j-1, 3) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i-1, 1, j-1, 3) + &
&                 tempd
                tempd = fourth*xcd(2)
                xcd(2) = 0.0_8
                flowdomsd(nn, groundlevel, sps)%x(i, 1, j, 2) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i, 1, j, 2) + tempd
                flowdomsd(nn, groundlevel, sps)%x(i, 1, j-1, 2) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i, 1, j-1, 2) + &
&                 tempd
                flowdomsd(nn, groundlevel, sps)%x(i-1, 1, j, 2) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i-1, 1, j, 2) + &
&                 tempd
                flowdomsd(nn, groundlevel, sps)%x(i-1, 1, j-1, 2) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i-1, 1, j-1, 2) + &
&                 tempd
                tempd = fourth*xcd(1)
                xcd(1) = 0.0_8
                flowdomsd(nn, groundlevel, sps)%x(i, 1, j, 1) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i, 1, j, 1) + tempd
                flowdomsd(nn, groundlevel, sps)%x(i, 1, j-1, 1) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i, 1, j-1, 1) + &
&                 tempd
                flowdomsd(nn, groundlevel, sps)%x(i-1, 1, j, 1) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i-1, 1, j, 1) + &
&                 tempd
                flowdomsd(nn, groundlevel, sps)%x(i-1, 1, j-1, 1) = &
&                 flowdomsd(nn, groundlevel, sps)%x(i-1, 1, j-1, 1) + &
&                 tempd
              end do
            end do
          end if
        else if (branch .eq. 5) then
          rotrated = 0.0_8
          velygridd = 0.0_8
          velzgridd = 0.0_8
          velxgridd = 0.0_8
          call popinteger4(ad_from2)
          call popinteger4(ad_to2)
          do j=ad_to2,ad_from2,-1
            call popinteger4(ad_from1)
            call popinteger4(ad_to1)
            do i=ad_to1,ad_from1,-1
              velzgridd = velzgridd + bcdatad(mm)%uslip(i, j, 3)
              derivrotationmatrixd(3, 1) = derivrotationmatrixd(3, 1) + &
&               xxc(1)*bcdatad(mm)%uslip(i, j, 3)
              xxcd(1) = xxcd(1) + derivrotationmatrix(3, 1)*bcdatad(mm)%&
&               uslip(i, j, 3) + derivrotationmatrix(2, 1)*bcdatad(mm)%&
&               uslip(i, j, 2) + derivrotationmatrix(1, 1)*bcdatad(mm)%&
&               uslip(i, j, 1)
              derivrotationmatrixd(3, 2) = derivrotationmatrixd(3, 2) + &
&               xxc(2)*bcdatad(mm)%uslip(i, j, 3)
              xxcd(2) = xxcd(2) + derivrotationmatrix(3, 2)*bcdatad(mm)%&
&               uslip(i, j, 3) + derivrotationmatrix(2, 2)*bcdatad(mm)%&
&               uslip(i, j, 2) + derivrotationmatrix(1, 2)*bcdatad(mm)%&
&               uslip(i, j, 1)
              derivrotationmatrixd(3, 3) = derivrotationmatrixd(3, 3) + &
&               xxc(3)*bcdatad(mm)%uslip(i, j, 3)
              xxcd(3) = xxcd(3) + derivrotationmatrix(3, 3)*bcdatad(mm)%&
&               uslip(i, j, 3) + derivrotationmatrix(2, 3)*bcdatad(mm)%&
&               uslip(i, j, 2) + derivrotationmatrix(1, 3)*bcdatad(mm)%&
&               uslip(i, j, 1)
              velygridd = velygridd + bcdatad(mm)%uslip(i, j, 2)
              derivrotationmatrixd(2, 1) = derivrotationmatrixd(2, 1) + &
&               xxc(1)*bcdatad(mm)%uslip(i, j, 2)
              derivrotationmatrixd(2, 2) = derivrotationmatrixd(2, 2) + &
&               xxc(2)*bcdatad(mm)%uslip(i, j, 2)
              derivrotationmatrixd(2, 3) = derivrotationmatrixd(2, 3) + &
&               xxc(3)*bcdatad(mm)%uslip(i, j, 2)
              velxgridd = velxgridd + bcdatad(mm)%uslip(i, j, 1)
              derivrotationmatrixd(1, 1) = derivrotationmatrixd(1, 1) + &
&               xxc(1)*bcdatad(mm)%uslip(i, j, 1)
              derivrotationmatrixd(1, 2) = derivrotationmatrixd(1, 2) + &
&               xxc(2)*bcdatad(mm)%uslip(i, j, 1)
              derivrotationmatrixd(1, 3) = derivrotationmatrixd(1, 3) + &
&               xxc(3)*bcdatad(mm)%uslip(i, j, 1)
              call popreal8(xxc(3))
              xcd(3) = xcd(3) + xxcd(3)
              call popreal8(xxc(2))
              xcd(2) = xcd(2) + xxcd(2)
              xxcd(2) = rotrate(1)*bcdatad(mm)%uslip(i, j, 3)
              call popreal8(xxc(1))
              xcd(1) = xcd(1) + xxcd(1)
              xxcd(1) = -(rotrate(2)*bcdatad(mm)%uslip(i, j, 3))
              rotrated(1) = rotrated(1) + xxc(2)*bcdatad(mm)%uslip(i, j&
&               , 3)
              rotrated(2) = rotrated(2) - xxc(1)*bcdatad(mm)%uslip(i, j&
&               , 3)
              bcdatad(mm)%uslip(i, j, 3) = 0.0_8
              xxcd(3) = -(rotrate(1)*bcdatad(mm)%uslip(i, j, 2))
              rotrated(3) = rotrated(3) + xxc(1)*bcdatad(mm)%uslip(i, j&
&               , 2)
              xxcd(1) = xxcd(1) + rotrate(3)*bcdatad(mm)%uslip(i, j, 2)
              rotrated(1) = rotrated(1) - xxc(3)*bcdatad(mm)%uslip(i, j&
&               , 2)
              bcdatad(mm)%uslip(i, j, 2) = 0.0_8
              rotrated(2) = rotrated(2) + xxc(3)*bcdatad(mm)%uslip(i, j&
&               , 1)
              xxcd(3) = xxcd(3) + rotrate(2)*bcdatad(mm)%uslip(i, j, 1)
              rotrated(3) = rotrated(3) - xxc(2)*bcdatad(mm)%uslip(i, j&
&               , 1)
              xxcd(2) = xxcd(2) - rotrate(3)*bcdatad(mm)%uslip(i, j, 1)
              bcdatad(mm)%uslip(i, j, 1) = 0.0_8
              call popreal8(xxc(3))
              xcd(3) = xcd(3) + xxcd(3)
              xxcd(3) = 0.0_8
              call popreal8(xxc(2))
              xcd(2) = xcd(2) + xxcd(2)
              xxcd(2) = 0.0_8
              call popreal8(xxc(1))
              xcd(1) = xcd(1) + xxcd(1)
              xxcd(1) = 0.0_8
              tempd = fourth*xcd(3)
              xcd(3) = 0.0_8
              flowdomsd(nn, groundlevel, sps)%x(il, i, j, 3) = flowdomsd&
&               (nn, groundlevel, sps)%x(il, i, j, 3) + tempd
              flowdomsd(nn, groundlevel, sps)%x(il, i, j-1, 3) = &
&               flowdomsd(nn, groundlevel, sps)%x(il, i, j-1, 3) + tempd
              flowdomsd(nn, groundlevel, sps)%x(il, i-1, j, 3) = &
&               flowdomsd(nn, groundlevel, sps)%x(il, i-1, j, 3) + tempd
              flowdomsd(nn, groundlevel, sps)%x(il, i-1, j-1, 3) = &
&               flowdomsd(nn, groundlevel, sps)%x(il, i-1, j-1, 3) + &
&               tempd
              tempd = fourth*xcd(2)
              xcd(2) = 0.0_8
              flowdomsd(nn, groundlevel, sps)%x(il, i, j, 2) = flowdomsd&
&               (nn, groundlevel, sps)%x(il, i, j, 2) + tempd
              flowdomsd(nn, groundlevel, sps)%x(il, i, j-1, 2) = &
&               flowdomsd(nn, groundlevel, sps)%x(il, i, j-1, 2) + tempd
              flowdomsd(nn, groundlevel, sps)%x(il, i-1, j, 2) = &
&               flowdomsd(nn, groundlevel, sps)%x(il, i-1, j, 2) + tempd
              flowdomsd(nn, groundlevel, sps)%x(il, i-1, j-1, 2) = &
&               flowdomsd(nn, groundlevel, sps)%x(il, i-1, j-1, 2) + &
&               tempd
              tempd = fourth*xcd(1)
              xcd(1) = 0.0_8
              flowdomsd(nn, groundlevel, sps)%x(il, i, j, 1) = flowdomsd&
&               (nn, groundlevel, sps)%x(il, i, j, 1) + tempd
              flowdomsd(nn, groundlevel, sps)%x(il, i, j-1, 1) = &
&               flowdomsd(nn, groundlevel, sps)%x(il, i, j-1, 1) + tempd
              flowdomsd(nn, groundlevel, sps)%x(il, i-1, j, 1) = &
&               flowdomsd(nn, groundlevel, sps)%x(il, i-1, j, 1) + tempd
              flowdomsd(nn, groundlevel, sps)%x(il, i-1, j-1, 1) = &
&               flowdomsd(nn, groundlevel, sps)%x(il, i-1, j-1, 1) + &
&               tempd
            end do
          end do
        else
          rotrated = 0.0_8
          velygridd = 0.0_8
          velzgridd = 0.0_8
          velxgridd = 0.0_8
          call popinteger4(ad_from0)
          call popinteger4(ad_to0)
          do j=ad_to0,ad_from0,-1
            call popinteger4(ad_from)
            call popinteger4(ad_to)
            do i=ad_to,ad_from,-1
              velzgridd = velzgridd + bcdatad(mm)%uslip(i, j, 3)
              derivrotationmatrixd(3, 1) = derivrotationmatrixd(3, 1) + &
&               xxc(1)*bcdatad(mm)%uslip(i, j, 3)
              xxcd(1) = xxcd(1) + derivrotationmatrix(3, 1)*bcdatad(mm)%&
&               uslip(i, j, 3) + derivrotationmatrix(2, 1)*bcdatad(mm)%&
&               uslip(i, j, 2) + derivrotationmatrix(1, 1)*bcdatad(mm)%&
&               uslip(i, j, 1)
              derivrotationmatrixd(3, 2) = derivrotationmatrixd(3, 2) + &
&               xxc(2)*bcdatad(mm)%uslip(i, j, 3)
              xxcd(2) = xxcd(2) + derivrotationmatrix(3, 2)*bcdatad(mm)%&
&               uslip(i, j, 3) + derivrotationmatrix(2, 2)*bcdatad(mm)%&
&               uslip(i, j, 2) + derivrotationmatrix(1, 2)*bcdatad(mm)%&
&               uslip(i, j, 1)
              derivrotationmatrixd(3, 3) = derivrotationmatrixd(3, 3) + &
&               xxc(3)*bcdatad(mm)%uslip(i, j, 3)
              xxcd(3) = xxcd(3) + derivrotationmatrix(3, 3)*bcdatad(mm)%&
&               uslip(i, j, 3) + derivrotationmatrix(2, 3)*bcdatad(mm)%&
&               uslip(i, j, 2) + derivrotationmatrix(1, 3)*bcdatad(mm)%&
&               uslip(i, j, 1)
              velygridd = velygridd + bcdatad(mm)%uslip(i, j, 2)
              derivrotationmatrixd(2, 1) = derivrotationmatrixd(2, 1) + &
&               xxc(1)*bcdatad(mm)%uslip(i, j, 2)
              derivrotationmatrixd(2, 2) = derivrotationmatrixd(2, 2) + &
&               xxc(2)*bcdatad(mm)%uslip(i, j, 2)
              derivrotationmatrixd(2, 3) = derivrotationmatrixd(2, 3) + &
&               xxc(3)*bcdatad(mm)%uslip(i, j, 2)
              velxgridd = velxgridd + bcdatad(mm)%uslip(i, j, 1)
              derivrotationmatrixd(1, 1) = derivrotationmatrixd(1, 1) + &
&               xxc(1)*bcdatad(mm)%uslip(i, j, 1)
              derivrotationmatrixd(1, 2) = derivrotationmatrixd(1, 2) + &
&               xxc(2)*bcdatad(mm)%uslip(i, j, 1)
              derivrotationmatrixd(1, 3) = derivrotationmatrixd(1, 3) + &
&               xxc(3)*bcdatad(mm)%uslip(i, j, 1)
              call popreal8(xxc(3))
              xcd(3) = xcd(3) + xxcd(3)
              call popreal8(xxc(2))
              xcd(2) = xcd(2) + xxcd(2)
              xxcd(2) = rotrate(1)*bcdatad(mm)%uslip(i, j, 3)
              call popreal8(xxc(1))
              xcd(1) = xcd(1) + xxcd(1)
              xxcd(1) = -(rotrate(2)*bcdatad(mm)%uslip(i, j, 3))
              rotrated(1) = rotrated(1) + xxc(2)*bcdatad(mm)%uslip(i, j&
&               , 3)
              rotrated(2) = rotrated(2) - xxc(1)*bcdatad(mm)%uslip(i, j&
&               , 3)
              bcdatad(mm)%uslip(i, j, 3) = 0.0_8
              xxcd(3) = -(rotrate(1)*bcdatad(mm)%uslip(i, j, 2))
              rotrated(3) = rotrated(3) + xxc(1)*bcdatad(mm)%uslip(i, j&
&               , 2)
              xxcd(1) = xxcd(1) + rotrate(3)*bcdatad(mm)%uslip(i, j, 2)
              rotrated(1) = rotrated(1) - xxc(3)*bcdatad(mm)%uslip(i, j&
&               , 2)
              bcdatad(mm)%uslip(i, j, 2) = 0.0_8
              rotrated(2) = rotrated(2) + xxc(3)*bcdatad(mm)%uslip(i, j&
&               , 1)
              xxcd(3) = xxcd(3) + rotrate(2)*bcdatad(mm)%uslip(i, j, 1)
              rotrated(3) = rotrated(3) - xxc(2)*bcdatad(mm)%uslip(i, j&
&               , 1)
              xxcd(2) = xxcd(2) - rotrate(3)*bcdatad(mm)%uslip(i, j, 1)
              bcdatad(mm)%uslip(i, j, 1) = 0.0_8
              call popreal8(xxc(3))
              xcd(3) = xcd(3) + xxcd(3)
              xxcd(3) = 0.0_8
              call popreal8(xxc(2))
              xcd(2) = xcd(2) + xxcd(2)
              xxcd(2) = 0.0_8
              call popreal8(xxc(1))
              xcd(1) = xcd(1) + xxcd(1)
              xxcd(1) = 0.0_8
              tempd = fourth*xcd(3)
              xcd(3) = 0.0_8
              flowdomsd(nn, groundlevel, sps)%x(1, i, j, 3) = flowdomsd(&
&               nn, groundlevel, sps)%x(1, i, j, 3) + tempd
              flowdomsd(nn, groundlevel, sps)%x(1, i, j-1, 3) = &
&               flowdomsd(nn, groundlevel, sps)%x(1, i, j-1, 3) + tempd
              flowdomsd(nn, groundlevel, sps)%x(1, i-1, j, 3) = &
&               flowdomsd(nn, groundlevel, sps)%x(1, i-1, j, 3) + tempd
              flowdomsd(nn, groundlevel, sps)%x(1, i-1, j-1, 3) = &
&               flowdomsd(nn, groundlevel, sps)%x(1, i-1, j-1, 3) + &
&               tempd
              tempd = fourth*xcd(2)
              xcd(2) = 0.0_8
              flowdomsd(nn, groundlevel, sps)%x(1, i, j, 2) = flowdomsd(&
&               nn, groundlevel, sps)%x(1, i, j, 2) + tempd
              flowdomsd(nn, groundlevel, sps)%x(1, i, j-1, 2) = &
&               flowdomsd(nn, groundlevel, sps)%x(1, i, j-1, 2) + tempd
              flowdomsd(nn, groundlevel, sps)%x(1, i-1, j, 2) = &
&               flowdomsd(nn, groundlevel, sps)%x(1, i-1, j, 2) + tempd
              flowdomsd(nn, groundlevel, sps)%x(1, i-1, j-1, 2) = &
&               flowdomsd(nn, groundlevel, sps)%x(1, i-1, j-1, 2) + &
&               tempd
              tempd = fourth*xcd(1)
              xcd(1) = 0.0_8
              flowdomsd(nn, groundlevel, sps)%x(1, i, j, 1) = flowdomsd(&
&               nn, groundlevel, sps)%x(1, i, j, 1) + tempd
              flowdomsd(nn, groundlevel, sps)%x(1, i, j-1, 1) = &
&               flowdomsd(nn, groundlevel, sps)%x(1, i, j-1, 1) + tempd
              flowdomsd(nn, groundlevel, sps)%x(1, i-1, j, 1) = &
&               flowdomsd(nn, groundlevel, sps)%x(1, i-1, j, 1) + tempd
              flowdomsd(nn, groundlevel, sps)%x(1, i-1, j-1, 1) = &
&               flowdomsd(nn, groundlevel, sps)%x(1, i-1, j-1, 1) + &
&               tempd
            end do
          end do
        end if
        velzgrid0d = velzgrid0d + velzgridd
        velygrid0d = velygrid0d + velygridd
        velxgrid0d = velxgrid0d + velxgridd
        call popreal8array(rotrate, 3)
        timerefd = timerefd + sum(cgnsdoms(nbkglobal)%bocoinfo(ii)%&
&         rotrate*rotrated)
        call popinteger4(ii)
      end do
      call derivativerotmatrixrigid_b(derivrotationmatrix, &
&                               derivrotationmatrixd, rotationpoint, t(1&
&                               ))
      ainfd = -(machgrid*veldirfreestream(3)*velzgrid0d) - machgrid*&
&       veldirfreestream(2)*velygrid0d - machgrid*veldirfreestream(1)*&
&       velxgrid0d
      machgridd = -(ainf*veldirfreestream(3)*velzgrid0d) - ainf*&
&       veldirfreestream(2)*velygrid0d - ainf*veldirfreestream(1)*&
&       velxgrid0d
      veldirfreestreamd(3) = veldirfreestreamd(3) - ainf*machgrid*&
&       velzgrid0d
      veldirfreestreamd(2) = veldirfreestreamd(2) - ainf*machgrid*&
&       velygrid0d
      veldirfreestreamd(1) = veldirfreestreamd(1) - ainf*machgrid*&
&       velxgrid0d
      if (gammainf*(pinf/rhoinf) .eq. 0.0_8) then
        tempd = 0.0_8
      else
        tempd = gammainf*ainfd/(rhoinf*2.0*sqrt(gammainf*(pinf/rhoinf)))
      end if
      pinfd = tempd
      rhoinfd = -(pinf*tempd/rhoinf)
    end if
  end subroutine slipvelocitiesfinelevel_block_b
 
  subroutine slipvelocitiesfinelevel_block(useoldcoor, t, sps, nn)
!
!       slipvelocitiesfinelevel computes the slip velocities for
!       viscous subfaces on all viscous boundaries on groundlevel for
!       the given spectral solution. if useoldcoor is .true. the
!       velocities are determined using the unsteady time integrator;
!       otherwise the analytic form is used.
!
    use constants
    use inputtimespectral
    use blockpointers
    use cgnsgrid
    use flowvarrefstate
    use inputmotion
    use inputunsteady
    use iteration
    use inputphysics
    use inputtsstabderiv
    use monitor
    use communication
    use flowutils_b, only : derivativerotmatrixrigid, getdirvector
    use utils_b, only : tsalpha, tsbeta, tsmach, terminate, &
&   rotmatrixrigidbody, setcoeftimeintegrator, getdirangle
    implicit none
!
!      subroutine arguments.
!
    integer(kind=inttype), intent(in) :: sps, nn
    logical, intent(in) :: useoldcoor
    real(kind=realtype), dimension(*), intent(in) :: t
!
!      local variables.
!
    integer(kind=inttype) :: mm, i, j, level, ii
    real(kind=realtype) :: oneover4dt
    real(kind=realtype) :: velxgrid, velygrid, velzgrid, ainf
    real(kind=realtype) :: velxgrid0, velygrid0, velzgrid0
    real(kind=realtype), dimension(3) :: xc, xxc
    real(kind=realtype), dimension(3) :: rotcenter, rotrate
    real(kind=realtype), dimension(3) :: rotationpoint
    real(kind=realtype), dimension(3, 3) :: rotationmatrix, &
&   derivrotationmatrix
    real(kind=realtype) :: tnew, told
    real(kind=realtype), dimension(:, :, :), pointer :: uslip
    real(kind=realtype), dimension(:, :, :), pointer :: xface
    real(kind=realtype), dimension(:, :, :, :), pointer :: xfaceold
    real(kind=realtype) :: intervalmach, alphats, alphaincrement, betats&
&   , betaincrement
    real(kind=realtype), dimension(3) :: veldir
    real(kind=realtype), dimension(3) :: refdirection
    intrinsic sqrt
! determine the situation we are having here.
    if (.not.useoldcoor) then
! the velocities must be determined analytically.
! compute the mesh velocity from the given mesh mach number.
!  ainf = sqrt(gammainf*pinf/rhoinf)
!  velxgrid = ainf*machgrid(1)
!  velygrid = ainf*machgrid(2)
!  velzgrid = ainf*machgrid(3)
      ainf = sqrt(gammainf*pinf/rhoinf)
      velxgrid0 = ainf*machgrid*(-veldirfreestream(1))
      velygrid0 = ainf*machgrid*(-veldirfreestream(2))
      velzgrid0 = ainf*machgrid*(-veldirfreestream(3))
! compute the derivative of the rotation matrix and the rotation
! point; needed for velocity due to the rigid body rotation of
! the entire grid. it is assumed that the rigid body motion of
! the grid is only specified if there is only 1 section present.
      call derivativerotmatrixrigid(derivrotationmatrix, rotationpoint, &
&                             t(1))
!compute the rotation matrix to update the velocities for the time
!spectral stability derivative case...
! loop over the number of viscous subfaces.
bocoloop2:do mm=1,nviscbocos
! store the rotation center and the rotation rate
! for this subface.
        ii = cgnssubface(mm)
        rotcenter = cgnsdoms(nbkglobal)%bocoinfo(ii)%rotcenter
        rotrate = timeref*cgnsdoms(nbkglobal)%bocoinfo(ii)%rotrate
! usewindaxis should go back here!
        velxgrid = velxgrid0
        velygrid = velygrid0
        velzgrid = velzgrid0
! loop over the quadrilateral faces of the viscous
! subface.
! the new procedure is less elegant as the previous one.
! but the new stands up to tapenade.
        if (bcfaceid(mm) .eq. imin) then
          do j=bcdata(mm)%jcbeg,bcdata(mm)%jcend
            do i=bcdata(mm)%icbeg,bcdata(mm)%icend
! compute the coordinates of the centroid of the face.
! normally this would be an average of i-1 and i, but
! due to the usage of the pointer xface and the fact
! that x starts at index 0 this is shifted 1 index.
              xc(1) = fourth*(flowdoms(nn, groundlevel, sps)%x(1, i, j, &
&               1)+flowdoms(nn, groundlevel, sps)%x(1, i, j-1, 1)+&
&               flowdoms(nn, groundlevel, sps)%x(1, i-1, j, 1)+flowdoms(&
&               nn, groundlevel, sps)%x(1, i-1, j-1, 1))
              xc(2) = fourth*(flowdoms(nn, groundlevel, sps)%x(1, i, j, &
&               2)+flowdoms(nn, groundlevel, sps)%x(1, i, j-1, 2)+&
&               flowdoms(nn, groundlevel, sps)%x(1, i-1, j, 2)+flowdoms(&
&               nn, groundlevel, sps)%x(1, i-1, j-1, 2))
              xc(3) = fourth*(flowdoms(nn, groundlevel, sps)%x(1, i, j, &
&               3)+flowdoms(nn, groundlevel, sps)%x(1, i, j-1, 3)+&
&               flowdoms(nn, groundlevel, sps)%x(1, i-1, j, 3)+flowdoms(&
&               nn, groundlevel, sps)%x(1, i-1, j-1, 3))
! determine the coordinates relative to the center
! of rotation.
              xxc(1) = xc(1) - rotcenter(1)
              xxc(2) = xc(2) - rotcenter(2)
              xxc(3) = xc(3) - rotcenter(3)
! compute the velocity, which is the cross product
! of rotrate and xc.
              bcdata(mm)%uslip(i, j, 1) = rotrate(2)*xxc(3) - rotrate(3)&
&               *xxc(2)
              bcdata(mm)%uslip(i, j, 2) = rotrate(3)*xxc(1) - rotrate(1)&
&               *xxc(3)
              bcdata(mm)%uslip(i, j, 3) = rotrate(1)*xxc(2) - rotrate(2)&
&               *xxc(1)
! determine the coordinates relative to the
! rigid body rotation point.
              xxc(1) = xc(1) - rotationpoint(1)
              xxc(2) = xc(2) - rotationpoint(2)
              xxc(3) = xc(3) - rotationpoint(3)
! determine the total velocity of the cell center.
! this is a combination of rotation speed of this
! block and the entire rigid body rotation.
              bcdata(mm)%uslip(i, j, 1) = bcdata(mm)%uslip(i, j, 1) + &
&               velxgrid + derivrotationmatrix(1, 1)*xxc(1) + &
&               derivrotationmatrix(1, 2)*xxc(2) + derivrotationmatrix(1&
&               , 3)*xxc(3)
              bcdata(mm)%uslip(i, j, 2) = bcdata(mm)%uslip(i, j, 2) + &
&               velygrid + derivrotationmatrix(2, 1)*xxc(1) + &
&               derivrotationmatrix(2, 2)*xxc(2) + derivrotationmatrix(2&
&               , 3)*xxc(3)
              bcdata(mm)%uslip(i, j, 3) = bcdata(mm)%uslip(i, j, 3) + &
&               velzgrid + derivrotationmatrix(3, 1)*xxc(1) + &
&               derivrotationmatrix(3, 2)*xxc(2) + derivrotationmatrix(3&
&               , 3)*xxc(3)
            end do
          end do
        else if (bcfaceid(mm) .eq. imax) then
          do j=bcdata(mm)%jcbeg,bcdata(mm)%jcend
            do i=bcdata(mm)%icbeg,bcdata(mm)%icend
! compute the coordinates of the centroid of the face.
! normally this would be an average of i-1 and i, but
! due to the usage of the pointer xface and the fact
! that x starts at index 0 this is shifted 1 index.
              xc(1) = fourth*(flowdoms(nn, groundlevel, sps)%x(il, i, j&
&               , 1)+flowdoms(nn, groundlevel, sps)%x(il, i, j-1, 1)+&
&               flowdoms(nn, groundlevel, sps)%x(il, i-1, j, 1)+flowdoms&
&               (nn, groundlevel, sps)%x(il, i-1, j-1, 1))
              xc(2) = fourth*(flowdoms(nn, groundlevel, sps)%x(il, i, j&
&               , 2)+flowdoms(nn, groundlevel, sps)%x(il, i, j-1, 2)+&
&               flowdoms(nn, groundlevel, sps)%x(il, i-1, j, 2)+flowdoms&
&               (nn, groundlevel, sps)%x(il, i-1, j-1, 2))
              xc(3) = fourth*(flowdoms(nn, groundlevel, sps)%x(il, i, j&
&               , 3)+flowdoms(nn, groundlevel, sps)%x(il, i, j-1, 3)+&
&               flowdoms(nn, groundlevel, sps)%x(il, i-1, j, 3)+flowdoms&
&               (nn, groundlevel, sps)%x(il, i-1, j-1, 3))
! determine the coordinates relative to the center
! of rotation.
              xxc(1) = xc(1) - rotcenter(1)
              xxc(2) = xc(2) - rotcenter(2)
              xxc(3) = xc(3) - rotcenter(3)
! compute the velocity, which is the cross product
! of rotrate and xc.
              bcdata(mm)%uslip(i, j, 1) = rotrate(2)*xxc(3) - rotrate(3)&
&               *xxc(2)
              bcdata(mm)%uslip(i, j, 2) = rotrate(3)*xxc(1) - rotrate(1)&
&               *xxc(3)
              bcdata(mm)%uslip(i, j, 3) = rotrate(1)*xxc(2) - rotrate(2)&
&               *xxc(1)
! determine the coordinates relative to the
! rigid body rotation point.
              xxc(1) = xc(1) - rotationpoint(1)
              xxc(2) = xc(2) - rotationpoint(2)
              xxc(3) = xc(3) - rotationpoint(3)
! determine the total velocity of the cell center.
! this is a combination of rotation speed of this
! block and the entire rigid body rotation.
              bcdata(mm)%uslip(i, j, 1) = bcdata(mm)%uslip(i, j, 1) + &
&               velxgrid + derivrotationmatrix(1, 1)*xxc(1) + &
&               derivrotationmatrix(1, 2)*xxc(2) + derivrotationmatrix(1&
&               , 3)*xxc(3)
              bcdata(mm)%uslip(i, j, 2) = bcdata(mm)%uslip(i, j, 2) + &
&               velygrid + derivrotationmatrix(2, 1)*xxc(1) + &
&               derivrotationmatrix(2, 2)*xxc(2) + derivrotationmatrix(2&
&               , 3)*xxc(3)
              bcdata(mm)%uslip(i, j, 3) = bcdata(mm)%uslip(i, j, 3) + &
&               velzgrid + derivrotationmatrix(3, 1)*xxc(1) + &
&               derivrotationmatrix(3, 2)*xxc(2) + derivrotationmatrix(3&
&               , 3)*xxc(3)
            end do
          end do
        else if (bcfaceid(mm) .eq. jmin) then
          do j=bcdata(mm)%jcbeg,bcdata(mm)%jcend
            do i=bcdata(mm)%icbeg,bcdata(mm)%icend
! compute the coordinates of the centroid of the face.
! normally this would be an average of i-1 and i, but
! due to the usage of the pointer xface and the fact
! that x starts at index 0 this is shifted 1 index.
              xc(1) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, 1, j, &
&               1)+flowdoms(nn, groundlevel, sps)%x(i, 1, j-1, 1)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, 1, j, 1)+flowdoms(&
&               nn, groundlevel, sps)%x(i-1, 1, j-1, 1))
              xc(2) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, 1, j, &
&               2)+flowdoms(nn, groundlevel, sps)%x(i, 1, j-1, 2)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, 1, j, 2)+flowdoms(&
&               nn, groundlevel, sps)%x(i-1, 1, j-1, 2))
              xc(3) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, 1, j, &
&               3)+flowdoms(nn, groundlevel, sps)%x(i, 1, j-1, 3)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, 1, j, 3)+flowdoms(&
&               nn, groundlevel, sps)%x(i-1, 1, j-1, 3))
! determine the coordinates relative to the center
! of rotation.
              xxc(1) = xc(1) - rotcenter(1)
              xxc(2) = xc(2) - rotcenter(2)
              xxc(3) = xc(3) - rotcenter(3)
! compute the velocity, which is the cross product
! of rotrate and xc.
              bcdata(mm)%uslip(i, j, 1) = rotrate(2)*xxc(3) - rotrate(3)&
&               *xxc(2)
              bcdata(mm)%uslip(i, j, 2) = rotrate(3)*xxc(1) - rotrate(1)&
&               *xxc(3)
              bcdata(mm)%uslip(i, j, 3) = rotrate(1)*xxc(2) - rotrate(2)&
&               *xxc(1)
! determine the coordinates relative to the
! rigid body rotation point.
              xxc(1) = xc(1) - rotationpoint(1)
              xxc(2) = xc(2) - rotationpoint(2)
              xxc(3) = xc(3) - rotationpoint(3)
! determine the total velocity of the cell center.
! this is a combination of rotation speed of this
! block and the entire rigid body rotation.
              bcdata(mm)%uslip(i, j, 1) = bcdata(mm)%uslip(i, j, 1) + &
&               velxgrid + derivrotationmatrix(1, 1)*xxc(1) + &
&               derivrotationmatrix(1, 2)*xxc(2) + derivrotationmatrix(1&
&               , 3)*xxc(3)
              bcdata(mm)%uslip(i, j, 2) = bcdata(mm)%uslip(i, j, 2) + &
&               velygrid + derivrotationmatrix(2, 1)*xxc(1) + &
&               derivrotationmatrix(2, 2)*xxc(2) + derivrotationmatrix(2&
&               , 3)*xxc(3)
              bcdata(mm)%uslip(i, j, 3) = bcdata(mm)%uslip(i, j, 3) + &
&               velzgrid + derivrotationmatrix(3, 1)*xxc(1) + &
&               derivrotationmatrix(3, 2)*xxc(2) + derivrotationmatrix(3&
&               , 3)*xxc(3)
            end do
          end do
        else if (bcfaceid(mm) .eq. jmax) then
          do j=bcdata(mm)%jcbeg,bcdata(mm)%jcend
            do i=bcdata(mm)%icbeg,bcdata(mm)%icend
! compute the coordinates of the centroid of the face.
! normally this would be an average of i-1 and i, but
! due to the usage of the pointer xface and the fact
! that x starts at index 0 this is shifted 1 index.
              xc(1) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, jl, j&
&               , 1)+flowdoms(nn, groundlevel, sps)%x(i, jl, j-1, 1)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, jl, j, 1)+flowdoms&
&               (nn, groundlevel, sps)%x(i-1, jl, j-1, 1))
              xc(2) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, jl, j&
&               , 2)+flowdoms(nn, groundlevel, sps)%x(i, jl, j-1, 2)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, jl, j, 2)+flowdoms&
&               (nn, groundlevel, sps)%x(i-1, jl, j-1, 2))
              xc(3) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, jl, j&
&               , 3)+flowdoms(nn, groundlevel, sps)%x(i, jl, j-1, 3)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, jl, j, 3)+flowdoms&
&               (nn, groundlevel, sps)%x(i-1, jl, j-1, 3))
! determine the coordinates relative to the center
! of rotation.
              xxc(1) = xc(1) - rotcenter(1)
              xxc(2) = xc(2) - rotcenter(2)
              xxc(3) = xc(3) - rotcenter(3)
! compute the velocity, which is the cross product
! of rotrate and xc.
              bcdata(mm)%uslip(i, j, 1) = rotrate(2)*xxc(3) - rotrate(3)&
&               *xxc(2)
              bcdata(mm)%uslip(i, j, 2) = rotrate(3)*xxc(1) - rotrate(1)&
&               *xxc(3)
              bcdata(mm)%uslip(i, j, 3) = rotrate(1)*xxc(2) - rotrate(2)&
&               *xxc(1)
! determine the coordinates relative to the
! rigid body rotation point.
              xxc(1) = xc(1) - rotationpoint(1)
              xxc(2) = xc(2) - rotationpoint(2)
              xxc(3) = xc(3) - rotationpoint(3)
! determine the total velocity of the cell center.
! this is a combination of rotation speed of this
! block and the entire rigid body rotation.
              bcdata(mm)%uslip(i, j, 1) = bcdata(mm)%uslip(i, j, 1) + &
&               velxgrid + derivrotationmatrix(1, 1)*xxc(1) + &
&               derivrotationmatrix(1, 2)*xxc(2) + derivrotationmatrix(1&
&               , 3)*xxc(3)
              bcdata(mm)%uslip(i, j, 2) = bcdata(mm)%uslip(i, j, 2) + &
&               velygrid + derivrotationmatrix(2, 1)*xxc(1) + &
&               derivrotationmatrix(2, 2)*xxc(2) + derivrotationmatrix(2&
&               , 3)*xxc(3)
              bcdata(mm)%uslip(i, j, 3) = bcdata(mm)%uslip(i, j, 3) + &
&               velzgrid + derivrotationmatrix(3, 1)*xxc(1) + &
&               derivrotationmatrix(3, 2)*xxc(2) + derivrotationmatrix(3&
&               , 3)*xxc(3)
            end do
          end do
        else if (bcfaceid(mm) .eq. kmin) then
          do j=bcdata(mm)%jcbeg,bcdata(mm)%jcend
            do i=bcdata(mm)%icbeg,bcdata(mm)%icend
! compute the coordinates of the centroid of the face.
! normally this would be an average of i-1 and i, but
! due to the usage of the pointer xface and the fact
! that x starts at index 0 this is shifted 1 index.
              xc(1) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, j, 1, &
&               1)+flowdoms(nn, groundlevel, sps)%x(i, j-1, 1, 1)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, j, 1, 1)+flowdoms(&
&               nn, groundlevel, sps)%x(i-1, j-1, 1, 1))
              xc(2) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, j, 1, &
&               2)+flowdoms(nn, groundlevel, sps)%x(i, j-1, 1, 2)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, j, 1, 2)+flowdoms(&
&               nn, groundlevel, sps)%x(i-1, j-1, 1, 2))
              xc(3) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, j, 1, &
&               3)+flowdoms(nn, groundlevel, sps)%x(i, j-1, 1, 3)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, j, 1, 3)+flowdoms(&
&               nn, groundlevel, sps)%x(i-1, j-1, 1, 3))
! determine the coordinates relative to the center
! of rotation.
              xxc(1) = xc(1) - rotcenter(1)
              xxc(2) = xc(2) - rotcenter(2)
              xxc(3) = xc(3) - rotcenter(3)
! compute the velocity, which is the cross product
! of rotrate and xc.
              bcdata(mm)%uslip(i, j, 1) = rotrate(2)*xxc(3) - rotrate(3)&
&               *xxc(2)
              bcdata(mm)%uslip(i, j, 2) = rotrate(3)*xxc(1) - rotrate(1)&
&               *xxc(3)
              bcdata(mm)%uslip(i, j, 3) = rotrate(1)*xxc(2) - rotrate(2)&
&               *xxc(1)
! determine the coordinates relative to the
! rigid body rotation point.
              xxc(1) = xc(1) - rotationpoint(1)
              xxc(2) = xc(2) - rotationpoint(2)
              xxc(3) = xc(3) - rotationpoint(3)
! determine the total velocity of the cell center.
! this is a combination of rotation speed of this
! block and the entire rigid body rotation.
              bcdata(mm)%uslip(i, j, 1) = bcdata(mm)%uslip(i, j, 1) + &
&               velxgrid + derivrotationmatrix(1, 1)*xxc(1) + &
&               derivrotationmatrix(1, 2)*xxc(2) + derivrotationmatrix(1&
&               , 3)*xxc(3)
              bcdata(mm)%uslip(i, j, 2) = bcdata(mm)%uslip(i, j, 2) + &
&               velygrid + derivrotationmatrix(2, 1)*xxc(1) + &
&               derivrotationmatrix(2, 2)*xxc(2) + derivrotationmatrix(2&
&               , 3)*xxc(3)
              bcdata(mm)%uslip(i, j, 3) = bcdata(mm)%uslip(i, j, 3) + &
&               velzgrid + derivrotationmatrix(3, 1)*xxc(1) + &
&               derivrotationmatrix(3, 2)*xxc(2) + derivrotationmatrix(3&
&               , 3)*xxc(3)
            end do
          end do
        else if (bcfaceid(mm) .eq. kmax) then
          do j=bcdata(mm)%jcbeg,bcdata(mm)%jcend
            do i=bcdata(mm)%icbeg,bcdata(mm)%icend
! compute the coordinates of the centroid of the face.
! normally this would be an average of i-1 and i, but
! due to the usage of the pointer xface and the fact
! that x starts at index 0 this is shifted 1 index.
              xc(1) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, j, kl&
&               , 1)+flowdoms(nn, groundlevel, sps)%x(i, j-1, kl, 1)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, j, kl, 1)+flowdoms&
&               (nn, groundlevel, sps)%x(i-1, j-1, kl, 1))
              xc(2) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, j, kl&
&               , 2)+flowdoms(nn, groundlevel, sps)%x(i, j-1, kl, 2)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, j, kl, 2)+flowdoms&
&               (nn, groundlevel, sps)%x(i-1, j-1, kl, 2))
              xc(3) = fourth*(flowdoms(nn, groundlevel, sps)%x(i, j, kl&
&               , 3)+flowdoms(nn, groundlevel, sps)%x(i, j-1, kl, 3)+&
&               flowdoms(nn, groundlevel, sps)%x(i-1, j, kl, 3)+flowdoms&
&               (nn, groundlevel, sps)%x(i-1, j-1, kl, 3))
! determine the coordinates relative to the center
! of rotation.
              xxc(1) = xc(1) - rotcenter(1)
              xxc(2) = xc(2) - rotcenter(2)
              xxc(3) = xc(3) - rotcenter(3)
! compute the velocity, which is the cross product
! of rotrate and xc.
              bcdata(mm)%uslip(i, j, 1) = rotrate(2)*xxc(3) - rotrate(3)&
&               *xxc(2)
              bcdata(mm)%uslip(i, j, 2) = rotrate(3)*xxc(1) - rotrate(1)&
&               *xxc(3)
              bcdata(mm)%uslip(i, j, 3) = rotrate(1)*xxc(2) - rotrate(2)&
&               *xxc(1)
! determine the coordinates relative to the
! rigid body rotation point.
              xxc(1) = xc(1) - rotationpoint(1)
              xxc(2) = xc(2) - rotationpoint(2)
              xxc(3) = xc(3) - rotationpoint(3)
! determine the total velocity of the cell center.
! this is a combination of rotation speed of this
! block and the entire rigid body rotation.
              bcdata(mm)%uslip(i, j, 1) = bcdata(mm)%uslip(i, j, 1) + &
&               velxgrid + derivrotationmatrix(1, 1)*xxc(1) + &
&               derivrotationmatrix(1, 2)*xxc(2) + derivrotationmatrix(1&
&               , 3)*xxc(3)
              bcdata(mm)%uslip(i, j, 2) = bcdata(mm)%uslip(i, j, 2) + &
&               velygrid + derivrotationmatrix(2, 1)*xxc(1) + &
&               derivrotationmatrix(2, 2)*xxc(2) + derivrotationmatrix(2&
&               , 3)*xxc(3)
              bcdata(mm)%uslip(i, j, 3) = bcdata(mm)%uslip(i, j, 3) + &
&               velzgrid + derivrotationmatrix(3, 1)*xxc(1) + &
&               derivrotationmatrix(3, 2)*xxc(2) + derivrotationmatrix(3&
&               , 3)*xxc(3)
            end do
          end do
        end if
      end do bocoloop2
    end if
  end subroutine slipvelocitiesfinelevel_block
 
!  differentiation of normalvelocities_block in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: *sfacei *sfacej *sfacek *si
!                *sj *sk *(*bcdata.rface)
!   with respect to varying inputs: *sfacei *sfacej *sfacek *si
!                *sj *sk *(*bcdata.rface)
!   rw status of diff variables: *sfacei:incr *sfacej:incr *sfacek:incr
!                *si:incr *sj:incr *sk:incr *(*bcdata.rface):in-out
!   plus diff mem management of: sfacei:in sfacej:in sfacek:in
!                si:in sj:in sk:in bcdata:in *bcdata.rface:in
  subroutine normalvelocities_block_b(sps)
!
!       normalvelocitiesalllevels computes the normal grid
!       velocities of some boundary faces of the moving blocks for
!       spectral mode sps. all grid levels from ground level to the
!       coarsest level are considered.
!
    use constants
    use blockpointers, only : il, jl, kl, addgridvelocities, nbocos, &
&   bcdata, bcdatad, sfacei, sfaceid, sfacej, sfacejd, sfacek, sfacekd, &
&   bcfaceid, si, sid, sj, sjd, sk, skd
    implicit none
!
!      subroutine arguments.
!
    integer(kind=inttype), intent(in) :: sps
!
!      local variables.
!
    integer(kind=inttype) :: mm
    integer(kind=inttype) :: i, j
    real(kind=realtype) :: weight, mult
    real(kind=realtype) :: weightd
    real(kind=realtype), dimension(:, :), pointer :: sface
    real(kind=realtype), dimension(:, :, :), pointer :: ss
    intrinsic associated
    intrinsic sqrt
    real(kind=realtype) :: temp
    real(kind=realtype) :: tempd
    real(kind=realtype) :: temp0
    real(kind=realtype) :: temp1
    real(kind=realtype) :: temp2
    real(kind=realtype) :: tempd0
    integer :: branch
    integer :: ad_from
    integer :: ad_to
    integer :: ad_from0
    integer :: ad_to0
    integer :: ad_from1
    integer :: ad_to1
    integer :: ad_from2
    integer :: ad_to2
    integer :: ad_from3
    integer :: ad_to3
    integer :: ad_from4
    integer :: ad_to4
    integer :: ad_from5
    integer :: ad_to5
    integer :: ad_from6
    integer :: ad_to6
    integer :: ad_from7
    integer :: ad_to7
    integer :: ad_from8
    integer :: ad_to8
    integer :: ad_from9
    integer :: ad_to9
    integer :: ad_from10
    integer :: ad_to10
! check for a moving block. as it is possible that in a
! multidisicplinary environment additional grid velocities
! are set, the test should be done on addgridvelocities
! and not on blockismoving.
    if (addgridvelocities) then
!
!             determine the normal grid velocities of the boundaries.
!             as these values are based on the unit normal. a division
!             by the length of the normal is needed.
!             furthermore the boundary unit normals are per definition
!             outward pointing, while on the imin, jmin and kmin
!             boundaries the face normals are inward pointing. this
!             is taken into account by the factor mult.
!
! loop over the boundary subfaces.
bocoloop:do mm=1,nbocos
! check whether rface is allocated.
        if (associated(bcdata(mm)%rface)) then
! determine the block face on which the subface is
! located and set some variables accordingly.
! the new procedure is less elegant as the previous one.
! but the new stands up to tapenade.
          if (bcfaceid(mm) .eq. imin) then
            call pushreal8(mult)
            mult = -one
            ad_from0 = bcdata(mm)%jcbeg
            do j=ad_from0,bcdata(mm)%jcend
              ad_from = bcdata(mm)%icbeg
              do i=ad_from,bcdata(mm)%icend
! compute the inverse of the length of the normal
! vector and possibly correct for inward pointing.
                call pushreal8(weight)
                weight = sqrt(si(1, i, j, 1)**2 + si(1, i, j, 2)**2 + si&
&                 (1, i, j, 3)**2)
                if (weight .gt. zero) then
                  call pushreal8(weight)
                  weight = mult/weight
                  call pushcontrol1b(0)
                else
                  call pushcontrol1b(1)
                end if
              end do
              call pushinteger4(i - 1)
              call pushinteger4(ad_from)
            end do
            call pushinteger4(j - 1)
            call pushinteger4(ad_from0)
            call pushcontrol3b(7)
          else if (bcfaceid(mm) .eq. imax) then
            call pushreal8(mult)
            mult = one
            ad_from2 = bcdata(mm)%jcbeg
            do j=ad_from2,bcdata(mm)%jcend
              ad_from1 = bcdata(mm)%icbeg
              do i=ad_from1,bcdata(mm)%icend
! compute the inverse of the length of the normal
! vector and possibly correct for inward pointing.
                call pushreal8(weight)
                weight = sqrt(si(il, i, j, 1)**2 + si(il, i, j, 2)**2 + &
&                 si(il, i, j, 3)**2)
                if (weight .gt. zero) then
                  call pushreal8(weight)
                  weight = mult/weight
                  call pushcontrol1b(0)
                else
                  call pushcontrol1b(1)
                end if
              end do
              call pushinteger4(i - 1)
              call pushinteger4(ad_from1)
            end do
            call pushinteger4(j - 1)
            call pushinteger4(ad_from2)
            call pushcontrol3b(6)
          else if (bcfaceid(mm) .eq. jmin) then
            call pushreal8(mult)
            mult = -one
            ad_from4 = bcdata(mm)%jcbeg
            do j=ad_from4,bcdata(mm)%jcend
              ad_from3 = bcdata(mm)%icbeg
              do i=ad_from3,bcdata(mm)%icend
! compute the inverse of the length of the normal
! vector and possibly correct for inward pointing.
                call pushreal8(weight)
                weight = sqrt(sj(i, 1, j, 1)**2 + sj(i, 1, j, 2)**2 + sj&
&                 (i, 1, j, 3)**2)
                if (weight .gt. zero) then
                  call pushreal8(weight)
                  weight = mult/weight
                  call pushcontrol1b(0)
                else
                  call pushcontrol1b(1)
                end if
              end do
              call pushinteger4(i - 1)
              call pushinteger4(ad_from3)
            end do
            call pushinteger4(j - 1)
            call pushinteger4(ad_from4)
            call pushcontrol3b(5)
          else if (bcfaceid(mm) .eq. jmax) then
            call pushreal8(mult)
            mult = one
            ad_from6 = bcdata(mm)%jcbeg
            do j=ad_from6,bcdata(mm)%jcend
              ad_from5 = bcdata(mm)%icbeg
              do i=ad_from5,bcdata(mm)%icend
! compute the inverse of the length of the normal
! vector and possibly correct for inward pointing.
                call pushreal8(weight)
                weight = sqrt(sj(i, jl, j, 1)**2 + sj(i, jl, j, 2)**2 + &
&                 sj(i, jl, j, 3)**2)
                if (weight .gt. zero) then
                  call pushreal8(weight)
                  weight = mult/weight
                  call pushcontrol1b(0)
                else
                  call pushcontrol1b(1)
                end if
              end do
              call pushinteger4(i - 1)
              call pushinteger4(ad_from5)
            end do
            call pushinteger4(j - 1)
            call pushinteger4(ad_from6)
            call pushcontrol3b(4)
          else if (bcfaceid(mm) .eq. kmin) then
            call pushreal8(mult)
            mult = -one
            ad_from8 = bcdata(mm)%jcbeg
            do j=ad_from8,bcdata(mm)%jcend
              ad_from7 = bcdata(mm)%icbeg
              do i=ad_from7,bcdata(mm)%icend
! compute the inverse of the length of the normal
! vector and possibly correct for inward pointing.
                call pushreal8(weight)
                weight = sqrt(sk(i, j, 1, 1)**2 + sk(i, j, 1, 2)**2 + sk&
&                 (i, j, 1, 3)**2)
                if (weight .gt. zero) then
                  call pushreal8(weight)
                  weight = mult/weight
                  call pushcontrol1b(0)
                else
                  call pushcontrol1b(1)
                end if
              end do
              call pushinteger4(i - 1)
              call pushinteger4(ad_from7)
            end do
            call pushinteger4(j - 1)
            call pushinteger4(ad_from8)
            call pushcontrol3b(3)
          else if (bcfaceid(mm) .eq. kmax) then
            call pushreal8(mult)
            mult = one
            ad_from10 = bcdata(mm)%jcbeg
            do j=ad_from10,bcdata(mm)%jcend
              ad_from9 = bcdata(mm)%icbeg
              do i=ad_from9,bcdata(mm)%icend
! compute the inverse of the length of the normal
! vector and possibly correct for inward pointing.
                call pushreal8(weight)
                weight = sqrt(sk(i, j, kl, 1)**2 + sk(i, j, kl, 2)**2 + &
&                 sk(i, j, kl, 3)**2)
                if (weight .gt. zero) then
                  call pushreal8(weight)
                  weight = mult/weight
                  call pushcontrol1b(0)
                else
                  call pushcontrol1b(1)
                end if
              end do
              call pushinteger4(i - 1)
              call pushinteger4(ad_from9)
            end do
            call pushinteger4(j - 1)
            call pushinteger4(ad_from10)
            call pushcontrol3b(2)
          else
            call pushcontrol3b(1)
          end if
        else
          call pushcontrol3b(0)
        end if
      end do bocoloop
      do mm=nbocos,1,-1
        call popcontrol3b(branch)
        if (branch .lt. 4) then
          if (branch .ge. 2) then
            if (branch .eq. 2) then
              call popinteger4(ad_from10)
              call popinteger4(ad_to10)
              do j=ad_to10,ad_from10,-1
                call popinteger4(ad_from9)
                call popinteger4(ad_to9)
                do i=ad_to9,ad_from9,-1
                  weightd = sfacek(i, j, kl)*bcdatad(mm)%rface(i, j)
                  sfacekd(i, j, kl) = sfacekd(i, j, kl) + weight*bcdatad&
&                   (mm)%rface(i, j)
                  bcdatad(mm)%rface(i, j) = 0.0_8
                  call popcontrol1b(branch)
                  if (branch .eq. 0) then
                    call popreal8(weight)
                    weightd = -(mult*weightd/weight**2)
                  end if
                  call popreal8(weight)
                  temp2 = sk(i, j, kl, 3)
                  temp1 = sk(i, j, kl, 2)
                  temp0 = sk(i, j, kl, 1)
                  if (temp0**2 + temp1**2 + temp2**2 .eq. 0.0_8) then
                    tempd = 0.0_8
                  else
                    tempd = weightd/(2.0*sqrt(temp0**2+temp1**2+temp2**2&
&                     ))
                  end if
                  skd(i, j, kl, 1) = skd(i, j, kl, 1) + 2*temp0*tempd
                  skd(i, j, kl, 2) = skd(i, j, kl, 2) + 2*temp1*tempd
                  skd(i, j, kl, 3) = skd(i, j, kl, 3) + 2*temp2*tempd
                end do
              end do
              call popreal8(mult)
            else
              call popinteger4(ad_from8)
              call popinteger4(ad_to8)
              do j=ad_to8,ad_from8,-1
                call popinteger4(ad_from7)
                call popinteger4(ad_to7)
                do i=ad_to7,ad_from7,-1
                  weightd = sfacek(i, j, 1)*bcdatad(mm)%rface(i, j)
                  sfacekd(i, j, 1) = sfacekd(i, j, 1) + weight*bcdatad(&
&                   mm)%rface(i, j)
                  bcdatad(mm)%rface(i, j) = 0.0_8
                  call popcontrol1b(branch)
                  if (branch .eq. 0) then
                    call popreal8(weight)
                    weightd = -(mult*weightd/weight**2)
                  end if
                  call popreal8(weight)
                  if (sk(i, j, 1, 1)**2 + sk(i, j, 1, 2)**2 + sk(i, j, 1&
&                     , 3)**2 .eq. 0.0_8) then
                    tempd0 = 0.0_8
                  else
                    tempd0 = weightd/(2.0*sqrt(sk(i, j, 1, 1)**2+sk(i, j&
&                     , 1, 2)**2+sk(i, j, 1, 3)**2))
                  end if
                  skd(i, j, 1, 1) = skd(i, j, 1, 1) + 2*sk(i, j, 1, 1)*&
&                   tempd0
                  skd(i, j, 1, 2) = skd(i, j, 1, 2) + 2*sk(i, j, 1, 2)*&
&                   tempd0
                  skd(i, j, 1, 3) = skd(i, j, 1, 3) + 2*sk(i, j, 1, 3)*&
&                   tempd0
                end do
              end do
              call popreal8(mult)
            end if
          end if
        else if (branch .lt. 6) then
          if (branch .eq. 4) then
            call popinteger4(ad_from6)
            call popinteger4(ad_to6)
            do j=ad_to6,ad_from6,-1
              call popinteger4(ad_from5)
              call popinteger4(ad_to5)
              do i=ad_to5,ad_from5,-1
                weightd = sfacej(i, jl, j)*bcdatad(mm)%rface(i, j)
                sfacejd(i, jl, j) = sfacejd(i, jl, j) + weight*bcdatad(&
&                 mm)%rface(i, j)
                bcdatad(mm)%rface(i, j) = 0.0_8
                call popcontrol1b(branch)
                if (branch .eq. 0) then
                  call popreal8(weight)
                  weightd = -(mult*weightd/weight**2)
                end if
                call popreal8(weight)
                temp2 = sj(i, jl, j, 3)
                temp1 = sj(i, jl, j, 2)
                temp0 = sj(i, jl, j, 1)
                if (temp0**2 + temp1**2 + temp2**2 .eq. 0.0_8) then
                  tempd = 0.0_8
                else
                  tempd = weightd/(2.0*sqrt(temp0**2+temp1**2+temp2**2))
                end if
                sjd(i, jl, j, 1) = sjd(i, jl, j, 1) + 2*temp0*tempd
                sjd(i, jl, j, 2) = sjd(i, jl, j, 2) + 2*temp1*tempd
                sjd(i, jl, j, 3) = sjd(i, jl, j, 3) + 2*temp2*tempd
              end do
            end do
            call popreal8(mult)
          else
            call popinteger4(ad_from4)
            call popinteger4(ad_to4)
            do j=ad_to4,ad_from4,-1
              call popinteger4(ad_from3)
              call popinteger4(ad_to3)
              do i=ad_to3,ad_from3,-1
                weightd = sfacej(i, 1, j)*bcdatad(mm)%rface(i, j)
                sfacejd(i, 1, j) = sfacejd(i, 1, j) + weight*bcdatad(mm)&
&                 %rface(i, j)
                bcdatad(mm)%rface(i, j) = 0.0_8
                call popcontrol1b(branch)
                if (branch .eq. 0) then
                  call popreal8(weight)
                  weightd = -(mult*weightd/weight**2)
                end if
                call popreal8(weight)
                if (sj(i, 1, j, 1)**2 + sj(i, 1, j, 2)**2 + sj(i, 1, j, &
&                   3)**2 .eq. 0.0_8) then
                  tempd0 = 0.0_8
                else
                  tempd0 = weightd/(2.0*sqrt(sj(i, 1, j, 1)**2+sj(i, 1, &
&                   j, 2)**2+sj(i, 1, j, 3)**2))
                end if
                sjd(i, 1, j, 1) = sjd(i, 1, j, 1) + 2*sj(i, 1, j, 1)*&
&                 tempd0
                sjd(i, 1, j, 2) = sjd(i, 1, j, 2) + 2*sj(i, 1, j, 2)*&
&                 tempd0
                sjd(i, 1, j, 3) = sjd(i, 1, j, 3) + 2*sj(i, 1, j, 3)*&
&                 tempd0
              end do
            end do
            call popreal8(mult)
          end if
        else if (branch .eq. 6) then
          call popinteger4(ad_from2)
          call popinteger4(ad_to2)
          do j=ad_to2,ad_from2,-1
            call popinteger4(ad_from1)
            call popinteger4(ad_to1)
            do i=ad_to1,ad_from1,-1
              weightd = sfacei(il, i, j)*bcdatad(mm)%rface(i, j)
              sfaceid(il, i, j) = sfaceid(il, i, j) + weight*bcdatad(mm)&
&               %rface(i, j)
              bcdatad(mm)%rface(i, j) = 0.0_8
              call popcontrol1b(branch)
              if (branch .eq. 0) then
                call popreal8(weight)
                weightd = -(mult*weightd/weight**2)
              end if
              call popreal8(weight)
              temp = si(il, i, j, 3)
              temp0 = si(il, i, j, 2)
              temp1 = si(il, i, j, 1)
              if (temp1**2 + temp0**2 + temp**2 .eq. 0.0_8) then
                tempd0 = 0.0_8
              else
                tempd0 = weightd/(2.0*sqrt(temp1**2+temp0**2+temp**2))
              end if
              sid(il, i, j, 1) = sid(il, i, j, 1) + 2*temp1*tempd0
              sid(il, i, j, 2) = sid(il, i, j, 2) + 2*temp0*tempd0
              sid(il, i, j, 3) = sid(il, i, j, 3) + 2*temp*tempd0
            end do
          end do
          call popreal8(mult)
        else
          call popinteger4(ad_from0)
          call popinteger4(ad_to0)
          do j=ad_to0,ad_from0,-1
            call popinteger4(ad_from)
            call popinteger4(ad_to)
            do i=ad_to,ad_from,-1
              weightd = sfacei(1, i, j)*bcdatad(mm)%rface(i, j)
              sfaceid(1, i, j) = sfaceid(1, i, j) + weight*bcdatad(mm)%&
&               rface(i, j)
              bcdatad(mm)%rface(i, j) = 0.0_8
              call popcontrol1b(branch)
              if (branch .eq. 0) then
                call popreal8(weight)
                weightd = -(mult*weightd/weight**2)
              end if
              call popreal8(weight)
              if (si(1, i, j, 1)**2 + si(1, i, j, 2)**2 + si(1, i, j, 3)&
&                 **2 .eq. 0.0_8) then
                tempd = 0.0_8
              else
                tempd = weightd/(2.0*sqrt(si(1, i, j, 1)**2+si(1, i, j, &
&                 2)**2+si(1, i, j, 3)**2))
              end if
              sid(1, i, j, 1) = sid(1, i, j, 1) + 2*si(1, i, j, 1)*tempd
              sid(1, i, j, 2) = sid(1, i, j, 2) + 2*si(1, i, j, 2)*tempd
              sid(1, i, j, 3) = sid(1, i, j, 3) + 2*si(1, i, j, 3)*tempd
            end do
          end do
          call popreal8(mult)
        end if
      end do
    else
! block is not moving. loop over the boundary faces and set
! the normal grid velocity to zero if allocated.
      do mm=1,nbocos
        if (associated(bcdata(mm)%rface)) then
          call pushcontrol1b(1)
        else
          call pushcontrol1b(0)
        end if
      end do
      do mm=nbocos,1,-1
        call popcontrol1b(branch)
        if (branch .ne. 0) bcdatad(mm)%rface = 0.0_8
      end do
    end if
  end subroutine normalvelocities_block_b
 
  subroutine normalvelocities_block(sps)
!
!       normalvelocitiesalllevels computes the normal grid
!       velocities of some boundary faces of the moving blocks for
!       spectral mode sps. all grid levels from ground level to the
!       coarsest level are considered.
!
    use constants
    use blockpointers, only : il, jl, kl, addgridvelocities, nbocos, &
&   bcdata, sfacei, sfacej, sfacek, bcfaceid, si, sj, sk
    implicit none
!
!      subroutine arguments.
!
    integer(kind=inttype), intent(in) :: sps
!
!      local variables.
!
    integer(kind=inttype) :: mm
    integer(kind=inttype) :: i, j
    real(kind=realtype) :: weight, mult
    real(kind=realtype), dimension(:, :), pointer :: sface
    real(kind=realtype), dimension(:, :, :), pointer :: ss
    intrinsic associated
    intrinsic sqrt
! check for a moving block. as it is possible that in a
! multidisicplinary environment additional grid velocities
! are set, the test should be done on addgridvelocities
! and not on blockismoving.
    if (addgridvelocities) then
!
!             determine the normal grid velocities of the boundaries.
!             as these values are based on the unit normal. a division
!             by the length of the normal is needed.
!             furthermore the boundary unit normals are per definition
!             outward pointing, while on the imin, jmin and kmin
!             boundaries the face normals are inward pointing. this
!             is taken into account by the factor mult.
!
! loop over the boundary subfaces.
bocoloop:do mm=1,nbocos
! check whether rface is allocated.
        if (associated(bcdata(mm)%rface)) then
! determine the block face on which the subface is
! located and set some variables accordingly.
! the new procedure is less elegant as the previous one.
! but the new stands up to tapenade.
          if (bcfaceid(mm) .eq. imin) then
            mult = -one
            do j=bcdata(mm)%jcbeg,bcdata(mm)%jcend
              do i=bcdata(mm)%icbeg,bcdata(mm)%icend
! compute the inverse of the length of the normal
! vector and possibly correct for inward pointing.
                weight = sqrt(si(1, i, j, 1)**2 + si(1, i, j, 2)**2 + si&
&                 (1, i, j, 3)**2)
                if (weight .gt. zero) weight = mult/weight
! compute the normal velocity based on the outward
! pointing unit normal.
                bcdata(mm)%rface(i, j) = weight*sfacei(1, i, j)
              end do
            end do
          else if (bcfaceid(mm) .eq. imax) then
            mult = one
            do j=bcdata(mm)%jcbeg,bcdata(mm)%jcend
              do i=bcdata(mm)%icbeg,bcdata(mm)%icend
! compute the inverse of the length of the normal
! vector and possibly correct for inward pointing.
                weight = sqrt(si(il, i, j, 1)**2 + si(il, i, j, 2)**2 + &
&                 si(il, i, j, 3)**2)
                if (weight .gt. zero) weight = mult/weight
! compute the normal velocity based on the outward
! pointing unit normal.
                bcdata(mm)%rface(i, j) = weight*sfacei(il, i, j)
              end do
            end do
          else if (bcfaceid(mm) .eq. jmin) then
            mult = -one
            do j=bcdata(mm)%jcbeg,bcdata(mm)%jcend
              do i=bcdata(mm)%icbeg,bcdata(mm)%icend
! compute the inverse of the length of the normal
! vector and possibly correct for inward pointing.
                weight = sqrt(sj(i, 1, j, 1)**2 + sj(i, 1, j, 2)**2 + sj&
&                 (i, 1, j, 3)**2)
                if (weight .gt. zero) weight = mult/weight
! compute the normal velocity based on the outward
! pointing unit normal.
                bcdata(mm)%rface(i, j) = weight*sfacej(i, 1, j)
              end do
            end do
          else if (bcfaceid(mm) .eq. jmax) then
            mult = one
            do j=bcdata(mm)%jcbeg,bcdata(mm)%jcend
              do i=bcdata(mm)%icbeg,bcdata(mm)%icend
! compute the inverse of the length of the normal
! vector and possibly correct for inward pointing.
                weight = sqrt(sj(i, jl, j, 1)**2 + sj(i, jl, j, 2)**2 + &
&                 sj(i, jl, j, 3)**2)
                if (weight .gt. zero) weight = mult/weight
! compute the normal velocity based on the outward
! pointing unit normal.
                bcdata(mm)%rface(i, j) = weight*sfacej(i, jl, j)
              end do
            end do
          else if (bcfaceid(mm) .eq. kmin) then
            mult = -one
            do j=bcdata(mm)%jcbeg,bcdata(mm)%jcend
              do i=bcdata(mm)%icbeg,bcdata(mm)%icend
! compute the inverse of the length of the normal
! vector and possibly correct for inward pointing.
                weight = sqrt(sk(i, j, 1, 1)**2 + sk(i, j, 1, 2)**2 + sk&
&                 (i, j, 1, 3)**2)
                if (weight .gt. zero) weight = mult/weight
! compute the normal velocity based on the outward
! pointing unit normal.
                bcdata(mm)%rface(i, j) = weight*sfacek(i, j, 1)
              end do
            end do
          else if (bcfaceid(mm) .eq. kmax) then
            mult = one
            do j=bcdata(mm)%jcbeg,bcdata(mm)%jcend
              do i=bcdata(mm)%icbeg,bcdata(mm)%icend
! compute the inverse of the length of the normal
! vector and possibly correct for inward pointing.
                weight = sqrt(sk(i, j, kl, 1)**2 + sk(i, j, kl, 2)**2 + &
&                 sk(i, j, kl, 3)**2)
                if (weight .gt. zero) weight = mult/weight
! compute the normal velocity based on the outward
! pointing unit normal.
                bcdata(mm)%rface(i, j) = weight*sfacek(i, j, kl)
              end do
            end do
          end if
        end if
      end do bocoloop
    else
! block is not moving. loop over the boundary faces and set
! the normal grid velocity to zero if allocated.
      do mm=1,nbocos
        if (associated(bcdata(mm)%rface)) bcdata(mm)%rface = zero
      end do
    end if
  end subroutine normalvelocities_block
! ----------------------------------------------------------------------
!                                                                      |
!                    no tapenade routine below this line               |
!                                                                      |
! ----------------------------------------------------------------------
 
end module solverutils_b