solverUtils_d.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_d
  implicit none
 
contains
!  differentiation of timestep_block in forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: *dtl *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:in pinfcorr:in *rev:in
!                *dtl:in-out *p:in *sfacei:in *sfacej:in *sfacek:in
!                *w:in *rlv:in *vol:in *si:in *sj:in *sk:in *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_d(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_d, 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 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) :: arg1
    real(kind=realtype) :: arg1d
    real(kind=realtype) :: result1
    real(kind=realtype) :: result1d
    real(kind=realtype) :: temp
    real(kind=realtype) :: temp0
! 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.
      plimd = 0.001_realtype*pinfcorrd
      plim = 0.001_realtype*pinfcorr
      rlim = 0.001_realtype*rhoinf
      clim2d = gammainf*0.000001_realtype*(pinfcorrd-pinfcorr*rhoinfd/&
&       rhoinf)/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) 
        sfaced = 0.0_8
! no preconditioner. simply the standard spectral radius.
! loop over the cells, including the first level halo.
        do k=1,ke
          do j=1,je
            do i=1,ie
! compute the velocities and speed of sound squared.
              uuxd = wd(i, j, k, ivx)
              uux = w(i, j, k, ivx)
              uuyd = wd(i, j, k, ivy)
              uuy = w(i, j, k, ivy)
              uuzd = wd(i, j, k, ivz)
              uuz = w(i, j, k, ivz)
              temp = w(i, j, k, irho)
              temp0 = p(i, j, k)/temp
              cc2d = gamma(i, j, k)*(pd(i, j, k)-temp0*wd(i, j, k, irho)&
&               )/temp
              cc2 = gamma(i, j, k)*temp0
              if (cc2 .lt. clim2) then
                cc2d = clim2d
                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) then
                sfaced = sfaceid(i-1, j, k) + sfaceid(i, j, k)
                sface = sfacei(i-1, j, k) + sfacei(i, j, k)
              end if
! spectral radius in i-direction.
              sxd = sid(i-1, j, k, 1) + sid(i, j, k, 1)
              sx = si(i-1, j, k, 1) + si(i, j, k, 1)
              syd = sid(i-1, j, k, 2) + sid(i, j, k, 2)
              sy = si(i-1, j, k, 2) + si(i, j, k, 2)
              szd = sid(i-1, j, k, 3) + sid(i, j, k, 3)
              sz = si(i-1, j, k, 3) + si(i, j, k, 3)
              qsid = sx*uuxd + uux*sxd + sy*uuyd + uuy*syd + sz*uuzd + &
&               uuz*szd - sfaced
              qsi = uux*sx + uuy*sy + uuz*sz - sface
              if (qsi .ge. 0.) then
                abs0d = qsid
                abs0 = qsi
              else
                abs0d = -qsid
                abs0 = -qsi
              end if
              temp0 = sx*sx + sy*sy + sz*sz
              arg1d = temp0*cc2d + cc2*(2*sx*sxd+2*sy*syd+2*sz*szd)
              arg1 = cc2*temp0
              temp0 = sqrt(arg1)
              if (arg1 .eq. 0.0_8) then
                result1d = 0.0_8
              else
                result1d = arg1d/(2.0*temp0)
              end if
              result1 = temp0
              rid = half*(abs0d+acousticscalefactor*result1d)
              ri = half*(abs0+acousticscalefactor*result1)
! the grid velocity in j-direction.
              if (addgridvelocities) then
                sfaced = sfacejd(i, j-1, k) + sfacejd(i, j, k)
                sface = sfacej(i, j-1, k) + sfacej(i, j, k)
              end if
! spectral radius in j-direction.
              sxd = sjd(i, j-1, k, 1) + sjd(i, j, k, 1)
              sx = sj(i, j-1, k, 1) + sj(i, j, k, 1)
              syd = sjd(i, j-1, k, 2) + sjd(i, j, k, 2)
              sy = sj(i, j-1, k, 2) + sj(i, j, k, 2)
              szd = sjd(i, j-1, k, 3) + sjd(i, j, k, 3)
              sz = sj(i, j-1, k, 3) + sj(i, j, k, 3)
              qsjd = sx*uuxd + uux*sxd + sy*uuyd + uuy*syd + sz*uuzd + &
&               uuz*szd - sfaced
              qsj = uux*sx + uuy*sy + uuz*sz - sface
              if (qsj .ge. 0.) then
                abs1d = qsjd
                abs1 = qsj
              else
                abs1d = -qsjd
                abs1 = -qsj
              end if
              temp0 = sx*sx + sy*sy + sz*sz
              arg1d = temp0*cc2d + cc2*(2*sx*sxd+2*sy*syd+2*sz*szd)
              arg1 = cc2*temp0
              temp0 = sqrt(arg1)
              if (arg1 .eq. 0.0_8) then
                result1d = 0.0_8
              else
                result1d = arg1d/(2.0*temp0)
              end if
              result1 = temp0
              rjd = half*(abs1d+acousticscalefactor*result1d)
              rj = half*(abs1+acousticscalefactor*result1)
! the grid velocity in k-direction.
              if (addgridvelocities) then
                sfaced = sfacekd(i, j, k-1) + sfacekd(i, j, k)
                sface = sfacek(i, j, k-1) + sfacek(i, j, k)
              end if
! spectral radius in k-direction.
              sxd = skd(i, j, k-1, 1) + skd(i, j, k, 1)
              sx = sk(i, j, k-1, 1) + sk(i, j, k, 1)
              syd = skd(i, j, k-1, 2) + skd(i, j, k, 2)
              sy = sk(i, j, k-1, 2) + sk(i, j, k, 2)
              szd = skd(i, j, k-1, 3) + skd(i, j, k, 3)
              sz = sk(i, j, k-1, 3) + sk(i, j, k, 3)
              qskd = sx*uuxd + uux*sxd + sy*uuyd + uuy*syd + sz*uuzd + &
&               uuz*szd - sfaced
              qsk = uux*sx + uuy*sy + uuz*sz - sface
              if (qsk .ge. 0.) then
                abs2d = qskd
                abs2 = qsk
              else
                abs2d = -qskd
                abs2 = -qsk
              end if
              temp0 = sx*sx + sy*sy + sz*sz
              arg1d = temp0*cc2d + cc2*(2*sx*sxd+2*sy*syd+2*sz*szd)
              arg1 = cc2*temp0
              temp0 = sqrt(arg1)
              if (arg1 .eq. 0.0_8) then
                result1d = 0.0_8
              else
                result1d = arg1d/(2.0*temp0)
              end if
              result1 = temp0
              rkd = half*(abs2d+acousticscalefactor*result1d)
              rk = half*(abs2+acousticscalefactor*result1)
! compute the inviscid contribution to the time step.
              if (.not.onlyradii) then
                dtld(i, j, k) = rid + rjd + rkd
                dtl(i, j, k) = ri + rj + rk
              end if
!
!           adapt the spectral radii if directional scaling must be
!           applied.
!
              if (doscaling) then
                if (ri .lt. eps) then
                  ri = eps
                  rid = 0.0_8
                else
                  ri = ri
                end if
                if (rj .lt. eps) then
                  rj = eps
                  rjd = 0.0_8
                else
                  rj = rj
                end if
                if (rk .lt. eps) then
                  rk = eps
                  rkd = 0.0_8
                else
                  rk = rk
                end if
! compute the scaling in the three coordinate
! directions.
                temp0 = ri/rj
                if (temp0 .le. 0.0_8 .and. (adis .eq. 0.0_8 .or. adis &
&                   .ne. int(adis))) then
                  rijd = 0.0_8
                else
                  rijd = adis*temp0**(adis-1)*(rid-temp0*rjd)/rj
                end if
                rij = temp0**adis
                temp0 = rj/rk
                if (temp0 .le. 0.0_8 .and. (adis .eq. 0.0_8 .or. adis &
&                   .ne. int(adis))) then
                  rjkd = 0.0_8
                else
                  rjkd = adis*temp0**(adis-1)*(rjd-temp0*rkd)/rk
                end if
                rjk = temp0**adis
                temp0 = rk/ri
                if (temp0 .le. 0.0_8 .and. (adis .eq. 0.0_8 .or. adis &
&                   .ne. int(adis))) then
                  rkid = 0.0_8
                else
                  rkid = adis*temp0**(adis-1)*(rkd-temp0*rid)/ri
                end if
                rki = temp0**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.
                temp0 = one/rij
                radid(i, j, k) = (one+temp0+rki)*rid + ri*(rkid-temp0*&
&                 rijd/rij)
                radi(i, j, k) = ri*(one+temp0+rki)
                temp0 = one/rjk
                radjd(i, j, k) = (one+temp0+rij)*rjd + rj*(rijd-temp0*&
&                 rjkd/rjk)
                radj(i, j, k) = rj*(one+temp0+rij)
                temp0 = one/rki
                radkd(i, j, k) = (one+temp0+rjk)*rkd + rk*(rjkd-temp0*&
&                 rkid/rki)
                radk(i, j, k) = rk*(one+temp0+rjk)
              else
                radid(i, j, k) = rid
                radi(i, j, k) = ri
                radjd(i, j, k) = rjd
                radj(i, j, k) = rj
                radkd(i, j, k) = rkd
                radk(i, j, k) = rk
              end if
            end do
          end do
        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.
                rmud = rlvd(i, j, k)
                rmu = rlv(i, j, k)
                if (eddymodel) then
                  rmud = rmud + revd(i, j, k)
                  rmu = rmu + rev(i, j, k)
                end if
                temp0 = w(i, j, k, irho)
                temp = temp0*vol(i, j, k)
                rmud = half*(rmud-rmu*(vol(i, j, k)*wd(i, j, k, irho)+&
&                 temp0*vold(i, j, k))/temp)/temp
                rmu = half*(rmu/temp)
