ets_math.f90

Go to the documentation of this file.

! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
!-------------------------------------------------------!
!                                                       !
!____________ ETS  MATHEMATICAL SUBROUTINES: ___________!
!                                                       !
!-------------------------------------------------------!
! These subroutines are aimed to provide the transport  !
! solver with the simple mathematics                    !
!-------------------------------------------------------!
MODULE ets_math
CONTAINS




!-------------------------------------------------------!
!-------------------------------------------------------!
SUBROUTINE deriv_fun(N, X, Y, DY)
!-------------------------------------------------------!
! These subroutine calculate the derivative, DY         !
! of the function Y with respect to argument X          !
!-------------------------------------------------------!

  USE itm_types

  IMPLICIT NONE

  INTEGER        :: n               ! number of radial points   (input)
  INTEGER        :: i               ! radial index

  REAL (R8)      :: x(n)            ! argument array            (input)
  REAL (R8)      :: y(n)            ! function array            (input)
  REAL (R8)      :: dy(n)           ! function derivative array (output)
  REAL (R8)      :: h(n)


  DO i=1,n-1
     h(i)=x(i+1)-x(i)
  END DO

  DO i=2,n-1
     dy(i) = ((y(i+1)-y(i))*h(i-1)/h(i) + (y(i)-y(i-1))*h(i)/h(i-1)) / (h(i)+h(i-1))
  END DO

  dy(1)    =  dy(2)  - 2.e0_r8*((y(1)-y(2))/h(1)+(y(3)-y(2))/h(2))/(h(2)+h(1)) * h(1)

  dy(n)    =  dy(n-1)+ 2.e0_r8*((y(n-2)-y(n-1))/h(n-2)+(y(n)-y(n-1))/h(n-1))/(h(n-1)+h(n-2)) * h(n-1)
! DY(N)    =  (Y(N)-Y(N-1))/H(N-1)

  RETURN

END SUBROUTINE deriv_fun
!-------------------------------------------------------!
!-------------------------------------------------------!





!-------------------------------------------------------!
!-------------------------------------------------------!
SUBROUTINE integr_fun(N,X,Y,INTY) 
!-------------------------------------------------------!
!  This subroutine calculates integral of function      !
!  Y(X) from X=0 until X=X(N)                           !
!-------------------------------------------------------!

  USE itm_types

  IMPLICIT NONE
  INTEGER        :: n               ! number of radial points (input)
  INTEGER        :: i

  REAL (R8)      :: x(n)            ! argument array (input)
  REAL (R8)      :: y(n)            ! function array (input)
  REAL (R8)      :: inty(n)         ! function integral array (output)


  inty(1)    = y(1)*x(1)/2.e0_r8
  DO i=2,n
     inty(i) = inty(i-1)+(y(i-1)+y(i))*(x(i)-x(i-1))/2.e0_r8
  END DO

  RETURN

END SUBROUTINE integr_fun
!-------------------------------------------------------!
!-------------------------------------------------------!





!-------------------------------------------------------!
!-------------------------------------------------------!
SUBROUTINE integr_xfun(N,X,Y,INTXY)
!-------------------------------------------------------!
!  This subroutine calculates integral of function      !
!  Y(X)*X from X=0 until X=X(N)                         !
!-------------------------------------------------------!

  USE itm_types

  IMPLICIT NONE
  INTEGER        :: n               ! number of radial points (input)
  INTEGER        :: i

  REAL (R8)      :: x(n)            ! argument array (input)
  REAL (R8)      :: y(n)            ! function array (input)
  REAL (R8)      :: intxy(n)         ! function integral array (output)

  intxy(1)    = y(1)*x(1)**2/2.e0_r8
  DO i=2,n
     intxy(i) = intxy(i-1)+(y(i-1)*x(i-1)+y(i)*x(i))*(x(i)-x(i-1))/2.e0_r8
  END DO

  RETURN
END SUBROUTINE integr_xfun
!-------------------------------------------------------!
!-------------------------------------------------------!





!-------------------------------------------------------!
!-------------------------------------------------------!
SUBROUTINE fix_axis_value(GEOMETRY,PROFILE,N)
!-------------------------------------------------------!
!                                                       !
! This subroutine shall fix central points of PROFILE   !
! from (1) to (N) using parabolic approximation         !
!                                                       !
!-------------------------------------------------------!
  USE itm_types
  USE itm_constants
  USE ets_plasma

  IMPLICIT NONE


! +++ Plasma parameters:
  TYPE (magnetic_geometry)      :: geometry               
  REAL (R8)                     :: profile(geometry%nrho)    
  REAL (R8)                     :: a, b, c 
  REAL (R8)                     :: x(n+2)
  INTEGER                       :: i,n

  DO i = 1,n+2
     x(i) = geometry%RHO(i)
  END DO

  a = (profile(n)-profile(n+2))/(x(n)-x(n+2))/x(n+1)/2.0_r8
  c = profile(n+1) - a*x(n+1)**2.0_r8

  DO i = 1,n
     profile(i)  = a*x(i)**2.0_r8 + c
  END DO

  RETURN
END SUBROUTINE fix_axis_value
!-------------------------------------------------------!
!-------------------------------------------------------!



END MODULE ets_math