utils_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 utils_d
  implicit none
 
contains
  function char2str(chararray, n)
    use constants
    implicit none
!
!      function arguments.
!
    character, dimension(maxcgnsnamelen), intent(in) :: chararray
    integer(kind=inttype), intent(in) :: n
!
!      function type
!
    character(len=n) :: char2str
!
!      local variables.
!
    integer(kind=inttype) :: i
    do i=1,n
      char2str(i:i) = chararray(i)
    end do
  end function char2str
 
  function tsbeta(degreepolbeta, coefpolbeta, degreefourbeta, &
&   omegafourbeta, coscoeffourbeta, sincoeffourbeta, t)
!
!       tsbeta computes the angle of attack for a given time interval
!       in a time spectral solution.
!
    use constants
    use inputphysics, only : equationmode
    implicit none
!
!      function type
!
    real(kind=realtype) :: tsbeta
!
!      function arguments.
!
    integer(kind=inttype), intent(in) :: degreepolbeta
    integer(kind=inttype), intent(in) :: degreefourbeta
    real(kind=realtype), intent(in) :: omegafourbeta, t
    real(kind=realtype), dimension(0:*), intent(in) :: coefpolbeta
    real(kind=realtype), dimension(0:*), intent(in) :: coscoeffourbeta
    real(kind=realtype), dimension(*), intent(in) :: sincoeffourbeta
!
!      local variables.
!
    integer(kind=inttype) :: nn
    real(kind=realtype) :: beta, val
    intrinsic cos
    intrinsic sin
! return immediately if this is a steady computation.
    if (equationmode .eq. steady) then
      tsbeta = zero
      return
    else
! compute the polynomial contribution. if no polynomial was
! specified, the value of index 0 is set to zero automatically.
      beta = coefpolbeta(0)
      do nn=1,degreepolbeta
        beta = beta + coefpolbeta(nn)*t**nn
      end do
! compute the fourier contribution. again the cosine coefficient
! of index 0 is defaulted to zero if not specified.
      beta = beta + coscoeffourbeta(0)
      do nn=1,degreefourbeta
        val = nn*omegafourbeta*t
        beta = beta + coscoeffourbeta(nn)*cos(val) + sincoeffourbeta(nn)&
&         *sin(val)
      end do
! set tsbeta to phi.
      tsbeta = beta
    end if
  end function tsbeta
 
  function tsbetadot(degreepolbeta, coefpolbeta, degreefourbeta, &
&   omegafourbeta, coscoeffourbeta, sincoeffourbeta, t)
!
!       tsbeta computes the angle of attack for a given time interval
!       in a time spectral solution.
!
    use constants
    use inputphysics, only : equationmode
    implicit none
!
!      function type
!
    real(kind=realtype) :: tsbetadot
!
!      function arguments.
!
    integer(kind=inttype), intent(in) :: degreepolbeta
    integer(kind=inttype), intent(in) :: degreefourbeta
    real(kind=realtype), intent(in) :: omegafourbeta, t
    real(kind=realtype), dimension(0:*), intent(in) :: coefpolbeta
    real(kind=realtype), dimension(0:*), intent(in) :: coscoeffourbeta
    real(kind=realtype), dimension(*), intent(in) :: sincoeffourbeta
!
!      local variables.
!
    integer(kind=inttype) :: nn
    real(kind=realtype) :: betadot, val
    intrinsic sin
    intrinsic cos
! return immediately if this is a steady computation.
    if (equationmode .eq. steady) then
      tsbetadot = zero
      return
    else
! compute the polynomial contribution. if no polynomial was
! specified, the value of index 0 is set to zero automatically.
      betadot = zero
      do nn=1,degreepolbeta
        betadot = betadot + nn*coefpolbeta(nn)*t**(nn-1)
      end do
! compute the fourier contribution. again the cosine coefficient
! of index 0 is defaulted to zero if not specified.
      do nn=1,degreefourbeta
        val = nn*omegafourbeta
        betadot = betadot - val*coscoeffourbeta(nn)*sin(val*t) + val*&
&         sincoeffourbeta(nn)*cos(val*t)
      end do
! set tsbeta to phi.
      tsbetadot = betadot
    end if
  end function tsbetadot
 
  function tsmach(degreepolmach, coefpolmach, degreefourmach, &
&   omegafourmach, coscoeffourmach, sincoeffourmach, t)
!
!       tsmach computes the mach number for a given time interval
!       in a time spectral solution.
!
    use constants
    use inputphysics, only : equationmode
    implicit none
!
!      function type
!
    real(kind=realtype) :: tsmach
!
!      function arguments.
!
    integer(kind=inttype), intent(in) :: degreepolmach
    integer(kind=inttype), intent(in) :: degreefourmach
    real(kind=realtype), intent(in) :: omegafourmach, t
    real(kind=realtype), dimension(0:*), intent(in) :: coefpolmach
    real(kind=realtype), dimension(0:*), intent(in) :: coscoeffourmach
    real(kind=realtype), dimension(*), intent(in) :: sincoeffourmach
!
!      local variables.
!
    integer(kind=inttype) :: nn
    real(kind=realtype) :: intervalmach, val
    intrinsic cos
    intrinsic sin
! return immediately if this is a steady computation.
    if (equationmode .eq. steady) then
      tsmach = zero
      return
    else
! compute the polynomial contribution. if no polynomial was
! specified, the value of index 0 is set to zero automatically.
      intervalmach = coefpolmach(0)
      do nn=1,degreepolmach
        intervalmach = intervalmach + coefpolmach(nn)*t**nn
      end do
! compute the fourier contribution. again the cosine coefficient
! of index 0 is defaulted to zero if not specified.
      intervalmach = intervalmach + coscoeffourmach(0)
      do nn=1,degreefourmach
        val = nn*omegafourmach*t
        intervalmach = intervalmach + coscoeffourmach(nn)*cos(val) + &
&         sincoeffourmach(nn)*sin(val)
      end do
      print*, 'intsmach', intervalmach, nn, val, t
! set tsmach to phi.
      tsmach = intervalmach
    end if
  end function tsmach
 
  function tsmachdot(degreepolmach, coefpolmach, degreefourmach, &
&   omegafourmach, coscoeffourmach, sincoeffourmach, t)
!
!       tsmach computes the angle of attack for a given time interval
!       in a time spectral solution.
!
    use constants
    use inputphysics, only : equationmode
    implicit none
!
!      function type
!
    real(kind=realtype) :: tsmachdot
!
!      function arguments.
!
    integer(kind=inttype), intent(in) :: degreepolmach
    integer(kind=inttype), intent(in) :: degreefourmach
    real(kind=realtype), intent(in) :: omegafourmach, t
    real(kind=realtype), dimension(0:*), intent(in) :: coefpolmach
    real(kind=realtype), dimension(0:*), intent(in) :: coscoeffourmach
    real(kind=realtype), dimension(*), intent(in) :: sincoeffourmach
!
!      local variables.
!
    integer(kind=inttype) :: nn
    real(kind=realtype) :: machdot, val
    intrinsic sin
    intrinsic cos
! return immediately if this is a steady computation.
    if (equationmode .eq. steady) then
      tsmachdot = zero
      return
    else
! compute the polynomial contribution. if no polynomial was
! specified, the value of index 0 is set to zero automatically.
      machdot = zero
      do nn=1,degreepolmach
        machdot = machdot + nn*coefpolmach(nn)*t**(nn-1)
      end do
! compute the fourier contribution. again the cosine coefficient
! of index 0 is defaulted to zero if not specified.
      do nn=1,degreefourmach
        val = nn*omegafourmach
        machdot = machdot - val*coscoeffourmach(nn)*sin(val*t) + val*&
&         sincoeffourmach(nn)*cos(val*t)
      end do
! set tsmach to phi.
      tsmachdot = machdot
    end if
  end function tsmachdot
 
  function tsalpha(degreepolalpha, coefpolalpha, degreefouralpha, &
&   omegafouralpha, coscoeffouralpha, sincoeffouralpha, t)
!
!       tsalpha computes the angle of attack for a given time interval
!       in a time spectral solution.
!
    use constants
    use inputphysics, only : equationmode
    implicit none
!
!      function type
!
    real(kind=realtype) :: tsalpha
