#include "cfx5ext.h"
C=======================================================================
C
C Purpose: User Fortran for importing/interpolating external
C          velocity profile.
C
C Junction Box Routines:
C          USER_INPUT -> USER_INPUT_TABLE
C
C User CEL Functions:
C          USER_UPROFILE -> USER_UPROFILE_CAL -> USER_TABINT
C                                             -> USER_PRINT_PROFILE
C
C Command File data areas used:
C          /USER/INPUT_FILE
C          /USER/USER_PRINTING
C
C User Input data areas used:
C          /USER_DATA/NCOL
C          /USER_DATA/NLIN
C          /USER_DATA/TABLE(NCOL,NLIN)
C
C=======================================================================
dllexport(user_input)
      SUBROUTINE USER_INPUT ( CZ,DZ,IZ,LZ,RZ )
CC
CC --------------------
CC        Input
CC --------------------
CC  none
CC
CC --------------------
CC      Modified
CC --------------------
CC
CC  Stacks possibly.
CC
CC --------------------
CC        Output
CC --------------------
CC  none
CC
CC --------------------
CC       Details
CC --------------------
CC
CC  Load a velocity profile from filename given in /USER/INPUT_FILE.
CC
CC  This routine is defined with 'Junction Box Location = User Input'.
CC  It reads data into the directory /USER_DATA.
CC
CC  In a parallel run, this routine is only called from the master
CC  process, but any /USER_DATA data will then be automatically 
CC  copied to the /USER_DATA directory on all slave processes.
CC======================================================================
C ------------------------------
C     Preprocessor includes
C ------------------------------
#include "MMS.h"
#include "stack_point.h"
C
C ------------------------------
C        Argument list
C ------------------------------
      INTEGER IZ(*)
      CHARACTER CZ(*)*(1)
      DOUBLE PRECISION DZ(*)
      LOGICAL LZ(*)
      REAL RZ(*)
C
C ------------------------------
C        Local Parameters
C ------------------------------
C Fortran input channel
      INTEGER NCHANNEL
      PARAMETER (NCHANNEL=59)
C
C ------------------------------
C        Local Variables
C ------------------------------
C Result flag
      CHARACTER*4 CRESLT
C
C System info
      CHARACTER*(MXPNAM) WHICH_CALL
C
C Command file info
      CHARACTER*(MXPNAM) INPUT_FILE
      CHARACTER*(MXPNAM) USER_PRINTING
C
      INTEGER NCOL, NLIN
C
C ------------------------------
C        Stack pointers
C ------------------------------
C
      __stack_point__ pOINT, pTABLE
C
C=======================================================================
C ---------------------------
C    Executable Statements
C ---------------------------
C
C Look up the system info for junction box calls: WHICH_CALL
      WHICH_CALL = 'Unknown'
      CALL PEEKCS( '/USER/WHICH_CALL', WHICH_CALL, 'SKIP', CRESLT, CZ )
C
C Look up command file info: INPUT_FILE, USER_PRINTING.
      INPUT_FILE = 'Unknown'
      CALL PEEKCS( '/USER/INPUT_FILE', INPUT_FILE, 'STOP', CRESLT, CZ )
C
      USER_PRINTING = 'No'
      CALL PEEKCS( '/USER/USER_PRINTING', USER_PRINTING, 'SKIP',
     &            CRESLT, CZ )