! add the viscous contribution in i-direction to the
! (inverse) of the time step.
                sxd = sid(i, j, k, 1) + sid(i-1, j, k, 1)
                sx = si(i, j, k, 1) + si(i-1, j, k, 1)
                syd = sid(i, j, k, 2) + sid(i-1, j, k, 2)
                sy = si(i, j, k, 2) + si(i-1, j, k, 2)
                szd = sid(i, j, k, 3) + sid(i-1, j, k, 3)
                sz = si(i, j, k, 3) + si(i-1, j, k, 3)
                temp0 = sx*sx + sy*sy + sz*sz
                vsid = temp0*rmud + rmu*(2*sx*sxd+2*sy*syd+2*sz*szd)
                vsi = rmu*temp0
                dtld(i, j, k) = dtld(i, j, k) + vsid
                dtl(i, j, k) = dtl(i, j, k) + vsi
! add the viscous contribution in j-direction to the
! (inverse) of the time step.
                sxd = sjd(i, j, k, 1) + sjd(i, j-1, k, 1)
                sx = sj(i, j, k, 1) + sj(i, j-1, k, 1)
                syd = sjd(i, j, k, 2) + sjd(i, j-1, k, 2)
                sy = sj(i, j, k, 2) + sj(i, j-1, k, 2)
                szd = sjd(i, j, k, 3) + sjd(i, j-1, k, 3)
                sz = sj(i, j, k, 3) + sj(i, j-1, k, 3)
                temp0 = sx*sx + sy*sy + sz*sz
                vsjd = temp0*rmud + rmu*(2*sx*sxd+2*sy*syd+2*sz*szd)
                vsj = rmu*temp0
                dtld(i, j, k) = dtld(i, j, k) + vsjd
                dtl(i, j, k) = dtl(i, j, k) + vsj
! add the viscous contribution in k-direction to the
! (inverse) of the time step.
                sxd = skd(i, j, k, 1) + skd(i, j, k-1, 1)
                sx = sk(i, j, k, 1) + sk(i, j, k-1, 1)
                syd = skd(i, j, k, 2) + skd(i, j, k-1, 2)
                sy = sk(i, j, k, 2) + sk(i, j, k-1, 2)
                szd = skd(i, j, k, 3) + skd(i, j, k-1, 3)
                sz = sk(i, j, k, 3) + sk(i, j, k-1, 3)
                temp0 = sx*sx + sy*sy + sz*sz
                vskd = temp0*rmud + rmu*(2*sx*sxd+2*sy*syd+2*sz*szd)
                vsk = rmu*temp0
                dtld(i, j, k) = dtld(i, j, k) + vskd
                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
                dtld(i, j, k) = dtld(i, j, k) + tmp*vold(i, j, k)
                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
                abs3d = pd(i+1, j, k) - two*pd(i, j, k) + pd(i-1, j, k)
                abs3 = p(i+1, j, k) - two*p(i, j, k) + p(i-1, j, k)
              else
                abs3d = two*pd(i, j, k) - pd(i+1, j, k) - pd(i-1, j, k)
                abs3 = -(p(i+1, j, k)-two*p(i, j, k)+p(i-1, j, k))
              end if
              temp0 = p(i+1, j, k) + two*p(i, j, k) + p(i-1, j, k) + &
&               plim
              dpid = (abs3d-abs3*(pd(i+1, j, k)+two*pd(i, j, k)+pd(i-1, &
&               j, k)+plimd)/temp0)/temp0
              dpi = abs3/temp0
              if (p(i, j+1, k) - two*p(i, j, k) + p(i, j-1, k) .ge. 0.) &
&             then
                abs4d = pd(i, j+1, k) - two*pd(i, j, k) + pd(i, j-1, k)
                abs4 = p(i, j+1, k) - two*p(i, j, k) + p(i, j-1, k)
              else
                abs4d = two*pd(i, j, k) - pd(i, j+1, k) - pd(i, j-1, k)
                abs4 = -(p(i, j+1, k)-two*p(i, j, k)+p(i, j-1, k))
              end if
              temp0 = p(i, j+1, k) + two*p(i, j, k) + p(i, j-1, k) + &
&               plim
              dpjd = (abs4d-abs4*(pd(i, j+1, k)+two*pd(i, j, k)+pd(i, j-&
&               1, k)+plimd)/temp0)/temp0
              dpj = abs4/temp0
              if (p(i, j, k+1) - two*p(i, j, k) + p(i, j, k-1) .ge. 0.) &
&             then
                abs5d = pd(i, j, k+1) - two*pd(i, j, k) + pd(i, j, k-1)
                abs5 = p(i, j, k+1) - two*p(i, j, k) + p(i, j, k-1)
              else
                abs5d = two*pd(i, j, k) - pd(i, j, k+1) - pd(i, j, k-1)
                abs5 = -(p(i, j, k+1)-two*p(i, j, k)+p(i, j, k-1))
              end if
              temp0 = p(i, j, k+1) + two*p(i, j, k) + p(i, j, k-1) + &
&               plim
              dpkd = (abs5d-abs5*(pd(i, j, k+1)+two*pd(i, j, k)+pd(i, j&
&               , k-1)+plimd)/temp0)/temp0
              dpk = abs5/temp0
              temp0 = one/(one+b*(dpi+dpj+dpk))
              rfld = -(temp0*b*(dpid+dpjd+dpkd)/(one+b*(dpi+dpj+dpk)))
              rfl = temp0
              temp0 = rfl/dtl(i, j, k)
              dtld(i, j, k) = (rfld-temp0*dtld(i, j, k))/dtl(i, j, k)
              dtl(i, j, k) = temp0
            end do
          end do
        end do
      end if
    end if
  end subroutine timestep_block_d
 
  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_d, 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 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
    real(kind=realtype) :: arg1
    real(kind=realtype) :: result1
! 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) 
! no preconditioner. simply the standard spectral radius.
! loop over the cells, including the first level halo.
        do k=1,ke
          do j=1,je
            do i=1,ie
! 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
              arg1 = cc2*(sx**2+sy**2+sz**2)
              result1 = sqrt(arg1)
              ri = half*(abs0+acousticscalefactor*result1)
! 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
              arg1 = cc2*(sx**2+sy**2+sz**2)
              result1 = sqrt(arg1)
              rj = half*(abs1+acousticscalefactor*result1)
! 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
              arg1 = cc2*(sx**2+sy**2+sz**2)
              result1 = sqrt(arg1)
              rk = half*(abs2+acousticscalefactor*result1)
! 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
          end do
        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 forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: *sfacei *sfacej *sfacek *s
!   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:in timeref:in rhoinf:in *(flowdoms.x):in
!                veldirfreestream:in machgrid:in *sfacei:in-out
!                *sfacej:in-out *sfacek:in-out *s:in-out *si:in
!                *sj:in *sk:in
!   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_d(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_d, only : derivativerotmatrixrigid, &
&   derivativerotmatrixrigid_d, getdirvector
    use utils_d, 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) :: arg1
    real(kind=realtype) :: arg1d
    real(kind=realtype) :: temp
    real(kind=realtype) :: temp0
    real(kind=realtype) :: temp1