!
!      function arguments.
!
    integer(kind=inttype), intent(in) :: degreepolalpha
    integer(kind=inttype), intent(in) :: degreefouralpha
    real(kind=realtype), intent(in) :: omegafouralpha, t
    real(kind=realtype), dimension(0:*), intent(in) :: coefpolalpha
    real(kind=realtype), dimension(0:*), intent(in) :: coscoeffouralpha
    real(kind=realtype), dimension(*), intent(in) :: sincoeffouralpha
!
!      local variables.
!
    integer(kind=inttype) :: nn
    real(kind=realtype) :: alpha, val
    intrinsic cos
    intrinsic sin
! return immediately if this is a steady computation.
    if (equationmode .eq. steady) then
      tsalpha = zero
      return
    else
! compute the polynomial contribution. if no polynomial was
! specified, the value of index 0 is set to zero automatically.
      alpha = coefpolalpha(0)
      do nn=1,degreepolalpha
        alpha = alpha + coefpolalpha(nn)*t**nn
      end do
! compute the fourier contribution. again the cosine coefficient
! of index 0 is defaulted to zero if not specified.
      alpha = alpha + coscoeffouralpha(0)
      do nn=1,degreefouralpha
        val = nn*omegafouralpha*t
        alpha = alpha + coscoeffouralpha(nn)*cos(val) + sincoeffouralpha&
&         (nn)*sin(val)
      end do
!print *,'intsalpha',alpha,nn,val,t
! set tsalpha to phi.
      tsalpha = alpha
    end if
  end function tsalpha
 
  function tsalphadot(degreepolalpha, coefpolalpha, degreefouralpha, &
&   omegafouralpha, coscoeffouralpha, sincoeffouralpha, t)
!
!       tsalpha computes the angle of attack for a given time interval
!       in a time spectral solution.
!
    use constants
    use inputphysics, only : equationmode
    implicit none
!
!      function type
!
    real(kind=realtype) :: tsalphadot
!
!      function arguments.
!
    integer(kind=inttype), intent(in) :: degreepolalpha
    integer(kind=inttype), intent(in) :: degreefouralpha
    real(kind=realtype), intent(in) :: omegafouralpha, t
    real(kind=realtype), dimension(0:*), intent(in) :: coefpolalpha
    real(kind=realtype), dimension(0:*), intent(in) :: coscoeffouralpha
    real(kind=realtype), dimension(*), intent(in) :: sincoeffouralpha
!
!      local variables.
!
    integer(kind=inttype) :: nn
    real(kind=realtype) :: alphadot, val
    intrinsic sin
    intrinsic cos
! return immediately if this is a steady computation.
    if (equationmode .eq. steady) then
      tsalphadot = zero
      return
    else
! compute the polynomial contribution. if no polynomial was
! specified, the value of index 0 is set to zero automatically.
      alphadot = zero
      do nn=1,degreepolalpha
        alphadot = alphadot + nn*coefpolalpha(nn)*t**(nn-1)
      end do
! compute the fourier contribution. again the cosine coefficient
! of index 0 is defaulted to zero if not specified.
      do nn=1,degreefouralpha
        val = nn*omegafouralpha
        alphadot = alphadot - val*coscoeffouralpha(nn)*sin(val*t) + val*&
&         sincoeffouralpha(nn)*cos(val*t)
      end do
! set tsalpha to phi.
      tsalphadot = alphadot
    end if
  end function tsalphadot
 
!  differentiation of derivativerigidrotangle in forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: derivativerigidrotangle
!   with respect to varying inputs: timeref
  real(kind=realtype) function derivativerigidrotangle_d(degreepolrot, &
&   coefpolrot, degreefourrot, omegafourrot, coscoeffourrot, &
&   sincoeffourrot, t, derivativerigidrotangle)
!
!       derivativerigidrotangle computes the time derivative of the
!       rigid body rotation angle at the given time for the given
!       arguments. the angle is described by a combination of a
!       polynomial and fourier series.
!
    use constants
    use inputphysics, only : equationmode
    use flowvarrefstate, only : timeref, timerefd
    implicit none
!
!      function type
!
    real(kind=realtype) :: derivativerigidrotangle
!
!      function arguments.
!
    integer(kind=inttype), intent(in) :: degreepolrot
    integer(kind=inttype), intent(in) :: degreefourrot
    real(kind=realtype), intent(in) :: omegafourrot, t
    real(kind=realtype), dimension(0:*), intent(in) :: coefpolrot
    real(kind=realtype), dimension(0:*), intent(in) :: coscoeffourrot
    real(kind=realtype), dimension(*), intent(in) :: sincoeffourrot
!
!      local variables.
!
    integer(kind=inttype) :: nn
    real(kind=realtype) :: dphi, val
    intrinsic sin
    intrinsic cos
! return immediately if this is a steady computation.
    if (equationmode .eq. steady) then
      derivativerigidrotangle = zero
      derivativerigidrotangle_d = 0.0_8
      return
    else
! compute the polynomial contribution.
      dphi = zero
      do nn=1,degreepolrot
        dphi = dphi + nn*coefpolrot(nn)*t**(nn-1)
      end do
! compute the fourier contribution.
      do nn=1,degreefourrot
        val = nn*omegafourrot
        dphi = dphi - val*coscoeffourrot(nn)*sin(val*t)
        dphi = dphi + val*sincoeffourrot(nn)*cos(val*t)
      end do
! set derivativerigidrotangle to dphi. multiply by timeref
! to obtain the correct non-dimensional value.
      derivativerigidrotangle_d = dphi*timerefd
      derivativerigidrotangle = timeref*dphi
    end if
  end function derivativerigidrotangle_d
 
  function derivativerigidrotangle(degreepolrot, coefpolrot, &
&   degreefourrot, omegafourrot, coscoeffourrot, sincoeffourrot, t)
!
!       derivativerigidrotangle computes the time derivative of the
!       rigid body rotation angle at the given time for the given
!       arguments. the angle is described by a combination of a
!       polynomial and fourier series.
!
    use constants
    use inputphysics, only : equationmode
    use flowvarrefstate, only : timeref
    implicit none
!
!      function type
!
    real(kind=realtype) :: derivativerigidrotangle
!
!      function arguments.
!
    integer(kind=inttype), intent(in) :: degreepolrot
    integer(kind=inttype), intent(in) :: degreefourrot
    real(kind=realtype), intent(in) :: omegafourrot, t
    real(kind=realtype), dimension(0:*), intent(in) :: coefpolrot
    real(kind=realtype), dimension(0:*), intent(in) :: coscoeffourrot
    real(kind=realtype), dimension(*), intent(in) :: sincoeffourrot
!
!      local variables.
!
    integer(kind=inttype) :: nn
    real(kind=realtype) :: dphi, val
    intrinsic sin
    intrinsic cos
! return immediately if this is a steady computation.
    if (equationmode .eq. steady) then
      derivativerigidrotangle = zero
      return
    else
! compute the polynomial contribution.
      dphi = zero
      do nn=1,degreepolrot
        dphi = dphi + nn*coefpolrot(nn)*t**(nn-1)
      end do
! compute the fourier contribution.
      do nn=1,degreefourrot
        val = nn*omegafourrot
        dphi = dphi - val*coscoeffourrot(nn)*sin(val*t)
        dphi = dphi + val*sincoeffourrot(nn)*cos(val*t)
      end do
! set derivativerigidrotangle to dphi. multiply by timeref
! to obtain the correct non-dimensional value.
      derivativerigidrotangle = timeref*dphi
    end if
  end function derivativerigidrotangle
 
!  differentiation of mydim in forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: mydim
!   with respect to varying inputs: x y
  real(kind=realtype) function mydim_d(x, xd, y, yd, mydim)
    use constants
    implicit none
    real(kind=realtype) :: x, y
    real(kind=realtype) :: xd, yd
    real(kind=realtype) :: mydim
    mydim_d = xd - yd
    mydim = x - y
    if (mydim .lt. 0.0) then
      mydim = 0.0
      mydim_d = 0.0_8
    end if
  end function mydim_d
 
  function mydim(x, y)
    use constants
    implicit none
    real(kind=realtype) :: x, y
    real(kind=realtype) :: mydim
    mydim = x - y
    if (mydim .lt. 0.0) mydim = 0.0
  end function mydim
 
  function getcorrectfork()
    use constants
    use flowvarrefstate, only : kpresent
    use iteration, only : currentlevel, groundlevel
    implicit none
    logical :: getcorrectfork
    if (kpresent .and. currentlevel .le. groundlevel) then
      getcorrectfork = .true.
    else
      getcorrectfork = .false.
    end if
  end function getcorrectfork
 
  subroutine terminate(routinename, errormessage)
