bivariate.f

Go to the documentation of this file.

!*******************************************************************************
!
!  Note separating the Doxygen comment block here so detailed decription is
!  found in the Module not the file.
!
!*******************************************************************************
      MODULE bivariate
 
      USE stel_kinds
      USE profiler
 
      IMPLICIT NONE
 
!*******************************************************************************
!  DERIVED-TYPE DECLARATIONS
!  1) bivariate_type
!
!*******************************************************************************
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      TYPE bivariate_type
         INTEGER                            :: nsu
         INTEGER                            :: nrz
         INTEGER, DIMENSION(:), POINTER     :: irz11_bi => null()
         INTEGER, DIMENSION(:), POINTER     :: irz12_bi => null()
         INTEGER, DIMENSION(:), POINTER     :: irz21_bi => null()
         INTEGER, DIMENSION(:), POINTER     :: irz22_bi => null()
         REAL(rprec), DIMENSION(:), POINTER :: w11_bi => null()
         REAL(rprec), DIMENSION(:), POINTER :: w12_bi => null()
         REAL(rprec), DIMENSION(:), POINTER :: w21_bi => null()
         REAL(rprec), DIMENSION(:), POINTER :: w22_bi => null()
      END TYPE
 
*******************************************************************************
!  INTERFACE BLOCKS
!*******************************************************************************
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      INTERFACE bivariate_set_grids
         MODULE PROCEDURE bivariate_set_grids_1d,                              &
     &                    bivariate_set_grids_2d
      END INTERFACE
 
      CONTAINS
!*******************************************************************************
!  CONSTRUCTION SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      FUNCTION bivariate_construct(ns, nu)
 
      IMPLICIT NONE
 
!  Declare Arguments
      TYPE (bivariate_type), POINTER :: bivariate_construct
      INTEGER, INTENT(in)            :: ns
      INTEGER, INTENT(in)            :: nu
 
!  local variables
      REAL (rprec)                   :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      ALLOCATE(bivariate_construct)
 
      bivariate_construct%nsu = ns*nu
 
      ALLOCATE (bivariate_construct%irz11_bi(bivariate_construct%nsu))
      ALLOCATE (bivariate_construct%irz12_bi(bivariate_construct%nsu))
      ALLOCATE (bivariate_construct%irz21_bi(bivariate_construct%nsu))
      ALLOCATE (bivariate_construct%irz22_bi(bivariate_construct%nsu))
      ALLOCATE (bivariate_construct%w11_bi(bivariate_construct%nsu))
      ALLOCATE (bivariate_construct%w12_bi(bivariate_construct%nsu))
      ALLOCATE (bivariate_construct%w21_bi(bivariate_construct%nsu))
      ALLOCATE (bivariate_construct%w22_bi(bivariate_construct%nsu))
 
      CALL profiler_set_stop_time('bivariate_construct', start_time)
 
      END FUNCTION
 
!*******************************************************************************
!  DESTRUCTION SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE bivariate_destruct(this)
 
      IMPLICIT NONE
 
!  Declare Arguments
      TYPE (bivariate_type), POINTER :: this
 
!  Start of executable code
      IF (ASSOCIATED(this%irz11_bi)) THEN
         DEALLOCATE (this%irz11_bi)
         this%irz11_bi => null()
      END IF
 
      IF (ASSOCIATED(this%irz12_bi)) THEN
         DEALLOCATE (this%irz12_bi)
         this%irz12_bi => null()
      END IF
 
      IF (ASSOCIATED(this%irz21_bi)) THEN
         DEALLOCATE (this%irz21_bi)
         this%irz21_bi => null()
      END IF
 
      IF (ASSOCIATED(this%irz22_bi)) THEN
         DEALLOCATE (this%irz22_bi)
         this%irz22_bi => null()
      END IF
 
      IF (ASSOCIATED(this%w11_bi)) THEN
         DEALLOCATE (this%w11_bi)
         this%w11_bi => null()
      END IF
 
      IF (ASSOCIATED(this%w12_bi)) THEN
         DEALLOCATE (this%w12_bi)
         this%w12_bi => null()
      END IF
 
      IF (ASSOCIATED(this%w21_bi)) THEN
         DEALLOCATE (this%w21_bi)
         this%w21_bi => null()
      END IF
 
      IF (ASSOCIATED(this%w22_bi)) THEN
         DEALLOCATE (this%w22_bi)
         this%w22_bi => null()
      END IF
 
      DEALLOCATE(this)
 
      END SUBROUTINE
 
