solverUtils_fast_b.f90

Go to the documentation of this file.

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