!
!       terminate writes an error message to standard output and
!       terminates the execution of the program.
!
    use constants
    use communication, only : adflow_comm_world, myid
    implicit none
!
!      subroutine arguments
!
    character(len=*), intent(in) :: routinename
    character(len=*), intent(in) :: errormessage
  end subroutine terminate
 
  subroutine rotmatrixrigidbody(tnew, told, rotationmatrix, &
&   rotationpoint)
!
!       rotmatrixrigidbody determines the rotation matrix and the
!       rotation point to determine the coordinates of the new time
!       level starting from the coordinates of the old time level.
!
    use constants
    use inputmotion
    use flowvarrefstate, only : lref
    implicit none
!
!      subroutine arguments.
!
    real(kind=realtype), intent(in) :: tnew, told
    real(kind=realtype), dimension(3), intent(out) :: rotationpoint
    real(kind=realtype), dimension(3, 3), intent(out) :: rotationmatrix
!
!      local variables.
!
    integer(kind=inttype) :: i, j
    real(kind=realtype) :: phi
    real(kind=realtype) :: cosx, cosy, cosz, sinx, siny, sinz
    real(kind=realtype), dimension(3, 3) :: mnew, mold
    intrinsic sin
    intrinsic cos
! determine the rotation angle around the x-axis for the new
! time level and the corresponding values of the sine and cosine.
    phi = rigidrotangle(degreepolxrot, coefpolxrot, degreefourxrot, &
&     omegafourxrot, coscoeffourxrot, sincoeffourxrot, tnew)
    sinx = sin(phi)
    cosx = cos(phi)
! idem for the y-axis.
    phi = rigidrotangle(degreepolyrot, coefpolyrot, degreefouryrot, &
&     omegafouryrot, coscoeffouryrot, sincoeffouryrot, tnew)
    siny = sin(phi)
    cosy = cos(phi)
! idem for the z-axis.
    phi = rigidrotangle(degreepolzrot, coefpolzrot, degreefourzrot, &
&     omegafourzrot, coscoeffourzrot, sincoeffourzrot, tnew)
    sinz = sin(phi)
    cosz = cos(phi)
! construct the transformation matrix at the new time level.
! it is assumed that the sequence of rotation is first around the
! x-axis then around the y-axis and finally around the z-axis.
    mnew(1, 1) = cosy*cosz
    mnew(2, 1) = cosy*sinz
    mnew(3, 1) = -siny
    mnew(1, 2) = sinx*siny*cosz - cosx*sinz
    mnew(2, 2) = sinx*siny*sinz + cosx*cosz
    mnew(3, 2) = sinx*cosy
    mnew(1, 3) = cosx*siny*cosz + sinx*sinz
    mnew(2, 3) = cosx*siny*sinz - sinx*cosz
    mnew(3, 3) = cosx*cosy
! determine the rotation angle around the x-axis for the old
! time level and the corresponding values of the sine and cosine.
    phi = rigidrotangle(degreepolxrot, coefpolxrot, degreefourxrot, &
&     omegafourxrot, coscoeffourxrot, sincoeffourxrot, told)
    sinx = sin(phi)
    cosx = cos(phi)
! idem for the y-axis.
    phi = rigidrotangle(degreepolyrot, coefpolyrot, degreefouryrot, &
&     omegafouryrot, coscoeffouryrot, sincoeffouryrot, told)
    siny = sin(phi)
    cosy = cos(phi)
! idem for the z-axis.
    phi = rigidrotangle(degreepolzrot, coefpolzrot, degreefourzrot, &
&     omegafourzrot, coscoeffourzrot, sincoeffourzrot, told)
    sinz = sin(phi)
    cosz = cos(phi)
! construct the transformation matrix at the old time level.
    mold(1, 1) = cosy*cosz
    mold(2, 1) = cosy*sinz
    mold(3, 1) = -siny
    mold(1, 2) = sinx*siny*cosz - cosx*sinz
    mold(2, 2) = sinx*siny*sinz + cosx*cosz
    mold(3, 2) = sinx*cosy
    mold(1, 3) = cosx*siny*cosz + sinx*sinz
    mold(2, 3) = cosx*siny*sinz - sinx*cosz
    mold(3, 3) = cosx*cosy
! construct the transformation matrix between the new and the
! old time level. this is mnew*inverse(mold). however the
! inverse of mold is the transpose.
    do j=1,3
      do i=1,3
        rotationmatrix(i, j) = mnew(i, 1)*mold(j, 1) + mnew(i, 2)*mold(j&
&         , 2) + mnew(i, 3)*mold(j, 3)
      end do
    end do
! determine the rotation point at the old time level; it is
! possible that this value changes due to translation of the grid.
!  ainf = sqrt(gammainf*pinf/rhoinf)
!  rotationpoint(1) = lref*rotpoint(1) &
!                   + machgrid(1)*ainf*told/timeref
!  rotationpoint(2) = lref*rotpoint(2) &
!                   + machgrid(2)*ainf*told/timeref
!  rotationpoint(3) = lref*rotpoint(3) &
!                   + machgrid(3)*ainf*told/timeref
    rotationpoint(1) = lref*rotpoint(1)
    rotationpoint(2) = lref*rotpoint(2)
    rotationpoint(3) = lref*rotpoint(3)
  end subroutine rotmatrixrigidbody
 
  function secondderivativerigidrotangle(degreepolrot, coefpolrot, &
&   degreefourrot, omegafourrot, coscoeffourrot, sincoeffourrot, t)
!
!       2ndderivativerigidrotangle computes the 2nd time derivative of
!       the rigid body rotation angle at the given time for the given
!       arguments. the angle is described by a combination of a
!       polynomial and fourier series.
!
    use constants
    use flowvarrefstate, only : timeref
    use inputphysics, only : equationmode
    implicit none
!
!      function type
!
    real(kind=realtype) :: secondderivativerigidrotangle
!
!      function arguments.
!
    integer(kind=inttype), intent(in) :: degreepolrot
    integer(kind=inttype), intent(in) :: degreefourrot
    real(kind=realtype), intent(in) :: omegafourrot, t
    real(kind=realtype), dimension(0:*), intent(in) :: coefpolrot
    real(kind=realtype), dimension(0:*), intent(in) :: coscoeffourrot
    real(kind=realtype), dimension(*), intent(in) :: sincoeffourrot
!
!      local variables.
!
    integer(kind=inttype) :: nn
    real(kind=realtype) :: dphi, val
    intrinsic sin
    intrinsic cos
! return immediately if this is a steady computation.
    if (equationmode .eq. steady) then
      secondderivativerigidrotangle = zero
      return
    else
! compute the polynomial contribution.
      dphi = zero
      do nn=2,degreepolrot
        dphi = dphi + (nn-1)*nn*coefpolrot(nn)*t**(nn-2)
      end do
! compute the fourier contribution.
      do nn=1,degreefourrot
        val = nn*omegafourrot
        dphi = dphi - val**2*sincoeffourrot(nn)*sin(val*t)
        dphi = dphi - val**2*coscoeffourrot(nn)*cos(val*t)
      end do
! set derivativerigidrotangle to dphi. multiply by timeref
! to obtain the correct non-dimensional value.
      secondderivativerigidrotangle = timeref**2*dphi
    end if
  end function secondderivativerigidrotangle
 
  function rigidrotangle(degreepolrot, coefpolrot, degreefourrot, &
&   omegafourrot, coscoeffourrot, sincoeffourrot, t)
!
!       rigidrotangle computes the rigid body rotation angle at the
!       given time for the given arguments. the angle is described by
!       a combination of a polynomial and fourier series.
!
    use constants
    use inputphysics, only : equationmode
    implicit none
!
!      function type
!
    real(kind=realtype) :: rigidrotangle
!
!      function arguments.
!
    integer(kind=inttype), intent(in) :: degreepolrot
    integer(kind=inttype), intent(in) :: degreefourrot
    real(kind=realtype), intent(in) :: omegafourrot, t
    real(kind=realtype), dimension(0:*), intent(in) :: coefpolrot
    real(kind=realtype), dimension(0:*), intent(in) :: coscoeffourrot
    real(kind=realtype), dimension(*), intent(in) :: sincoeffourrot
!
!      local variables.
!
    integer(kind=inttype) :: nn
    real(kind=realtype) :: phi, val
    intrinsic cos
    intrinsic sin