!*******************************************************************************
!  SETTER SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE bivariate_set_grids_2d(this, rsu, zsu, rgrid, zgrid)
 
      IMPLICIT NONE
 
!  Declare Arguments
      TYPE (bivariate_type), INTENT(inout)    :: this
      REAL(rprec), DIMENSION(:,:), INTENT(in) :: rsu
      REAL(rprec), DIMENSION(:,:), INTENT(in) :: zsu
      REAL(rprec), DIMENSION(:), INTENT(in)   :: rgrid
      REAL(rprec), DIMENSION(:), INTENT(in)   :: zgrid
 
!  local variables
      INTEGER                                 :: ku, js
      INTEGER                                 :: ns, nu
      INTEGER                                 :: nr, nz
      INTEGER                                 :: index1d
      INTEGER                                 :: ir, jz
      INTEGER                                 :: ir1, jz1
      REAL(rprec)                             :: rad0, zee0
      REAL(rprec)                             :: ri, zj
      REAL(rprec)                             :: pr, qz
      REAL(rprec)                             :: temp
      REAL(rprec)                             :: delr, delz
      REAL (rprec)                            :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      nr = SIZE(rgrid, 1)
      nz = SIZE(zgrid, 1)
      this%nrz = nr*nz
      ns = SIZE(rsu, 1)
      nu = SIZE(rsu, 2)
 
      delr = rgrid(2) - rgrid(1)
      delz = zgrid(2) - zgrid(1)
 
      index1d = 0
 
      DO ku = 1, nu
         DO js = 1, ns
 
            index1d = index1d + 1
!
!       CHECK THAT BOUNDARY POINTS ARE INSIDE GRID.  IF NOT, STOP!
!
            rad0 = rsu(js, ku)
            zee0 = zsu(js, ku)
 
            IF (rad0.lt.rgrid(1) .or. rad0.gt.rgrid(nr) .or.                   &
     &          zee0.lt.zgrid(1) .or. zee0.gt.zgrid(nz)) THEN
               stop 'Plasma point outside response function grid!'
            END IF
!
!       DETERMINE INTEGER INDICES (IR,JZ) FOR LOWER LEFT R, Z CORNER GRID POINT
!
            ir = int((rad0 - rgrid(1))/delr) + 1
            jz = int((zee0 - zgrid(1))/delz) + 1
            ir1 = min(nr, ir + 1)
            jz1 = min(nz, jz + 1)
 
!
!        STORE INDICES IN 1D ARRAYS
!
            this%irz11_bi(index1d) = ir  + nr*(jz - 1)
            this%irz22_bi(index1d) = ir1 + nr*(jz1 - 1)
            this%irz12_bi(index1d) = ir  + nr*(jz1 - 1)
            this%irz21_bi(index1d) = ir1 + nr*(jz - 1)
!
!       COMPUTE RI, ZJ AND PR , QZ AT GRID POINT (IR , JZ)
!       ALSO, COMPUTE WEIGHTS WIJ FOR 4 CORNER GRID POINTS
!
            ri = rgrid(ir)
            zj = zgrid(jz)
            pr = (rad0 - ri)/delr
            qz = (zee0 - zj)/delz
            temp = pr*qz
            this%w22_bi(index1d) = temp                   !p*q
            this%w21_bi(index1d) = pr - temp              !p*(1-q)
            this%w12_bi(index1d) = qz - temp              !q*(1-p)
            this%w11_bi(index1d) = 1 + temp - (pr + qz)   !(1-p)*(1-q)
 
         END DO
      END DO
 
      CALL profiler_set_stop_time('bivariate_set_grids_2d', start_time)
 
      END SUBROUTINE
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE bivariate_set_grids_1d(this, rsu, zsu, rgrid, zgrid)
 
      IMPLICIT NONE
 
!  Declare Arguments
      TYPE (bivariate_type), INTENT(inout)  :: this
      REAL(rprec), DIMENSION(:), INTENT(in) :: rsu
      REAL(rprec), DIMENSION(:), INTENT(in) :: zsu
      REAL(rprec), DIMENSION(:), INTENT(in) :: rgrid
      REAL(rprec), DIMENSION(:), INTENT(in) :: zgrid
 
