time_advancing Subroutine

    subroutine time_advancing(nr, dt)
      integer, intent(in) :: nr
      real(r8k), intent(in)  :: dt

      real(r8k), dimension(:), allocatable :: a, b, c, d
      real(r8k) t, dr_m,  rf

      real(r8k), parameter :: h=230      !! heat transfer coefficient
      real(r8k), parameter :: Tb=293.15  !! ambient temperature

      t=0.0

      allocate(a(nr),b(nr),c(nr),d(nr))

      do while(t<1000)

        call this%set_rho()
        call this%set_cp()

        i=1
        associate( &
          e    => 1.0/(this%rho(i)*this%cp(i)), &
          dr_f => this%v(i+1,1) - this%v(i,1), &
          rf   => 0.5*(this%v(i+1,1) + this%v(i,1)) &
        )
          b(i) =  6.0*e*kc(rf)/(dr_f**2) + 1.0/dt
          c(i) = -6.0*e*kc(rf)/(dr_f**2)
          d(i) = this%v(i,2)/dt
        end associate

        do concurrent( i=2:nr-1 )
          associate( &
            e    => 1.0/(this%rho(i)*this%cp(i)), &
            dr_f => this%v(i+1,1) - this%v(i  ,1), &
            dr_b => this%v(i  ,1) - this%v(i-1,1), &
            dr_m => 0.5*(this%v(i+1,1) - this%v(i-1,1)), &
            rf   => 0.5*(this%v(i+1,1) + this%v(i,1)), &
            rb   => 0.5*(this%v(i,1)   + this%v(i-1,1)) &
          )
            a(i) = -e*kc(rb)*rb**2/(dr_b*dr_m*this%v(i,1)**2)
            b(i) =  e*kc(rf)*rf**2/(dr_f*dr_m*this%v(i,1)**2) + &
                    e*kc(rb)*rb**2/(dr_b*dr_m*this%v(i,1)**2) + 1.0/dt
            c(i) = -e*kc(rf)*rf**2/(dr_f*dr_m*this%v(i,1)**2)
            d(i) = this%v(i,2)/dt
          end associate
        end do

        i=nr
        associate( e => 1.0/(this%rho(i)*this%cp(i)), dr_b => this%v(i,1) - this%v(i-1,1), rb => 0.5*(this%v(i,1) + this%v(i-1,1)))
          dr_m = dr_b                  !! ifort 18 bug precludes associate( dr_m => dr_b)
          rf = this%v(i,1) + 0.5*dr_b  !! ifort 18 bug precludes associate( rf => this%v(i,1) + 0.5*dr_b)
            a(i)= -1.0*e*kc(rb)*rb**2/(dr_b*dr_m*this%v(i,1)**2)
            b(i)= e*h*rf**2/(dr_m*this%v(i,1)**2) + e*kc(rb)*rb**2/(dr_b*dr_m*this%v(i,1)**2) + 1.0/dt
            d(i)= this%v(i,2)/dt+e*h*Tb*rf**2/(dr_m*this%v(i,1)**2)
        end associate

        this%v(:,2) = tridiagonal_matrix_algorithm(a,b,c,d)
        t=t+dt
      end do
    end subroutine time_advancing