text:
!
! CRTM_Interpolation
!
! Module containing routines for data interpolation in the CRTM
!
!
! CREATION HISTORY:
!       Written by:     Yong Han, 27-July-2005
!                       Yong.Han@noaa.gov
!

MODULE CRTM_Profile_Interpolation

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module use
  USE Type_Kinds,      ONLY: fp=>fp_kind
  USE CRTM_Parameters, ONLY: ZERO, ONE, TOLERANCE
  ! Disable all implicit typing
  IMPLICIT NONE

  ! ------------
  ! Visibilities
  ! ------------
  ! Everything private by default
  PRIVATE
  ! Routines in this module
  PUBLIC  :: Interpolate_Profile
  PUBLIC  :: Interpolate_Profile_TL
  PUBLIC  :: Interpolate_Profile_AD
  PUBLIC  :: Compute_Interp_Index

  ! ---------------------
  ! Procedure overloading
  ! ---------------------
  INTERFACE Interpolate_Profile
    MODULE PROCEDURE Interpolate_Profile_F1
    MODULE PROCEDURE Interpolate_Profile_F2
  END INTERFACE Interpolate_Profile

  INTERFACE Interpolate_Profile_TL
    MODULE PROCEDURE Interpolate_Profile_F1_TL
    MODULE PROCEDURE Interpolate_Profile_F2_TL
  END INTERFACE Interpolate_Profile_TL

  INTERFACE Interpolate_Profile_AD
    MODULE PROCEDURE Interpolate_Profile_F1_AD
    MODULE PROCEDURE Interpolate_Profile_F2_AD
  END INTERFACE Interpolate_Profile_AD

CONTAINS

!---------------------------------------------------------------------------------------------
! NAME: Interpolate_Profile
!
! PURPOSE:
!    Given x and u that are ascending arrays, it interpolates y with the abscissa x
!    on the abscissa u using the following algorithm:
!       y_int(i) = y(1)  if u(i) < x(1)
!       y_int(i) = y(nx) if u(i) > x(nx)
!       y_int(i) = y(ix1) + (y(ix2)-y(ix1))*(u(i) - x(ix1))/(x(ix2)-x(ix1))
!                        if x(ix1) <= u(i) <= x(ix2)
!
!    IThe index array interp_index contains the following content
!
!      interp_index(i, 1) = 1 and interp_index(i, 2) = 1, if u(i) < x(1)
!      interp_index(i, 1) = nx and interp_index(i, 2) = nx, if u(i) > x(nx),
!                                                          where nx = SIZE(x)
!      x(interp_index(i, 1)) <= u(i) <= x(interp_index(i, 2)) if x(1) <= u(i) <= x(nx)
!
! CALLING SEQUENCE:
!            CALL Interpolate_Profile(y, x, u, y_int)
!      or    CALL Interpolate_Profile(interp_index, y, x, u, y_int)
!
! INPUT ARGUMENTS:
!       y:            The data array to be interpolated.
!                     UNITS:      N/A
!                     TYPE:       fp
!                     DIMENSION:  rank-1
!                     ATTRIBUTES: INTENT(IN)
!
!       x:            The abscissa values for y and they must be monotonically
!                     ascending.
!                     UNITS:      N/A
!                     TYPE:       fp
!                     DIMENSION:  rank-1
!                     ATTRIBUTES: INTENT(IN)
!
!       u:            The abscissa values for the results
!                     and they must be monotonically ascending
!                     UNITS:      N/A
!                     TYPE:       fp
!                     DIMENSION:  rank-1
!                     ATTRIBUTES: INTENT(IN)
!
!    interp_index:    The index array of dimension (nu x 2), where nu = SIZE(u)
!                     UNITS:      N/A
!                     TYPE:       Integer
!                     DIMENSION:  rank-2
!                     ATTRIBUTES: INTENT(IN)
!
! OUTPUT ARGUMENTS:
!       y_int:        The array contains the results
!                     UNITS:      N/A
!                     TYPE:       Integer
!                     DIMENSION:  rank-1
!                     ATTRIBUTES: INTENT(OUT)
!
! RESTRICTIONS:
!     To be efficient, this routine does not check that x and u are both
!     monotonically ascending and the index bounds.
!
! CREATION HISTORY:
!       Written by:     Yong Han, 07-May-2004
!-----------------------------------------------------------------------------------

  !----------------------------------------------------------------------------
  !  Interpolation routine with interperlation index array already calculated.
  !----------------------------------------------------------------------------
  SUBROUTINE Interpolate_Profile_F1(interp_index, y, x, u, y_int)
    ! Arguments
    INTEGER,  DIMENSION(:,:), INTENT(IN)  :: interp_index
    REAL(fp), DIMENSION(:),   INTENT(IN)  :: y
    REAL(fp), DIMENSION(:),   INTENT(IN)  :: x
    REAL(fp), DIMENSION(:),   INTENT(IN)  :: u
    REAL(fp), DIMENSION(:),   INTENT(OUT) :: y_int
    ! Local variables
    INTEGER :: i, k1, k2

    DO i = 1, SIZE(u)
      k1 = interp_index(i, 1)
      k2 = interp_index(i, 2)
      IF( k1 == k2)THEN
        y_int(i) = y(k1)
      ELSE
        CALL Interp_linear(y(k1), x(k1), y(k2), x(k2), u(i), y_int(i))
      END IF
    END DO

  END SUBROUTINE Interpolate_Profile_F1

  !----------------------------------------------------------------------------
  !  Interpolation routine with the interperlation index array not supplied
  !----------------------------------------------------------------------------

  SUBROUTINE Interpolate_Profile_F2(y, x, u, y_int)
    ! Arguments
    REAL(fp), DIMENSION(:),   INTENT(IN)  :: y
    REAL(fp), DIMENSION(:),   INTENT(IN)  :: x
    REAL(fp), DIMENSION(:),   INTENT(IN)  :: u
    REAL(fp), DIMENSION(:),   INTENT(OUT) :: y_int
    ! Local variables
    INTEGER :: interp_index(SIZE(u), 2)

    ! Compute the index array for indexing the two interpolation points
    CALL Compute_Interp_Index(x, u, interp_index)

    ! Perform the interpolation
    CALL Interpolate_Profile_F1(interp_index, y, x, u, y_int)

  END SUBROUTINE Interpolate_Profile_F2

!---------------------------------------------------------------------------------------------
! NAME: Interpolate_Profile_TL
!
! PURPOSE:
!     The Tangent_Linear routine of Interpolate_Profile
! CALLING SEQUENCE:
!            CALL Interpolate_Profile_TL(y, x, u, y_TL, y_int_TL)
!      or    CALL Interpolate_Profile_TL(interp_index, y, x, u, y_TL, y_int_TL)
!
! INPUT ARGUMENTS:
!       y:            The data array to be interpolated.
!                     UNITS:      N/A
!                     TYPE:       fp
!                     DIMENSION:  rank-1
!                     ATTRIBUTES: INTENT(IN)
!
!       x:            The abscissa values for y and they must be monotonically
!                     ascending.
!                     UNITS:      N/A
!                     TYPE:       fp
!                     DIMENSION:  rank-1
!                     ATTRIBUTES: INTENT(IN)
!
!       u:            The abscissa values for the results
!                     and they must be monotonically ascending
!                     UNITS:      N/A
!                     TYPE:       fp
!                     DIMENSION:  rank-1
!                     ATTRIBUTES: INTENT(IN)
!
!       y_TL:         The Tangent-linear data array of y
!                     UNITS:      N/A
!                     TYPE:       fp
!                     DIMENSION:  rank-1
!                     ATTRIBUTES: INTENT(IN)
!
!    interp_index:    The index array of dimension (nu x 2), where nu = SIZE(u)
!                     UNITS:      N/A
!                     TYPE:       Integer
!                     DIMENSION:  rank-2
!                     ATTRIBUTES: INTENT(IN)
!
! OUTPUT ARGUMENTS:
!       y_int_TL:     The Tangent-linear array of y_int
!                     UNITS:      N/A
!                     TYPE:       Integer
!                     DIMENSION:  rank-1
!                     ATTRIBUTES: INTENT(OUT)
!
! RESTRICTIONS:
!     To be efficient, this routine does not check that x and u are both
!     monotonically ascending and the index bounds.
!
! CREATION HISTORY:
!       Written by:     Yong Han, 07-May-2004
!-----------------------------------------------------------------------------------

  SUBROUTINE Interpolate_Profile_F1_TL(interp_index, y, x, u, y_TL, y_int_TL)
    ! Arguments
    INTEGER,  DIMENSION(:,:), INTENT(IN)  :: interp_index
    REAL(fp), DIMENSION(:),   INTENT(IN)  :: y
    REAL(fp), DIMENSION(:),   INTENT(IN)  :: x
    REAL(fp), DIMENSION(:),   INTENT(IN)  :: u
    REAL(fp), DIMENSION(:),   INTENT(IN)  :: y_TL
    REAL(fp), DIMENSION(:),   INTENT(OUT) :: y_int_TL
    ! Local variables
    INTEGER :: i, k1, k2

    DO i = 1, SIZE(u)
      k1 = interp_index(i, 1)
      k2 = interp_index(i, 2)
      IF( k1 == k2)THEN
        y_int_TL(i) = y_TL(k1)
      ELSE
        CALL Interp_linear_TL(y(k1), x(k1), y(k2), x(k2), u(i), y_TL(k1), y_TL(k2), y_int_TL(i))
      END IF
    END DO

  END SUBROUTINE Interpolate_Profile_F1_TL

  SUBROUTINE Interpolate_Profile_F2_TL(y, x, u, y_TL, y_int_TL)
    ! Arguments
    REAL(fp), DIMENSION(:), INTENT(IN)  :: y
    REAL(fp), DIMENSION(:), INTENT(IN)  :: x
    REAL(fp), DIMENSION(:), INTENT(IN)  :: u
    REAL(fp), DIMENSION(:), INTENT(IN)  :: y_TL
    REAL(fp), DIMENSION(:), INTENT(OUT) :: y_int_TL
    ! Local variables
    INTEGER :: interp_index(SIZE(u), 2)

    ! Compute the index array for indexing the two interpolation points
    CALL Compute_Interp_Index(x, u, interp_index)

    ! Perform the TL interpolation
    CALL Interpolate_Profile_F1_TL(interp_index, y, x, u, y_TL, y_int_TL)

  END SUBROUTINE Interpolate_Profile_F2_TL