! return immediately if this is a steady computation.
    if (equationmode .eq. steady) then
      rigidrotangle = zero
      return
    else
! compute the polynomial contribution. if no polynomial was
! specified, the value of index 0 is set to zero automatically.
      phi = coefpolrot(0)
      do nn=1,degreepolrot
        phi = phi + coefpolrot(nn)*t**nn
      end do
! compute the fourier contribution. again the cosine coefficient
! of index 0 is defaulted to zero if not specified.
      phi = phi + coscoeffourrot(0)
      do nn=1,degreefourrot
        val = nn*omegafourrot*t
        phi = phi + coscoeffourrot(nn)*cos(val) + sincoeffourrot(nn)*sin&
&         (val)
      end do
! set rigidrotangle to phi.
      rigidrotangle = phi
    end if
  end function rigidrotangle
 
  subroutine setbcpointers(nn, spatialpointers)
!
!       setbcpointers sets the pointers needed for the boundary
!       condition treatment on a general face, such that the boundary
!       routines are only implemented once instead of 6 times.
!
    use constants
    use blockpointers, only : w, p, rlv, rev, gamma, x, d2wall, si, sj&
&   , sk, s, globalcell, bcdata, nx, il, ie, ib, ny, jl, je, jb, nz, kl,&
&   ke, kb, bcfaceid, addgridvelocities, sfacei, sfacej, sfacek, &
&   addgridvelocities
    use bcpointers_d, only : ww0, ww1, ww2, ww3, pp0, pp1, pp2, pp3, &
&   rlv0, rlv1, rlv2, rlv3, rev0, rev1, rev2, rev3, gamma0, gamma1, &
&   gamma2, gamma3, gcp, xx, ss, ssi, ssj, ssk, dd2wall, sface, istart, &
&   iend, jstart, jend, isize, jsize
    use inputphysics, only : cpmodel, equations
    implicit none
! subroutine arguments.
    integer(kind=inttype), intent(in) :: nn
    logical, intent(in) :: spatialpointers
! determine the sizes of each face and point to just the range we
! need on each face.
    istart = bcdata(nn)%icbeg
    iend = bcdata(nn)%icend
    jstart = bcdata(nn)%jcbeg
    jend = bcdata(nn)%jcend
! set the size of the subface
    isize = iend - istart + 1
    jsize = jend - jstart + 1
! determine the face id on which the subface is located and set
! the pointers accordinly.
    select case  (bcfaceid(nn)) 
    case (imin) 
!---------------------------------------------------------------------------
      ww3 => w(3, 1:, 1:, :)
      ww2 => w(2, 1:, 1:, :)
      ww1 => w(1, 1:, 1:, :)
      ww0 => w(0, 1:, 1:, :)
      pp3 => p(3, 1:, 1:)
      pp2 => p(2, 1:, 1:)
      pp1 => p(1, 1:, 1:)
      pp0 => p(0, 1:, 1:)
      rlv3 => rlv(3, 1:, 1:)
      rlv2 => rlv(2, 1:, 1:)
      rlv1 => rlv(1, 1:, 1:)
      rlv0 => rlv(0, 1:, 1:)
      rev3 => rev(3, 1:, 1:)
      rev2 => rev(2, 1:, 1:)
      rev1 => rev(1, 1:, 1:)
      rev0 => rev(0, 1:, 1:)
      gamma3 => gamma(3, 1:, 1:)
      gamma2 => gamma(2, 1:, 1:)
      gamma1 => gamma(1, 1:, 1:)
      gamma0 => gamma(0, 1:, 1:)
      gcp => globalcell(2, 1:, 1:)
!---------------------------------------------------------------------------
    case (imax) 
      ww3 => w(nx, 1:, 1:, :)
      ww2 => w(il, 1:, 1:, :)
      ww1 => w(ie, 1:, 1:, :)
      ww0 => w(ib, 1:, 1:, :)
      pp3 => p(nx, 1:, 1:)
      pp2 => p(il, 1:, 1:)
      pp1 => p(ie, 1:, 1:)
      pp0 => p(ib, 1:, 1:)
      rlv3 => rlv(nx, 1:, 1:)
      rlv2 => rlv(il, 1:, 1:)
      rlv1 => rlv(ie, 1:, 1:)
      rlv0 => rlv(ib, 1:, 1:)
      rev3 => rev(nx, 1:, 1:)
      rev2 => rev(il, 1:, 1:)
      rev1 => rev(ie, 1:, 1:)
      rev0 => rev(ib, 1:, 1:)
      gamma3 => gamma(nx, 1:, 1:)
      gamma2 => gamma(il, 1:, 1:)
      gamma1 => gamma(ie, 1:, 1:)
      gamma0 => gamma(ib, 1:, 1:)
      gcp => globalcell(il, 1:, 1:)
!---------------------------------------------------------------------------
    case (jmin) 
      ww3 => w(1:, 3, 1:, :)
      ww2 => w(1:, 2, 1:, :)
      ww1 => w(1:, 1, 1:, :)
      ww0 => w(1:, 0, 1:, :)
      pp3 => p(1:, 3, 1:)
      pp2 => p(1:, 2, 1:)
      pp1 => p(1:, 1, 1:)
      pp0 => p(1:, 0, 1:)
      rlv3 => rlv(1:, 3, 1:)
      rlv2 => rlv(1:, 2, 1:)
      rlv1 => rlv(1:, 1, 1:)
      rlv0 => rlv(1:, 0, 1:)
      rev3 => rev(1:, 3, 1:)
      rev2 => rev(1:, 2, 1:)
      rev1 => rev(1:, 1, 1:)
      rev0 => rev(1:, 0, 1:)
      gamma3 => gamma(1:, 3, 1:)
      gamma2 => gamma(1:, 2, 1:)
      gamma1 => gamma(1:, 1, 1:)
      gamma0 => gamma(1:, 0, 1:)
      gcp => globalcell(1:, 2, 1:)
!---------------------------------------------------------------------------
    case (jmax) 
      ww3 => w(1:, ny, 1:, :)
      ww2 => w(1:, jl, 1:, :)
      ww1 => w(1:, je, 1:, :)
      ww0 => w(1:, jb, 1:, :)
      pp3 => p(1:, ny, 1:)
      pp2 => p(1:, jl, 1:)
      pp1 => p(1:, je, 1:)
      pp0 => p(1:, jb, 1:)
      rlv3 => rlv(1:, ny, 1:)
      rlv2 => rlv(1:, jl, 1:)
      rlv1 => rlv(1:, je, 1:)
      rlv0 => rlv(1:, jb, 1:)
      rev3 => rev(1:, ny, 1:)
      rev2 => rev(1:, jl, 1:)
      rev1 => rev(1:, je, 1:)
      rev0 => rev(1:, jb, 1:)
      gamma3 => gamma(1:, ny, 1:)
      gamma2 => gamma(1:, jl, 1:)
      gamma1 => gamma(1:, je, 1:)
      gamma0 => gamma(1:, jb, 1:)
      gcp => globalcell(1:, jl, 1:)
!---------------------------------------------------------------------------
    case (kmin) 
      ww3 => w(1:, 1:, 3, :)
      ww2 => w(1:, 1:, 2, :)
      ww1 => w(1:, 1:, 1, :)
      ww0 => w(1:, 1:, 0, :)
      pp3 => p(1:, 1:, 3)
      pp2 => p(1:, 1:, 2)
      pp1 => p(1:, 1:, 1)
      pp0 => p(1:, 1:, 0)
      rlv3 => rlv(1:, 1:, 3)
      rlv2 => rlv(1:, 1:, 2)
      rlv1 => rlv(1:, 1:, 1)
      rlv0 => rlv(1:, 1:, 0)
      rev3 => rev(1:, 1:, 3)
      rev2 => rev(1:, 1:, 2)
      rev1 => rev(1:, 1:, 1)
      rev0 => rev(1:, 1:, 0)
      gamma3 => gamma(1:, 1:, 3)
      gamma2 => gamma(1:, 1:, 2)
      gamma1 => gamma(1:, 1:, 1)
      gamma0 => gamma(1:, 1:, 0)
      gcp => globalcell(1:, 1:, 2)