! compute the mesh velocity from the given mesh mach number.
! vel{x,y,z}grid0 is the actual velocity you want at the
! geometry.
    arg1d = gammainf*(pinfd-pinf*rhoinfd/rhoinf)/rhoinf
    arg1 = gammainf*pinf/rhoinf
    temp = sqrt(arg1)
    if (arg1 .eq. 0.0_8) then
      ainfd = 0.0_8
    else
      ainfd = arg1d/(2.0*temp)
    end if
    ainf = temp
    velxgrid0d = -(veldirfreestream(1)*(machgrid*ainfd+ainf*machgridd)+&
&     ainf*machgrid*veldirfreestreamd(1))
    velxgrid0 = ainf*machgrid*(-veldirfreestream(1))
    velygrid0d = -(veldirfreestream(2)*(machgrid*ainfd+ainf*machgridd)+&
&     ainf*machgrid*veldirfreestreamd(2))
    velygrid0 = ainf*machgrid*(-veldirfreestream(2))
    velzgrid0d = -(veldirfreestream(3)*(machgrid*ainfd+ainf*machgridd)+&
&     ainf*machgrid*veldirfreestreamd(3))
    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.
    derivrotationmatrixd = 0.0_8
    call derivativerotmatrixrigid_d(derivrotationmatrix, &
&                             derivrotationmatrixd, 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
        rotrated = cgnsdoms(j)%rotrate*timerefd
        rotrate = timeref*cgnsdoms(j)%rotrate
        velxgridd = velxgrid0d
        velxgrid = velxgrid0
        velygridd = velygrid0d
        velygrid = velygrid0
        velzgridd = velzgrid0d
        velzgrid = velzgrid0
        xcd = 0.0_8
        xxcd = 0.0_8
        scd = 0.0_8
!
!             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.
              xcd(1) = eighth*(flowdomsd(nn, groundlevel, sps)%x(i-1, j-&
&               1, k-1, 1)+flowdomsd(nn, groundlevel, sps)%x(i, j-1, k-1&
&               , 1)+flowdomsd(nn, groundlevel, sps)%x(i-1, j, k-1, 1)+&
&               flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 1)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, k, 1)+&
&               flowdomsd(nn, groundlevel, sps)%x(i, j-1, k, 1)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, k, 1)+&
&               flowdomsd(nn, groundlevel, sps)%x(i, j, k, 1))
              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))
              xcd(2) = eighth*(flowdomsd(nn, groundlevel, sps)%x(i-1, j-&
&               1, k-1, 2)+flowdomsd(nn, groundlevel, sps)%x(i, j-1, k-1&
&               , 2)+flowdomsd(nn, groundlevel, sps)%x(i-1, j, k-1, 2)+&
&               flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 2)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, k, 2)+&
&               flowdomsd(nn, groundlevel, sps)%x(i, j-1, k, 2)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, k, 2)+&
&               flowdomsd(nn, groundlevel, sps)%x(i, j, k, 2))
              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))
              xcd(3) = eighth*(flowdomsd(nn, groundlevel, sps)%x(i-1, j-&
&               1, k-1, 3)+flowdomsd(nn, groundlevel, sps)%x(i, j-1, k-1&
&               , 3)+flowdomsd(nn, groundlevel, sps)%x(i-1, j, k-1, 3)+&
&               flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 3)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, k, 3)+&
&               flowdomsd(nn, groundlevel, sps)%x(i, j-1, k, 3)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, k, 3)+&
&               flowdomsd(nn, groundlevel, sps)%x(i, j, k, 3))
              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.
              xxcd(1) = xcd(1)
              xxc(1) = xc(1) - rotcenter(1)
              xxcd(2) = xcd(2)
              xxc(2) = xc(2) - rotcenter(2)
              xxcd(3) = xcd(3)
              xxc(3) = xc(3) - rotcenter(3)
! determine the rotation speed of the cell center,
! which is omega*r.
              scd(1) = xxc(3)*rotrated(2) + rotrate(2)*xxcd(3) - xxc(2)*&
&               rotrated(3) - rotrate(3)*xxcd(2)
              sc(1) = rotrate(2)*xxc(3) - rotrate(3)*xxc(2)
              scd(2) = xxc(1)*rotrated(3) + rotrate(3)*xxcd(1) - xxc(3)*&
&               rotrated(1) - rotrate(1)*xxcd(3)
              sc(2) = rotrate(3)*xxc(1) - rotrate(1)*xxc(3)
              scd(3) = xxc(2)*rotrated(1) + rotrate(1)*xxcd(2) - xxc(1)*&
&               rotrated(2) - rotrate(2)*xxcd(1)
              sc(3) = rotrate(1)*xxc(2) - rotrate(2)*xxc(1)
! determine the coordinates relative to the
! rigid body rotation point.
              xxcd(1) = xcd(1)
              xxc(1) = xc(1) - rotationpoint(1)
              xxcd(2) = xcd(2)
              xxc(2) = xc(2) - rotationpoint(2)
              xxcd(3) = xcd(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.
              sd(i, j, k, 1) = scd(1) + velxgridd + xxc(1)*&
&               derivrotationmatrixd(1, 1) + derivrotationmatrix(1, 1)*&
&               xxcd(1) + xxc(2)*derivrotationmatrixd(1, 2) + &
&               derivrotationmatrix(1, 2)*xxcd(2) + xxc(3)*&
&               derivrotationmatrixd(1, 3) + derivrotationmatrix(1, 3)*&
&               xxcd(3)
              s(i, j, k, 1) = sc(1) + velxgrid + derivrotationmatrix(1, &
&               1)*xxc(1) + derivrotationmatrix(1, 2)*xxc(2) + &
&               derivrotationmatrix(1, 3)*xxc(3)
              sd(i, j, k, 2) = scd(2) + velygridd + xxc(1)*&
&               derivrotationmatrixd(2, 1) + derivrotationmatrix(2, 1)*&
&               xxcd(1) + xxc(2)*derivrotationmatrixd(2, 2) + &
&               derivrotationmatrix(2, 2)*xxcd(2) + xxc(3)*&
&               derivrotationmatrixd(2, 3) + derivrotationmatrix(2, 3)*&
&               xxcd(3)
              s(i, j, k, 2) = sc(2) + velygrid + derivrotationmatrix(2, &
&               1)*xxc(1) + derivrotationmatrix(2, 2)*xxc(2) + &
&               derivrotationmatrix(2, 3)*xxc(3)
              sd(i, j, k, 3) = scd(3) + velzgridd + xxc(1)*&
&               derivrotationmatrixd(3, 1) + derivrotationmatrix(3, 1)*&
&               xxcd(1) + xxc(2)*derivrotationmatrixd(3, 2) + &
&               derivrotationmatrix(3, 2)*xxcd(2) + xxc(3)*&
&               derivrotationmatrixd(3, 3) + derivrotationmatrix(3, 3)*&
&               xxcd(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.
              xcd(1) = fourth*(flowdomsd(nn, groundlevel, sps)%x(i, j-1&
&               , k-1, 1)+flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 1&
&               )+flowdomsd(nn, groundlevel, sps)%x(i, j-1, k, 1)+&
&               flowdomsd(nn, groundlevel, sps)%x(i, j, k, 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))
              xcd(2) = fourth*(flowdomsd(nn, groundlevel, sps)%x(i, j-1&
&               , k-1, 2)+flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 2&
&               )+flowdomsd(nn, groundlevel, sps)%x(i, j-1, k, 2)+&
&               flowdomsd(nn, groundlevel, sps)%x(i, j, k, 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))
              xcd(3) = fourth*(flowdomsd(nn, groundlevel, sps)%x(i, j-1&
&               , k-1, 3)+flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 3&
&               )+flowdomsd(nn, groundlevel, sps)%x(i, j-1, k, 3)+&
&               flowdomsd(nn, groundlevel, sps)%x(i, j, k, 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))
              rotcenterd = 0.0_8
              call cellfacevelocities_d(xc, xcd, rotcenter, rotcenterd, &
&                                 rotrate, rotrated, velxgrid, velxgridd&
&                                 , velygrid, velygridd, velzgrid, &
&                                 velzgridd, derivrotationmatrix, &
&                                 derivrotationmatrixd, sc, scd)
! store the dot product of grid velocity sc and
! the normal ss in sface.
              temp = si(i, j, k, 1)
              temp0 = si(i, j, k, 2)
              temp1 = si(i, j, k, 3)
              sfaceid(i, j, k) = temp*scd(1) + sc(1)*sid(i, j, k, 1) + &
&               temp0*scd(2) + sc(2)*sid(i, j, k, 2) + temp1*scd(3) + sc&
&               (3)*sid(i, j, k, 3)
              sfacei(i, j, k) = sc(1)*temp + sc(2)*temp0 + sc(3)*temp1
            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.
              xcd(1) = fourth*(flowdomsd(nn, groundlevel, sps)%x(i-1, j&
&               , k, 1)+flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 1)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, k-1, 1)+&
&               flowdomsd(nn, groundlevel, sps)%x(i, j, k, 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))
              xcd(2) = fourth*(flowdomsd(nn, groundlevel, sps)%x(i-1, j&
&               , k, 2)+flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 2)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, k-1, 2)+&
&               flowdomsd(nn, groundlevel, sps)%x(i, j, k, 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))
              xcd(3) = fourth*(flowdomsd(nn, groundlevel, sps)%x(i-1, j&
&               , k, 3)+flowdomsd(nn, groundlevel, sps)%x(i, j, k-1, 3)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, k-1, 3)+&
&               flowdomsd(nn, groundlevel, sps)%x(i, j, k, 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))
              rotcenterd = 0.0_8
              call cellfacevelocities_d(xc, xcd, rotcenter, rotcenterd, &
&                                 rotrate, rotrated, velxgrid, velxgridd&
&                                 , velygrid, velygridd, velzgrid, &
&                                 velzgridd, derivrotationmatrix, &
&                                 derivrotationmatrixd, sc, scd)
! store the dot product of grid velocity sc and
! the normal ss in sface.
              temp1 = sj(i, j, k, 1)
              temp0 = sj(i, j, k, 2)
              temp = sj(i, j, k, 3)
              sfacejd(i, j, k) = temp1*scd(1) + sc(1)*sjd(i, j, k, 1) + &
