#include "cfx5ext.h"
dllexport(userpart)
      SUBROUTINE USERPART(CZ,DZ,IZ,LZ,RZ)
C
C ---------------------------
C    Preprocessor includes
C ---------------------------
C
#include "stack_point.h"
#include "MMS.h"
#include "common_sizes.h"
C
C ------------------------------
C        Argument list
C ------------------------------
C
      INTEGER IZ(*)
      CHARACTER CZ(*)*(1)
      DOUBLE PRECISION DZ(*)
      LOGICAL LZ(*)
      REAL RZ(*)
C
C ------------------------------
C        External routines
C ------------------------------
C
      LOGICAL          LFIL
      EXTERNAL         LFIL
C
C ------------------------------
C        Local variables
C ------------------------------
C
      INTEGER          NPART, NVx
C
      LOGICAL          LPartWgt, LVxWgt
C
      CHARACTER*(4)    CRESLT, CDTYPE
      CHARACTER*(MXLEN_XALIAS) CZONE_ALIAS
C
C ---------------------------
C        Stack pointers
C ---------------------------
C
      __stack_point__    pKPartVx, pCrdVx, pPartWgt, pVxWgt,
     &                   pWORK_INTR, pWORK_REAL
C
C ------------------------------
C      Executable Statements
C ------------------------------
C
C---- Switch to user partitioning directory
C
      CALL PSHDIR('/USER_PARTITIONING','STOP',CRESLT)
C
C-----------------------------------------------------------------------
C     Get information from MMS
C-----------------------------------------------------------------------
C
        CALL PEEKCS('CZONE_ALIAS',CZONE_ALIAS,'STOP',CRESLT,CZ)
        CALL PEEKI('NPART',1,NPART,'STOP',CRESLT,IZ)
        CALL LOCDAT('KPartVx',CDTYPE,'STOP',NVx,pKPartVx,CRESLT)
        CALL LOCDAT('CrdVx',CDTYPE,'STOP',NSIZE,pCrdVx,CRESLT)
C
        IF (LFIL('PART_WEIGHTS')) THEN
          CALL LOCDAT('PART_WEIGHTS',
     &                CDTYPE,'STOP',NSIZE,pPartWgt,CRESLT)
          LPartWgt = .TRUE.
        ELSE
          pPartWgt = __dummy_stack_point_rz__
          LPartWgt = .FALSE.
        ENDIF
C
        IF (LFIL('VxWgt')) THEN
          CALL LOCDAT('VxWgt',
     &                CDTYPE,'STOP',NSIZE,pVxWgt,CRESLT)
          LVxWgt = .TRUE.
        ELSE
          pVxWgt = __dummy_stack_point_iz__
          LVxWgt = .FALSE.
        ENDIF
C
C-----------------------------------------------------------------------
C     Create working arrays
C-----------------------------------------------------------------------
C
        CALL MAKDAT('WORK_INTR','INTR','STOP',NVx,pWORK_INTR,CRESLT)
        CALL MAKDAT('WORK_REAL','REAL','STOP',NVx,pWORK_REAL,CRESLT)
C
C-----------------------------------------------------------------------
C     Call low-level partitioning routine
C-----------------------------------------------------------------------
C
        CALL USERPART_SUB(CZONE_ALIAS,NPART,NVx,
     &                    IZ(pKPartVx),RZ(pCrdVx),
     &                    LPartWgt,RZ(pPartWgt),LVxWgt,IZ(pVxWgt),
     &                    IZ(pWORK_INTR),RZ(pWORK_REAL))
C
      CALL POPDIR('STOP',CRESLT)
C
      END

C=======================================================================
C=======================================================================
C=======================================================================
C=======================================================================
C=======================================================================

      SUBROUTINE USERPART_SUB(CZONE_ALIAS,NPART,NVx,
     &                        KPartVx,CrdVx,
     &                        LPartWeight,PartWeight,LVxWeight,VxWeight,
     &                        WORK_INTR,WORK_REAL)