!---------------------------------------------------------------------------
    case (kmax) 
      ww3 => w(1:, 1:, nz, :)
      ww2 => w(1:, 1:, kl, :)
      ww1 => w(1:, 1:, ke, :)
      ww0 => w(1:, 1:, kb, :)
      pp3 => p(1:, 1:, nz)
      pp2 => p(1:, 1:, kl)
      pp1 => p(1:, 1:, ke)
      pp0 => p(1:, 1:, kb)
      rlv3 => rlv(1:, 1:, nz)
      rlv2 => rlv(1:, 1:, kl)
      rlv1 => rlv(1:, 1:, ke)
      rlv0 => rlv(1:, 1:, kb)
      rev3 => rev(1:, 1:, nz)
      rev2 => rev(1:, 1:, kl)
      rev1 => rev(1:, 1:, ke)
      rev0 => rev(1:, 1:, kb)
      gamma3 => gamma(1:, 1:, nz)
      gamma2 => gamma(1:, 1:, kl)
      gamma1 => gamma(1:, 1:, ke)
      gamma0 => gamma(1:, 1:, kb)
      gcp => globalcell(1:, 1:, kl)
    end select
    if (spatialpointers) then
      select case  (bcfaceid(nn)) 
      case (imin) 
        xx => x(1, :, :, :)
        ssi => si(1, :, :, :)
        ssj => sj(2, :, :, :)
        ssk => sk(2, :, :, :)
        ss => s(2, :, :, :)
      case (imax) 
        xx => x(il, :, :, :)
        ssi => si(il, :, :, :)
        ssj => sj(il, :, :, :)
        ssk => sk(il, :, :, :)
        ss => s(il, :, :, :)
      case (jmin) 
        xx => x(:, 1, :, :)
        ssi => sj(:, 1, :, :)
        ssj => si(:, 2, :, :)
        ssk => sk(:, 2, :, :)
        ss => s(:, 2, :, :)
      case (jmax) 
        xx => x(:, jl, :, :)
        ssi => sj(:, jl, :, :)
        ssj => si(:, jl, :, :)
        ssk => sk(:, jl, :, :)
        ss => s(:, jl, :, :)
      case (kmin) 
        xx => x(:, :, 1, :)
        ssi => sk(:, :, 1, :)
        ssj => si(:, :, 2, :)
        ssk => sj(:, :, 2, :)
        ss => s(:, :, 2, :)
      case (kmax) 
        xx => x(:, :, kl, :)
        ssi => sk(:, :, kl, :)
        ssj => si(:, :, kl, :)
        ssk => sj(:, :, kl, :)
        ss => s(:, :, kl, :)
      end select
      if (addgridvelocities) then
        select case  (bcfaceid(nn)) 
        case (imin) 
          sface => sfacei(1, :, :)
        case (imax) 
          sface => sfacei(il, :, :)
        case (jmin) 
          sface => sfacej(:, 1, :)
        case (jmax) 
          sface => sfacej(:, jl, :)
        case (kmin) 
          sface => sfacek(:, :, 1)
        case (kmax) 
          sface => sfacek(:, :, kl)
        end select
      end if
      if (equations .eq. ransequations) then
        select case  (bcfaceid(nn)) 
        case (imin) 
          dd2wall => d2wall(2, :, :)
        case (imax) 
          dd2wall => d2wall(il, :, :)
        case (jmin) 
          dd2wall => d2wall(:, 2, :)
        case (jmax) 
          dd2wall => d2wall(:, jl, :)
        case (kmin) 
          dd2wall => d2wall(:, :, 2)
        case (kmax) 
          dd2wall => d2wall(:, :, kl)
        end select
      end if
    end if
  end subroutine setbcpointers
 
  subroutine computerootbendingmoment(cf, cm, bendingmoment)
!                                                      *
! compute a normalized bending moment coefficient from *
! the force and moment coefficient. at the moment this *
! routine only works for a half body. additional logic *
! would be needed for a full body.                     *
!                                                      *
    use constants
    use inputphysics, only : lengthref, pointref, pointrefec, &
&   liftindex
    implicit none
!input/output variables
    real(kind=realtype), dimension(3), intent(in) :: cf, cm
    real(kind=realtype), intent(out) :: bendingmoment
!subroutine variables
    real(kind=realtype) :: elasticmomentx, elasticmomenty, &
&   elasticmomentz
    intrinsic sqrt
    real(kind=realtype) :: arg1
    bendingmoment = zero
    if (liftindex .eq. 2) then
!z out wing sum momentx,momentz
      elasticmomentx = cm(1) + cf(2)*(pointrefec(3)-pointref(3))/&
&       lengthref - cf(3)*(pointrefec(2)-pointref(2))/lengthref
      elasticmomentz = cm(3) - cf(2)*(pointrefec(1)-pointref(1))/&
&       lengthref + cf(1)*(pointrefec(2)-pointref(2))/lengthref
      arg1 = elasticmomentx**2 + elasticmomentz**2
      bendingmoment = sqrt(arg1)
    else if (liftindex .eq. 3) then
!y out wing sum momentx,momenty
      elasticmomentx = cm(1) + cf(3)*(pointrefec(2)-pointref(2))/&
&       lengthref + cf(3)*(pointrefec(3)-pointref(3))/lengthref
      elasticmomenty = cm(2) + cf(3)*(pointrefec(1)-pointref(1))/&
&       lengthref + cf(1)*(pointrefec(3)-pointref(3))/lengthref
      arg1 = elasticmomentx**2 + elasticmomenty**2
      bendingmoment = sqrt(arg1)
    end if
  end subroutine computerootbendingmoment
 
  subroutine computeleastsquaresregression(y, x, npts, m, b)
!
!       computes the slope of best fit for a set of x,y data of length
!       npts
!
    use constants
    implicit none
!subroutine arguments
    integer(kind=inttype) :: npts
    real(kind=realtype), dimension(npts) :: x, y
    real(kind=realtype) :: m, b
!local variables
    real(kind=realtype) :: sumx, sumy, sumx2, sumxy
    integer(kind=inttype) :: i
!begin execution
    sumx = 0.0
    sumy = 0.0
    sumx2 = 0.0
    sumxy = 0.0
    do i=1,npts
      sumx = sumx + x(i)
      sumy = sumy + y(i)
      sumx2 = sumx2 + x(i)*x(i)
      sumxy = sumxy + x(i)*y(i)
    end do
    m = (npts*sumxy-sumy*sumx)/(npts*sumx2-sumx**2)
    b = (sumy*sumx2-sumx*sumxy)/(npts*sumx2-sumx**2)
  end subroutine computeleastsquaresregression
 
  subroutine computetsderivatives(force, moment, coef0, dcdalpha, &
&   dcdalphadot, dcdq, dcdqdot)
!
!      computes the stability derivatives based on the time spectral
!      solution of a given mesh. takes in the force coefficients at
!      all time instantces and computes the agregate parameters
!
    use constants
    use communication
    use inputphysics
    use inputtimespectral
    use inputtsstabderiv
    use flowvarrefstate
    use monitor
    use section
    use inputmotion
    implicit none
!
!     subroutine arguments.
!
    real(kind=realtype), dimension(3, ntimeintervalsspectral) :: force, &
&   moment
    real(kind=realtype), dimension(8) :: dcdq, dcdqdot
    real(kind=realtype), dimension(8) :: dcdalpha, dcdalphadot
    real(kind=realtype), dimension(8) :: coef0
! working variables
    real(kind=realtype), dimension(ntimeintervalsspectral, 8) :: &
&   basecoef
    real(kind=realtype), dimension(8) :: coef0dot
    real(kind=realtype), dimension(ntimeintervalsspectral, 8) :: &
&   resbasecoef
    real(kind=realtype), dimension(ntimeintervalsspectral) :: &
&   intervalalpha, intervalalphadot
    real(kind=realtype), dimension(ntimeintervalsspectral) :: &
&   intervalmach, intervalmachdot
    real(kind=realtype), dimension(nsections) :: t
    integer(kind=inttype) :: i, sps, nn
!speed of sound: for normalization of q derivatives
    real(kind=realtype) :: a
    real(kind=realtype) :: fact, factmoment
! functions
    real(kind=realtype), dimension(ntimeintervalsspectral) :: dphix, &
&   dphiy, dphiz
    real(kind=realtype), dimension(ntimeintervalsspectral) :: dphixdot, &
&   dphiydot, dphizdot
    real(kind=realtype) :: derivativerigidrotangle, &
&   secondderivativerigidrotangle
    intrinsic sqrt
    real(kind=realtype) :: arg1
    fact = two/(gammainf*pinf*machcoef**2*surfaceref*lref**2)
    factmoment = fact/(lengthref*lref)
    if (tsqmode) then
      print*, &