&               temp0*scd(2) + sc(2)*sjd(i, j, k, 2) + temp*scd(3) + sc(&
&               3)*sjd(i, j, k, 3)
              sfacej(i, j, k) = sc(1)*temp1 + sc(2)*temp0 + sc(3)*temp
            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.
              xcd(1) = fourth*(flowdomsd(nn, groundlevel, sps)%x(i, j-1&
&               , k, 1)+flowdomsd(nn, groundlevel, sps)%x(i-1, j, k, 1)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, k, 1)+&
&               flowdomsd(nn, groundlevel, sps)%x(i, j, k, 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))
              xcd(2) = fourth*(flowdomsd(nn, groundlevel, sps)%x(i, j-1&
&               , k, 2)+flowdomsd(nn, groundlevel, sps)%x(i-1, j, k, 2)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, k, 2)+&
&               flowdomsd(nn, groundlevel, sps)%x(i, j, k, 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))
              xcd(3) = fourth*(flowdomsd(nn, groundlevel, sps)%x(i, j-1&
&               , k, 3)+flowdomsd(nn, groundlevel, sps)%x(i-1, j, k, 3)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, k, 3)+&
&               flowdomsd(nn, groundlevel, sps)%x(i, j, k, 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))
              rotcenterd = 0.0_8
              call cellfacevelocities_d(xc, xcd, rotcenter, rotcenterd, &
&                                 rotrate, rotrated, velxgrid, velxgridd&
&                                 , velygrid, velygridd, velzgrid, &
&                                 velzgridd, derivrotationmatrix, &
&                                 derivrotationmatrixd, sc, scd)
! store the dot product of grid velocity sc and
! the normal ss in sface.
              temp1 = sk(i, j, k, 1)
              temp0 = sk(i, j, k, 2)
              temp = sk(i, j, k, 3)
              sfacekd(i, j, k) = temp1*scd(1) + sc(1)*skd(i, j, k, 1) + &
&               temp0*scd(2) + sc(2)*skd(i, j, k, 2) + temp*scd(3) + sc(&
&               3)*skd(i, j, k, 3)
              sfacek(i, j, k) = sc(1)*temp1 + sc(2)*temp0 + sc(3)*temp
            end do
          end do
        end do
      end if
    end if
  end subroutine gridvelocitiesfinelevel_block_d
 
  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_d, only : derivativerotmatrixrigid, getdirvector
    use utils_d, 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
    real(kind=realtype) :: arg1
! compute the mesh velocity from the given mesh mach number.
! vel{x,y,z}grid0 is the actual velocity you want at the
! geometry.
    arg1 = gammainf*pinf/rhoinf
    ainf = sqrt(arg1)
    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 forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: sc
!   with respect to varying inputs: rotrate velygrid xc velzgrid
!                derivrotationmatrix sc rotcenter velxgrid
!   rw status of diff variables: rotrate:in velygrid:in xc:in velzgrid:in
!                derivrotationmatrix:in sc:in-out rotcenter:in
!                velxgrid:in
  subroutine cellfacevelocities_d(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), intent(in) :: xcd, rotcenterd, &
&   rotrated
    real(kind=realtype), intent(in) :: velxgrid, velygrid, velzgrid
    real(kind=realtype), intent(in) :: velxgridd, velygridd, velzgridd
    real(kind=realtype), dimension(3, 3), intent(in) :: &
&   derivrotationmatrix
    real(kind=realtype), dimension(3, 3), intent(in) :: &
&   derivrotationmatrixd
    real(kind=realtype), dimension(3), intent(out) :: sc
    real(kind=realtype), dimension(3), intent(out) :: 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.
    xxcd = 0.0_8
    xxcd(1) = xcd(1) - rotcenterd(1)
    xxc(1) = xc(1) - rotcenter(1)
    xxcd(2) = xcd(2) - rotcenterd(2)
    xxc(2) = xc(2) - rotcenter(2)
    xxcd(3) = xcd(3) - rotcenterd(3)
    xxc(3) = xc(3) - rotcenter(3)
! determine the rotation speed of the face center,
! which is omega*r.
    scd(1) = xxc(3)*rotrated(2) + rotrate(2)*xxcd(3) - xxc(2)*rotrated(3&
&     ) - rotrate(3)*xxcd(2)
    sc(1) = rotrate(2)*xxc(3) - rotrate(3)*xxc(2)
    scd(2) = xxc(1)*rotrated(3) + rotrate(3)*xxcd(1) - xxc(3)*rotrated(1&
&     ) - rotrate(1)*xxcd(3)
    sc(2) = rotrate(3)*xxc(1) - rotrate(1)*xxc(3)
    scd(3) = xxc(2)*rotrated(1) + rotrate(1)*xxcd(2) - xxc(1)*rotrated(2&
&     ) - rotrate(2)*xxcd(1)
    sc(3) = rotrate(1)*xxc(2) - rotrate(2)*xxc(1)
! determine the coordinates relative to the
! rigid body rotation point.
    xxcd(1) = xcd(1)
    xxc(1) = xc(1) - rotationpoint(1)
    xxcd(2) = xcd(2)
    xxc(2) = xc(2) - rotationpoint(2)
    xxcd(3) = xcd(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.
    scd(1) = scd(1) + velxgridd + xxc(1)*derivrotationmatrixd(1, 1) + &
&     derivrotationmatrix(1, 1)*xxcd(1) + xxc(2)*derivrotationmatrixd(1&
&     , 2) + derivrotationmatrix(1, 2)*xxcd(2) + xxc(3)*&
&     derivrotationmatrixd(1, 3) + derivrotationmatrix(1, 3)*xxcd(3)
    sc(1) = sc(1) + velxgrid + derivrotationmatrix(1, 1)*xxc(1) + &
&     derivrotationmatrix(1, 2)*xxc(2) + derivrotationmatrix(1, 3)*xxc(3&
&     )
    scd(2) = scd(2) + velygridd + xxc(1)*derivrotationmatrixd(2, 1) + &
&     derivrotationmatrix(2, 1)*xxcd(1) + xxc(2)*derivrotationmatrixd(2&
&     , 2) + derivrotationmatrix(2, 2)*xxcd(2) + xxc(3)*&
&     derivrotationmatrixd(2, 3) + derivrotationmatrix(2, 3)*xxcd(3)
    sc(2) = sc(2) + velygrid + derivrotationmatrix(2, 1)*xxc(1) + &
&     derivrotationmatrix(2, 2)*xxc(2) + derivrotationmatrix(2, 3)*xxc(3&
&     )
    scd(3) = scd(3) + velzgridd + xxc(1)*derivrotationmatrixd(3, 1) + &
&     derivrotationmatrix(3, 1)*xxcd(1) + xxc(2)*derivrotationmatrixd(3&
&     , 2) + derivrotationmatrix(3, 2)*xxcd(2) + xxc(3)*&
&     derivrotationmatrixd(3, 3) + derivrotationmatrix(3, 3)*xxcd(3)
    sc(3) = sc(3) + velzgrid + derivrotationmatrix(3, 1)*xxc(1) + &
&     derivrotationmatrix(3, 2)*xxc(2) + derivrotationmatrix(3, 3)*xxc(3&
&     )
  end subroutine cellfacevelocities_d
 
  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 forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: *(*bcdata.uslip)
!   with respect to varying inputs: pinf timeref rhoinf *(flowdoms.x)
!                veldirfreestream machgrid *(*bcdata.uslip)
!   rw status of diff variables: pinf:in timeref:in rhoinf:in *(flowdoms.x):in
!                veldirfreestream:in machgrid:in *(*bcdata.uslip):in-out
!   plus diff mem management of: flowdoms.x:in bcdata:in *bcdata.uslip:in
  subroutine slipvelocitiesfinelevel_block_d(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_d, only : derivativerotmatrixrigid, &
&   derivativerotmatrixrigid_d, getdirvector
    use utils_d, 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) :: arg1
    real(kind=realtype) :: arg1d
    real(kind=realtype) :: temp
! 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)
      arg1d = gammainf*(pinfd-pinf*rhoinfd/rhoinf)/rhoinf
      arg1 = gammainf*pinf/rhoinf
      temp = sqrt(arg1)
      if (arg1 .eq. 0.0_8) then
        ainfd = 0.0_8
      else
        ainfd = arg1d/(2.0*temp)
      end if
      ainf = temp
      velxgrid0d = -(veldirfreestream(1)*(machgrid*ainfd+ainf*machgridd)&
&       +ainf*machgrid*veldirfreestreamd(1))
      velxgrid0 = ainf*machgrid*(-veldirfreestream(1))
      velygrid0d = -(veldirfreestream(2)*(machgrid*ainfd+ainf*machgridd)&
&       +ainf*machgrid*veldirfreestreamd(2))
      velygrid0 = ainf*machgrid*(-veldirfreestream(2))
      velzgrid0d = -(veldirfreestream(3)*(machgrid*ainfd+ainf*machgridd)&
&       +ainf*machgrid*veldirfreestreamd(3))
      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.
      derivrotationmatrixd = 0.0_8
      call derivativerotmatrixrigid_d(derivrotationmatrix, &
&                               derivrotationmatrixd, rotationpoint, t(1&
&                               ))
      xcd = 0.0_8
      xxcd = 0.0_8
