radial_mesh.f90

Go to the documentation of this file.

 subroutine radial_mesh(rpsi, zpsi, dzpsi, ddzpsi)
!-------------------------------------------------------------------
! function to determine a non-equidistant grid in sqrt(psi).
! grid will be equidistant in rpsi.
! rpsi must be a monotonic function of psi between 0. and 1.
!-------------------------------------------------------------------

  use itm_types
  use mod_meshacc
  
  implicit none
  
  real(r8), intent(in) :: rpsi
  
  real(r8), intent(out) :: zpsi, dzpsi, ddzpsi
  
  real(r8) :: a1, a2
  
  zpsi = rpsi**2
  dzpsi = 2._r8 * rpsi
  ddzpsi = 2._r8
!TODO: remove this redundant re-accumulation of the grid
  if (imesh == 1) then
    a1 = amesh
    a2 = bmesh
  
    zpsi = rpsi**2 + a1 * rpsi**3 + a2 * rpsi**4  - (a1 + a2) &
     * rpsi**5
    dzpsi = 2._r8 * rpsi + 3._r8 * a1 * rpsi**2 + 4._r8 * a2 &
     * rpsi**3 - 5._r8 * (a1 + a2) * rpsi**4
    ddzpsi = 2._r8 + 6._r8 * a1 * rpsi + 12._r8 * a2 * rpsi**2 &
     - 20._r8 * (a1 + a2) * rpsi**3
  end if
  
  return
end subroutine radial_mesh