&     'ts q mode code needs to be updated in computetsderivatives!'
      stop
! !q is pitch
! do sps =1,ntimeintervalsspectral
!    !compute the time of this intervavc
!    t = timeunsteadyrestart
!    if(equationmode == timespectral) then
!       do nn=1,nsections
!          t(nn) = t(nn) + (sps-1)*sections(nn)%timeperiod &
!               /         (ntimeintervalsspectral*1.0)
!       enddo
!    endif
!    ! compute the time derivative of the rotation angles around the
!    ! z-axis. i.e. compute q
!    dphiz(sps) = derivativerigidrotangle(degreepolzrot,   &
!         coefpolzrot,     &
!         degreefourzrot,  &
!         omegafourzrot,   &
!         coscoeffourzrot, &
!         sincoeffourzrot, t)
!    ! add in q_dot computation
!    dphizdot(sps) = secondderivativerigidrotangle(degreepolzrot,   &
!         coefpolzrot,     &
!         degreefourzrot,  &
!         omegafourzrot,   &
!         coscoeffourzrot, &
!         sincoeffourzrot, t)
! end do
! !now compute dcl/dq
! do i =1,8
!    call computeleastsquaresregression(basecoef(:,i),dphiz,ntimeintervalsspectral,dcdq(i),coef0(i))
! end do
! ! now subtract off estimated cl,cmz and use remainder to compute
! ! clqdot and cmzqdot.
! do i = 1,8
!    do sps = 1,ntimeintervalsspectral
!       resbasecoef(sps,i) = basecoef(sps,i)-(dcdq(i)*dphiz(sps)+coef0(i))
!    enddo
! enddo
! !now normalize the results...
! a  = sqrt(gammainf*pinfdim/rhoinfdim)
! dcdq = dcdq*timeref*2*(machgrid*a)/lengthref
! !now compute dcl/dpdot
! do i = 1,8
!    call computeleastsquaresregression(resbasecoef(:,i),dphizdot,ntimeintervalsspectral,dcdqdot(i),coef0dot(i))
! enddo
    else if (tsalphamode) then
      do sps=1,ntimeintervalsspectral
!compute the time of this interval
        t = timeunsteadyrestart
        if (equationmode .eq. timespectral) then
          do nn=1,nsections
            t(nn) = t(nn) + (sps-1)*sections(nn)%timeperiod/(&
&             ntimeintervalsspectral*1.0)
          end do
        end if
        intervalalpha(sps) = tsalpha(degreepolalpha, coefpolalpha, &
&         degreefouralpha, omegafouralpha, coscoeffouralpha, &
&         sincoeffouralpha, t(1))
        intervalalphadot(sps) = tsalphadot(degreepolalpha, coefpolalpha&
&         , degreefouralpha, omegafouralpha, coscoeffouralpha, &
&         sincoeffouralpha, t(1))
! this call is wrong!!!!
!call getdirangle(veldirfreestream,liftdirection,liftindex,alpha+intervalalpha(sps), beta)
        basecoef(sps, 1) = fact*(force(1, sps)*liftdirection(1)+force(2&
&         , sps)*liftdirection(2)+force(3, sps)*liftdirection(3))
        basecoef(sps, 2) = fact*(force(1, sps)*dragdirection(1)+force(2&
&         , sps)*dragdirection(2)+force(3, sps)*dragdirection(3))
        basecoef(sps, 3) = force(1, sps)*fact
        basecoef(sps, 4) = force(2, sps)*fact
        basecoef(sps, 5) = force(3, sps)*fact
        basecoef(sps, 6) = moment(1, sps)*factmoment
        basecoef(sps, 7) = moment(2, sps)*factmoment
        basecoef(sps, 8) = moment(3, sps)*factmoment
      end do
!now compute dcl/dalpha
      do i=1,8
        call computeleastsquaresregression(basecoef(:, i), intervalalpha&
&                                    , ntimeintervalsspectral, dcdalpha(&
&                                    i), coef0(i))
      end do
! now subtract off estimated cl,cmz and use remainder to compute
! clalphadot and cmzalphadot.
      do i=1,8
        do sps=1,ntimeintervalsspectral
          resbasecoef(sps, i) = basecoef(sps, i) - (dcdalpha(i)*&
&           intervalalpha(sps)+coef0(i))
        end do
      end do
!now compute dci/dalphadot
      do i=1,8
        call computeleastsquaresregression(resbasecoef(:, i), &
&                                    intervalalphadot, &
&                                    ntimeintervalsspectral, dcdalphadot&
&                                    (i), coef0dot(i))
      end do
      arg1 = gammainf*pinfdim/rhoinfdim
      a = sqrt(arg1)
      dcdalphadot = dcdalphadot*2*(machgrid*a)/lengthref
    else
      call terminate('computetsderivatives', &
&              'not a valid stability motion')
    end if
  end subroutine computetsderivatives
 
  subroutine getdirangle(freestreamaxis, liftaxis, liftindex, alpha, &
&   beta)
!
!      convert the wind axes to angle of attack and side slip angle.
!      the direction angles alpha and beta are computed given the
!      components of the wind direction vector (freestreamaxis), the
!      lift direction vector (liftaxis) and assuming that the
!      body direction (xb,yb,zb) is in the default ijk coordinate
!      system. the rotations are determined by first determining
!      whether the lift is primarily in the j or k direction and then
!      determining the angles accordingly.
!      direction vector:
!        1) rotation about the zb or yb -axis: alpha clockwise (cw)
!           (xb,yb,zb) -> (x1,y1,z1)
!        2) rotation about the yl or z1 -axis: beta counter-clockwise
!           (ccw) (x1,y1,z1) -> (xw,yw,zw)
!         input arguments:
!            freestreamaxis = wind vector in body axes
!            liftaxis       = lift direction vector in body axis
!         output arguments:
!            alpha    = angle of attack in radians
!            beta     = side slip angle in radians
!
    use constants
    implicit none
!
!     subroutine arguments.
!
!      real(kind=realtype), intent(in)  :: xw, yw, zw
    real(kind=realtype), dimension(3), intent(in) :: freestreamaxis
    real(kind=realtype), dimension(3), intent(in) :: liftaxis
    real(kind=realtype), intent(out) :: alpha, beta
    integer(kind=inttype), intent(out) :: liftindex
!
!     local variables.
!
    real(kind=realtype) :: rnorm
    integer(kind=inttype) :: flowindex, i
    real(kind=realtype), dimension(3) :: freestreamaxisnorm
    integer(kind=inttype) :: temp
    intrinsic abs
    intrinsic sqrt
    intrinsic asin
    intrinsic atan2
    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) :: abs6
    real(kind=realtype) :: abs7
    real(kind=realtype) :: arg1
! assume domoniate flow is x
    flowindex = 1
    if (liftaxis(1) .ge. 0.) then
      abs0 = liftaxis(1)
    else
      abs0 = -liftaxis(1)
    end if
    if (liftaxis(2) .ge. 0.) then
      abs2 = liftaxis(2)
    else
      abs2 = -liftaxis(2)
    end if
    if (liftaxis(1) .ge. 0.) then
      abs4 = liftaxis(1)
    else
      abs4 = -liftaxis(1)
    end if
    if (liftaxis(3) .ge. 0.) then
      abs6 = liftaxis(3)
    else
      abs6 = -liftaxis(3)
    end if
! determine the dominant lift direction
    if (abs0 .gt. abs2 .and. abs4 .gt. abs6) then
      temp = 1
    else
      if (liftaxis(2) .ge. 0.) then
        abs1 = liftaxis(2)
      else
        abs1 = -liftaxis(2)
      end if
      if (liftaxis(1) .ge. 0.) then
        abs3 = liftaxis(1)
      else
        abs3 = -liftaxis(1)
      end if
      if (liftaxis(2) .ge. 0.) then
        abs5 = liftaxis(2)
      else
        abs5 = -liftaxis(2)
      end if
      if (liftaxis(3) .ge. 0.) then
        abs7 = liftaxis(3)
      else
        abs7 = -liftaxis(3)
      end if
      if (abs1 .gt. abs3 .and. abs5 .gt. abs7) then
        temp = 2
      else
        temp = 3
      end if
    end if
    liftindex = temp
! normalize the freestreamdirection vector.
    arg1 = freestreamaxis(1)**2 + freestreamaxis(2)**2 + freestreamaxis(&
&     3)**2
    rnorm = sqrt(arg1)
    do i=1,3
      freestreamaxisnorm(i) = freestreamaxis(i)/rnorm
    end do
    if (liftindex .eq. 2) then
