#include "cfx5ext.h"
dllexport(pt_tpd1)
      SUBROUTINE PT_TPD1(CZ,DZ,IZ,LZ,RZ)
C
C=======================================================================
C     Calculate the total particle mass within a sphere, defined
C     by a centre (SPHERE_CENTRE) and radius (SPHERE_RADIUS)
C=======================================================================
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
C ------------------------------
C        Local variables
C ------------------------------
C
      INTEGER            NPART, NSIZE
C
      REAL               TOTAL_MASS(3), SPHERE_CENTRE(3), 
     &                   SPHERE_RADIUS(3)
C
      CHARACTER*(MXDNAM) CVAR, COPER
      CHARACTER*(4)      CRESLT, CDTYPE
      CHARACTER*(MXLEN_XALIAS) CVAR_USER
C
C ---------------------------
C        Stack pointers
C ---------------------------
C
      __stack_point__    pCRD, pMASST, pRATE
C
C ------------------------------
C      Executable Statements
C ------------------------------
C
C-----------------------------------------------------------------------
C     Sphere definition
C-----------------------------------------------------------------------
C
C---- Centre of spheres
C
      SPHERE_CENTRE(1) = 0.5
      SPHERE_CENTRE(2) = 0.5
      SPHERE_CENTRE(3) = 0.5
C
C---- Sphere radii
C
      SPHERE_RADIUS(1) = 0.25
      SPHERE_RADIUS(2) = 0.35
      SPHERE_RADIUS(3) = 0.45
C
C-----------------------------------------------------------------------
C     Get number of particles and particle type from the MMS
C-----------------------------------------------------------------------
C
      CALL PEEKI('NPART',1,NPART,'STOP',CRESLT,IZ)
C
C-----------------------------------------------------------------------
C     To compute the mass of particles inside the spheres we need to 
C     provide the following data:
C       1. Particle position    (to decide whether particle is inside
C                                or outside of sphere)
C       2. Particle mass        (this is the mass of the representative
C                                particle)
C       3. Particle number rate (this is the number of 'real' particles
C                                that are represented by the 
C                                representative particle)
C
C     The following steps are required to get this data:
C       1. Convert user names into solver names
C       2. Lookup pointers for particle coordinates, mass and number
C          rate
C-----------------------------------------------------------------------
C
      CVAR_USER = 'Particle Position'
      CALL CONVERT_USER_NAMES(' ',CVAR_USER,CVAR,COPER,CRESLT,
     &                        CZ,DZ,IZ,LZ,RZ)
      CALL LOCDAT(CVAR,CDTYPE,'STOP',NSIZE,pCRD  ,CRESLT)
C
      CVAR_USER = 'Total Particle Mass'
      CALL CONVERT_USER_NAMES(' ',CVAR_USER,CVAR,COPER,CRESLT,
     &                        CZ,DZ,IZ,LZ,RZ)
      CALL LOCDAT(CVAR,CDTYPE,'STOP',NSIZE,pMASST,CRESLT)
C
      CVAR_USER = 'Particle Number Rate'
      CALL CONVERT_USER_NAMES(' ',CVAR_USER,CVAR,COPER,CRESLT,
     &                        CZ,DZ,IZ,LZ,RZ)
      CALL LOCDAT(CVAR,CDTYPE,'STOP',NSIZE,pRATE ,CRESLT)
C
C-----------------------------------------------------------------------
C     Pass down data to CALC_TPD1 and calculate the total mass within 
C     the three spheres
C-----------------------------------------------------------------------
C
      CALL CALC_TPD1(TOTAL_MASS,SPHERE_CENTRE,SPHERE_RADIUS,
     &               NPART,RZ(pCRD),RZ(pMASST),RZ(pRATE))
C
C-----------------------------------------------------------------------
C     NOTE: Transient particle diagnostics does not return results
C           via a return variable but requires the data to be stored
C           in a data area called TPD_VALUE within the current MMS
C           directory 
C           The solver then picks up the data from there and allows to 
C           monitor the quantities
C-----------------------------------------------------------------------
C
      CALL POKER('TPD_VALUE',1,TOTAL_MASS(1),'STOP',CRESLT,RZ)
      CALL POKER('TPD_VALUE',2,TOTAL_MASS(2),'STOP',CRESLT,RZ)
      CALL POKER('TPD_VALUE',3,TOTAL_MASS(3),'STOP',CRESLT,RZ)
C
      END

      SUBROUTINE CALC_TPD1(TOTAL_MASS,SPHERE_CENTRE,
     &                     SPHERE_RADIUS,NPART,CRD,MASST,RATE)
C
C=======================================================================
C     Subroutine for TPD1 which does the real calculation
C=======================================================================
C
      INTEGER NPART
      REAL    TOTAL_MASS(3), SPHERE_CENTRE(3), SPHERE_RADIUS(3),
     &        CRD(3,NPART), MASST(NPART), RATE(NPART)
C
      INTEGER IPART, I
      REAL    RADIUS
C
C---- Initialize return data
C
      DO I=1,3
        TOTAL_MASS(I) = 0.0
      ENDDO
C
C---- Loop over all representative particles
C
      DO IPART=1,NPART
C
C---- Compute distance of representative particle from sphere centre
C
        RADIUS = SQRT( (CRD(1,IPART)-SPHERE_CENTRE(1))**2
     &               + (CRD(2,IPART)-SPHERE_CENTRE(2))**2
     &               + (CRD(3,IPART)-SPHERE_CENTRE(3))**2 )
C
C---- Loop over all three spheres and add particle mass, if inside
C     sphere
C
        DO I=1,3
          IF (RADIUS.LE.SPHERE_RADIUS(I)) THEN
            TOTAL_MASS(I) = TOTAL_MASS(I) + MASST(IPART)*RATE(IPART)
          ENDIF
        ENDDO
      ENDDO
C
      END