!---------------------------------------------------------------------------------------------
! NAME: Interpolate_Profile_AD
!
! PURPOSE:
!     The Adjoint routine of Interpolate_Profile
!
! CALLING SEQUENCE:
!            CALL Interpolate_Profile_AD(y, x, u, y_int_AD, y_AD)
!      or    CALL Interpolate_Profile_AD(interp_index, y, x, u, y_int_AD, y_AD)
!
! INPUT ARGUMENTS:
!       y:            The data array to be interpolated.
!                     UNITS:      N/A
!                     TYPE:       fp
!                     DIMENSION:  rank-1
!                     ATTRIBUTES: INTENT(IN)
!
!       x:            The abscissa values for y and they must be monotonically
!                     ascending.
!                     UNITS:      N/A
!                     TYPE:       fp
!                     DIMENSION:  rank-1
!                     ATTRIBUTES: INTENT(IN)
!
!       u:            The abscissa values for the results
!                     and they must be monotonically ascending
!                     UNITS:      N/A
!                     TYPE:       fp
!                     DIMENSION:  rank-1
!                     ATTRIBUTES: INTENT(IN)
!
!       y_int_AD:     The Adjoint array of y_int
!                     UNITS:      N/A
!                     TYPE:       Integer
!                     DIMENSION:  rank-1
!                     ATTRIBUTES: INTENT(IN)
!
!    interp_index:    The index array of dimension (nu x 2), where nu = SIZE(u)
!                     UNITS:      N/A
!                     TYPE:       Integer
!                     DIMENSION:  rank-2
!                     ATTRIBUTES: INTENT(IN)
!
! IN/OUTPUT ARGUMENTS:
!       y_AD:         The Adjoint data array of y
!                     UNITS:      N/A
!                     TYPE:       fp
!                     DIMENSION:  rank-1
!                     ATTRIBUTES: INTENT(IN)
!
! RESTRICTIONS:
!     To be efficient, this routine does not check that x and u are both
!     monotonically ascending and the index bounds.
!
! CREATION HISTORY:
!       Written by:     Yong Han, 07-May-2004
!-----------------------------------------------------------------------------------

  SUBROUTINE Interpolate_Profile_F1_AD(interp_index, y, x, u, y_int_AD, y_AD)
    ! Arguments
    INTEGER,  DIMENSION(:,:), INTENT(IN)      :: interp_index
    REAL(fp), DIMENSION(:),   INTENT(IN)      :: y
    REAL(fp), DIMENSION(:),   INTENT(IN)      :: x
    REAL(fp), DIMENSION(:),   INTENT(IN)      :: u
    REAL(fp), DIMENSION(:),   INTENT(IN OUT)  :: y_int_AD
    REAL(fp), DIMENSION(:),   INTENT(IN OUT)  :: y_AD
    ! Local variables
    INTEGER :: i, k1, k2

    DO i = SIZE(u), 1, -1
      k1 = interp_index(i, 1)
      k2 = interp_index(i, 2)
      IF( k1 == k2)THEN
        y_AD(k1) = y_AD(k1) + y_int_AD(i)
        y_int_AD(i) = ZERO
      ELSE
        CALL Interp_linear_AD(y(k1), x(k1), y(k2), x(k2), u(i), y_int_AD(i), y_AD(k1), y_AD(k2))
      END IF
    END DO

  END SUBROUTINE Interpolate_Profile_F1_AD

  SUBROUTINE Interpolate_Profile_F2_AD(y, x, u, y_int_AD, y_AD)
    ! Arguments
    REAL(fp), DIMENSION(:),   INTENT(IN)      :: y
    REAL(fp), DIMENSION(:),   INTENT(IN)      :: x
    REAL(fp), DIMENSION(:),   INTENT(IN)      :: u
    REAL(fp), DIMENSION(:),   INTENT(IN OUT)  :: y_int_AD
    REAL(fp), DIMENSION(:),   INTENT(IN OUT)  :: y_AD
    ! Local variables
    INTEGER :: interp_index(SIZE(u), 2)

    ! Compute the index array for indexing the two interpolation points
    CALL Compute_Interp_Index(x, u, interp_index)

    ! Perform the AD interpolation
    CALL Interpolate_Profile_F1_AD(interp_index, y, x, u, y_int_AD, y_AD)

  END SUBROUTINE Interpolate_Profile_F2_AD

!---------------------------------------------------------------------------------------------
! NAME: Compute_Interp_Index
!
! PURPOSE:
!    Given x and u that are ascending arrays, it computes an index array, interp_index,
!    such that
!
!      interp_index(i, 1) = 1 and interp_index(i, 2) = 1, if u(i) < x(1)
!      interp_index(i, 1) = nx and interp_index(i, 2) = nx, if u(i) > x(nx),
!                                                          where nx = SIZE(x)
!      x(interp_index(i, 1)) <= u(i) <= x(interp_index(i, 2)) if x(1) <= u(i) <= x(nx)
!
! CALLING SEQUENCE:
!            CALL Compute_Interp_Index(x, u, interp_index)
!
! INPUT ARGUMENTS:
!       x:            The abscissa values for the data to be interpolated and
!                     they must be monotonically ascending.
!                     UNITS:      N/A
!                     TYPE:       fp
!                     DIMENSION:  rank-1
!                     ATTRIBUTES: INTENT(IN)
!
!       u:            The abscissa values on which the data are interpolated
!                     they must be monotonically ascending
!                     the elements of array x.
!                     UNITS:      N/A
!                     TYPE:       fp
!                     DIMENSION:  rank-1
!                     ATTRIBUTES: INTENT(IN)
!
! OUTPUT ARGUMENTS:
!       interp_index: The index array of dimension (nu x 2), where nu = SIZE(u)
!                     UNITS:      N/A
!                     TYPE:       Integer
!                     DIMENSION:  rank-2
!                     ATTRIBUTES: INTENT(OUT)
!
! RESTRICTIONS:
!     To be efficient, this routine does not check that x and u are both
!     monotonically ascending and the index bounds.
!
! CREATION HISTORY:
!       Written by:     Yong Han, 07-May-2004
!-----------------------------------------------------------------------------------

  SUBROUTINE Compute_Interp_Index(x, u, interp_index)
    ! Arguments
    REAL(fp), DIMENSION(:),   INTENT(IN)  :: x
    REAL(fp), DIMENSION(:),   INTENT(IN)  :: u
    INTEGER,  DIMENSION(:,:), INTENT(OUT) :: interp_index
    ! Local variables
    INTEGER :: nx, nu, ix, iu, j, k1, k2

    nx = SIZE(x)
    nu = SIZE(u)

    ! Set the indexes to 1 for the elements in u that are smaller than x(1)
    k1 = nu + 1
    LessThan_Loop: DO iu = 1, nu
      IF(u(iu) < x(1))THEN
        interp_index(iu, 1) = 1
        interp_index(iu, 2) = 1
      ELSE
        k1 = iu
        EXIT LessThan_Loop
      END IF
    END DO LessThan_Loop

    ! Set the indexes to nx for the elements in u that are larger than x(nx)
    k2 = 0
    GreaterThan_Loop: DO iu = nu, k1, -1
      IF(u(iu) > x(nx))THEN
        interp_index(iu, 1) = nx
        interp_index(iu, 2) = nx
      ELSE
        k2 = iu
        EXIT GreaterThan_Loop
      END IF
    END DO GreaterThan_Loop

    ! Set the indexes for the elements in u that are in the range
    ! between x1(1) and x(nx)
    j = 1
    Outer_Loop: DO iu = k1, k2
      Inner_Loop: DO ix = j, nx-1
        IF(u(iu) >= x(ix) .AND. u(iu) <= x(ix+1))THEN
          interp_index(iu, 1) = ix
          interp_index(iu, 2) = ix+1
          j = ix
          EXIT Inner_Loop
        ENDIF
      END DO Inner_Loop
    END DO Outer_Loop

  END SUBROUTINE Compute_Interp_Index

  !---------------------------------------------
  ! Function for two points linear interpolation
  !---------------------------------------------
  SUBROUTINE Interp_linear(y1, x1, y2, x2, x, y)
    REAL(fp), INTENT(IN)  :: y1, x1, y2, x2, x
    REAL(fp), INTENT(OUT) :: y
    y = y1 + (y2-y1)*(x - x1)/(x2 - x1)
  END SUBROUTINE Interp_linear

  SUBROUTINE Interp_linear_TL(y1, x1, y2, x2, x, y1_TL, y2_TL, y_TL)
    REAL(fp), INTENT(IN)  :: y1, x1, y2, x2, x, y1_TL, y2_TL
    REAL(fp), INTENT(OUT) :: y_TL
    y_TL = y1_TL + (y2_TL-y1_TL)*(x - x1)/(x2 - x1)
  END SUBROUTINE Interp_linear_TL

  SUBROUTINE Interp_linear_AD(y1, x1, y2, x2, x, y_AD, y1_AD, y2_AD)
    REAL(fp), INTENT(IN)     :: y1, x1, y2, x2, x
    REAL(fp), INTENT(IN OUT) :: y_AD
    REAL(fp), INTENT(IN OUT) :: y1_AD, y2_AD
    ! Local variables
    REAL(fp) :: fac
    fac = (x - x1)/(x2 - x1)
    y1_AD = y1_AD + y_AD
    y1_AD = y1_AD - fac*y_AD
    y2_AD = y2_AD + fac*y_AD
    y_AD = ZERO
  END SUBROUTINE Interp_linear_AD

  !----------------------------------------------
  ! Function for two points log_log interpolation
  !----------------------------------------------
  SUBROUTINE Interp_loglog(y1, x1, y2, x2, x, y)
    REAL(fp), INTENT(IN)  :: y1, x1, y2, x2, x
    REAL(fp), INTENT(OUT) :: y
    ! Local variables
    REAL(fp) :: xx

    IF ( x1 <= TOLERANCE .OR. x2 <= TOLERANCE .OR. x <= TOLERANCE .OR. &
         y1 <= TOLERANCE .OR. y2 <= TOLERANCE                          ) THEN
      ! Linear interpolation
      CALL Interp_linear(y1, x1, y2, x2, x, y)
    ELSE
      ! Log interpolation
      xx = LOG(x/x1) / LOG(x2/x1)
      y = y1 * (y2/y1)**xx
    END IF

  END SUBROUTINE Interp_loglog

  SUBROUTINE Interp_loglog_TL(y1, x1, y2, x2, x, y1_TL, y2_TL, y_TL)
    REAL(fp), INTENT(IN)  :: y1, x1, y2, x2, x, y1_TL, y2_TL
    REAL(fp), INTENT(OUT) :: y_TL
    ! Local variables
    REAL(fp) :: xx, r, rxx

    IF( x1 <= TOLERANCE .OR. x2 <= TOLERANCE .OR. x <= TOLERANCE .OR. &
        y1 <= TOLERANCE .OR. y2 <= TOLERANCE                          ) THEN
      ! Linear interpolation
      CALL Interp_linear_TL(y1, x1, y2, x2, x, y1_TL, y2_TL, y_TL)
    ELSE
      ! Log interpolation
      xx  = Log(x/x1) / Log(x2/x1)
      r   = (y2/y1)
      rxx = r**xx
      y_TL = (ONE - xx)*rxx*y1_TL + xx*(rxx/r)*y2_TL
    END IF

  END SUBROUTINE Interp_loglog_TL

  SUBROUTINE Interp_loglog_AD(y1, x1, y2, x2, x, y_AD, y1_AD, y2_AD)
    REAL(fp), INTENT(IN)     :: y1, x1, y2, x2, x  ! Input
    REAL(fp), INTENT(IN OUT) :: y_AD               ! Input
    REAL(fp), INTENT(IN OUT) :: y1_AD, y2_AD       ! Output
    ! Local variables
    REAL(fp) :: xx, r, rxx

    IF( x1 <= TOLERANCE .OR. x2 <= TOLERANCE .OR. x <= TOLERANCE .OR. &
        y1 <= TOLERANCE .OR. y2 <= TOLERANCE                          ) THEN
      ! Linear interpolation
      CALL Interp_linear_AD(y1, x1, y2, x2, x, y_AD, y1_AD, y2_AD)
    ELSE
      ! Log interpolation
      xx  = Log(x/x1) / Log(x2/x1)
      r   = (y2/y1)
      rxx = r**xx
      y1_AD = y1_AD + (ONE - xx)*rxx*y_AD
      y2_AD = y2_AD + xx*(rxx/r)*y_AD
      y_AD = ZERO
    END IF

  END SUBROUTINE Interp_loglog_AD