! different coordinate system for aerosurf
! wing is in z- direction
! compute angle of attack alpha.
      alpha = asin(freestreamaxisnorm(2))
! compute side-slip angle beta.
      beta = -atan2(freestreamaxisnorm(3), freestreamaxisnorm(1))
    else if (liftindex .eq. 3) then
! wing is in y- direction
! compute angle of attack alpha.
      alpha = asin(freestreamaxisnorm(3))
! compute side-slip angle beta.
      beta = atan2(freestreamaxisnorm(2), freestreamaxisnorm(1))
    else
      call terminate('getdirangle', 'invalid lift direction')
    end if
  end subroutine getdirangle
 
  subroutine stabilityderivativedriver()
!
!      runs the time spectral stability derivative routines from the
!      main program file
!
    use precision
    implicit none
!
!     local variables.
!
    real(kind=realtype), dimension(8) :: dcdalpha, dcdalphadot, dcdbeta&
&   , dcdbetadot, dcdmach, dcdmachdot
    real(kind=realtype), dimension(8) :: dcdp, dcdpdot, dcdq, dcdqdot, &
&   dcdr, dcdrdot
    real(kind=realtype), dimension(8) :: coef0, coef0dot
!call computetsderivatives(coef0,dcdalpha,dcdalphadot,dcdq,dcdqdot)
  end subroutine stabilityderivativedriver
 
  subroutine setcoeftimeintegrator()
!
!       setcoeftimeintegrator determines the coefficients of the
!       time integration scheme in unsteady mode. normally these are
!       equal to the coefficients corresponding to the specified
!       accuracy. however during the initial phase there are not
!       enough states in the past and the accuracy is reduced.
!
    use constants
    use inputunsteady
    use inputphysics
    use iteration
    use monitor
    implicit none
!
!      local variables.
!
    integer(kind=inttype) :: nn, nlevelsset
! determine which time integrator must be used.
! modified by hdn
    select case  (timeaccuracy) 
    case (firstorder) 
! 1st order. no need to check the number of available
! states in the past. set the two coefficients and
! nlevelsset to 2.
      coeftime(0) = 1.0_realtype
      coeftime(1) = -1.0_realtype
      if (useale .and. equationmode .eq. unsteady) then
        coeftimeale(1) = 1.0_realtype
        coefmeshale(1, 1) = half
        coefmeshale(1, 2) = half
      end if
      nlevelsset = 2
!--------------------------------------------------
    case (secondorder) 
! second order time integrator. determine the amount of
! available states and set the coefficients accordingly.
      select case  (noldsolavail) 
      case (1_inttype) 
        coeftime(0) = 1.0_realtype
        coeftime(1) = -1.0_realtype
        if (useale .and. equationmode .eq. unsteady) then
          coeftimeale(1) = half
          coeftimeale(2) = half
          coeftimeale(3) = zero
          coeftimeale(4) = zero
          coefmeshale(1, 1) = half
          coefmeshale(1, 2) = half
          coefmeshale(2, 1) = half
          coefmeshale(2, 2) = half
        end if
        nlevelsset = 2
      case default
! 2 or bigger.
        coeftime(0) = 1.5_realtype
        coeftime(1) = -2.0_realtype
        coeftime(2) = 0.5_realtype
        if (useale .and. equationmode .eq. unsteady) then
          coeftimeale(1) = threefourth
          coeftimeale(2) = threefourth
          coeftimeale(3) = -fourth
          coeftimeale(4) = -fourth
          coefmeshale(1, 1) = half*(1.0_realtype+1.0_realtype/sqrtthree)
          coefmeshale(1, 2) = half*(1.0_realtype-1.0_realtype/sqrtthree)
          coefmeshale(2, 1) = coefmeshale(1, 2)
          coefmeshale(2, 2) = coefmeshale(1, 1)
        end if
        nlevelsset = 3
      end select
!--------------------------------------------------
    case (thirdorder) 
! third order time integrator.  determine the amount of
! available states and set the coefficients accordingly.
      select case  (noldsolavail) 
      case (1_inttype) 
        coeftime(0) = 1.0_realtype
        coeftime(1) = -1.0_realtype
        if (useale .and. equationmode .eq. unsteady) then
          coeftimeale(1) = 1.0_realtype
          coefmeshale(1, 1) = half
          coefmeshale(1, 2) = half
        end if
        nlevelsset = 2
      case (2_inttype) 
        coeftime(0) = 1.5_realtype
        coeftime(1) = -2.0_realtype
        coeftime(2) = 0.5_realtype
        if (useale .and. equationmode .eq. unsteady) then
          coeftimeale(1) = threefourth
          coeftimeale(2) = -fourth
          coefmeshale(1, 1) = half*(1.0_realtype+1.0_realtype/sqrtthree)
          coefmeshale(1, 2) = half*(1.0_realtype-1.0_realtype/sqrtthree)
          coefmeshale(2, 1) = coefmeshale(1, 2)
          coefmeshale(2, 2) = coefmeshale(1, 1)
        end if
        nlevelsset = 3
      case default
! 3 or bigger.
        coeftime(0) = 11.0_realtype/6.0_realtype
        coeftime(1) = -3.0_realtype
        coeftime(2) = 1.5_realtype
        coeftime(3) = -(1.0_realtype/3.0_realtype)
! these numbers are not correct
! do not use 3rd order ale for now
        if (useale .and. equationmode .eq. unsteady) then
          print*, 'third-order ale not implemented yet.'
          coeftimeale(1) = threefourth
          coeftimeale(2) = threefourth
          coeftimeale(3) = -fourth
          coeftimeale(4) = -fourth
          coefmeshale(1, 1) = half*(1.0_realtype+1.0_realtype/sqrtthree)
          coefmeshale(1, 2) = half*(1.0_realtype-1.0_realtype/sqrtthree)
          coefmeshale(2, 1) = coefmeshale(1, 2)
          coefmeshale(2, 2) = coefmeshale(1, 1)
          coefmeshale(3, 1) = coefmeshale(1, 2)
          coefmeshale(3, 2) = coefmeshale(1, 1)
        end if
        nlevelsset = 4
      end select
    end select
! set the rest of the coefficients to 0 if not enough states
! in the past are available.
    do nn=nlevelsset,noldlevels
      coeftime(nn) = zero
    end do
  end subroutine setcoeftimeintegrator
 
!  differentiation of mynorm2 in forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: mynorm2
!   with respect to varying inputs: x
  real(kind=realtype) function mynorm2_d(x, xd, mynorm2)
    use constants
    implicit none
    real(kind=realtype), dimension(3), intent(in) :: x
    real(kind=realtype), dimension(3), intent(in) :: xd
    real(kind=realtype) :: mynorm2
    intrinsic sqrt
    real(kind=realtype) :: arg1
    real(kind=realtype) :: arg1d
    real(kind=realtype) :: temp
    arg1d = 2*x(1)*xd(1) + 2*x(2)*xd(2) + 2*x(3)*xd(3)
    arg1 = x(1)**2 + x(2)**2 + x(3)**2
    temp = sqrt(arg1)
    if (arg1 .eq. 0.0_8) then
      mynorm2_d = 0.0_8
    else
      mynorm2_d = arg1d/(2.0*temp)
    end if
    mynorm2 = temp
  end function mynorm2_d
 
  function mynorm2(x)
    use constants
    implicit none
    real(kind=realtype), dimension(3), intent(in) :: x
    real(kind=realtype) :: mynorm2
    intrinsic sqrt
    real(kind=realtype) :: arg1
    arg1 = x(1)**2 + x(2)**2 + x(3)**2
    mynorm2 = sqrt(arg1)
  end function mynorm2
 
  function iswalltype(btype)
    use constants
    implicit none
    integer(kind=inttype) :: btype
    logical :: iswalltype
    iswalltype = .false.
    if ((btype .eq. nswalladiabatic .or. btype .eq. nswallisothermal) &
&       .or. btype .eq. eulerwall) iswalltype = .true.
  end function iswalltype
 
!  differentiation of cross_prod in forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: c
!   with respect to varying inputs: a b c
  subroutine cross_prod_d(a, ad, b, bd, c, cd)
    use precision
    implicit none
! inputs
    real(kind=realtype), dimension(3), intent(in) :: a, b
    real(kind=realtype), dimension(3), intent(in) :: ad, bd
