#include "cfx5ext.h"
dllexport(pt_wall_splash)
      SUBROUTINE PT_WALL_SPLASH (NLOC,NRET,NARG,RET,ARG,CRESLT,
     &                           CZ,DZ,IZ,LZ,RZ)
CC
CD User routine: template for particle user routine
CC
CC --------------------
CC        Input
CC --------------------
CC
CC  NLOC   - number of entities
CC  NRET   - length of return stack
CC  NARG   - length of argument stack
CC  ARG    - argument values
CC
CC --------------------
CC      Modified
CC --------------------
CC
CC --------------------
CC        Output
CC --------------------
CC
CC  RET    - return values
CC
CC --------------------
CC       Details
CC --------------------
CC
CC======================================================================
C
C ------------------------------
C     Preprocessor includes
C ------------------------------
C
C ------------------------------
C        Argument list
C ------------------------------
C
      INTEGER NARG, NRET, NLOC
C
      CHARACTER*(4) CRESLT
C
      REAL ARG(NLOC,NARG), RET(NLOC,NRET)
C
      INTEGER IZ(*)
      CHARACTER CZ(*)*(1)
      DOUBLE PRECISION DZ(*)
      LOGICAL LZ(*)
      REAL RZ(*)
C
C=======================================================================
C
C ---------------------------
C    Executable Statements
C ---------------------------
C
C=======================================================================
C
C     Return variables:
C     -----------------
C
C     Child droplet breakup indicator                       : RET(1,1)
C     Child droplet diameter                                : RET(1,5)
C     Child droplet number rate                             : RET(1,9)
C     Child droplet normal coef of restitution              : RET(1,13)
C     Child droplet parallel coef of restitution            : RET(1,17)
C     Child droplet mode                                    : RET(1,21)
C
C     Argument variables 
C     -------------------
C
C     Particle impact angle                                : ARG(1,1)
C     Particle velocity                                    : ARG(1,2)
C     Particle diameter                                    : ARG(1,5)
C     Particle number rate                                 : ARG(1,6)
C     Fluid density                                        : ARG(1,7)
C
C     We know that NLOC is 1 for the particle user source routines!!!!
C=======================================================================
C
C-----------------------------------------------------------------------
C     Calculate the return variables
C-----------------------------------------------------------------------
C
      CALL SPLASH (RET(1,1),RET(1,5),RET(1,9),RET(1,13),RET(1,17),
     &             RET(1,21),ARG(1,1),ARG(1,2),ARG(1,5),ARG(1,6),
     &             ARG(1,7))
C
      END
C
      SUBROUTINE SPLASH (BRKUP_IND,CHILD_DIAM,CHILD_NRATE,
     &                   CHILD_COEF_N,CHILD_COEF_P,CHILD_MODE,
     &                   ANGLE,VEL_PT,DIAM_PT,NRATE_PT,DENS_FL)
C
C ---------------------------
C    Preprocessor includes
C ---------------------------
C
#include "cfd_constants.h"
C
C---- Do not change these defines, they have to be consistent with what
C     the solver uses
C
#define   __regular_particle__ 1
#define      __wall_particle__ 20
#define     __pt_uf_get_data__ 1
#define    __pt_uf_save_data__ 2
#define        __pt_uf_iseed__ 1
#define        __pt_uf_pmode__ 2
C
C ------------------------------
C        Argument list
C ------------------------------
C
      REAL   ANGLE, VEL_PT(3), DIAM_PT, NRATE_PT,
     &       BRKUP_IND(3), CHILD_DIAM(3), CHILD_NRATE(3),
     &       CHILD_COEF_N(3), CHILD_COEF_P(3), CHILD_MODE(3), DENS_FL
C
C ------------------------------
C        Local variables
C ------------------------------
C
      REAL    SIGMA, VISC_PT, R, SPEED_PT, ANGLE_DEG, ALFA, BETA,
     &        GAMMA, VEL_NORM, OH, RE, K, M1M0
      INTEGER ICHILD, ISEED, ACTION
C
C ------------------------------
C    Executable statements
C ------------------------------
C
C=======================================================================
C     Prologue
C=======================================================================
C
C---- Initialization
C
      ICHILD = 0
