analytical_solution Subroutine

    subroutine analytical_solution(nr)
      integer, intent(in) :: nr
      integer(i4k) i, n
      real(r8k) t, r_0, R, T0, T_inf
      real(r8k), parameter :: pi = 4.0_r8k*atan(1.0_r8k)
      real(r8k), parameter, dimension(*) :: mu = [(( 2*i-1)*pi/2.0, i = 1, 20 )]

      T0=1073.15
      T_inf=293.15
      t=1000.0
      r_0=0.2

      i=1
      R=this%v(i,1)/r_0
      this%T_analytical = &
        sum( [( 2.0*(sin(mu(n))-mu(n)*cos(mu(n)))/(mu(n)-sin(mu(n))*cos(mu(n)))*exp(-1.0e-5*(mu(n)/r_0)**2*t), n = 1, size(mu))] )
      this%T_analytical(i) = (T0-T_inf)*this%T_analytical(i) + T_inf

      do concurrent(i=2:nr)
        R=this%v(i,1)/r_0
        this%T_analytical(i) = &
          sum( [( 2.0*(sin(mu(n))-mu(n)*cos(mu(n)))/(mu(n)-sin(mu(n))*cos(mu(n)))* &
          sin(mu(n)*R)/(mu(n)*R)*exp(-1.0e-5*(mu(n)/r_0)**2*t), n = 1, 6 )] )
        this%T_analytical(i) = (T0-T_inf)*this%T_analytical(i) + T_inf
      end do
    end subroutine analytical_solution