suydam_pq.f90

Go to the documentation of this file.

subroutine suydam_pq(ias)
!-----------------------------------------------------------------------
! subroutine to calculate the p and q coefficents in the ballooning
! equation in all grid points.
! the terms depending on the extended ballooning angle are seperated.
! the ballooning equation (poguste and yurchenko) reads:
!       dt((p0 + t * p1 + t^2 p2) dt(f)) - (q0 + t * q1) f = 0 
!
! note : the final arrays cp and cq have a different index from the
! gij_hel arrays. the j index runs from 0 to 2pi (including 2pi), the
! number of points in poloidal direction is nchi+1 for ias=1 and 
! 2*nchi-1 for ias=0.
!-----------------------------------------------------------------------      

  use itm_types
  use mod_map
  use mod_spline
  use mod_suydam

  implicit none

  integer(itm_i4), intent(in) :: ias

  integer(itm_i4) :: i, j, ij, ij1, ij2
  real(r8) :: sps2, zq, zt, zdp, zdq
  real(r8) :: beta

  if (ias == 1) then
    ncpq = nchi + 1
  else
    ncpq = 2 * nchi - 1
  end if

  allocate(cp0(npsi * ncpq))
  allocate(cp1(npsi * ncpq))
  allocate(cp2(npsi * ncpq))
  allocate(cq0(npsi * ncpq))
  allocate(cq1(npsi * ncpq))

  do i = 2, npsi
    sps2 = 2._r8 * cpsurf * cs(i)
    zq = qs(i)
    zt = rbphi(i)
    zdp = p_spline%sp2(i)
    zdq = dqs(i)
    do j = 1, nchi
      ij = (i - 1) * nchi + j
      ij1 = (i - 1) * ncpq + j
!----------------------------------- lower index geometric coeff. -------     
!      g33 = 1._r8 / g33_hel(i, j)
!      g11 = sps2**2 * (1._r8 + zq**2 / zt**2 / g33_hel(i, j) &
!       * g12_hel(i, j)**2 ) / g11_hel(i, j)
!      g12 = - sps2 * zq**2 / zt**2 * g12_hel(i, j) / g33_hel(i, j)
!      g22 = zq**2 / zt**2 * g11_hel(i, j) / g33_hel(i, j)
!      zj  = sps2 * zq / g33_hel(i, j) / zt
!----------------------------------------------------------------------- 
      beta =  -g12_hel(i, j) / g11_hel(i, j) 
      cp0(ij1) = 1._r8 * g33_hel(i, j) / g11_hel(i, j) + zq**2 & 
       * g11_hel(i, j) * beta**2 * g33_hel(i, j) / b02(ij)
      cp1(ij1) = 2._r8 * beta * g11_hel(i, j) * (zdq * zq) / (sps2 &
       / g33_hel(i, j) * b02(ij))
      cp2(ij1) = g11_hel(i, j) * zdq**2 / (sps2**2 / g33_hel(i, j) &
       * b02(ij))
      cq0(ij1) = -zdp * zq**2 / g33_hel(i, j) / (sps2 * zt * b02(ij))**2 &
       * ((2._r8 * zdp + ds_b02(ij)) * b02(ij) + sps2 * beta &
       * g11_hel(i, j) * g33_hel(i, j) * dt_b02(ij))     
      cq1(ij1) = zdp * zdq * zq / (sps2 * b02(ij))**2 * dt_b02(ij)   
!----------------------------------------------------------------------- 
      if (ias == 0) then
        ij2 = i * (2 * nchi - 1) - j + 1 
        cp0(ij2) = cp0(ij1)
        cp1(ij2) = -cp1(ij1)
        cp2(ij2) = cp2(ij1)
        cq0(ij2) = cq0(ij1)
        cq1(ij2) = -cq1(ij1)
      else if (ias == 1 .and. j == 1) then
        ij2 = i * (nchi + 1)
        cp0(ij2) = cp0(ij1)
        cp1(ij2) = cp1(ij1)
        cp2(ij2) = cp2(ij1)
        cq0(ij2) = cq0(ij1)
        cq1(ij2) = cq1(ij1)
      end if    
    end do         
  end do

  return

end subroutine suydam_pq