!  local variables
      INTEGER                               :: ku
      INTEGER                               :: ns, nu
      INTEGER                               :: nr, nz
      INTEGER                               :: index1d
      INTEGER                               :: ir, jz
      INTEGER                               :: ir1, jz1
      REAL(rprec)                           :: rad0, zee0
      REAL(rprec)                           :: ri, zj
      REAL(rprec)                           :: pr, qz
      REAL(rprec)                           :: temp
      REAL(rprec)                           :: delr, delz
      REAL (rprec)                          :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      nr = SIZE(rgrid, 1)
      nz = SIZE(zgrid, 1)
      this%nrz = nr*nz
      ns = 1
      nu = SIZE(rsu)
 
      delr = rgrid(2) - rgrid(1)
      delz = zgrid(2) - zgrid(1)
 
      index1d = 0
 
      DO ku = 1, nu
         index1d = index1d + 1
!
!       CHECK THAT BOUNDARY POINTS ARE INSIDE GRID.  IF NOT, STOP!
!
         rad0 = rsu(ku)
         zee0 = zsu(ku)
 
         IF (rad0.lt.rgrid(1) .or. rad0.gt.rgrid(nr) .or.                      &
     &       zee0.lt.zgrid(1) .or. zee0.gt.zgrid(nz)) THEN
            stop 'Plasma point outside response function grid!'
         END IF
!
!       DETERMINE INTEGER INDICES (IR,JZ) FOR LOWER LEFT R, Z CORNER GRID POINT
!
         ir = int((rad0 - rgrid(1))/delr) + 1
         jz = int((zee0 - zgrid(1))/delz) + 1
         ir1 = min(nr, ir + 1)
         jz1 = min(nz, jz + 1)
 
!
!        STORE INDICES IN 1D ARRAYS
!
         this%irz11_bi(index1d) = ir  + nr*(jz - 1)
         this%irz22_bi(index1d) = ir1 + nr*(jz1 - 1)
         this%irz12_bi(index1d) = ir  + nr*(jz1 - 1)
         this%irz21_bi(index1d) = ir1 + nr*(jz - 1)
!
!       COMPUTE RI, ZJ AND PR , QZ AT GRID POINT (IR , JZ)
!       ALSO, COMPUTE WEIGHTS WIJ FOR 4 CORNER GRID POINTS
!
         ri = rgrid(ir)
         zj = zgrid(jz)
         pr = (rad0 - ri)/delr
         qz = (zee0 - zj)/delz
         temp = pr*qz
         this%w22_bi(index1d) = temp                   !p*q
         this%w21_bi(index1d) = pr - temp              !p*(1-q)
         this%w12_bi(index1d) = qz - temp              !q*(1-p)
         this%w11_bi(index1d) = 1 + temp - (pr + qz)   !(1-p)*(1-q)
      END DO
 
      CALL profiler_set_stop_time('bivariate_set_grids_2d', start_time)
 
      END SUBROUTINE
 
!*******************************************************************************
!  GETTER SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE bivariate_get_4pt(this, resp_rz, resp_su)
      IMPLICIT NONE
      TYPE (bivariate_type), INTENT(inout)          :: this
      REAL(rprec), DIMENSION(this%nrz), INTENT(in)  :: resp_rz
      REAL(rprec), DIMENSION(this%nsu), INTENT(out) :: resp_su
 
!  local variables
      INTEGER                                       :: jsu
      REAL (rprec)                                  :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      DO jsu = 1, this%nsu
!
!       COMPUTE RESPONSE AT R, PHI (fixed kin), Z BY 4-PT INTERPOLATION
!
         resp_su(jsu) = this%w11_bi(jsu)*resp_rz(this%irz11_bi(jsu))           &
     &                + this%w22_bi(jsu)*resp_rz(this%irz22_bi(jsu))           &
     &                + this%w21_bi(jsu)*resp_rz(this%irz21_bi(jsu))           &
     &                + this%w12_bi(jsu)*resp_rz(this%irz12_bi(jsu))
 
      END DO
 
      CALL profiler_set_stop_time('bivariate_get_4pt', start_time)
 
      END SUBROUTINE
 
      END MODULE