!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
        rotrated = cgnsdoms(nbkglobal)%bocoinfo(ii)%rotrate*timerefd
        rotrate = timeref*cgnsdoms(nbkglobal)%bocoinfo(ii)%rotrate
! usewindaxis should go back here!
        velxgridd = velxgrid0d
        velxgrid = velxgrid0
        velygridd = velygrid0d
        velygrid = velygrid0
        velzgridd = velzgrid0d
        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.
              xcd(1) = fourth*(flowdomsd(nn, groundlevel, sps)%x(1, i, j&
&               , 1)+flowdomsd(nn, groundlevel, sps)%x(1, i, j-1, 1)+&
&               flowdomsd(nn, groundlevel, sps)%x(1, i-1, j, 1)+&
&               flowdomsd(nn, groundlevel, sps)%x(1, i-1, j-1, 1))
              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))
              xcd(2) = fourth*(flowdomsd(nn, groundlevel, sps)%x(1, i, j&
&               , 2)+flowdomsd(nn, groundlevel, sps)%x(1, i, j-1, 2)+&
&               flowdomsd(nn, groundlevel, sps)%x(1, i-1, j, 2)+&
&               flowdomsd(nn, groundlevel, sps)%x(1, i-1, j-1, 2))
              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))
              xcd(3) = fourth*(flowdomsd(nn, groundlevel, sps)%x(1, i, j&
&               , 3)+flowdomsd(nn, groundlevel, sps)%x(1, i, j-1, 3)+&
&               flowdomsd(nn, groundlevel, sps)%x(1, i-1, j, 3)+&
&               flowdomsd(nn, groundlevel, sps)%x(1, i-1, j-1, 3))
              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.
              xxcd(1) = xcd(1)
              xxc(1) = xc(1) - rotcenter(1)
              xxcd(2) = xcd(2)
              xxc(2) = xc(2) - rotcenter(2)
              xxcd(3) = xcd(3)
              xxc(3) = xc(3) - rotcenter(3)
! compute the velocity, which is the cross product
! of rotrate and xc.
              bcdatad(mm)%uslip(i, j, 1) = xxc(3)*rotrated(2) + rotrate(&
&               2)*xxcd(3) - xxc(2)*rotrated(3) - rotrate(3)*xxcd(2)
              bcdata(mm)%uslip(i, j, 1) = rotrate(2)*xxc(3) - rotrate(3)&
&               *xxc(2)
              bcdatad(mm)%uslip(i, j, 2) = xxc(1)*rotrated(3) + rotrate(&
&               3)*xxcd(1) - xxc(3)*rotrated(1) - rotrate(1)*xxcd(3)
              bcdata(mm)%uslip(i, j, 2) = rotrate(3)*xxc(1) - rotrate(1)&
&               *xxc(3)
              bcdatad(mm)%uslip(i, j, 3) = xxc(2)*rotrated(1) + rotrate(&
&               1)*xxcd(2) - xxc(1)*rotrated(2) - rotrate(2)*xxcd(1)
              bcdata(mm)%uslip(i, j, 3) = rotrate(1)*xxc(2) - rotrate(2)&
&               *xxc(1)
! determine the coordinates relative to the
! rigid body rotation point.
              xxcd(1) = xcd(1)
              xxc(1) = xc(1) - rotationpoint(1)
              xxcd(2) = xcd(2)
              xxc(2) = xc(2) - rotationpoint(2)
              xxcd(3) = xcd(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.
              bcdatad(mm)%uslip(i, j, 1) = bcdatad(mm)%uslip(i, j, 1) + &
&               velxgridd + xxc(1)*derivrotationmatrixd(1, 1) + &
&               derivrotationmatrix(1, 1)*xxcd(1) + xxc(2)*&
&               derivrotationmatrixd(1, 2) + derivrotationmatrix(1, 2)*&
&               xxcd(2) + xxc(3)*derivrotationmatrixd(1, 3) + &
&               derivrotationmatrix(1, 3)*xxcd(3)
              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)
              bcdatad(mm)%uslip(i, j, 2) = bcdatad(mm)%uslip(i, j, 2) + &
&               velygridd + xxc(1)*derivrotationmatrixd(2, 1) + &
&               derivrotationmatrix(2, 1)*xxcd(1) + xxc(2)*&
&               derivrotationmatrixd(2, 2) + derivrotationmatrix(2, 2)*&
&               xxcd(2) + xxc(3)*derivrotationmatrixd(2, 3) + &
&               derivrotationmatrix(2, 3)*xxcd(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)
              bcdatad(mm)%uslip(i, j, 3) = bcdatad(mm)%uslip(i, j, 3) + &
&               velzgridd + xxc(1)*derivrotationmatrixd(3, 1) + &
&               derivrotationmatrix(3, 1)*xxcd(1) + xxc(2)*&
&               derivrotationmatrixd(3, 2) + derivrotationmatrix(3, 2)*&
&               xxcd(2) + xxc(3)*derivrotationmatrixd(3, 3) + &
&               derivrotationmatrix(3, 3)*xxcd(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.
              xcd(1) = fourth*(flowdomsd(nn, groundlevel, sps)%x(il, i, &
&               j, 1)+flowdomsd(nn, groundlevel, sps)%x(il, i, j-1, 1)+&
&               flowdomsd(nn, groundlevel, sps)%x(il, i-1, j, 1)+&
&               flowdomsd(nn, groundlevel, sps)%x(il, i-1, j-1, 1))
              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))
              xcd(2) = fourth*(flowdomsd(nn, groundlevel, sps)%x(il, i, &
&               j, 2)+flowdomsd(nn, groundlevel, sps)%x(il, i, j-1, 2)+&
&               flowdomsd(nn, groundlevel, sps)%x(il, i-1, j, 2)+&
&               flowdomsd(nn, groundlevel, sps)%x(il, i-1, j-1, 2))
              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))
              xcd(3) = fourth*(flowdomsd(nn, groundlevel, sps)%x(il, i, &
&               j, 3)+flowdomsd(nn, groundlevel, sps)%x(il, i, j-1, 3)+&
&               flowdomsd(nn, groundlevel, sps)%x(il, i-1, j, 3)+&
&               flowdomsd(nn, groundlevel, sps)%x(il, i-1, j-1, 3))
              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.
              xxcd(1) = xcd(1)
              xxc(1) = xc(1) - rotcenter(1)
              xxcd(2) = xcd(2)
              xxc(2) = xc(2) - rotcenter(2)
              xxcd(3) = xcd(3)
              xxc(3) = xc(3) - rotcenter(3)
! compute the velocity, which is the cross product
! of rotrate and xc.
              bcdatad(mm)%uslip(i, j, 1) = xxc(3)*rotrated(2) + rotrate(&
&               2)*xxcd(3) - xxc(2)*rotrated(3) - rotrate(3)*xxcd(2)
              bcdata(mm)%uslip(i, j, 1) = rotrate(2)*xxc(3) - rotrate(3)&
&               *xxc(2)
              bcdatad(mm)%uslip(i, j, 2) = xxc(1)*rotrated(3) + rotrate(&
&               3)*xxcd(1) - xxc(3)*rotrated(1) - rotrate(1)*xxcd(3)
              bcdata(mm)%uslip(i, j, 2) = rotrate(3)*xxc(1) - rotrate(1)&
&               *xxc(3)
              bcdatad(mm)%uslip(i, j, 3) = xxc(2)*rotrated(1) + rotrate(&
&               1)*xxcd(2) - xxc(1)*rotrated(2) - rotrate(2)*xxcd(1)
              bcdata(mm)%uslip(i, j, 3) = rotrate(1)*xxc(2) - rotrate(2)&
