transition Subroutine

  subroutine transition ( rho, R, S, fT, fTp, fTpp )
  !! Routine to calculate the smooth transition function, \(fT\), and it's first
  !! \(fT'\) and second \(fT''\) derivative with respect to the computational
  !! coordinate, \(\rho\).

    implicit none

    real(wp), intent(in) :: rho
    !! The computational coordinate, \(\rho\), where the transition function
    !! should be calculated. Should be between \(R\) and \(S\).
    real(wp), intent(in) :: R
    !! The value of \(\rho\) where \(fT\) should be \(0\).
    real(wp), intent(in) :: S
    !! The value of \(\rho\) where \(fT\) should be \(1\).
    real(wp), intent(out) :: fT
    !! On return contains \(fT\).
    real(wp), intent(out) :: fTp
    !! On return contains \(fT'\).
    real(wp), intent(out) :: fTpp
    !! On return contains \(fT''\).
    real(wp), parameter :: es = 1.5_wp
    real(wp), parameter :: q2 = 1.3_wp
    real(wp), parameter :: pi = acos(-1.0_wp)
    real(wp), parameter :: half = 0.5_wp
    real(wp) :: x, tanx, cotx, fac, tanhfac, cscx2, secx2, sechfac2

    x = max ( half * pi * ( rho - R ) / ( S - R ), 0.0_wp )
    if ( x > 0.0_wp ) then
!      print*,'x = ', x
      tanx = tan ( x )
!      print*,'tanx = ', tanx
      cotx = 1.0_wp / tanx
!      print*,'cotx = ', cotx
      fac = es/pi * ( tanx - q2 * cotx )
!      print*,'fac = ', fac
      tanhfac = tanh ( fac)
!      print*,'tanhfac = ', tanhfac
      fT = half + half * tanhfac
!      print*,'fT = ', fT
      cscx2 = 1.0_wp / sin ( x )**2
      secx2 = 1.0_wp / cos ( x )**2
      sechfac2 = 1.0_wp / cosh ( fac )**2
      fTp = 0.25_wp * es * ( q2 * cscx2 + secx2 ) * sechfac2 / ( S - R )
      fTpp = -0.25_wp * es * sechfac2 * ( &
               pi * ( q2 * cotx * cscx2 - secx2 * tanx ) + &
               es * ( q2 * cscx2 + secx2 )**2 * tanhfac ) / ( S - R )**2
    else
      fT = 0.0_wp
      fTp = 0.0_wp
      fTpp = 0.0_wp
    end if
  end subroutine transition