suydam.f90

Go to the documentation of this file.

subroutine suydam(ipsi, t0, tbb, tbf, fact, bal)
!-----------------------------------------------------------------------
!  subroutine to verify the sign of the ballooning energy        
!  functional according to the suydam finite difference scheme.  
!-----------------------------------------------------------------------

  use itm_types
  use mod_parameter
  use mod_suydam

  implicit none

  integer(itm_i4), intent(in) :: ipsi
  real(r8), intent(in) :: t0, tbb, tbf, fact
  character (len=25), intent(out) :: bal

  integer(itm_i4) :: i, n
  real(r8) :: dbt, alp, a01, a11
  real(r8) :: t1, t2, ts
  real(r8) :: kp1, kp2, gp1, gp2
  real(r8) :: gp2temp

  bal = ' stable'
  dbt = 4._r8 * asin(1._r8) / dble(ncpq - 1)

  n  = int((tbf - tbb) / dbt)

  alp = 1._r8

  t2 = tbb - dbt / 2._r8

  kp2 = 0._r8
  gp2 = 0._r8

  do i = 1, n
    t1  = t2
    kp1 = kp2
    gp1 = gp2
    t2  = t1 + dbt
    call suydam_kgs(t2, t0, kp2, gp2temp, ipsi)
    gp2 = gp2temp * fact
    a11 = (kp2 + kp1) / dbt + 0.25_r8 * (gp1 + gp2) * dbt       
    a01 = - kp1 / dbt + 0.25_r8 * gp1 * dbt
    if (i == 1) a01 = 0._r8
    alp = a11 - a01**2 / alp
    if(alp <= 0._r8) then
      ts = (t1 + t2) / 2._r8
      write(bal, 11) ts
      exit     
    end if
  end do

 11 format(' unstable at t = ',f8.3)

  return

end subroutine suydam