&               *xxc(1)
! determine the coordinates relative to the
! rigid body rotation point.
              xxcd(1) = xcd(1)
              xxc(1) = xc(1) - rotationpoint(1)
              xxcd(2) = xcd(2)
              xxc(2) = xc(2) - rotationpoint(2)
              xxcd(3) = xcd(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.
              bcdatad(mm)%uslip(i, j, 1) = bcdatad(mm)%uslip(i, j, 1) + &
&               velxgridd + xxc(1)*derivrotationmatrixd(1, 1) + &
&               derivrotationmatrix(1, 1)*xxcd(1) + xxc(2)*&
&               derivrotationmatrixd(1, 2) + derivrotationmatrix(1, 2)*&
&               xxcd(2) + xxc(3)*derivrotationmatrixd(1, 3) + &
&               derivrotationmatrix(1, 3)*xxcd(3)
              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)
              bcdatad(mm)%uslip(i, j, 2) = bcdatad(mm)%uslip(i, j, 2) + &
&               velygridd + xxc(1)*derivrotationmatrixd(2, 1) + &
&               derivrotationmatrix(2, 1)*xxcd(1) + xxc(2)*&
&               derivrotationmatrixd(2, 2) + derivrotationmatrix(2, 2)*&
&               xxcd(2) + xxc(3)*derivrotationmatrixd(2, 3) + &
&               derivrotationmatrix(2, 3)*xxcd(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)
              bcdatad(mm)%uslip(i, j, 3) = bcdatad(mm)%uslip(i, j, 3) + &
&               velzgridd + xxc(1)*derivrotationmatrixd(3, 1) + &
&               derivrotationmatrix(3, 1)*xxcd(1) + xxc(2)*&
&               derivrotationmatrixd(3, 2) + derivrotationmatrix(3, 2)*&
&               xxcd(2) + xxc(3)*derivrotationmatrixd(3, 3) + &
&               derivrotationmatrix(3, 3)*xxcd(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.
              xcd(1) = fourth*(flowdomsd(nn, groundlevel, sps)%x(i, 1, j&
&               , 1)+flowdomsd(nn, groundlevel, sps)%x(i, 1, j-1, 1)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, 1, j, 1)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, 1, j-1, 1))
              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))
              xcd(2) = fourth*(flowdomsd(nn, groundlevel, sps)%x(i, 1, j&
&               , 2)+flowdomsd(nn, groundlevel, sps)%x(i, 1, j-1, 2)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, 1, j, 2)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, 1, j-1, 2))
              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))
              xcd(3) = fourth*(flowdomsd(nn, groundlevel, sps)%x(i, 1, j&
&               , 3)+flowdomsd(nn, groundlevel, sps)%x(i, 1, j-1, 3)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, 1, j, 3)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, 1, j-1, 3))
              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.
              xxcd(1) = xcd(1)
              xxc(1) = xc(1) - rotcenter(1)
              xxcd(2) = xcd(2)
              xxc(2) = xc(2) - rotcenter(2)
              xxcd(3) = xcd(3)
              xxc(3) = xc(3) - rotcenter(3)
! compute the velocity, which is the cross product
! of rotrate and xc.
              bcdatad(mm)%uslip(i, j, 1) = xxc(3)*rotrated(2) + rotrate(&
&               2)*xxcd(3) - xxc(2)*rotrated(3) - rotrate(3)*xxcd(2)
              bcdata(mm)%uslip(i, j, 1) = rotrate(2)*xxc(3) - rotrate(3)&
&               *xxc(2)
              bcdatad(mm)%uslip(i, j, 2) = xxc(1)*rotrated(3) + rotrate(&
&               3)*xxcd(1) - xxc(3)*rotrated(1) - rotrate(1)*xxcd(3)
              bcdata(mm)%uslip(i, j, 2) = rotrate(3)*xxc(1) - rotrate(1)&
&               *xxc(3)
              bcdatad(mm)%uslip(i, j, 3) = xxc(2)*rotrated(1) + rotrate(&
&               1)*xxcd(2) - xxc(1)*rotrated(2) - rotrate(2)*xxcd(1)
              bcdata(mm)%uslip(i, j, 3) = rotrate(1)*xxc(2) - rotrate(2)&
&               *xxc(1)
! determine the coordinates relative to the
! rigid body rotation point.
              xxcd(1) = xcd(1)
              xxc(1) = xc(1) - rotationpoint(1)
              xxcd(2) = xcd(2)
              xxc(2) = xc(2) - rotationpoint(2)
              xxcd(3) = xcd(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.
              bcdatad(mm)%uslip(i, j, 1) = bcdatad(mm)%uslip(i, j, 1) + &
&               velxgridd + xxc(1)*derivrotationmatrixd(1, 1) + &
&               derivrotationmatrix(1, 1)*xxcd(1) + xxc(2)*&
&               derivrotationmatrixd(1, 2) + derivrotationmatrix(1, 2)*&
&               xxcd(2) + xxc(3)*derivrotationmatrixd(1, 3) + &
&               derivrotationmatrix(1, 3)*xxcd(3)
              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)
              bcdatad(mm)%uslip(i, j, 2) = bcdatad(mm)%uslip(i, j, 2) + &
&               velygridd + xxc(1)*derivrotationmatrixd(2, 1) + &
&               derivrotationmatrix(2, 1)*xxcd(1) + xxc(2)*&
&               derivrotationmatrixd(2, 2) + derivrotationmatrix(2, 2)*&
&               xxcd(2) + xxc(3)*derivrotationmatrixd(2, 3) + &
&               derivrotationmatrix(2, 3)*xxcd(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)
              bcdatad(mm)%uslip(i, j, 3) = bcdatad(mm)%uslip(i, j, 3) + &
&               velzgridd + xxc(1)*derivrotationmatrixd(3, 1) + &
&               derivrotationmatrix(3, 1)*xxcd(1) + xxc(2)*&
&               derivrotationmatrixd(3, 2) + derivrotationmatrix(3, 2)*&
&               xxcd(2) + xxc(3)*derivrotationmatrixd(3, 3) + &
&               derivrotationmatrix(3, 3)*xxcd(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.
              xcd(1) = fourth*(flowdomsd(nn, groundlevel, sps)%x(i, jl, &
&               j, 1)+flowdomsd(nn, groundlevel, sps)%x(i, jl, j-1, 1)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, jl, j, 1)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, jl, j-1, 1))
              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))
              xcd(2) = fourth*(flowdomsd(nn, groundlevel, sps)%x(i, jl, &
&               j, 2)+flowdomsd(nn, groundlevel, sps)%x(i, jl, j-1, 2)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, jl, j, 2)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, jl, j-1, 2))
              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))
              xcd(3) = fourth*(flowdomsd(nn, groundlevel, sps)%x(i, jl, &
&               j, 3)+flowdomsd(nn, groundlevel, sps)%x(i, jl, j-1, 3)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, jl, j, 3)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, jl, j-1, 3))
              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.
              xxcd(1) = xcd(1)
              xxc(1) = xc(1) - rotcenter(1)
              xxcd(2) = xcd(2)
              xxc(2) = xc(2) - rotcenter(2)
              xxcd(3) = xcd(3)
              xxc(3) = xc(3) - rotcenter(3)
! compute the velocity, which is the cross product
! of rotrate and xc.
              bcdatad(mm)%uslip(i, j, 1) = xxc(3)*rotrated(2) + rotrate(&
&               2)*xxcd(3) - xxc(2)*rotrated(3) - rotrate(3)*xxcd(2)
              bcdata(mm)%uslip(i, j, 1) = rotrate(2)*xxc(3) - rotrate(3)&
&               *xxc(2)
              bcdatad(mm)%uslip(i, j, 2) = xxc(1)*rotrated(3) + rotrate(&
&               3)*xxcd(1) - xxc(3)*rotrated(1) - rotrate(1)*xxcd(3)
              bcdata(mm)%uslip(i, j, 2) = rotrate(3)*xxc(1) - rotrate(1)&
&               *xxc(3)
              bcdatad(mm)%uslip(i, j, 3) = xxc(2)*rotrated(1) + rotrate(&
&               1)*xxcd(2) - xxc(1)*rotrated(2) - rotrate(2)*xxcd(1)
              bcdata(mm)%uslip(i, j, 3) = rotrate(1)*xxc(2) - rotrate(2)&