END MODULE CRTM_Profile_Interpolation

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
!
!     (c) 2019 Guide Star Engineering, LLC
!     This Software was developed for the US Nuclear Regulatory Commission (US NRC)
!     under contract "Multi-Dimensional Physics Implementation into Fuel Analysis under
!     Steady-state and Transients (FAST)", contract # NRC-HQ-60-17-C-0007
!
!
!    NEMO - Numerical Engine (for) Multiphysics Operators
! Copyright (c) 2007, Stefano Toninel
!                     Gian Marco Bianchi  University of Bologna
!              David P. Schmidt    University of Massachusetts - Amherst
!              Salvatore Filippone University of Rome Tor Vergata
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification,
! are permitted provided that the following conditions are met:
!
!     1. Redistributions of source code must retain the above copyright notice,
!        this list of conditions and the following disclaimer.
!     2. Redistributions in binary form must reproduce the above copyright notice,
!        this list of conditions and the following disclaimer in the documentation
!        and/or other materials provided with the distribution.
!     3. Neither the name of the NEMO project nor the names of its contributors
!        may be used to endorse or promote products derived from this software
!        without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
! ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
! WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
! ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
!---------------------------------------------------------------------------------
!
! $Id: class_bc_wall.f90 8157 2014-10-09 13:02:44Z sfilippo $
!
! Description:
!    WALL boundary condition class.
!
SUBMODULE(class_bc_wall) class_bc_wall_procedures

    USE class_psblas
    USE class_bc_math
    USE class_material
    USE class_vector
    USE class_face

    IMPLICIT NONE

CONTAINS

    ! REMARK: the implementation of run-time polymorphism requires
    ! specific BC object as POINTERS!

    MODULE PROCEDURE nemo_bc_wall_sizeof
        USE psb_base_mod

        INTEGER(kind=nemo_int_long_)   :: val

        val = 2 * nemo_sizeof_int
        val = val + bc%temp%nemo_sizeof() + SUM(bc%vel%nemo_sizeof())
        nemo_bc_wall_sizeof = val

    END PROCEDURE nemo_bc_wall_sizeof

    ! ----- Constructor -----

    MODULE PROCEDURE create_bc_wall
        IMPLICIT NONE
        !
        INTEGER :: info

        ! Alloc bc_wall target on every process

        IF(ASSOCIATED(bc)) THEN
            WRITE(*,100)
            CALL abort_psblas
        ELSE
            ALLOCATE(bc,stat=info)
            IF(info /= 0) THEN
                WRITE(*,200)
                CALL abort_psblas
            END IF
        END IF

        ! Allocates and sets BC class members, according to the parameters
        ! get from the input file.
        CALL rd_inp_bc_wall(input_file,sec,nbf,bc%id,bc%temp,bc%conc,bc%vel,bc%stress)

100     FORMAT(' ERROR! Illegal call to CREATE_BC_WALL: pointer already associated')
200     FORMAT(' ERROR! Memory allocation failure in CREATE_BC_WALL')

    END PROCEDURE create_bc_wall

    ! ----- Destructor -----

    MODULE PROCEDURE free_bc_wall
        IMPLICIT NONE
        !
        INTEGER :: info

        IF(bc%temp%is_allocated()) CALL bc%temp%dealloc_bc_math()
        IF(ANY(bc%vel%is_allocated()))  THEN
            CALL bc%vel(1)%dealloc_bc_math()
            CALL bc%vel(2)%dealloc_bc_math()
            CALL bc%vel(3)%dealloc_bc_math()
        END IF

        DEALLOCATE(bc,stat=info)
        IF(info /= 0) THEN
            WRITE(*,100)
            CALL abort_psblas
        END IF

100     FORMAT(' ERROR! Memory allocation failure in FREE_BC_WALL')

    END PROCEDURE free_bc_wall

    ! ----- Getter -----

    MODULE PROCEDURE get_abc_wall_s
        USE class_dimensions
        IMPLICIT NONE

        IF(dim == temperature_) THEN
            CALL get_abc_math(bc%temp,id,a,b,c)
        ELSEIF(dim == density_) THEN
            CALL get_abc_math(bc%temp,id,a,b,c)
        ELSE

            WRITE(*,100)
            CALL abort_psblas
        END IF

100     FORMAT(' ERROR! Unsupported field dimensions in GET_ABC_WALL_S')

    END PROCEDURE get_abc_wall_s

    MODULE PROCEDURE get_abc_wall_v
        USE class_dimensions
        IMPLICIT NONE

        IF(dim == velocity_) THEN
            CALL get_abc_math(bc%vel,id,a,b,c)
        ELSE IF(dim == length_) THEN
            CALL get_abc_math(bc%vel,id,a,b,c)
        ELSE IF(dim == pressure_) THEN
            CALL get_abc_math(bc%stress,id,a,b,c)
        ELSE
            WRITE(*,100)
            CALL abort_psblas
        END IF

        ! REMARK: BC(:) elements are supposed to differ only in "C" term

100     FORMAT(' ERROR! Unsupported field dimensions in GET_ABC_WALL_V')

    END PROCEDURE get_abc_wall_v

    ! ----- Setter -----

    MODULE PROCEDURE set_bc_wall_map_s
        USE tools_bc
        USE class_bc_math
        IMPLICIT NONE

        SELECT CASE(bc%id(1))
        CASE(bc_temp_convection_map_)
            CALL bc%temp%set_bc_math_map(i,a,b,c)
        CASE DEFAULT
            WRITE(*,100)
            CALL abort_psblas
        END SELECT

100     FORMAT(' ERROR! Unsupported BC type in SET_BC_WALL_MAP')

    END PROCEDURE set_bc_wall_map_s

    MODULE PROCEDURE set_bc_wall_map_v
        USE tools_bc
        USE class_bc_math
        IMPLICIT NONE

        SELECT CASE(bc%id(2))
        CASE(bc_vel_free_sliding_)
            CALL bc%vel(1)%set_bc_math_map(i,a,b,c%x_())
            CALL bc%vel(2)%set_bc_math_map(i,a,b,c%y_())
            CALL bc%vel(3)%set_bc_math_map(i,a,b,c%x_())
        CASE DEFAULT
            WRITE(*,100)
            CALL abort_psblas
        END SELECT

100     FORMAT(' ERROR! Unsupported BC type in SET_BC_WALL_MAP')

    END PROCEDURE set_bc_wall_map_v

    ! ----- Boundary Values Updater -----

    MODULE PROCEDURE update_boundary_wall_s
        USE class_dimensions
        USE class_material
        USE class_mesh
        USE tools_bc
        IMPLICIT NONE
        !
        INTEGER :: i, id, IF, info, ib_offset, n
        REAL(psb_dpk_), ALLOCATABLE :: a(:), b(:), c(:)

        ! WARNING!
        ! - Use intent(inout) for BX(:).

        ! Number of boundary faces with flag < IB
        ib_offset = COUNT(msh%faces%flag_() > 0 .AND. msh%faces%flag_() < ib)

        n = COUNT(msh%faces%flag_() == ib)

        ALLOCATE(a(n),b(n),c(n),stat=info)
        IF(info /= 0) THEN
            WRITE(*,100)
            CALL abort_psblas
        END IF

        IF(dim == temperature_) THEN
            CALL get_abc_math(bc%temp,id,a,b,c)
            IF(  id == bc_neumann_flux_ .OR. &
                id == bc_robin_convection_) THEN
                DO i = 1, n
                    IF = ib_offset + i
                    CALL matlaw(mats,im(i),bx(IF),conductivity_,b(i))
                END DO
            END IF
        ELSEIF(dim == density_) THEN
            CALL get_abc_math(bc%temp,id,a,b,c)
            IF(  id == bc_neumann_flux_ .OR. &
                id == bc_robin_convection_) THEN
                DO i = 1, n
                    IF = ib_offset + i
                    CALL matlaw(mats,im(i),bx(IF),conductivity_,b(i))
                END DO
            END IF
        ELSE

            WRITE(*,200)
            CALL abort_psblas
        END IF

        CALL apply_abc_to_boundary(id,a,b,c,ib,msh,x,bx)

        DEALLOCATE(a,b,c,stat=info)
        IF(info /= 0) THEN
            WRITE(*,300)
            CALL abort_psblas
        END IF

100     FORMAT(' ERROR! Memory allocation failure in UPDATE_BOUNDARY_WALL_S')
200     FORMAT(' ERROR! Unsupported field dimensions in UPDATE_BOUNDARY_WALL_S')
300     FORMAT(' ERROR! Memory deallocation failure in UPDATE_BOUNDARY_WALL_S')

    END PROCEDURE update_boundary_wall_s

    MODULE PROCEDURE update_boundary_wall_v
        USE class_dimensions
        USE class_mesh
        USE tools_bc
        IMPLICIT NONE
        !
        INTEGER :: id, info, ib_offset, n
        REAL(psb_dpk_), ALLOCATABLE :: a(:), b(:)
        TYPE(vector),     ALLOCATABLE :: c(:)

        ! WARNING!
        ! - Use intent(inout) for BX(:).

        ! Number of boundary faces with flag < IB
        ib_offset = COUNT(msh%faces%flag_() > 0 .AND. msh%faces%flag_() < ib)

        n = COUNT(msh%faces%flag_() == ib)

        ALLOCATE(a(n),b(n),c(n),stat=info)
        IF(info /= 0) THEN
            WRITE(*,100)
            CALL abort_psblas
        END IF

        IF(dim == velocity_ .OR. dim == length_) THEN
            CALL get_abc_math(bc%vel,id,a,b,c)
            CALL apply_abc_to_boundary(id,a,b,c,ib,msh,x,bx)
        ELSE IF(dim == pressure_) THEN
            CALL get_abc_math(bc%stress,id,a,b,c)
            CALL apply_abc_to_boundary(id,a,b,c,ib,msh,x,bx)
        ELSE
            WRITE(*,200)
            CALL abort_psblas
        END IF

        DEALLOCATE(a,b,c,stat=info)
        IF(info /= 0) THEN
            WRITE(*,300)
            CALL abort_psblas
        END IF