! outputs
    real(kind=realtype), dimension(3), intent(out) :: c
    real(kind=realtype), dimension(3), intent(out) :: cd
    cd(1) = b(3)*ad(2) + a(2)*bd(3) - b(2)*ad(3) - a(3)*bd(2)
    c(1) = a(2)*b(3) - a(3)*b(2)
    cd(2) = b(1)*ad(3) + a(3)*bd(1) - b(3)*ad(1) - a(1)*bd(3)
    c(2) = a(3)*b(1) - a(1)*b(3)
    cd(3) = b(2)*ad(1) + a(1)*bd(2) - b(1)*ad(2) - a(2)*bd(1)
    c(3) = a(1)*b(2) - a(2)*b(1)
  end subroutine cross_prod_d
 
  subroutine cross_prod(a, b, c)
    use precision
    implicit none
! inputs
    real(kind=realtype), dimension(3), intent(in) :: a, b
! outputs
    real(kind=realtype), dimension(3), intent(out) :: c
    c(1) = a(2)*b(3) - a(3)*b(2)
    c(2) = a(3)*b(1) - a(1)*b(3)
    c(3) = a(1)*b(2) - a(2)*b(1)
  end subroutine cross_prod
 
  subroutine siangle(angle, mult, trans)
    use constants
    use su_cgns, only : radian, degree
    implicit none
!
!      subroutine arguments.
!
    integer, intent(in) :: angle
    real(kind=realtype), intent(out) :: mult, trans
! determine the situation we are having here.
    if (angle .eq. radian) then
! angle is already given in radians. no need for a conversion.
      mult = one
      trans = zero
    else if (angle .eq. degree) then
! angle is given in degrees. a multiplication must be performed.
      mult = pi/180.0_realtype
      trans = zero
    else
      call terminate('siangle', &
&              'no idea how to convert this to si units')
    end if
  end subroutine siangle
 
  subroutine sidensity(mass, len, mult, trans)
!
!       sidensity computes the conversion from the given density
!       unit, which can be constructed from mass and length, to the
!       si-unit kg/m^3. the conversion will look like:
!       density in kg/m^3 = mult*(density in ncu) + trans.
!       ncu means non-christian units, i.e. everything that is not si.
!
    use constants
    use su_cgns, only : kilogram, meter
    implicit none
!
!      subroutine arguments.
!
    integer, intent(in) :: mass, len
    real(kind=realtype), intent(out) :: mult, trans
! determine the situation we are having here.
    if (mass .eq. kilogram .and. len .eq. meter) then
! density is given in kg/m^3, i.e. no need for a conversion.
      mult = one
      trans = zero
    else
      call terminate('sidensity', &
&              'no idea how to convert this to si units')
    end if
  end subroutine sidensity
 
  subroutine silen(len, mult, trans)
!
!       silen computes the conversion from the given length unit to
!       the si-unit meter. the conversion will look like:
!       length in meter = mult*(length in ncu) + trans.
!       ncu means non-christian units, i.e. everything that is not si.
!
    use constants
    use su_cgns, only : meter, centimeter, millimeter, foot, inch
    implicit none
!
!      subroutine arguments.
!
    integer, intent(in) :: len
    real(kind=realtype), intent(out) :: mult, trans
! determine the situation we are having here.
    select case  (len) 
    case (meter) 
      mult = one
      trans = zero
    case (centimeter) 
      mult = 0.01_realtype
      trans = zero
    case (millimeter) 
      mult = 0.001_realtype
      trans = zero
    case (foot) 
      mult = 0.3048_realtype
      trans = zero
    case (inch) 
      mult = 0.0254_realtype
      trans = zero
    case default
      call terminate('silen', 'no idea how to convert this to si units')
    end select
  end subroutine silen
 
  subroutine sipressure(mass, len, time, mult, trans)
!
!       sipressure computes the conversion from the given pressure
!       unit, which can be constructed from mass, length and time, to
!       the si-unit pa. the conversion will look like:
!       pressure in pa = mult*(pressure in ncu) + trans.
!       ncu means non-christian units, i.e. everything that is not si.
!
    use constants
    use su_cgns, only : kilogram, meter, second
    implicit none
!
!      subroutine arguments.
!
    integer, intent(in) :: mass, len, time
    real(kind=realtype), intent(out) :: mult, trans
! determine the situation we are having here.
    if (mass .eq. kilogram .and. len .eq. meter .and. time .eq. second) &
&   then
! pressure is given in pa, i.e. no need for a conversion.
      mult = one
      trans = zero
    else
      call terminate('sipressure', &
&              'no idea how to convert this to si units')
    end if
  end subroutine sipressure
 
  subroutine sitemperature(temp, mult, trans)
!
!       sitemperature computes the conversion from the given
!       temperature unit to the si-unit kelvin. the conversion will
!       look like:
!       temperature in k = mult*(temperature in ncu) + trans.
!       ncu means non-christian units, i.e. everything that is not si.
!
    use constants
    use su_cgns, only : kelvin, celsius, rankine, fahrenheit
    implicit none
!
!      subroutine arguments.
!
    integer, intent(in) :: temp
    real(kind=realtype), intent(out) :: mult, trans
! determine the situation we are having here.
    select case  (temp) 
    case (kelvin) 
! temperature is already given in kelvin. no need to convert.
      mult = one
      trans = zero
    case (celsius) 
! is it celcius or celsius?
! temperature is in celsius. only an offset must be applied.
      mult = one
      trans = 273.16_realtype
    case (rankine) 
! temperature is in rankine. only a multiplication needs to
! be performed.
      mult = 5.0_realtype/9.0_realtype
      trans = zero
    case (fahrenheit) 
! temperature is in fahrenheit. both a multiplication and an
! offset must be applied.
      mult = 5.0_realtype/9.0_realtype
      trans = 255.382
    case default
! unknown temperature unit.
      call terminate('sitemperature', &
&              'no idea how to convert this to si units')
    end select
  end subroutine sitemperature
 
  subroutine siturb(mass, len, time, temp, turbname, mult, trans)
!
!       siturb computes the conversion from the given turbulence
!       unit, which can be constructed from mass, len, time and temp,
!       to the si-unit for the given variable. the conversion will
!       look like: var in si = mult*(var in ncu) + trans.
!       ncu means non-christian units, i.e. everything that is not si.
!
    use constants
    use su_cgns, only : kilogram, meter, second, kelvin
    implicit none
!
!      subroutine arguments.
!
    integer, intent(in) :: mass, len, time, temp
    character(len=*), intent(in) :: turbname
    real(kind=realtype), intent(out) :: mult, trans
! determine the situation we are having here.
    if (mass .eq. kilogram .and. len .eq. meter .and. time .eq. second &
&       .and. temp .eq. kelvin) then
! everthing is already in si units. no conversion needed.
      mult = one
      trans = zero
    else
      call terminate('siturb', 'no idea how to convert this to si units'&
&             )
    end if
  end subroutine siturb
 
  subroutine sivelocity(length, time, mult, trans)
!
!       sivelocity computes the conversion from the given velocity
!       unit, which can be constructed from length and time, to the
!       si-unit m/s. the conversion will look like:
!       velocity in m/s = mult*(velocity in ncu) + trans.
!       ncu means non-christian units, i.e. everything that is not si.
!
    use constants
    use su_cgns, only : meter, centimeter, millimeter, foot, inch, &
&   second
    implicit none
!
!      subroutine arguments.
!
    integer, intent(in) :: length, time
    real(kind=realtype), intent(out) :: mult, trans
! determine the situation we are having here.
! first the length.
    select case  (length) 
    case (meter) 
      mult = one
      trans = zero
    case (centimeter) 
      mult = 0.01_realtype
      trans = zero
    case (millimeter) 
      mult = 0.001_realtype
      trans = zero
    case (foot) 
      mult = 0.3048_realtype
      trans = zero
    case (inch) 
      mult = 0.0254_realtype
      trans = zero
    case default
      call terminate('sivelocity', &
&              'no idea how to convert this length to si units')
    end select
! and the time.
    select case  (time) 
    case (second) 
      mult = mult
    case default
      call terminate('sivelocity', &
&              'no idea how to convert this time to si units')
    end select
  end subroutine sivelocity
! ----------------------------------------------------------------------
!                                                                      |
!                    no tapenade routine below this line               |
!                                                                      |
! ----------------------------------------------------------------------
 
end module utils_d