CC
CD Low level partitioning routine
CC
CC --------------------
CC        Input
CC --------------------
CC
CI CZONE_ALIAS  : Domain alias name
CI NPART        : Number of partitions
CI NVx          : Number of vertices in domain (could be reduced due
CI                to periodic boundary conditions)
CI CrdVx        : Cartesian vertex coordinates
CI LPartWeight  : Logical flag indicating if partition weights exist
CI PartWeight   : Partition weighting factors (optional)
CI LVxWeight    : Logical flag indicating if vertex weights exist
CI VxWeight     : Vertex weighting factors (optional)
CI WORK_INTR    : INTEGER work array of length NVx
CI WORK_REAL    : REAL work array of length NVx
CC
CC --------------------
CC      Modified
CC --------------------
CC
CC --------------------
CC        Output
CC --------------------
CC
CO KPartVx       : Partition-vertex assignment
CC
CC --------------------
CC       Details
CC --------------------
CC
CC --------------------
CC Local variables & parameters
CC --------------------
CC
CC======================================================================
C
C ---------------------------
C    Preprocessor includes
C ---------------------------
C
#include "cfd_sysdep.h"
#include "cfd_constants.h"
C
C ------------------------------
C        Argument list
C ------------------------------
C
      INTEGER            NPART, NVX,
     &                   KPartVx(NVx), VxWeight(*), WORK_INTR(NVx)
C
      REAL               CrdVx(3,NVx), PartWeight(*), WORK_REAL(NVx)
C
      LOGICAL            LPartWeight, LVxWeight
C
      CHARACTER*(*)      CZONE_ALIAS
C
C ------------------------------
C        External routines
C ------------------------------
C
C ------------------------------
C        Local Parameters
C ------------------------------
C
      CHARACTER*(*)      ROUTIN
      PARAMETER         (ROUTIN = 'USERPART_SUB')
C
C ------------------------------
C        Local variables
C ------------------------------
C
      INTEGER            IPART, IVx, IVx_S, IVx_F
C
      REAL               PartWeightSum, PartWeightUni, XORIGIN, ALPHA
C
C ------------------------------
C    Executable statements
C ------------------------------
C
C---- Set partition weight for uniform partitioning and initialize
C     local sum
C
      PartWeightUni = ONE/FLOAT(NPART)
      PartWeightSum = ZERO
C
C=======================================================================
C     Calculate angle for each vertex and sort vertices with respect to
C     this angle
C     --> distinguish between two regions:
C         Region 1:   y < 0   -->   0   < ALPHA < PI
C         Region 2:   y > 0   -->   PI  < ALPHA < 2*PI
C=======================================================================
C
      DO IVx=1,NVx
        IF (CrdVx(2,IVx).LE.0.0 .AND. CrdVx(1,IVx).LT.5.0) THEN
          XORIGIN = 2.0
          ALPHA = ATAN2(-CrdVx(2,IVx),-(CrdVx(1,IVx)-XORIGIN))
        ELSE 
          XORIGIN = 5.0
          ALPHA = TWO*PI - ATAN2(CrdVx(2,IVx),(CrdVx(1,IVx)-XORIGIN))
        ENDIF
        WORK_REAL(IVx) = ALPHA
      ENDDO
C
      CALL INDEX_SORT_REAL(NVx,WORK_REAL,WORK_INTR)
C
C=======================================================================
C     Assign partition numbers ro vertices
C=======================================================================
C
      IVx_F = 0
C
      DO IPART=1,NPART
C
C-----------------------------------------------------------------------
C     Increase sum of partition weights
C-----------------------------------------------------------------------
C
C---- Weighted partitioning
C
        IF (LPartWeight) THEN
C
          PartWeightSum = PartWeightSum + PartWeight(IPART)
C
C---- Uniform partitioning
C
        ELSE
C
          PartWeightSum = PartWeightSum + PartWeightUni
C
        ENDIF
C
C-----------------------------------------------------------------------
C     Determine vertex range for this partition
C-----------------------------------------------------------------------
C
        IVx_S = IVx_F + 1
        IVx_F = NINT(PartWeightSum*FLOAT(NVx))
        IVx_F = MIN(IVx_F,NVx)
        IVx_F = MAX(IVx_F,IVx_S)
C
C-----------------------------------------------------------------------
C     Set partition number for index rabge
C-----------------------------------------------------------------------
C
        DO IVx=IVx_S,IVx_F
          KPartVx(WORK_INTR(IVx)) = IPART
        ENDDO
C
      ENDDO
C
      END