100     FORMAT(' ERROR! Memory allocation failure in UPDATE_BOUNDARY_WALL_V')
200     FORMAT(' ERROR! Unsupported field dimensions in UPDATE_BOUNDARY_WALL_V')
300     FORMAT(' ERROR! Memory deallocation failure in UPDATE_BOUNDARY_WALL_V')

    END PROCEDURE update_boundary_wall_v

SUBROUTINE rd_inp_bc_wall(input_file,sec,nbf,id,bc_temp,bc_conc,bc_vel,bc_stress)
    USE class_psblas
    USE json_module, ONLY : json_file
    USE class_bc_math
    USE tools_bc
    USE tools_input

    IMPLICIT NONE
    !
    TYPE(json_file) :: nemo_json
    CHARACTER(len=*), INTENT(IN) :: input_file
    CHARACTER(len=*), INTENT(IN) :: sec
    INTEGER, INTENT(IN) :: nbf
    INTEGER, INTENT(OUT) :: id(:)
    TYPE(bc_math), INTENT(OUT) :: bc_temp
    TYPE(bc_math), INTENT(OUT) :: bc_conc
    TYPE(bc_math), INTENT(OUT) :: bc_vel(3)
    TYPE(bc_math), INTENT(OUT) :: bc_stress(3)
    !
    LOGICAL, PARAMETER :: debug = .FALSE.
    LOGICAL :: found
    !
    INTEGER :: mypnum, icontxt
    INTEGER :: id_sec, inp
    REAL(psb_dpk_) :: work(3,5), wtemp
    CHARACTER(len=15) :: par

    icontxt = icontxt_()
    mypnum  = mypnum_()

    id   = -1 ! DEFAULT value
    work = 0.d0

    IF(mypnum == 0) THEN

        CALL open_file(input_file,nemo_json)
        WRITE(*,*) '- Reading ', TRIM(sec), ' section: type WALL'

        CALL nemo_json%get(trim(sec)//'.temperature.id', id_sec, found)

        IF (.NOT.found) THEN
            WRITE(*,*) 'Temperature BC not found in RD_INP_BC_WALL'
            CALL abort_psblas
        ELSE
            id(bc_temp_) = id_sec
            SELECT CASE(id_sec)
            CASE(bc_temp_fixed_)      ! Fixed temperature
                CALL nemo_json%get(trim(sec)//'.temperature.value', work(1,bc_temp_), found)
            CASE(bc_temp_adiabatic_)  ! Adiabatic wall
                !READ(inp,'()')

            CASE(bc_temp_flux_)       ! Fixed heat flux
                CALL nemo_json%get(trim(sec)//'.temperature.value', work(1,bc_temp_), found)

            CASE(bc_temp_convection_) ! Convection
                CALL nemo_json%get(trim(sec)//'.temperature.value', work(1,bc_temp_), found)
                CALL nemo_json%get(trim(sec)//'.temperature.value2', work(2,bc_temp_), found)

            CASE(bc_temp_convection_map_) ! Convection map
                !READ(inp,'()')         ! To be actually set at the first mapping

            CASE DEFAULT
                WRITE(*,210)
                CALL abort_psblas
            END SELECT
        END IF

        CALL nemo_json%get(trim(sec)//'.concentration.id', id_sec, found)

        IF (.NOT.found) THEN
            WRITE(*,*) 'Concentration BC not found in RD_INP_BC_WALL'
            ! CALL abort_psblas
        ELSE
            id(bc_conc_) = id_sec
            SELECT CASE(id_sec)
            CASE(bc_conc_fixed_)      ! Fixed temperature
                CALL nemo_json%get(trim(sec)//'.concetration.value', work(1,bc_conc_), found)

            CASE(bc_conc_adiabatic_)  ! Adiabatic wall
                !READ(inp,'()')

            CASE DEFAULT
                WRITE(*,210)
                CALL abort_psblas
            END SELECT
        END IF

        CALL nemo_json%get(trim(sec)//'.velocity.id', id_sec, found)

        IF (.NOT.found) THEN
            WRITE(*,*) 'Velocity BC not found in RD_INP_BC_WALL'
        ELSE
            id(bc_vel_) = id_sec

            SELECT CASE(id_sec)
            CASE(bc_vel_no_slip_, bc_vel_free_slip_, bc_vel_free_sliding_)
                !READ(inp,'()')
            CASE(bc_vel_sliding_)
                CALL nemo_json%get(trim(sec)//'.velocity.v1', work(1,bc_vel_), found)
                CALL nemo_json%get(trim(sec)//'.velocity.v2', work(2,bc_vel_), found)
                CALL nemo_json%get(trim(sec)//'.velocity.v3', work(3,bc_vel_), found)
                WRITE(0,*) 'rd_inp_bc_wall: Debug: work', work(1:3,bc_vel_)
            CASE(bc_vel_moving_)
                !READ(inp,'()')
            CASE DEFAULT
                WRITE(*,220)
                CALL abort_psblas
            END SELECT
        END IF

        CALL nemo_json%get(trim(sec)//'.stress.id', id_sec, found)

        IF (.NOT.found) THEN
            WRITE(*,*) 'Stress BC not found in RD_INP_BC_WALL'
        ELSE
            id(bc_stress_) = id_sec

            SELECT CASE(id_sec)
            CASE(bc_stress_free_)
                !READ(inp,'()')
            CASE(bc_stress_prescribed_)
                CALL nemo_json%get(trim(sec)//'.stress.v1', work(1,bc_stress_), found)
                CALL nemo_json%get(trim(sec)//'.stress.v2', work(2,bc_stress_), found)
                CALL nemo_json%get(trim(sec)//'.stress.v3', work(3,bc_stress_), found)
                WRITE(0,*) 'rd_inp_bc_wall: Debug: work', work(1:3,bc_stress_)
            CASE DEFAULT
                WRITE(*,230)
                CALL abort_psblas
            END SELECT
        END IF
    END IF

    ! Broadcast
    CALL psb_bcast(icontxt,id)
    CALL psb_bcast(icontxt,work)

    ! TEMPERATURE section
    SELECT CASE(id(bc_temp_))
    CASE(bc_temp_fixed_)      ! Fixed temperature
        CALL bc_temp%alloc_bc_math(bc_dirichlet_,nbf,&
            & a=(/1.d0/),b=(/0.d0/),c=(/work(1,bc_temp_)/))

    CASE(bc_temp_adiabatic_)  ! Adiabatic wall
        CALL bc_temp%alloc_bc_math(bc_neumann_,nbf,&
            & a=(/0.d0/),b=(/1.d0/),c=(/0.d0/))

    CASE(bc_temp_flux_)       ! Fixed heat flux
        CALL bc_temp%alloc_bc_math(bc_neumann_flux_,nbf,&
            & a=(/0.d0/),b=(/1.d0/),c=(/work(1,bc_temp_)/))

    CASE(bc_temp_convection_) ! Convection
        CALL bc_temp%alloc_bc_math(bc_robin_convection_,nbf,&
            & a=(/work(1,bc_temp_)/),b=(/1.d0/),&
            & c=(/work(1,bc_temp_)*work(2,bc_temp_)/))

    CASE(bc_temp_convection_map_) ! Convection map
        CALL bc_temp%alloc_bc_math(bc_robin_map_,nbf,&
            & a=(/(0.d0, inp=1,nbf)/),b=(/(1.d0, inp=1,nbf)/),&
            & c=(/(0.d0, inp=1,nbf)/))

    END SELECT

    ! VELOCITY section
    SELECT CASE(id(bc_vel_))
    CASE(bc_vel_no_slip_)

        CALL bc_vel(1)%alloc_bc_math(bc_dirichlet_,nbf,&
            & a=(/1.d0/),b=(/0.d0/),c=(/0.d0/))
        CALL bc_vel(2)%alloc_bc_math(bc_dirichlet_,nbf,&
            & a=(/1.d0/),b=(/0.d0/),c=(/0.d0/))
        CALL bc_vel(3)%alloc_bc_math(bc_dirichlet_,nbf,&
            & a=(/1.d0/),b=(/0.d0/),c=(/0.d0/))

    CASE(bc_vel_free_slip_)

        CALL bc_vel(1)%alloc_bc_math(bc_neumann_,nbf,&
            & a=(/0.d0/),b=(/1.d0/),c=(/0.d0/))
        CALL bc_vel(2)%alloc_bc_math(bc_neumann_,nbf,&
            & a=(/0.d0/),b=(/1.d0/),c=(/0.d0/))
        CALL bc_vel(3)%alloc_bc_math(bc_neumann_,nbf,&
            & a=(/0.d0/),b=(/1.d0/),c=(/0.d0/))

    CASE(bc_vel_sliding_)
        WRITE(0,*) 'rd_inp_bc_wall: Debug: setting', work(1:3,bc_vel_)
        CALL bc_vel(1)%alloc_bc_math(bc_dirichlet_,nbf,&
            & a=(/1.d0/),b=(/0.d0/),c=(/work(1,bc_vel_)/))
        CALL bc_vel(2)%alloc_bc_math(bc_dirichlet_,nbf,&
            & a=(/1.d0/),b=(/0.d0/),c=(/work(2,bc_vel_)/))
        CALL bc_vel(3)%alloc_bc_math(bc_dirichlet_,nbf,&
            & a=(/1.d0/),b=(/0.d0/),c=(/work(3,bc_vel_)/))
        WRITE(0,*) 'rd_inp_bc_wall: Debug: set'

    CASE(bc_vel_free_sliding_)

        CALL bc_vel(1)%alloc_bc_math(bc_dirichlet_map_,nbf,&
            & a=(/1.d0/),b=(/0.d0/),c=(/(0.d0, inp=1,nbf)/))
        CALL bc_vel(2)%alloc_bc_math(bc_dirichlet_map_,nbf,&
            & a=(/1.d0/),b=(/0.d0/),c=(/(0.d0, inp=1,nbf)/))
        CALL bc_vel(3)%alloc_bc_math(bc_dirichlet_map_,nbf,&
            & a=(/1.d0/),b=(/0.d0/),c=(/(0.d0, inp=1,nbf)/))
    CASE(-1)
    CASE DEFAULT
        WRITE(0,*) 'Fix here unimplemented BC setup !!'

    END SELECT

    ! Stress section
    SELECT CASE(id(bc_stress_))
    CASE(bc_stress_free_)
        CALL bc_stress(1)%alloc_bc_math(bc_neumann_,nbf,&
            & a=(/0.d0/),b=(/1.d0/),c=(/0.d0/))
        CALL bc_stress(2)%alloc_bc_math(bc_neumann_,nbf,&
            & a=(/0.d0/),b=(/1.d0/),c=(/0.d0/))
        CALL bc_stress(3)%alloc_bc_math(bc_neumann_,nbf,&
            & a=(/0.d0/),b=(/1.d0/),c=(/0.d0/))

    CASE(bc_stress_prescribed_)
        CALL bc_stress(1)%alloc_bc_math(bc_dirichlet_,nbf,&
            & a=(/1.d0/),b=(/0.d0/),c=(/work(1,bc_stress_)/))
        CALL bc_stress(2)%alloc_bc_math(bc_dirichlet_,nbf,&
            & a=(/1.d0/),b=(/0.d0/),c=(/work(2,bc_stress_)/))
        CALL bc_stress(3)%alloc_bc_math(bc_dirichlet_,nbf,&
            & a=(/1.d0/),b=(/0.d0/),c=(/work(3,bc_stress_)/))

    CASE DEFAULT
        WRITE(0,*) 'Fix here unimplemented BC setup !!', id(bc_stress_)

    END SELECT

    ! ***
    ! If ID = -1 (Default)  the corresponding BC will not be allocated
    ! ***

    IF(debug) THEN
        WRITE(*,*)
        WRITE(*,400) mypnum
        WRITE(*,500) TRIM(sec),' - Type: Wall'

        WRITE(*,600) ' * TEMPERATURE Section *'
        WRITE(*,700) '  BC%id(bc_temp_) = ', id(bc_temp_)
        IF (id(bc_temp_) > 0) CALL bc_temp%debug_bc_math()

        WRITE(*,600) ' * VELOCITY Section *'
        WRITE(*,700) '  BC%id(bc_vel_) = ', id(bc_vel_)
        IF (id(bc_vel_) > 0) THEN
            CALL bc_vel(1)%debug_bc_math()
            CALL bc_vel(2)%debug_bc_math()
            CALL bc_vel(3)%debug_bc_math()
            WRITE(*,*)
        END IF
    END IF

100 FORMAT(a,i1)
210 FORMAT(' ERROR! Unsupported ID(BC_TEMP_) in RD_INP_BC_WALL')
220 FORMAT(' ERROR! Unsupported ID(BC_VEL_) in RD_INP_BC_WALL')
230 FORMAT(' ERROR! Unsupported ID(BC_STRESS_) in RD_INP_BC_WALL')
999 FORMAT(' ERROR! Missing boundary condition in RC_INP_BC_WALL')

400 FORMAT(' ----- Process ID = ',i2,' -----')
500 FORMAT(1x,a,a)
600 FORMAT(1x,a)
700 FORMAT(1x,a,i2)

END SUBROUTINE rd_inp_bc_wall

END SUBMODULE class_bc_wall_procedures

max_stars_repo_name: sourceryinstitute/MORFEUS-Source

max_issues_repo_name: sourceryinstitute/MORFEUS-Source

max_forks_repo_name: sourceryinstitute/MORFEUS-Source


!----------------------------------------------------------------------------------------------
! MODULE: write_all
!> @author Marin Sapunar, Ruđer Bošković Institute
!> @date May, 2019
!
!> @brief Contains subroutines for writing data to a variety of quantum chemistry formats.
!----------------------------------------------------------------------------------------------
module write_all_mod
    use global_defs
    use molden_write_mod
    implicit none

    private
    public :: write_all

contains

    !----------------------------------------------------------------------------------------------
    ! SUBROUTINE: write_all
    !> @brief Write all available data to a given format.
    !< @todo make input variables optional
    !----------------------------------------------------------------------------------------------
    subroutine write_all(dformat, path, title, geom, atom_symbol, atom_number, basis, mos)
        use element_mod
        use basis_set_mod, only : basis_set
        use molecular_orbitals_mod, only : molecular_orbitals
        character(len=*), intent(in) :: dformat
        character(len=*), intent(in) :: path
        character(len=*), intent(in), optional :: title
        real(dp), intent(in) :: geom(:)
        character(len=2), intent(in) :: atom_symbol(:)
        integer, intent(in) :: atom_number(:)
        type(basis_set), intent(in) :: basis
        type(molecular_orbitals), intent(in) :: mos
        integer :: outunit

        select case(dformat)
        case('molden_cart')
            open(newunit=outunit, file=path, action='write')
            write(outunit, '(a)') '[Molden Format]'
            write(outunit, '(a)') '[Title]'
            if (present(title)) then
                write(outunit, '(a)') title
            else
                write(outunit, *) 'TheoChemUtil write'
            end if
            call molden_write_atoms(outunit, geom, atom_symbol, atom_number)
            call molden_write_gto(outunit, basis)
            call molden_write_mo(outunit, mos)
        case default
            write(stderr, *) 'Error in write_all subroutine.'
            write(stderr, *) 'Input format not implemented: '//trim(adjustl(dformat))
        end select
    end subroutine write_all

end module write_all_mod


   module universe
   use iso_fortran_env, only: real32, int64, output_unit
   implicit none
   private
   public      pr,  pp2d, stack, particle, pi, latt_const

      integer, parameter      :: pr = real32

      real(pr), parameter     :: pi = 3.14159265359_pr

      type particle
          integer             :: name
          real(pr)            :: pos(2)
      contains
          procedure           :: move
      end type particle

      interface cube
         procedure               init_cube
      end interface cube

      type cube
         integer              :: cap
         integer              :: occ
         integer, allocatable :: idx(:)
         type(particle),                 &
                  allocatable :: hld(:)
         integer              :: su(0:8)
         integer              :: cc(2)
         real(pr)             :: llim(2)
         real(pr)             :: ulim(2)
         type(cube), pointer  :: sur1, sur2,&
                                 sur3, sur4
         integer              :: tag
      contains
         procedure            :: add    => add_cube
         procedure            :: pop    => pop_cube
         procedure            :: ext    => ext_cube
         procedure            :: shrink => shrink_cube
         procedure            :: show   => show_cube
         procedure            :: rel    => relative_cube
         procedure            :: abs    => absolute_cube
      end type cube

      interface pp2d
         procedure               init_pp2d, load_pp2d, hex_pp2d, read_lammps_pp2d
      end interface

      type pp2d
         real(pr)             :: ll(2), lx, ly
         real(pr)             :: ww(2), wx, wy
         integer              :: nxy(2), nx, ny
         integer              :: noc
         integer              :: nop
         type(cube), &
          allocatable         :: cup(:)
         integer              :: stat
         integer, allocatable :: cidx(:)
         integer, allocatable :: kidx(:)
         integer              :: lnop
      contains
         procedure            :: assoc    => associate_pp2d
         procedure            :: add      => add_pp2d
         procedure            :: relative => relative_pp2d
         procedure            :: absolute => absolute_pp2d
         procedure            :: shrink   => shrink_pp2d
         procedure            :: write    => write_pp2d
         procedure            :: reserve  => reserve_pp2d
         procedure            :: stage    => active_pp2d
         procedure            :: original => original_pp2d
         procedure            :: random   => random_pp2d
         procedure            :: randcup
         procedure            :: move     => move_pp2d
         procedure            :: zoom_on  => scube_pp2d
         procedure            :: zoom_u   => scube_up_pp2d
         procedure            :: zoom_r   => scube_right_pp2d
         procedure            :: look_up
         procedure            :: metro    => metropolis_pp2d
      end type pp2d

      integer, parameter      :: su_cubes(2,8) = &
                                 reshape( [1,0,0,1,-1,0,0,-1, &    ! 7 2 6
                                          1,-1,1,1,-1,1,-1,-1], &  ! 3 0 1
                                          [2,8] )                  ! 8 4 5
      real(pr)                :: su_ww(2,8)

      integer, parameter      :: nstk=1000
      type stack
         integer              :: n
         integer              :: id(nstk)
         real(pr)             :: x(nstk)
         real(pr)             :: y(nstk)
         real(pr)             :: r(nstk)
         real(pr)             :: f(nstk)
         real(pr)             :: g(nstk)
         real(pr)             :: h(nstk)
      contains
         procedure            :: xx => x_stack
         procedure            :: yy => y_stack
         procedure            :: hexorder
         procedure            :: show => show_stack
         procedure            :: rotate => rotate_stack
      end type stack

   contains

      !//////////////////////////////////////
      !//////////////////////////////////////           particle
      !//////////////////////////////////////

      subroutine move(part,dir,delta)
      implicit none
      class(particle), intent(inout) :: part
      integer,           intent(in)  :: dir
      real(pr),          intent(in)  :: delta
      part%pos(dir) = part%pos(dir) + delta
      end subroutine move

      !//////////////////////////////////////
      !//////////////////////////////////////           cube
      !//////////////////////////////////////

      function init_cube(cap) result(cub)
      implicit none
      integer, intent(in) :: cap
      type(cube)          :: cub
      cub%cap = cap
      cub%occ = 0
      allocate( cub%idx(cap) )
      allocate( cub%hld(cap) )
      end function init_cube

      subroutine ext_cube(cub)
      implicit none
      class(cube), intent(inout)  :: cub
      integer,        allocatable :: idx(:)
      type(particle), allocatable :: hld(:)
      allocate(idx,source=[cub%idx,0])
      deallocate(cub%idx)
      call move_alloc(idx,cub%idx)
      allocate(hld,source=[cub%hld,cub%hld(1)])
      deallocate(cub%hld)
      call move_alloc(hld,cub%hld)
      cub%cap = cub%cap + 1
      end subroutine ext_cube

      subroutine shrink_cube(cub)
      implicit none
      class(cube), intent(inout)  :: cub
      integer,        allocatable :: idx(:)
      type(particle), allocatable :: hld(:)
      integer                     :: occ
      occ = cub%occ
      allocate(idx,source=cub%idx(1:occ))
      deallocate(cub%idx)
      call move_alloc(idx,cub%idx)
      allocate(hld,source=cub%hld(1:occ))
      deallocate(cub%hld)
      call move_alloc(hld,cub%hld)
      cub%cap = occ
      end subroutine shrink_cube

      subroutine add_cube(cub,idx,p)
      implicit none
      class(cube), intent(inout) :: cub
      integer,        intent(in) :: idx
      type(particle), intent(in) :: p
      integer                    :: k
      k = cub%occ + 1
      if( k>cub%cap ) call cub%ext()
      cub%hld(k) = p
      cub%idx(k) = idx
      cub%occ = k
      end subroutine add_cube

      subroutine pop_cube(cub,k,idx)
      implicit none
      class(cube), intent(inout) :: cub
      integer,       intent(in)  :: k
      integer,      intent(out)  :: idx
      idx = cub%idx(cub%occ)
      if( k<cub%occ ) then
         cub%idx(k) = idx
         cub%hld(k) = cub%hld(cub%occ)
      end if
      cub%occ = cub%occ - 1
      end subroutine pop_cube

      subroutine show_cube(cub)
      implicit none
      class(cube), intent(in) :: cub
      integer                 :: i
      write(*,*)
      write(*,'(i3,a,i3)') cub%occ, "/", cub%cap
      do i=1, cub%occ
         write(*,'(i7,6x,i7,2f8.3)') cub%idx(i), &
                   cub%hld(i)%name, &
                   cub%hld(i)%pos
      end do
      write(*,*)
      end subroutine show_cube

      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ rel/abs

      subroutine relative_cube(cub)
      implicit none
      class(cube), intent(inout) :: cub
      cub%hld(1:cub%occ)%pos(1) = &
      cub%hld(1:cub%occ)%pos(1) - cub%llim(1)
      cub%hld(1:cub%occ)%pos(2) = &
      cub%hld(1:cub%occ)%pos(2) - cub%llim(2)
      end subroutine relative_cube

      subroutine absolute_cube(cub)
      implicit none
      class(cube), intent(inout) :: cub
      cub%hld(1:cub%occ)%pos(1) = &
      cub%hld(1:cub%occ)%pos(1) + cub%llim(1)
      cub%hld(1:cub%occ)%pos(2) = &
      cub%hld(1:cub%occ)%pos(2) + cub%llim(2)
      end subroutine absolute_cube

      !///////////////////////////////////////////
      !///////////////////////////////////////////        pp2d
      !///////////////////////////////////////////

      function init_pp2d(ll,wth) result(pos)
      implicit none
      type(pp2d)           :: pos
      real(pr), intent(in) :: ll(2), wth
      integer              :: cx, cy, cl, i, j, k
      pos%ll = ll; pos%lx = ll(1); pos%ly = ll(2)
      pos%nxy = floor(ll/wth); pos%nx = pos%nxy(1); pos%ny = pos%nxy(2)
      pos%ww = ll/pos%nxy; pos%wx = pos%ww(1); pos%wy = pos%ww(2)
      pos%noc = pos%nx*pos%ny
      allocate(pos%cup(0:pos%noc-1))
      do cy = 0, pos%ny-1
         do cx = 0, pos%nx-1
            cl = cy*pos%nx + cx
            pos%cup(cl) = cube(1)
            pos%cup(cl)%cc = (/cx,cy/)
            pos%cup(cl)%llim = (/ cx   *pos%wx, cy   *pos%wy/)
            pos%cup(cl)%ulim = (/(cx+1)*pos%wx,(cy+1)*pos%wy/)
            pos%cup(cl)%su(0) = cl
            do k = 1, 8
               i = mod(cx+su_cubes(1,k),pos%nx)
               j = mod(cy+su_cubes(2,k),pos%ny)
               if(i<0) i = i + pos%nx
               if(j<0) j = j + pos%ny
               pos%cup(cl)%su(k) =  j*pos%nx + i
            end do
         end do
      end do
      pos%nop = 0
      do k = 1, 8
         su_ww(:,k) = su_cubes(:,k)*pos%ww
      end do
      pos%stat = 1
      call pos%assoc()
      end function init_pp2d

      subroutine associate_pp2d(pos)
      implicit none
      class(pp2d), intent(inout), target :: pos
      integer cl
      do cl = 0, pos%noc-1
         pos%cup(cl)%sur1 => pos%cup(pos%cup(cl)%su(1))
         pos%cup(cl)%sur2 => pos%cup(pos%cup(cl)%su(2))
         pos%cup(cl)%sur3 => pos%cup(pos%cup(cl)%su(3))
         pos%cup(cl)%sur4 => pos%cup(pos%cup(cl)%su(4))
      end do
      end subroutine associate_pp2d

      subroutine add_pp2d(pos,idx,p)
      implicit none
      class(pp2d), intent(inout)    :: pos
      integer,       intent(in)     :: idx
      type(particle), intent(inout) :: p
      integer                       :: cc(2), cl
      type(particle)                :: rp
      p%pos = mod(p%pos,pos%ll)
      where( p%pos<0.0_pr ) p%pos = p%pos+pos%ll
      cc = floor(p%pos/pos%ww)
      ! *** bug *** when r = l-eps; floor(r/w)=floor(l/w)=nxy exceeds dims
      ! fixed by this line:
      where( cc>=pos%nxy ) cc=pos%nxy-1
      ! end bug
      cl = cc(2)*pos%nx + cc(1)
      if( pos%stat==0 ) then
         call pos%cup(cl)%add(idx,p)
      else if( pos%stat==1 ) then
         rp = p
         rp%pos = rp%pos - pos%cup(cl)%llim
         call pos%cup(cl)%add(idx,rp)
      end if
      pos%nop = pos%nop + 1
      end subroutine add_pp2d

      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      function load_pp2d(handle, wth) result(pos)
      implicit none
      class(*),         intent(in)    :: handle
      real(pr),         intent(in)    :: wth
      type(pp2d)                      :: pos
      integer                         :: uin, i, num
      real(pr)                        :: lx, ly
      type(particle)                  :: p
      logical                         :: doclose

      doclose = .false.
      select type(handle)
      type is (integer)
         uin = handle
      type is (character(len=*))
         open( newunit=uin, file=handle )
         doclose = .true.
      class default
         stop
      end select

      read(uin,*)
      read(uin,*)
      read(uin,*)
      read(uin,*) num
      read(uin,*) lx, ly
      pos = pp2d([lx,ly],wth)
      do i = 0, num-1
          read(uin,*) p
          call pos%add(i,p)
      end do

      if( doclose ) close(uin)
      end function load_pp2d

      subroutine write_pp2d(pos,handle,string,ints,reals)
      implicit none
      class(pp2d), intent(in)        :: pos
      class(*),     intent(in)       :: handle
      character(len=*), intent(in), &
                         optional    :: string
      class(*), intent(in), optional :: ints(:)
      real(pr), intent(in), optional :: reals(:)
      integer                        :: i, j, uout
      !
      select type (handle)
      type is (integer)
         uout = handle
      type is (character(len=*))
         open(newunit=uout,file=handle)
      end select
      ! header
      if( present(string) ) then
         write(uout,*) string
      else
         write(uout,*)
      end if
      if( present(ints) ) then
         select type (ints)
         type is (integer(int64))
            write(uout,*) ints
         type is (integer)
            write(uout,*) ints
         end select
      else
         write(uout,*)
      end if
      if( present(reals) ) then
         write(uout,*) reals
      else
         write(uout,*)
      end if
      ! body
      write(uout,*) pos%nop
      write(uout,*) pos%ll
      if( pos%stat==0) then
         do i = 0, pos%noc-1
            do j = 1, pos%cup(i)%occ
               write(uout,*) pos%cup(i)%hld(j)
            end do
         end do
      else if( pos%stat==1 ) then
          do i = 0, pos%noc-1
            do j = 1, pos%cup(i)%occ
               write(uout,*) pos%cup(i)%hld(j)%name, &
               pos%cup(i)%llim + pos%cup(i)%hld(j)%pos
            end do
         end do
      end if
      !
      select type (handle)
      type is (character(len=*))
         close(uout)
      end select
      end subroutine write_pp2d

      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ lattices

      function hex_pp2d(nxy,rho,wth,f_distort) result(pos)
      implicit none
      integer,  intent(in) :: nxy(2)
      real(pr), intent(in) :: rho, wth
      real(pr), intent(in),   &
             optional      :: f_distort
      type(pp2d)           :: pos
      real(pr)             :: v1(2), v2(2), a, ax, ay, rnd(2), delta
      integer              :: i, j, idx
      type(particle)       :: part
      ax = 0.5_pr
      ay = sqrt(3.0_pr)/2
      a = sqrt(1.0_pr/(rho*ay))
      v1 = a*[1.0_pr,0.0_pr]
      v2 = a*[ax,ay]
      pos = pp2d(nxy*[v1(1),v2(2)],wth)
      if( present(f_distort) ) then
              delta = f_distort
      else
              delta = 0.0_pr
      end if
      idx = -1
      do j = 0, nxy(2)-1
         do i = 0, nxy(1)-1
            call random_number(rnd)
            idx = idx + 1
            part%name = idx
            part%pos  = i*v1 + j*v2 + ( rnd(1)*v1 + rnd(2)*v2 )*delta
            call pos%add(idx, part)
         end do
      end do
      end function hex_pp2d

      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ read lammps

      function read_lammps_pp2d(fname, wth,lammps) result(pos)
      implicit none
      character(len=*), intent(in)    :: fname,lammps
      real(pr),         intent(in)    :: wth
      type(pp2d)                      :: pos
      integer                         :: uin, i, num
      real(pr)                        :: lx, ly, skip
      type(particle)                  :: p
      if(lammps /= "lammps") then
         write(*,*)"not lammps";stop
      end if
      open( newunit=uin, file=fname )
         read(uin,*)
         read(uin,*)
         read(uin,*)
         read(uin,*) num
         read(uin,*)
         read(uin,*) skip, lx
         read(uin,*) skip, ly
         read(uin,*)
         read(uin,*)
         pos = pp2d([lx,ly],wth)
         do i = 0, num-1
             read(uin,*) p
             call pos%add(i,p)
         end do
      close(uin)
      end function read_lammps_pp2d

      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      subroutine relative_pp2d(pos)
      implicit none
      class(pp2d), intent(inout) :: pos
      integer                    :: cl
      if( pos%stat==0 ) then
         do cl = 0, pos%noc-1
            call pos%cup(cl)%rel()
         end do
         pos%stat = 1
      end if
      end subroutine relative_pp2d

      subroutine absolute_pp2d(pos)
      implicit none
      class(pp2d), intent(inout) :: pos
      integer                    :: cl
      if( pos%stat==1 ) then
         do cl = 0, pos%noc-1
            call pos%cup(cl)%abs()
         end do
         pos%stat = 0
      end if
      end subroutine absolute_pp2d

      subroutine shrink_pp2d(pos)
      implicit none
      class(pp2d), intent(inout) :: pos
      integer                    :: cl
      do cl = 0, pos%noc-1
         call pos%cup(cl)%shrink()
      end do
      end subroutine shrink_pp2d

      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      subroutine reserve_pp2d(pos,maxnum)
      implicit none
      class(pp2d), intent(inout)    :: pos
      integer, intent(in), optional :: maxnum
      integer                       :: ub
      if( present(maxnum) ) then
         ub = maxnum - 1
      else
         ub = pos%nop - 1
      end if
      if( .not.allocated(pos%cidx) ) then
          allocate(pos%cidx(0:ub))
          allocate(pos%kidx(0:ub))
      else if( size(pos%cidx)<ub+1 ) then
          deallocate(pos%cidx)
          deallocate(pos%kidx)
          allocate(pos%cidx(0:ub))
          allocate(pos%kidx(0:ub))
      end if
      end subroutine reserve_pp2d

      subroutine active_pp2d(pos,c,cw)
      implicit none
      class(pp2d), intent(inout) :: pos
      integer,      intent(in),&
                      optional   :: c(2), cw(2)
      integer                    :: c1(2), c2(2)
      integer                    :: cx, cy, cl
      integer                    :: k, idx, kx, ky
      if(present(c) .and. present(cw)) then
         c1 = c; c2 = c1+cw-1
      else
         c1 = [0,0]; c2 = pos%nxy-1
      end if
      ! first count lnop and reserve
      idx = 0
      do cy = c1(2), c2(2)
         ky = mod(cy,pos%ny)
         if( ky<0 ) ky = ky + pos%ny
         do cx = c1(1), c2(1)
            kx = mod(cx,pos%nx)
            if( kx<0 ) kx = kx + pos%nx
            cl = ky*pos%nx + kx
            idx = idx + pos%cup(cl)%occ
         end do
      end do
      call pos%reserve(idx)
      pos%lnop = idx
      ! assign tags then
      pos%cup%tag = 1
      idx = -1
      do cy = c1(2), c2(2)
         ky = mod(cy,pos%ny)
         if( ky<0 ) ky = ky + pos%ny
         do cx = c1(1), c2(1)
            kx = mod(cx,pos%nx)
            if( kx<0 ) kx = kx + pos%nx
            cl = ky*pos%nx + kx
            pos%cup(cl)%tag = 0
            do k = 1, pos%cup(cl)%occ
               idx = idx + 1
               pos%cup(cl)%idx(k) = idx
               pos%cidx(idx) = cl
               pos%kidx(idx) = k
            end do
         end do
      end do
      end subroutine active_pp2d

      subroutine original_pp2d(pos)
      implicit none
      class(pp2d), intent(inout) :: pos
      integer                    :: cl, k, idx
      call pos%reserve()
      pos%cup%tag = 0
      do cl = 0, pos%noc-1
         do k = 1, pos%cup(cl)%occ
            idx = pos%cup(cl)%hld(k)%name
            pos%cup(cl)%idx(k) = idx
            pos%cidx(idx) = cl
            pos%kidx(idx) = k
         end do
      end do
      pos%lnop = pos%nop
      end subroutine original_pp2d

      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      recursive &
      subroutine random_pp2d(pos,dmin,dmax,idx,dir,delta)
      implicit none
      class(pp2d), intent(in) :: pos
      real(pr),    intent(in) :: dmin, dmax
      integer,    intent(out) :: idx, dir
      real(pr),   intent(out) :: delta
      real(pr)                :: rand_pr
      integer                 :: ncl, cl, k
      call random_number(rand_pr)
      idx = floor(rand_pr*pos%lnop)
      cl = pos%cidx(idx)
      k = pos%kidx(idx)
      call random_number(rand_pr)
      if(rand_pr<0.5_pr) then
         dir = 1
      else
         dir = 2
      end if
      call random_number(rand_pr)
      delta = dmin + rand_pr*(dmax-dmin)
      ! check if forbiden
      rand_pr = pos%cup(cl)%hld(k)%pos(dir) + delta
      if( rand_pr>=pos%ww(dir) ) then
         ncl = pos%cup(cl)%su(dir)
         if( pos%cup(ncl)%tag==1 ) &
         call pos%random(dmin,dmax,idx,dir,delta)
      else if( rand_pr<0.0_pr ) then
         ncl = pos%cup(cl)%su(dir+2)
         if( pos%cup(ncl)%tag==1 ) &
         call pos%random(dmin,dmax,idx,dir,delta)
      end if
      end subroutine random_pp2d

      function randcup(pos) result(cc)
      implicit none
      class(pp2d), intent(in) :: pos
      integer                 :: cc(2)
      real(pr)                :: rnd(2)
      call random_number(rnd)
      cc = floor(rnd*pos%nxy)
      end function randcup

      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      subroutine look_up(pos,idx,xy,env)
      implicit none
      class(pp2d), intent(in), target :: pos
      integer,   intent(in), optional :: idx
      real(pr),      intent(in)       :: xy(2)
      type(stack), intent(out)        :: env
      integer                         :: i, j, cl, cc(2)
      real(pr)                        :: r(2)
      type(cube), pointer             :: ccub
      ! find r
      if( present(idx) ) then
         cl = pos%cidx(idx)
         j = pos%kidx(idx)
         r = pos%cup(cl)%llim + pos%cup(cl)%hld(j)%pos + xy
      else
         r = xy
      end if
      ! find cl of r
      r = mod(r,pos%ll)
      where( r<0.0_pr ) r = r+pos%ll
      cc = floor(r/pos%ww)
      where( cc>=pos%nxy ) cc=pos%nxy-1
      cl = cc(2)*pos%nx + cc(1)
      ! standardize r
      ccub => pos%cup(cl)
      r = r - ccub%llim
      ! look around then
      env%n = 0
      do i = 1, ccub%occ
          env%n = env%n + 1
          env%id(env%n) = ccub%idx(i)
          env%x(env%n) = ccub%hld(i)%pos(1)
          env%y(env%n) = ccub%hld(i)%pos(2)
      end do
      ccub => ccub%sur1
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(1)<r(1) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) + su_ww(1,1)
                   env%y(env%n) = ccub%hld(i)%pos(2) !+ su_ww(2,1)
                end if
             end do
      ccub => ccub%sur2
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(1)<r(1) .and. &
                     ccub%hld(i)%pos(2)<r(2) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) + su_ww(1,6)
                   env%y(env%n) = ccub%hld(i)%pos(2) + su_ww(2,6)
                end if
             end do
      ccub => ccub%sur3
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(2)<r(2) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) !+ su_ww(1,2)
                   env%y(env%n) = ccub%hld(i)%pos(2) + su_ww(2,2)
                end if
             end do
      ccub => ccub%sur3
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(1)>r(1) .and. &
                     ccub%hld(i)%pos(2)<r(2) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) + su_ww(1,7)
                   env%y(env%n) = ccub%hld(i)%pos(2) + su_ww(2,7)
                end if
             end do
      ccub => ccub%sur4
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(1)>r(1) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) + su_ww(1,3)
                   env%y(env%n) = ccub%hld(i)%pos(2) !+ su_ww(2,3)
                end if
             end do
      ccub => ccub%sur4
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(1)>r(1) .and. &
                     ccub%hld(i)%pos(2)>r(2) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) + su_ww(1,8)
                   env%y(env%n) = ccub%hld(i)%pos(2) + su_ww(2,8)
                end if
             end do
      ccub => ccub%sur1
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(2)>r(2) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) !+ su_ww(1,4)
                   env%y(env%n) = ccub%hld(i)%pos(2) + su_ww(2,4)
                end if
             end do
      ccub => ccub%sur1
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(1)<r(1) .and. &
                     ccub%hld(i)%pos(2)>r(2) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) + su_ww(1,5)
                   env%y(env%n) = ccub%hld(i)%pos(2) + su_ww(2,5)
                end if
             end do
      env%x(1:env%n) = env%x(1:env%n)-r(1)
      env%y(1:env%n) = env%y(1:env%n)-r(2)
      end subroutine look_up

      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      subroutine scube_pp2d(pos,idx,env)
      implicit none
      class(pp2d), intent(in), &
                       target   :: pos
      integer,       intent(in) :: idx
      type(stack), intent(out)  :: env
      integer                   :: i, cl, k
      real(pr)                  :: r(2)
      type(cube), pointer       :: ccub
      cl = pos%cidx(idx)
      k = pos%kidx(idx)
      ccub => pos%cup(cl)
      r = ccub%hld(k)%pos
      env%n = 0
      do i = 1, ccub%occ
         if(i/=k) then
                 env%n = env%n + 1
                 env%id(env%n) = ccub%idx(i)
                 env%x(env%n) = ccub%hld(i)%pos(1)
                 env%y(env%n) = ccub%hld(i)%pos(2)
         end if
      end do
      ccub => ccub%sur1
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(1)<r(1) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) + su_ww(1,1)
                   env%y(env%n) = ccub%hld(i)%pos(2) !+ su_ww(2,1)
                end if
             end do
      ccub => ccub%sur2
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(1)<r(1) .and. &
                     ccub%hld(i)%pos(2)<r(2) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) + su_ww(1,6)
                   env%y(env%n) = ccub%hld(i)%pos(2) + su_ww(2,6)
                end if
             end do
      ccub => ccub%sur3
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(2)<r(2) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) !+ su_ww(1,2)
                   env%y(env%n) = ccub%hld(i)%pos(2) + su_ww(2,2)
                end if
             end do
      ccub => ccub%sur3
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(1)>r(1) .and. &
                     ccub%hld(i)%pos(2)<r(2) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) + su_ww(1,7)
                   env%y(env%n) = ccub%hld(i)%pos(2) + su_ww(2,7)
                end if
             end do
      ccub => ccub%sur4
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(1)>r(1) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) + su_ww(1,3)
                   env%y(env%n) = ccub%hld(i)%pos(2) !+ su_ww(2,3)
                end if
             end do
      ccub => ccub%sur4
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(1)>r(1) .and. &
                     ccub%hld(i)%pos(2)>r(2) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) + su_ww(1,8)
                   env%y(env%n) = ccub%hld(i)%pos(2) + su_ww(2,8)
                end if
             end do
      ccub => ccub%sur1
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(2)>r(2) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) !+ su_ww(1,4)
                   env%y(env%n) = ccub%hld(i)%pos(2) + su_ww(2,4)
                end if
             end do
      ccub => ccub%sur1
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(1)<r(1) .and. &
                     ccub%hld(i)%pos(2)>r(2) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) + su_ww(1,5)
                   env%y(env%n) = ccub%hld(i)%pos(2) + su_ww(2,5)
                end if
             end do
      env%x(1:env%n) = env%x(1:env%n)-r(1)
      env%y(1:env%n) = env%y(1:env%n)-r(2)
      end subroutine scube_pp2d

      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      subroutine scube_up_pp2d(pos,idx,env)
      implicit none
      class(pp2d), intent(in), &
                       target   :: pos
      integer,       intent(in) :: idx
      type(stack), intent(out)  :: env
      integer                   :: i, cl, k
      real(pr)                  :: r(2)
      type(cube), pointer       :: ccub
      cl = pos%cidx(idx)
      k = pos%kidx(idx)
      ccub => pos%cup(cl)
      r = ccub%hld(k)%pos
      env%n = 0
      do i = 1, ccub%occ
         if( i/=k .and. ccub%hld(i)%pos(2)>r(2) ) then
                 env%n = env%n + 1
                 env%id(env%n) = ccub%idx(i)
                 env%x(env%n) = ccub%hld(i)%pos(1)
                 env%y(env%n) = ccub%hld(i)%pos(2)
         end if
      end do
      ccub => ccub%sur1
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(1)<r(1) .and. &
                     ccub%hld(i)%pos(2)>r(2) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) + su_ww(1,1)
                   env%y(env%n) = ccub%hld(i)%pos(2) !+ su_ww(2,1)
                end if
             end do
      ccub => ccub%sur2
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(1)<r(1) .and. &
                     ccub%hld(i)%pos(2)<r(2) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) + su_ww(1,6)
                   env%y(env%n) = ccub%hld(i)%pos(2) + su_ww(2,6)
                end if
             end do
      ccub => ccub%sur3
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(2)<r(2) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) !+ su_ww(1,2)
                   env%y(env%n) = ccub%hld(i)%pos(2) + su_ww(2,2)
                end if
             end do
      ccub => ccub%sur3
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(1)>r(1) .and. &
                     ccub%hld(i)%pos(2)<r(2) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) + su_ww(1,7)
                   env%y(env%n) = ccub%hld(i)%pos(2) + su_ww(2,7)
                end if
             end do
      ccub => ccub%sur4
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(1)>r(1) .and. &
                     ccub%hld(i)%pos(2)>r(2) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) + su_ww(1,3)
                   env%y(env%n) = ccub%hld(i)%pos(2) !+ su_ww(2,3)
                end if
             end do
      env%x(1:env%n) = env%x(1:env%n)-r(1)
      env%y(1:env%n) = env%y(1:env%n)-r(2)
      end subroutine scube_up_pp2d

      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      subroutine scube_right_pp2d(pos,idx,env)
      implicit none
      class(pp2d), intent(in), &
                       target   :: pos
      integer,       intent(in) :: idx
      type(stack), intent(out)  :: env
      integer                   :: i, cl, k
      real(pr)                  :: r(2)
      type(cube), pointer       :: ccub
      cl = pos%cidx(idx)
      k = pos%kidx(idx)
      ccub => pos%cup(cl)
      r = ccub%hld(k)%pos
      env%n = 0
      do i = 1, ccub%occ
         if(i/=k .and. ccub%hld(i)%pos(1)>r(1) ) then
                 env%n = env%n + 1
                 env%id(env%n) = ccub%idx(i)
                 env%x(env%n) = ccub%hld(i)%pos(1)
                 env%y(env%n) = ccub%hld(i)%pos(2)
         end if
      end do
      ccub => ccub%sur4
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(2)>r(2) .and. &
                     ccub%hld(i)%pos(1)>r(1) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) !+ su_ww(1,4)
                   env%y(env%n) = ccub%hld(i)%pos(2) + su_ww(2,4)
                end if
             end do
      ccub => ccub%sur1
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(1)<r(1) .and. &
                     ccub%hld(i)%pos(2)>r(2) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) + su_ww(1,5)
                   env%y(env%n) = ccub%hld(i)%pos(2) + su_ww(2,5)
                end if
             end do

      ccub => ccub%sur2
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(1)<r(1) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) + su_ww(1,1)
                   env%y(env%n) = ccub%hld(i)%pos(2) !+ su_ww(2,1)
                end if
             end do
      ccub => ccub%sur2
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(1)<r(1) .and. &
                     ccub%hld(i)%pos(2)<r(2) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) + su_ww(1,6)
                   env%y(env%n) = ccub%hld(i)%pos(2) + su_ww(2,6)
                end if
             end do
      ccub => ccub%sur3
             do i = 1, ccub%occ
                if( ccub%hld(i)%pos(2)<r(2) .and. &
                     ccub%hld(i)%pos(1)>r(1) ) then
                   env%n = env%n + 1
                   env%id(env%n) = ccub%idx(i)
                   env%x(env%n) = ccub%hld(i)%pos(1) !+ su_ww(1,2)
                   env%y(env%n) = ccub%hld(i)%pos(2) + su_ww(2,2)
                end if
             end do
      env%x(1:env%n) = env%x(1:env%n)-r(1)
      env%y(1:env%n) = env%y(1:env%n)-r(2)
      end subroutine scube_right_pp2d

      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      subroutine move_pp2d(pos,idx,dir,delta)
      implicit none
      class(pp2d), intent(inout), &
                          target :: pos
      integer,       intent(in)  :: idx, dir
      real(pr),      intent(in)  :: delta
      integer                    :: kk, side_eff, cl, k
      type(cube), pointer        :: cub
      cl = pos%cidx(idx)
      k = pos%kidx(idx)
      cub => pos%cup(cl)
      call cub%hld(k)%move(dir,delta)
      if( cub%hld(k)%pos(dir)>=pos%ww(dir) ) then
            call cub%hld(k)%move(dir,-pos%ww(dir))
            if( dir==1 ) then
               call cub%sur1%add(idx,cub%hld(k))
               kk = cub%sur1%occ
            else
               call cub%sur2%add(idx,cub%hld(k))
               kk = cub%sur2%occ
            end if
            call cub%pop( k, side_eff )
            pos%kidx( side_eff ) = k
            pos%cidx( idx ) = cub%su(dir)
            pos%kidx( idx ) = kk
      else  if( cub%hld(k)%pos(dir)<0.0_pr ) then
            call cub%hld(k)%move(dir,pos%ww(dir))
            if( dir==1 ) then
               call cub%sur3%add(idx,cub%hld(k))
               kk = cub%sur3%occ
            else
               call cub%sur4%add(idx,cub%hld(k))
               kk = cub%sur4%occ
            end if
            call cub%pop( k, side_eff )
            pos%kidx( side_eff ) = k
            pos%cidx( idx ) = cub%su(dir+2)
            pos%kidx( idx ) = kk
      end if
      end subroutine move_pp2d

      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      subroutine metropolis_pp2d(pos,efunc,eargs,tem,dmax,steps,nsuccess)
      implicit none
      class(pp2d), intent(inout) :: pos
      interface
         function efunc(x,y,args) result(en)
         import pr
         implicit none
         real(pr), intent(in) :: x(:), y(:), args(:)
         real(pr)             :: en(size(x))
         end function
      end interface
      real(pr),    intent(in)    :: eargs(:), tem, dmax
      integer,     intent(in)    :: steps
      integer,     intent(out)   :: nsuccess
      integer                    :: step, idx, dir, n
      real(pr)                   :: delta, de, rnd
      type(stack)                :: env
      nsuccess = 0
      do step = 1, steps
         call pos%random(-dmax,dmax,idx,dir,delta)
         call pos%zoom_on(idx,env)
         n = env%n
         env%f(1:n) = efunc(env%x(1:n),env%y(1:n),eargs)
         if( dir==1 ) then
            env%g(1:n) = efunc(env%x(1:n)-delta,env%y(1:n),eargs)
         else
            env%g(1:n) = efunc(env%x(1:n),env%y(1:n)-delta,eargs)
         end if
         de = sum(env%g(1:n)-env%f(1:n))
         if( de<=0.0_pr ) then
            call pos%move(idx,dir,delta) ; nsuccess = nsuccess + 1
         else
            call random_number(rnd)
            if(rnd<=exp(-de/tem)) then
               call pos%move(idx,dir,delta) ; nsuccess = nsuccess + 1
            end if
         end if
      end do
      end subroutine metropolis_pp2d

      !//////////////////////////////////////
      !//////////////////////////////////////           stack
      !//////////////////////////////////////

      pure &
      function x_stack(env) result(x)
      implicit none
      class(stack), intent(in) :: env
      real(pr)                 :: x(env%n)
      x = env%x(1:env%n)
      end function x_stack

      pure &
      function y_stack(env) result(y)
      implicit none
      class(stack), intent(in) :: env
      real(pr)                 :: y(env%n)
      y = env%x(1:env%n)
      end function y_stack

      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      pure &
      function hexorder(env) result(psi)
      implicit none
      class(stack), intent(in) :: env
      integer                  :: i, n, l(1)
      real(pr)                 :: psi(2), r, x, y, theta, &
                                  d(max(1,env%n))
      logical                  :: mask(max(1,env%n))
      n = env%n
      psi = 0.0_pr
      if( n>=6 ) then
         d = env%x(1:n)**2 + env%y(1:n)**2
         mask = .true.
         do i = 1, 6
            l = minloc(d,mask)
            r = sqrt(d(l(1)))
            x = env%x(l(1))
            y = env%y(l(1))
            theta = acos(x/r)
            if( y<0.0_pr ) theta = 2*pi - theta
            psi = psi + [cos(6*theta),sin(6*theta)]
            mask(l(1)) = .false.
         end do
         psi = psi/6
      end if
      end function hexorder

      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      subroutine rotate_stack(env,theta,matrix)
      implicit none
      class(stack), intent(inout)    :: env
      real(pr), intent(in), optional :: theta, matrix(2,2)
      real(pr)                       :: mat(2,2), r(2)
      integer                        :: i
      if( present(matrix) ) then
         mat = matrix
      elseif( present(theta) ) then
         mat = reshape( [cos(theta), sin(theta), -sin(theta), cos(theta)], [2,2] )
      else
         write(output_unit,*)"rotate stack"; stop
      end if
      do i = 1, env%n
         r = matmul(mat,[env%x(i),env%y(i)])
         env%x(i) = r(1)
         env%y(i) = r(2)
      end do
      end subroutine rotate_stack

      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      subroutine show_stack(env)
      implicit none
      class(stack), intent(in) :: env
      integer                  :: i
      do i = 1, env%n
         write(*,*) env%id(i), env%x(i), env%y(i)
      end do
      end subroutine show_stack

      !//////////////////////////////////////
      !//////////////////////////////////////           other
      !//////////////////////////////////////

      function latt_const(rho) result(a)
      implicit none
      real(pr),   intent(in)  :: rho
      real(pr)                :: a
      a = sqrt(2.0_pr/(rho*sqrt(3.0_pr)))
      end function latt_const

   end module universe