&               *xxc(1)
! determine the coordinates relative to the
! rigid body rotation point.
              xxcd(1) = xcd(1)
              xxc(1) = xc(1) - rotationpoint(1)
              xxcd(2) = xcd(2)
              xxc(2) = xc(2) - rotationpoint(2)
              xxcd(3) = xcd(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.
              bcdatad(mm)%uslip(i, j, 1) = bcdatad(mm)%uslip(i, j, 1) + &
&               velxgridd + xxc(1)*derivrotationmatrixd(1, 1) + &
&               derivrotationmatrix(1, 1)*xxcd(1) + xxc(2)*&
&               derivrotationmatrixd(1, 2) + derivrotationmatrix(1, 2)*&
&               xxcd(2) + xxc(3)*derivrotationmatrixd(1, 3) + &
&               derivrotationmatrix(1, 3)*xxcd(3)
              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)
              bcdatad(mm)%uslip(i, j, 2) = bcdatad(mm)%uslip(i, j, 2) + &
&               velygridd + xxc(1)*derivrotationmatrixd(2, 1) + &
&               derivrotationmatrix(2, 1)*xxcd(1) + xxc(2)*&
&               derivrotationmatrixd(2, 2) + derivrotationmatrix(2, 2)*&
&               xxcd(2) + xxc(3)*derivrotationmatrixd(2, 3) + &
&               derivrotationmatrix(2, 3)*xxcd(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)
              bcdatad(mm)%uslip(i, j, 3) = bcdatad(mm)%uslip(i, j, 3) + &
&               velzgridd + xxc(1)*derivrotationmatrixd(3, 1) + &
&               derivrotationmatrix(3, 1)*xxcd(1) + xxc(2)*&
&               derivrotationmatrixd(3, 2) + derivrotationmatrix(3, 2)*&
&               xxcd(2) + xxc(3)*derivrotationmatrixd(3, 3) + &
&               derivrotationmatrix(3, 3)*xxcd(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.
              xcd(1) = fourth*(flowdomsd(nn, groundlevel, sps)%x(i, j, 1&
&               , 1)+flowdomsd(nn, groundlevel, sps)%x(i, j-1, 1, 1)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, 1, 1)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, 1, 1))
              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))
              xcd(2) = fourth*(flowdomsd(nn, groundlevel, sps)%x(i, j, 1&
&               , 2)+flowdomsd(nn, groundlevel, sps)%x(i, j-1, 1, 2)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, 1, 2)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, 1, 2))
              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))
              xcd(3) = fourth*(flowdomsd(nn, groundlevel, sps)%x(i, j, 1&
&               , 3)+flowdomsd(nn, groundlevel, sps)%x(i, j-1, 1, 3)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, 1, 3)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, 1, 3))
              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.
              xxcd(1) = xcd(1)
              xxc(1) = xc(1) - rotcenter(1)
              xxcd(2) = xcd(2)
              xxc(2) = xc(2) - rotcenter(2)
              xxcd(3) = xcd(3)
              xxc(3) = xc(3) - rotcenter(3)
! compute the velocity, which is the cross product
! of rotrate and xc.
              bcdatad(mm)%uslip(i, j, 1) = xxc(3)*rotrated(2) + rotrate(&
&               2)*xxcd(3) - xxc(2)*rotrated(3) - rotrate(3)*xxcd(2)
              bcdata(mm)%uslip(i, j, 1) = rotrate(2)*xxc(3) - rotrate(3)&
&               *xxc(2)
              bcdatad(mm)%uslip(i, j, 2) = xxc(1)*rotrated(3) + rotrate(&
&               3)*xxcd(1) - xxc(3)*rotrated(1) - rotrate(1)*xxcd(3)
              bcdata(mm)%uslip(i, j, 2) = rotrate(3)*xxc(1) - rotrate(1)&
&               *xxc(3)
              bcdatad(mm)%uslip(i, j, 3) = xxc(2)*rotrated(1) + rotrate(&
&               1)*xxcd(2) - xxc(1)*rotrated(2) - rotrate(2)*xxcd(1)
              bcdata(mm)%uslip(i, j, 3) = rotrate(1)*xxc(2) - rotrate(2)&
&               *xxc(1)
! determine the coordinates relative to the
! rigid body rotation point.
              xxcd(1) = xcd(1)
              xxc(1) = xc(1) - rotationpoint(1)
              xxcd(2) = xcd(2)
              xxc(2) = xc(2) - rotationpoint(2)
              xxcd(3) = xcd(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.
              bcdatad(mm)%uslip(i, j, 1) = bcdatad(mm)%uslip(i, j, 1) + &
&               velxgridd + xxc(1)*derivrotationmatrixd(1, 1) + &
&               derivrotationmatrix(1, 1)*xxcd(1) + xxc(2)*&
&               derivrotationmatrixd(1, 2) + derivrotationmatrix(1, 2)*&
&               xxcd(2) + xxc(3)*derivrotationmatrixd(1, 3) + &
&               derivrotationmatrix(1, 3)*xxcd(3)
              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)
              bcdatad(mm)%uslip(i, j, 2) = bcdatad(mm)%uslip(i, j, 2) + &
&               velygridd + xxc(1)*derivrotationmatrixd(2, 1) + &
&               derivrotationmatrix(2, 1)*xxcd(1) + xxc(2)*&
&               derivrotationmatrixd(2, 2) + derivrotationmatrix(2, 2)*&
&               xxcd(2) + xxc(3)*derivrotationmatrixd(2, 3) + &
&               derivrotationmatrix(2, 3)*xxcd(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)
              bcdatad(mm)%uslip(i, j, 3) = bcdatad(mm)%uslip(i, j, 3) + &
&               velzgridd + xxc(1)*derivrotationmatrixd(3, 1) + &
&               derivrotationmatrix(3, 1)*xxcd(1) + xxc(2)*&
&               derivrotationmatrixd(3, 2) + derivrotationmatrix(3, 2)*&
&               xxcd(2) + xxc(3)*derivrotationmatrixd(3, 3) + &
&               derivrotationmatrix(3, 3)*xxcd(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.
              xcd(1) = fourth*(flowdomsd(nn, groundlevel, sps)%x(i, j, &
&               kl, 1)+flowdomsd(nn, groundlevel, sps)%x(i, j-1, kl, 1)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, kl, 1)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, kl, 1))
              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))
              xcd(2) = fourth*(flowdomsd(nn, groundlevel, sps)%x(i, j, &
&               kl, 2)+flowdomsd(nn, groundlevel, sps)%x(i, j-1, kl, 2)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, kl, 2)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, kl, 2))
              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))
              xcd(3) = fourth*(flowdomsd(nn, groundlevel, sps)%x(i, j, &
&               kl, 3)+flowdomsd(nn, groundlevel, sps)%x(i, j-1, kl, 3)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j, kl, 3)+&
&               flowdomsd(nn, groundlevel, sps)%x(i-1, j-1, kl, 3))
              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.
              xxcd(1) = xcd(1)
              xxc(1) = xc(1) - rotcenter(1)
              xxcd(2) = xcd(2)
              xxc(2) = xc(2) - rotcenter(2)
              xxcd(3) = xcd(3)
              xxc(3) = xc(3) - rotcenter(3)
! compute the velocity, which is the cross product
! of rotrate and xc.
              bcdatad(mm)%uslip(i, j, 1) = xxc(3)*rotrated(2) + rotrate(&
&               2)*xxcd(3) - xxc(2)*rotrated(3) - rotrate(3)*xxcd(2)
              bcdata(mm)%uslip(i, j, 1) = rotrate(2)*xxc(3) - rotrate(3)&
&               *xxc(2)
              bcdatad(mm)%uslip(i, j, 2) = xxc(1)*rotrated(3) + rotrate(&
&               3)*xxcd(1) - xxc(3)*rotrated(1) - rotrate(1)*xxcd(3)
              bcdata(mm)%uslip(i, j, 2) = rotrate(3)*xxc(1) - rotrate(1)&
&               *xxc(3)
              bcdatad(mm)%uslip(i, j, 3) = xxc(2)*rotrated(1) + rotrate(&
&               1)*xxcd(2) - xxc(1)*rotrated(2) - rotrate(2)*xxcd(1)
              bcdata(mm)%uslip(i, j, 3) = rotrate(1)*xxc(2) - rotrate(2)&
&               *xxc(1)
! determine the coordinates relative to the
! rigid body rotation point.
              xxcd(1) = xcd(1)
              xxc(1) = xc(1) - rotationpoint(1)
              xxcd(2) = xcd(2)
              xxc(2) = xc(2) - rotationpoint(2)
              xxcd(3) = xcd(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.
              bcdatad(mm)%uslip(i, j, 1) = bcdatad(mm)%uslip(i, j, 1) + &
&               velxgridd + xxc(1)*derivrotationmatrixd(1, 1) + &
&               derivrotationmatrix(1, 1)*xxcd(1) + xxc(2)*&
&               derivrotationmatrixd(1, 2) + derivrotationmatrix(1, 2)*&
&               xxcd(2) + xxc(3)*derivrotationmatrixd(1, 3) + &
&               derivrotationmatrix(1, 3)*xxcd(3)
              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)
              bcdatad(mm)%uslip(i, j, 2) = bcdatad(mm)%uslip(i, j, 2) + &
