soloviev.f90

Go to the documentation of this file.

subroutine soloviev(x, y, solpsi, psix, psiy, psixy)
!-----------------------------------------------------------------------
! soloviev equilibrium (not called by Helena)
!-----------------------------------------------------------------------
  
  use itm_types
  use mod_dat
  
  implicit none
  
  real(r8), intent(in) :: x
  real(r8), intent(inout) :: y
  real(r8), intent(out) :: solpsi, psix, psiy, psixy
  
  y = -y
  par4 = quad
  par3 = tria
  par2 = eps
  par1 = ellip
  
  solpsi = (1._r8 - 0.25_r8 * par2**2) * (1._r8 + par2 * x)**2 &
   * (y / par1)**2 + (x - 0.5_r8 * par2 * (1._r8 - x * x))**2 &
   + (4._r8 - par2**2) * (x - 0.5_r8 * par2 * (1._r8 - x**2)) &
   * par4 * y / (2._r8 * par2)    
  
  psix = 2._r8* (x - 0.5_r8 * par2 * (1._r8 - x * x)) * (1._r8 &
   + par2 * x) + (2._r8 * par2 * (1._r8 - 0.25_r8 * par2**2) &
   * (1._r8 + par2 * x)) * (y / par1)**2 + (4._r8 - par2**2) &
   * (1._r8 + 0.5_r8 * par2 * (2._r8 * x)) * par4 &
   * y / (2._r8 * par2)
  
  psiy = 2._r8 * y / par1**2 * (1._r8 - 0.25_r8 * par2**2) &
   * (1._r8 +  par2 * x)**2 + (4._r8 - par2**2) * (x - 0.5_r8 &
   * par2 * (1._r8 - x**2)) * par4 / (2._r8 * par2)      
  
  psiy = -psiy
  
  psixy = 2._r8 * y / par1**2 * (1._r8 - 0.25_r8 * par2**2) &
   * (1._r8 + par2 * x) * 2._r8 * par2 + (4._r8 - par2**2) &
   * (1._r8 + 0.5_r8 * par2 * (2._r8 * x)) * par4 &
   / (2._r8 * par2)
  
  psixy = -psixy
  
  return
end subroutine soloviev