C
C Send any diagnostic messages via master process.
      IF ( USER_PRINTING(1:1) .NE. 'N' ) THEN
         CALL MESAGE( 'WRITE', 'Start USER_INPUT')
         CALL MESAGE( 'WRITE', 'WHICH_CALL='//WHICH_CALL )
         CALL MESAGE( 'WRITE', 'INPUT_FILE='//INPUT_FILE )
      END IF
C
C Ensure that directory /USER_DATA exists and make this the current directory.
C
      CALL PSHDIR( '/', 'STOP', CRESLT )
      CALL CHGDIR( 'USER_DATA', 'STOP', CRESLT )
C
C Read data from file.
C
      OPEN( NCHANNEL, FILE=INPUT_FILE )
C
C Create space for array sizes NCOL and NLIN
      CALL MAKDAT( 'NCOL', 'INTR', 'STOP', 1, pOINT, CRESLT )
      CALL MAKDAT( 'NLIN', 'INTR', 'STOP', 1, pOINT, CRESLT )
C
C Read data in NCOL and NLIN
      READ ( NCHANNEL, * ) NCOL, NLIN
      CALL POKEI( 'NCOL', 1, NCOL, 'STOP', CRESLT, IZ )
      CALL POKEI( 'NLIN', 1, NLIN, 'STOP', CRESLT, IZ )
C
C Create space for array TABLE(NCOL,NLIN)
      CALL MAKDAT( 'TABLE', 'REAL', 'STOP', NCOL*NLIN, pTABLE, CRESLT )
C
C Read data into TABLE(NLIN,NCOL) located in RZ at pointer pTABLE
      CALL USER_INPUT_TABLE( NCHANNEL, RZ(pTABLE), NCOL, NLIN )
C
      CLOSE(NCHANNEL)
C
C Return to current directory on entry.
C
      CALL POPDIR( 'STOP', CRESLT )
C
C Send any diagnostic messages via master process.
      IF ( USER_PRINTING(1:1) .NE. 'N' ) THEN
         CALL MESAGE( 'WRITE','End USER_INPUT' )
      END IF
C
      END
C
C=======================================================================
      SUBROUTINE USER_INPUT_TABLE( NCHANNEL, TABLE, NCOL, NLIN )
C     ----------------------------------------------------------
C Read data into TABLE(NLIN,NCOL)
C
C Input arguments
      INTEGER NCHANNEL, NCOL, NLIN
C
C Output arguments
      REAL TABLE(NCOL,NLIN)
C
      READ( NCHANNEL, *, END=900, ERR=900 ) TABLE
      RETURN
C
C Error Handling via master process.
 900  CONTINUE
      CALL MESAGE( 'WRITE', 'Unexpected error/end on user file' )
      CALL CFXSTP( 'USER_INPUT_TABLE' )
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCcc
C=======================================================================
dllexport(user_uprofile)
      SUBROUTINE USER_UPROFILE ( 
     & NLOC, NRET, NARG, RET, ARGS, CRESLT, CZ,DZ,IZ,LZ,RZ )
CC
CD User routine: interpolate U profile from imported data.
CC
CC --------------------
CC        Input
CC --------------------
CC  NLOC   - size of current locale
CC  NRET   - number of components in result
CC  NARG   - number of arguments in call
CC  ARGS() - (NLOC,NARG) argument values
CC
CC --------------------
CC      Modified
CC --------------------
CC  Stacks possibly.
CC
CC --------------------
CC        Output
CC --------------------
CC  RET()  - (NLOC,NRET) return values
CC  CRESLT - 'GOOD' for success
CC
CC --------------------
CC       Details
CC --------------------
CC  An array of profiles TABLE(NCOL,NLIN) has been read from file
CC  by junction box routine USER_INPUT. This routine now uses
CC  the table to return U values for each y value given in the
CC  first argument.
CC
CC======================================================================
C ------------------------------
C     Preprocessor includes
C ------------------------------
#include "MMS.h"
#include "stack_point.h"
C
C ------------------------------
C        Argument list
C ------------------------------
      INTEGER NLOC,NARG,NRET
C
      CHARACTER*4 CRESLT
C
      REAL ARGS(NLOC,NARG), RET(NLOC,NRET)
C
      INTEGER IZ(*)
      CHARACTER CZ(*)*(1)
      DOUBLE PRECISION DZ(*)
      LOGICAL LZ(*)
      REAL RZ(*)
C
C ------------------------------
C        Local Variables
C ------------------------------
C Table info
      INTEGER NCOL, NLIN
C
C Command file info
      CHARACTER*(MXPNAM) USER_PRINTING
C
C Arguments for LOCDAT call
      CHARACTER*4 CDTYPE
      INTEGER ISIZE
C
C ------------------------------
C        Stack pointers
C ------------------------------
      __stack_point__ pTABLE
C
C=======================================================================
C ---------------------------
C    Executable Statements
C ---------------------------
C Look up command file info: USER_PRINTING.
      USER_PRINTING = 'No'
      CALL PEEKCS( '/USER/USER_PRINTING', USER_PRINTING, 'SKIP',
     & CRESLT, CZ )
C
C Send any diagnostic messages via master process.
      IF ( USER_PRINTING(1:1) .NE. 'N' ) THEN
         CALL MESAGE( 'WRITE','Start USER_UPROFILE' )
      END IF
C
C Initialise RET() to zero.
      CALL SET_A_0( RET, NLOC*NRET )
C
C Find TABLE(NCOL,NLIN)
      CALL PEEKI( '/USER_DATA/NCOL', 1, NCOL, 'STOP', CRESLT, IZ )
      CALL PEEKI( '/USER_DATA/NLIN', 1, NLIN, 'STOP', CRESLT, IZ )
      CALL LOCDAT( '/USER_DATA/TABLE', CDTYPE, 'STOP', ISIZE,
     &            pTABLE, CRESLT )
C
C Compute profile U(N) for Y(N) using TABLE(NCOL,NLIN)
      CALL USER_UPROFILE_CAL(
     & RET, ARGS, NLOC, RZ(pTABLE), NCOL, NLIN, USER_PRINTING )
C
C Send any diagnostic messages via master process.
      IF ( USER_PRINTING(1:1) .NE. 'N' ) THEN
         CALL MESAGE( 'WRITE','End USER_UPROFILE' )
      END IF
C
C Set success flag.
      CRESLT = 'GOOD'
C
      END
C=======================================================================
      SUBROUTINE USER_UPROFILE_CAL(
     & U, Y, N, TABLE, NCOL, NLIN, USER_PRINTING )
C     ---------------------------------------------
C Compute profile U(N) for Y(N) using TABLE(NCOL,NLIN)
C-----------------------
C       Input arguments
C-----------------------
C number of points
      INTEGER N
C
C y coordinate
      REAL Y(N)
C
C input table
C Y values are in the first column of the table.
C U values are in the second column of the table.
      INTEGER NCOL, NLIN
      REAL TABLE(NCOL,NLIN)
C
C printing control
      CHARACTER*(*) USER_PRINTING
C--------------------------
C       Output arguments
C--------------------------
      REAL U(N)
C----------------------
C        Locals
C----------------------
      INTEGER I
C
      DO I = 1, N
C
C interpolate U(I) for Y(I)
      CALL USER_TABINT( Y(I), U(I), 1, 2, TABLE, NCOL, NLIN )
C
C print I, U(I), Y(I)
C note that this data is unstructured, not ordered by Y value.
C
      IF ( USER_PRINTING(1:1) .NE. 'N' ) THEN
         CALL USER_PRINT_PROFILE( I, Y(I), U(I) )
      END IF
C
      END DO
C
      END
C=======================================================================
      SUBROUTINE USER_TABINT( X, Y, IX, IY, TABLE, NCOL, NLIN )
C     ---------------------------------------------------------
C Linear interpolation of table data Y in column IY of TABLE(),
C against an increasing co-ordinate X in column IX of TABLE().
C---------------
C        Inputs:
C---------------
      REAL X
      INTEGER IX, IY
      INTEGER NCOL, NLIN
      REAL TABLE(NCOL,NLIN)
C---------------
C       Outputs:
C---------------
      REAL Y
C---------------
C       Locals:
C---------------
      INTEGER I, J
      REAL S
C
C Find J such that TABLE(IX,J) <= X < TABLE(IX,J+1)
      J = 0
      DO I = 1, NLIN
         IF ( TABLE(IX,I) .GT. X ) GO TO 100
         J = I
      END DO
 100  CONTINUE
C
C Linearly interpolate Y
      IF( J .LE. 0 ) THEN
         Y = TABLE(IY,1)
      ELSE IF ( J .GE. NLIN ) THEN
         Y = TABLE(IY,NLIN)
      ELSE
         S = (X-TABLE(IX,J))/(TABLE(IX,J+1)-TABLE(IX,J))
         Y = TABLE(IY,J)*(1.0-S) + TABLE(IY,J+1)*S
      END IF
      END
C
C=======================================================================
      SUBROUTINE USER_PRINT_PROFILE( I, Y, U )
C     ----------------------------------------
C Print (I, Y, U) via master process.
C
      INTEGER I
      REAL Y, U
C
      CHARACTER CFROMI*5, CFROMR*15
      EXTERNAL CFROMI, CFROMR
C
      CALL MESAGE( 'BUFF-ASIS', CFROMI(I) )
      CALL MESAGE( 'BUFF-ASIS', ' Y(I) = ' )
      CALL MESAGE( 'BUFF-ASIS', CFROMR(Y) )
      CALL MESAGE( 'BUFF-ASIS', ' U(I) = ' )
      CALL MESAGE( 'BUFF-ASIS', CFROMR(U) )
      CALL MESAGE( 'BUFF-OUT', ' ' )
      END