&               velygridd + xxc(1)*derivrotationmatrixd(2, 1) + &
&               derivrotationmatrix(2, 1)*xxcd(1) + xxc(2)*&
&               derivrotationmatrixd(2, 2) + derivrotationmatrix(2, 2)*&
&               xxcd(2) + xxc(3)*derivrotationmatrixd(2, 3) + &
&               derivrotationmatrix(2, 3)*xxcd(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)
              bcdatad(mm)%uslip(i, j, 3) = bcdatad(mm)%uslip(i, j, 3) + &
&               velzgridd + xxc(1)*derivrotationmatrixd(3, 1) + &
&               derivrotationmatrix(3, 1)*xxcd(1) + xxc(2)*&
&               derivrotationmatrixd(3, 2) + derivrotationmatrix(3, 2)*&
&               xxcd(2) + xxc(3)*derivrotationmatrixd(3, 3) + &
&               derivrotationmatrix(3, 3)*xxcd(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_d
 
  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_d, only : derivativerotmatrixrigid, getdirvector
    use utils_d, 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
    real(kind=realtype) :: arg1
! 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)
      arg1 = gammainf*pinf/rhoinf
      ainf = sqrt(arg1)
      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 forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: *(*bcdata.rface)
!   with respect to varying inputs: *sfacei *sfacej *sfacek *si
!                *sj *sk *(*bcdata.rface)
!   rw status of diff variables: *sfacei:in *sfacej:in *sfacek:in
!                *si:in *sj:in *sk:in *(*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_d(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) :: arg1
    real(kind=realtype) :: arg1d
    real(kind=realtype) :: temp
    real(kind=realtype) :: temp0
    real(kind=realtype) :: temp1
! 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.
                arg1d = 2*si(1, i, j, 1)*sid(1, i, j, 1) + 2*si(1, i, j&
&                 , 2)*sid(1, i, j, 2) + 2*si(1, i, j, 3)*sid(1, i, j, 3&
&                 )
                arg1 = si(1, i, j, 1)**2 + si(1, i, j, 2)**2 + si(1, i, &
&                 j, 3)**2
                temp = sqrt(arg1)
                if (arg1 .eq. 0.0_8) then
                  weightd = 0.0_8
                else
                  weightd = arg1d/(2.0*temp)
                end if
                weight = temp
                if (weight .gt. zero) then
                  weightd = -(mult*weightd/weight**2)
                  weight = mult/weight
                end if
! compute the normal velocity based on the outward
! pointing unit normal.
                bcdatad(mm)%rface(i, j) = sfacei(1, i, j)*weightd + &
&                 weight*sfaceid(1, i, j)
                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.
                temp = si(il, i, j, 1)
                temp0 = si(il, i, j, 2)
                temp1 = si(il, i, j, 3)
                arg1d = 2*temp*sid(il, i, j, 1) + 2*temp0*sid(il, i, j, &
&                 2) + 2*temp1*sid(il, i, j, 3)
                arg1 = temp*temp + temp0*temp0 + temp1*temp1
                temp1 = sqrt(arg1)
                if (arg1 .eq. 0.0_8) then
                  weightd = 0.0_8
                else
                  weightd = arg1d/(2.0*temp1)
                end if
                weight = temp1
                if (weight .gt. zero) then
                  weightd = -(mult*weightd/weight**2)
                  weight = mult/weight
                end if
! compute the normal velocity based on the outward
! pointing unit normal.
                bcdatad(mm)%rface(i, j) = sfacei(il, i, j)*weightd + &
&                 weight*sfaceid(il, i, j)
                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.
                arg1d = 2*sj(i, 1, j, 1)*sjd(i, 1, j, 1) + 2*sj(i, 1, j&
&                 , 2)*sjd(i, 1, j, 2) + 2*sj(i, 1, j, 3)*sjd(i, 1, j, 3&
&                 )
                arg1 = sj(i, 1, j, 1)**2 + sj(i, 1, j, 2)**2 + sj(i, 1, &
&                 j, 3)**2
                temp1 = sqrt(arg1)
                if (arg1 .eq. 0.0_8) then
                  weightd = 0.0_8
                else
                  weightd = arg1d/(2.0*temp1)
                end if
                weight = temp1
                if (weight .gt. zero) then
                  weightd = -(mult*weightd/weight**2)
                  weight = mult/weight
                end if
! compute the normal velocity based on the outward
! pointing unit normal.
                bcdatad(mm)%rface(i, j) = sfacej(i, 1, j)*weightd + &
&                 weight*sfacejd(i, 1, j)
                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.
                temp1 = sj(i, jl, j, 1)
                temp0 = sj(i, jl, j, 2)
                temp = sj(i, jl, j, 3)
                arg1d = 2*temp1*sjd(i, jl, j, 1) + 2*temp0*sjd(i, jl, j&
&                 , 2) + 2*temp*sjd(i, jl, j, 3)
                arg1 = temp1*temp1 + temp0*temp0 + temp*temp
                temp1 = sqrt(arg1)
                if (arg1 .eq. 0.0_8) then
                  weightd = 0.0_8
                else
                  weightd = arg1d/(2.0*temp1)
                end if
                weight = temp1
                if (weight .gt. zero) then
                  weightd = -(mult*weightd/weight**2)
                  weight = mult/weight
                end if
! compute the normal velocity based on the outward
! pointing unit normal.
                bcdatad(mm)%rface(i, j) = sfacej(i, jl, j)*weightd + &
&                 weight*sfacejd(i, jl, j)
                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.
                arg1d = 2*sk(i, j, 1, 1)*skd(i, j, 1, 1) + 2*sk(i, j, 1&
&                 , 2)*skd(i, j, 1, 2) + 2*sk(i, j, 1, 3)*skd(i, j, 1, 3&
&                 )
                arg1 = sk(i, j, 1, 1)**2 + sk(i, j, 1, 2)**2 + sk(i, j, &
&                 1, 3)**2
                temp1 = sqrt(arg1)
                if (arg1 .eq. 0.0_8) then
                  weightd = 0.0_8
                else
                  weightd = arg1d/(2.0*temp1)
                end if
                weight = temp1
                if (weight .gt. zero) then
                  weightd = -(mult*weightd/weight**2)
                  weight = mult/weight
                end if
! compute the normal velocity based on the outward
! pointing unit normal.
                bcdatad(mm)%rface(i, j) = sfacek(i, j, 1)*weightd + &
&                 weight*sfacekd(i, j, 1)
                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.
                temp1 = sk(i, j, kl, 1)
                temp0 = sk(i, j, kl, 2)
                temp = sk(i, j, kl, 3)
                arg1d = 2*temp1*skd(i, j, kl, 1) + 2*temp0*skd(i, j, kl&
&                 , 2) + 2*temp*skd(i, j, kl, 3)
                arg1 = temp1*temp1 + temp0*temp0 + temp*temp
                temp1 = sqrt(arg1)
                if (arg1 .eq. 0.0_8) then
                  weightd = 0.0_8
                else
                  weightd = arg1d/(2.0*temp1)
                end if
                weight = temp1
                if (weight .gt. zero) then
                  weightd = -(mult*weightd/weight**2)
                  weight = mult/weight
                end if
! compute the normal velocity based on the outward
! pointing unit normal.
                bcdatad(mm)%rface(i, j) = sfacek(i, j, kl)*weightd + &
&                 weight*sfacekd(i, j, kl)
                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)) then
          bcdatad(mm)%rface = 0.0_8
          bcdata(mm)%rface = zero
        end if
      end do
    end if
  end subroutine normalvelocities_block_d
 
  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
    real(kind=realtype) :: arg1
! 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.
                arg1 = si(1, i, j, 1)**2 + si(1, i, j, 2)**2 + si(1, i, &
&                 j, 3)**2
                weight = sqrt(arg1)
                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.
                arg1 = si(il, i, j, 1)**2 + si(il, i, j, 2)**2 + si(il, &
&                 i, j, 3)**2
                weight = sqrt(arg1)
                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.
                arg1 = sj(i, 1, j, 1)**2 + sj(i, 1, j, 2)**2 + sj(i, 1, &
&                 j, 3)**2
                weight = sqrt(arg1)
                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.
                arg1 = sj(i, jl, j, 1)**2 + sj(i, jl, j, 2)**2 + sj(i, &
&                 jl, j, 3)**2
                weight = sqrt(arg1)
                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.
                arg1 = sk(i, j, 1, 1)**2 + sk(i, j, 1, 2)**2 + sk(i, j, &
&                 1, 3)**2
                weight = sqrt(arg1)
                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.
                arg1 = sk(i, j, kl, 1)**2 + sk(i, j, kl, 2)**2 + sk(i, j&
&                 , kl, 3)**2
                weight = sqrt(arg1)
                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_d