C
C---- Get particle seed info (needed for random numer generator)
C
      ACTION = __pt_uf_save_data__
      CALL USER_PARTICLE_INFO(ACTION,__pt_uf_iseed__,ISEED)
C
C---- Some material properties (water)
C     --> Surface tension coefficient
C     --> Particle viscosity
C
      SIGMA   = 0.0201
      VISC_PT = 0.001
C
C---- Random number (required for deviation angle)
C
      CALL GET_RANDOM(R,1,ISEED)
C
C=======================================================================
C     Computation Section
C=======================================================================
C
C---- Particle velocity magnitude and impact angle
C
      SPEED_PT  = SQRT(VEL_PT(1)**2 + VEL_PT(2)**2 + VEL_PT(3)**2)
      ANGLE_DEG = ANGLE*180./PI
C
C---- Ohnesorge and Reynolds number 
C     --> Re, based on wall normal impact velocity!
C
      VEL_NORM = SIN(ANGLE)*SPEED_PT
      OH = VISC_PT/SQRT(DENS_FL*DIAM_PT*SIGMA)
      RE = DENS_FL*DIAM_PT*VEL_NORM/VISC_PT
C
C-----------------------------------------------------------------------
C     Precompute reflection angle
C-----------------------------------------------------------------------
C
C---- Impact angle, ALFA, average reflection angle, BETA. Add random 
C     deviation angle, GAMMA, if ALFA < 15 deg.
C
      ALFA  = 90.   - ANGLE_DEG
      BETA  = 61.88 + 0.326*ALFA
      GAMMA = 17.6  - 0.18*ALFA
C
      IF (ALFA .LT. 15.) BETA = BETA + (R-0.5)*GAMMA
C
C---- Compute K-value
C     --> Clip K to avoid overflow during m1/m0 calculation
C         (221^9.2133 = 3.9773E+21)
C
      K = MIN(OH*RE**1.25,221.)
C
C-----------------------------------------------------------------------
C     K < 57 -> Particle deposited at wall
C-----------------------------------------------------------------------
C
      IF (K .LE. 57) THEN
C
        ICHILD = ICHILD + 1
        BRKUP_IND(ICHILD)    = 1.
        CHILD_DIAM(ICHILD)   = DIAM_PT
        CHILD_NRATE(ICHILD)  = NRATE_PT
        CHILD_COEF_N(ICHILD) = 0.0
        CHILD_COEF_P(ICHILD) = 0.0
        CHILD_MODE(ICHILD)   = __wall_particle__
C
C-----------------------------------------------------------------------
C     ... else partially reflected
C-----------------------------------------------------------------------
C
      ELSE
C
C---- Mass fraction across reflection
C
        M1M0 = 3.9869E-21*K**9.2133
C
C---- Reflected portion of particle
C
        ICHILD = ICHILD + 1
        BRKUP_IND(ICHILD)    = 1.
        CHILD_DIAM(ICHILD)   = DIAM_PT*M1M0
        CHILD_NRATE(ICHILD)  = NRATE_PT
        CHILD_COEF_N(ICHILD) = COS(BETA*PI/180.)*SPEED_PT
        CHILD_COEF_P(ICHILD) = SIN(BETA*PI/180.)*SPEED_PT
        CHILD_MODE(ICHILD)   = __regular_particle__
C
C---- Deposited portion of particle
C
        ICHILD = ICHILD + 1
        BRKUP_IND(ICHILD)    = 1.
        CHILD_DIAM(ICHILD)   = DIAM_PT*(1. - M1M0)
        CHILD_NRATE(ICHILD)  = NRATE_PT
        CHILD_COEF_N(ICHILD) = 0.0
        CHILD_COEF_P(ICHILD) = 0.0
        CHILD_MODE(ICHILD)   = __wall_particle__
C
      ENDIF
C
C=======================================================================
C     Epiloge
C=======================================================================
C
C---- Update ISEED in particle database with value as computed in 
C     this routine
C
      ACTION = __pt_uf_get_data__
      CALL USER_PARTICLE_INFO(ACTION,__pt_uf_iseed__,ISEED)
C
      END