text:
! =========================================================
subroutine tfluct1(maxnx,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,amdq2)
! =========================================================
!
!	dummy tfluct function (see example euler_1d for more information)

    implicit double precision (a-h,o-z)

    write(*,*) "Error: solver.tfluct_solver=True, but you have not defined"
    write(*,*) "a tfluct solver, have you imported your solver and set solver.tfluct?"
    write(*,*) "see example/euler_1d for more information"

    stop

    return
end subroutine tfluct1

! =========================================================
subroutine tfluct2(ixyz,maxnx,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,amdq2)
! =========================================================
!
!	dummy tfluct function (see example euler_1d for more information)

    implicit double precision (a-h,o-z)

    write(*,*) "Error: solver.tfluct_solver=True, but you have not defined"
    write(*,*) "a tfluct solver, have you imported your solver and set solver.tfluct?"
    write(*,*) "see example/euler_1d for more information"

    stop

    return
end subroutine tfluct2

! =========================================================
subroutine tfluct3(ixyz,maxnx,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,amdq2)
! =========================================================
!
!	dummy tfluct function (see example euler_1d for more information)

    implicit double precision (a-h,o-z)

    write(*,*) "Error: solver.tfluct_solver=True, but you have not defined"
    write(*,*) "a tfluct solver, have you imported your solver and set solver.tfluct?"
    write(*,*) "see example/euler_1d for more information"

    stop

    return
end subroutine tfluct3

max_stars_repo_name: gyanz/pyclaw

max_issues_repo_name: gyanz/pyclaw

max_forks_repo_name: gyanz/pyclaw