MODULE TRIVARIATEVAR
!     Global variables used in calculation of TRIVARIATE
!     normal and TRIVARIATE student T probabilties
      INTEGER :: NU
      DOUBLE PRECISION :: H1, H2, H3, R23, RUA, RUB, AR, RUC
      END MODULE TRIVARIATEVAR
!
! FIMOD contains functions for calculating 1D, 2D and 3D Normal and student T probabilites
!       and  1D expectations
      MODULE FIMOD
!      USE ERFCOREMOD
      IMPLICIT NONE
      PRIVATE
      PUBLIC :: NORMPRB, FI, FIINV, MVNLIMITS, MVNLMS, BVU,BVNMVN
      PUBLIC :: GAUSINT, GAUSINT2, EXLMS, EXINV
      PUBLIC :: STUDNT, BVTL, TVTL, TVNMVN

      INTERFACE NORMPRB
      MODULE PROCEDURE NORMPRB
      END INTERFACE

      INTERFACE FI
      MODULE PROCEDURE FI
      END INTERFACE

      INTERFACE FI2
      MODULE PROCEDURE FI2
      END INTERFACE

      INTERFACE FIINV
      MODULE PROCEDURE FIINV
      END INTERFACE

      INTERFACE  MVNLIMITS
      MODULE PROCEDURE  MVNLIMITS
      END INTERFACE

      INTERFACE  MVNLMS
      MODULE PROCEDURE  MVNLMS
      END INTERFACE

      INTERFACE BVU
      MODULE PROCEDURE BVU
      END INTERFACE

      INTERFACE BVNMVN
      MODULE PROCEDURE BVNMVN
      END INTERFACE

      INTERFACE STUDNT
      MODULE PROCEDURE STUDNT
      END INTERFACE

      INTERFACE BVTL
      MODULE PROCEDURE BVTL
      END INTERFACE

      INTERFACE TVTL
      MODULE PROCEDURE TVTL
      END INTERFACE

      INTERFACE GAUSINT
      MODULE PROCEDURE GAUSINT
      END INTERFACE

      INTERFACE GAUSINT2
      MODULE PROCEDURE GAUSINT2
      END INTERFACE

      INTERFACE EXLMS
      MODULE PROCEDURE EXLMS
      END INTERFACE

      INTERFACE EXINV
      MODULE PROCEDURE EXINV
      END INTERFACE

      CONTAINS
      FUNCTION FIINV(P) RESULT (VAL)
      IMPLICIT NONE
*
*     ALGORITHM AS241  APPL. STATIST. (1988) VOL. 37, NO. 3
*
*     Produces the normal deviate Z corresponding to a given lower
*     tail area of P.
*       Absolute error less than 1e-13
*       Relative error less than 1e-15 for abs(VAL)>0.1
*
*     The hash sums below are the sums of the mantissas of the
*     coefficients.   They are included for use in checking
*     transcription.
*
      DOUBLE PRECISION, INTENT(in) :: P
      DOUBLE PRECISION :: VAL
!local variables
      DOUBLE PRECISION SPLIT1, SPLIT2, CONST1, CONST2, ONE, ZERO, HALF,
     &     A0, A1, A2, A3, A4, A5, A6, A7, B1, B2, B3, B4, B5, B6, B7,
     &     C0, C1, C2, C3, C4, C5, C6, C7, D1, D2, D3, D4, D5, D6, D7,
     &     E0, E1, E2, E3, E4, E5, E6, E7, F1, F2, F3, F4, F5, F6, F7,
     &     Q, R
      PARAMETER ( SPLIT1 = 0.425D0, SPLIT2 = 5.D0,
     &     CONST1 = 0.180625D0, CONST2 = 1.6D0,
     &     ONE = 1.D0, ZERO = 0.D0, HALF = 0.5D0 )
*
*     Coefficients for P close to 0.5
*
      PARAMETER (
     *     A0 = 3.38713 28727 96366 6080D0,
     *     A1 = 1.33141 66789 17843 7745D+2,
     *     A2 = 1.97159 09503 06551 4427D+3,
     *     A3 = 1.37316 93765 50946 1125D+4,
     *     A4 = 4.59219 53931 54987 1457D+4,
     *     A5 = 6.72657 70927 00870 0853D+4,
     *     A6 = 3.34305 75583 58812 8105D+4,
     *     A7 = 2.50908 09287 30122 6727D+3,
     *     B1 = 4.23133 30701 60091 1252D+1,
     *     B2 = 6.87187 00749 20579 0830D+2,
     *     B3 = 5.39419 60214 24751 1077D+3,
     *     B4 = 2.12137 94301 58659 5867D+4,
     *     B5 = 3.93078 95800 09271 0610D+4,
     *     B6 = 2.87290 85735 72194 2674D+4,
     *     B7 = 5.22649 52788 52854 5610D+3 )
*     HASH SUM AB    55.88319 28806 14901 4439
*
*     Coefficients for P not close to 0, 0.5 or 1.
*
      PARAMETER (
     *     C0 = 1.42343 71107 49683 57734D0,
     *     C1 = 4.63033 78461 56545 29590D0,
     *     C2 = 5.76949 72214 60691 40550D0,
     *     C3 = 3.64784 83247 63204 60504D0,
     *     C4 = 1.27045 82524 52368 38258D0,
     *     C5 = 2.41780 72517 74506 11770D-1,
     *     C6 = 2.27238 44989 26918 45833D-2,
     *     C7 = 7.74545 01427 83414 07640D-4,
     *     D1 = 2.05319 16266 37758 82187D0,
     *     D2 = 1.67638 48301 83803 84940D0,
     *     D3 = 6.89767 33498 51000 04550D-1,
     *     D4 = 1.48103 97642 74800 74590D-1,
     *     D5 = 1.51986 66563 61645 71966D-2,
     *     D6 = 5.47593 80849 95344 94600D-4,
     *     D7 = 1.05075 00716 44416 84324D-9 )
*     HASH SUM CD    49.33206 50330 16102 89036
*
*	Coefficients for P near 0 or 1.
*
      PARAMETER (
     *     E0 = 6.65790 46435 01103 77720D0,
     *     E1 = 5.46378 49111 64114 36990D0,
     *     E2 = 1.78482 65399 17291 33580D0,
     *     E3 = 2.96560 57182 85048 91230D-1,
     *     E4 = 2.65321 89526 57612 30930D-2,
     *     E5 = 1.24266 09473 88078 43860D-3,
     *     E6 = 2.71155 55687 43487 57815D-5,
     *     E7 = 2.01033 43992 92288 13265D-7,
     *     F1 = 5.99832 20655 58879 37690D-1,
     *     F2 = 1.36929 88092 27358 05310D-1,
     *     F3 = 1.48753 61290 85061 48525D-2,
     *     F4 = 7.86869 13114 56132 59100D-4,
     *     F5 = 1.84631 83175 10054 68180D-5,
     *     F6 = 1.42151 17583 16445 88870D-7,
     *     F7 = 2.04426 31033 89939 78564D-15 )
*     HASH SUM EF    47.52583 31754 92896 71629
*
      Q = ( P - HALF)
      IF ( ABS(Q) .LE. SPLIT1 ) THEN ! Central range.
         R = CONST1 - Q*Q
         VAL = Q*( ( ( ((((A7*R + A6)*R + A5)*R + A4)*R + A3)
     *                  *R + A2 )*R + A1 )*R + A0 )
     *            /( ( ( ((((B7*R + B6)*R + B5)*R + B4)*R + B3)
     *                  *R + B2 )*R + B1 )*R + ONE)
      ELSE ! near the endpoints
         R = MIN( P, ONE - P )
         IF  (R .GT.ZERO) THEN ! ( 2.d0*R .GT. CFxCutOff) THEN ! R .GT.0.d0
            R = SQRT( -LOG(R) )
            IF ( R .LE. SPLIT2 ) THEN
               R = R - CONST2
               VAL = ( ( ( ((((C7*R + C6)*R + C5)*R + C4)*R + C3)
     *                      *R + C2 )*R + C1 )*R + C0 )
     *                /( ( ( ((((D7*R + D6)*R + D5)*R + D4)*R + D3)
     *                      *R + D2 )*R + D1 )*R + ONE )
            ELSE
               R = R - SPLIT2
               VAL = ( ( ( ((((E7*R + E6)*R + E5)*R + E4)*R + E3)
     *                      *R + E2 )*R + E1 )*R + E0 )
     *                /( ( ( ((((F7*R + F6)*R + F5)*R + F4)*R + F3)
     *                      *R + F2 )*R + F1 )*R + ONE )
            END IF
         ELSE
            VAL = 37.0d0 !9.D0 !XMAX 9.d0
         END IF
         IF ( Q  <  ZERO ) VAL = - VAL
      END IF
      RETURN
      END FUNCTION FIINV
                                ! *********************************
      SUBROUTINE NORMPRB(Z, P, Q)
!      USE ERFCOREMOD
!      USE GLOBALDATA, ONLY : XMAX
! Normal distribution probabilities accurate to 18 digits between
! -XMAX and XMAX
!
! Z    = no. of standard deviations from the mean.
! P, Q = probabilities to the left & right of Z.   P + Q = 1.
!
! by pab 23.03.2003
!
      IMPLICIT NONE
      DOUBLE PRECISION, INTENT(IN)            :: Z
      DOUBLE PRECISION, INTENT(OUT)           :: P
      DOUBLE PRECISION, INTENT(OUT), OPTIONAL ::  Q
!Local variables
      DOUBLE PRECISION            :: PP, ZABS
      DOUBLE PRECISION, PARAMETER :: ZERO  = 0.0D0
      DOUBLE PRECISION, PARAMETER :: ONE   = 1.0D0
      DOUBLE PRECISION, PARAMETER :: SQ2M1 = 0.70710678118655D0 !     1/SQRT(2)
      DOUBLE PRECISION, PARAMETER :: HALF  = 0.5D0
      DOUBLE PRECISION, PARAMETER :: XMAX  = 37D0
      ZABS = ABS(Z)
!
!     |Z| > 37  (or XMAX)
!
      IF ( ZABS .GT. XMAX ) THEN
         IF (Z > ZERO) THEN
            P = ONE
            IF (PRESENT(Q)) Q = ZERO
         ELSE
            P = ZERO
            IF (PRESENT(Q)) Q = ONE
         END IF
      ELSE
!
!     |Z| <= 37
!
         PP = DERFC(ZABS*SQ2M1)*HALF

         IF (Z  <  ZERO) THEN
            P = PP
            IF (PRESENT(Q)) Q = ONE - PP
         ELSE
            P = ONE - PP
            IF (PRESENT(Q)) Q = PP
         END IF
      END IF

      RETURN
      END SUBROUTINE NORMPRB
      FUNCTION FI( Z ) RESULT (VALUE)
!      USE ERFCOREMOD
!      USE GLOBALDATA, ONLY : XMAX
      IMPLICIT NONE
      DOUBLE PRECISION, INTENT(in) :: Z
      DOUBLE PRECISION :: VALUE
! Local variables
      DOUBLE PRECISION :: ZABS
      DOUBLE PRECISION, PARAMETER:: SQ2M1 = 0.70710678118655D0 !     1/SQRT(2)
      DOUBLE PRECISION, PARAMETER:: HALF = 0.5D0
      DOUBLE PRECISION, PARAMETER:: XMAX = 37.D0
      ZABS = ABS(Z)
*
*     |Z| > 37  (or XMAX)
*
      IF ( ZABS .GT. XMAX ) THEN
         IF (Z  <  0.0D0) THEN
            VALUE = 0.0D0
         ELSE
            VALUE = 1.0D0
         ENDIF
      ELSE
         VALUE = DERFC(-Z*SQ2M1)*HALF
      ENDIF
      RETURN
      END FUNCTION FI

      FUNCTION FI2( Z ) RESULT (VALUE)
!      USE GLOBALDATA, ONLY : XMAX
      IMPLICIT NONE
      DOUBLE PRECISION, INTENT(in) :: Z
      DOUBLE PRECISION :: VALUE
*
*     Normal distribution probabilities accurate to 1.e-15.
*     relative error less than 1e-8;
*     Z = no. of standard deviations from the mean.
*
*     Based upon algorithm 5666 for the error function, from:
*     Hart, J.F. et al, 'Computer Approximations', Wiley 1968
*
*     Programmer: Alan Miller
*
*     Latest revision - 30 March 1986
*
      DOUBLE PRECISION :: P0, P1, P2, P3, P4, P5, P6,
     *     Q0, Q1, Q2, Q3, Q4, Q5, Q6, Q7,XMAX,
     *     P, EXPNTL, CUTOFF, ROOTPI, ZABS, Z2
      PARAMETER(
     *     P0 = 220.20 68679 12376 1D0,
     *     P1 = 221.21 35961 69931 1D0,
     *     P2 = 112.07 92914 97870 9D0,
     *     P3 = 33.912 86607 83830 0D0,
     *     P4 = 6.3739 62203 53165 0D0,
     *     P5 = 0.70038 30644 43688 1D0,
     *     P6 = 0.035262 49659 98910 9D0 )
      PARAMETER(
     *     Q0 = 440.41 37358 24752 2D0,
     *     Q1 = 793.82 65125 19948 4D0,
     *     Q2 = 637.33 36333 78831 1D0,
     *     Q3 = 296.56 42487 79673 7D0,
     *     Q4 = 86.780 73220 29460 8D0,
     *     Q5 = 16.064 17757 92069 5D0,
     *     Q6 = 1.7556 67163 18264 2D0,
     *     Q7 = 0.088388 34764 83184 4D0 )
      PARAMETER( ROOTPI = 2.5066 28274 63100 1D0 )
      PARAMETER( CUTOFF = 7.0710 67811 86547 5D0 )
      PARAMETER( XMAX   = 37.D0 )
*
      ZABS = ABS(Z)
*
*     |Z| > 37  (or XMAX)
*
      IF ( ZABS .GT. XMAX ) THEN
         P = 0.d0
      ELSE
*
*     |Z| <= 37
*
         Z2 = ZABS * ZABS
         EXPNTL = EXP( -Z2 * 0.5D0 )
*
*     |Z| < CUTOFF = 10/SQRT(2)
*
         IF ( ZABS  <  CUTOFF ) THEN
            P = EXPNTL*( (((((P6*ZABS + P5)*ZABS + P4)*ZABS + P3)*ZABS
     *           + P2)*ZABS + P1)*ZABS + P0)/(((((((Q7*ZABS + Q6)*ZABS
     *           + Q5)*ZABS + Q4)*ZABS + Q3)*ZABS + Q2)*ZABS + Q1)*ZABS
     *           + Q0 )
*
*     |Z| >= CUTOFF.
*
         ELSE
            P = EXPNTL/( ZABS + 1.d0/( ZABS + 2.d0/( ZABS + 3.d0/( ZABS
     *                        + 4.d0/( ZABS + 0.65D0 ) ) ) ) )/ROOTPI
         END IF
      END IF
      IF ( Z .GT. 0.d0 ) P = 1.d0 - P
      VALUE = P
      RETURN
      END FUNCTION FI2

      SUBROUTINE MVNLIMITS( A, B, INFIN, AP, PRB, AQ)
! RETURN probabilities for being between A and B
!  WHERE
!  AP = FI(A), AQ = 1 - FI(A)
!  BP = FI(B), BQ = 1 - FI(B)
!  PRB = BP-AP IF BP+AP<1
!      = AQ-BQ OTHERWISE
!
      IMPLICIT NONE
      DOUBLE PRECISION, INTENT(in) :: A, B
      DOUBLE PRECISION, INTENT(out) :: AP
      DOUBLE PRECISION, INTENT(out),OPTIONAL :: PRB,AQ
      INTEGER,INTENT(in) :: INFIN
!     LOCAL VARIABLES
      DOUBLE PRECISION :: BP,AQQ, BQQ
      DOUBLE PRECISION, PARAMETER :: ONE=1.D0, ZERO = 0.D0

      SELECT CASE (infin)
      CASE (:-1)
         AP  = ZERO
!     BP  = ONE
         IF (PRESENT(PRB)) PRB = ONE
         IF (PRESENT(AQ)) AQ = ONE
!     IF (PRESENT(BQ)) BQ = ZERO
      CASE (0)
         AP  = ZERO
         CALL NORMPRB(B,BP)     !,BQQ)
         IF (PRESENT(PRB)) PRB = BP
         IF (PRESENT(AQ)) AQ = ONE
!     IF (PRESENT(BQ)) BQ = BQQ
      CASE (1)
!     BP = ONE
         CALL NORMPRB(A,AP,AQQ)
         IF (PRESENT(PRB)) PRB = AQQ
         IF (PRESENT(AQ)) AQ = AQQ
!     IF (PRESENT(BQ)) BQ = ZERO
      CASE (2:)
         CALL NORMPRB(A,AP,AQQ)
         CALL NORMPRB(B,BP,BQQ)
         IF (PRESENT(PRB)) THEN
            IF (AP+BP  <  ONE) THEN
               PRB = BP - AP
            ELSE
               PRB = AQQ - BQQ
            END IF
         ENDIF
         IF (PRESENT(AQ)) AQ = AQQ
!     IF (PRESENT(BQ)) BQ = BQQ
      END SELECT
      RETURN
      END SUBROUTINE MVNLIMITS


      SUBROUTINE MVNLMS( A, B, INFIN, LOWER, UPPER )
      IMPLICIT NONE
      DOUBLE PRECISION, INTENT(in) :: A, B
      DOUBLE PRECISION, INTENT(out) :: LOWER, UPPER
      INTEGER,INTENT(in) :: INFIN

      LOWER = 0.0D0
      UPPER = 1.0D0
      IF ( INFIN  <  0 ) RETURN
      IF ( INFIN .NE. 0 ) LOWER = FI(A)
      IF ( INFIN .NE. 1 ) UPPER = FI(B)
      RETURN
      END SUBROUTINE MVNLMS

      FUNCTION TVNMVN(A, B, INFIN, R, EPSI ) RESULT (VAL)
      IMPLICIT NONE
*
*     A function for computing trivariate normal probabilities.
*
*  Parameters
*
*     A  REAL, array of lower integration limits.
*     B  REAL, array of upper integration limits.
!     R  REAL, array of correlation coefficents
!         R = [r12 r13 r23]
!    EPSI = REAL tolerance
*     INFIN  INTEGER, array of integration limits flags:
*            if INFIN(I) = 0, Ith limits are (-infinity, B(I)];
*            if INFIN(I) = 1, Ith limits are [A(I), infinity);
*            if INFIN(I) = 2, Ith limits are [A(I), B(I)].
      DOUBLE PRECISION, DIMENSION(:), INTENT (IN) :: A, B , R
      DOUBLE PRECISION,               INTENT (IN) :: EPSI
      INTEGER,          DIMENSION(:), INTENT (IN) :: INFIN
      DOUBLE PRECISION :: VAL

      IF ( INFIN(1) .EQ. 2 ) THEN
         IF (  INFIN(2) .EQ. 2 ) THEN
            IF (INFIN(3) .EQ. 2 ) THEN !OK
               VAL =  TVNL( B(1), B(2), B(3),R(1),R(2),R(3),EPSI )
     &              - TVNL( A(1), B(2), B(3),R(1),R(2),R(3),EPSI )
     &              - TVNL( B(1), A(2), B(3),R(1),R(2),R(3),EPSI )
     &              - TVNL( B(1), B(2), A(3),R(1),R(2),R(3),EPSI )
     &              + TVNL( A(1), A(2), B(3),R(1),R(2),R(3),EPSI )
     &              + TVNL( A(1), B(2), A(3),R(1),R(2),R(3),EPSI )
     $              + TVNL( B(1), A(2), A(3),R(1),R(2),R(3),EPSI )
     &              - TVNL( A(1), A(2), A(3),R(1),R(2),R(3),EPSI )
            ELSE IF (INFIN(3) .EQ. 1 ) THEN ! B(3) = inf ok
               VAL = TVNL( B(1), B(2), -A(3),R(1),-R(2),-R(3),EPSI )
     &              - TVNL( A(1), B(2), -A(3),R(1),-R(2),-R(3),EPSI )
     &              - TVNL( B(1), A(2), -A(3),R(1),-R(2),-R(3),EPSI )
     &              + TVNL( A(1), A(2), -A(3),R(1),-R(2),-R(3),EPSI )
            ELSE IF (INFIN(3) .EQ. 0 ) THEN !OK A(3) = -inf
               VAL =  TVNL( B(1), B(2), B(3),R(1),R(2),R(3),EPSI )
     &              - TVNL( A(1), B(2), B(3),R(1),R(2),R(3),EPSI )
     &              - TVNL( B(1), A(2), B(3),R(1),R(2),R(3),EPSI )
     &              + TVNL( A(1), A(2), B(3),R(1),R(2),R(3),EPSI )
            ELSE ! INFIN(1:2)=2
               VAL = BVNMVN( A ,B, INFIN, R(1) )
            ENDIF
         ELSE IF (INFIN(2) .EQ. 1 ) THEN ! B(2) = inf
            IF (INFIN(3) .EQ. 2 ) THEN
               VAL = TVNL( B(1), -A(2), B(3),-R(1),R(2),-R(3),EPSI )
     &              - TVNL( A(1), -A(2), B(3),-R(1),R(2),-R(3),EPSI )
     &              - TVNL( B(1), -A(2), A(3),-R(1),R(2),-R(3),EPSI )
     &              + TVNL( A(1), -A(2), A(3),-R(1),R(2),-R(3),EPSI )
            ELSE IF (INFIN(3) .EQ. 1 ) THEN
               VAL = TVNL( B(1), -A(2), -A(3),-R(1),-R(2),R(3),EPSI )
     $              - TVNL( A(1), -A(2), -A(3),-R(1),-R(2),R(3),EPSI )
            ELSE IF (INFIN(3) .EQ. 0 ) THEN
               VAL = TVNL( B(1), -A(2), B(3),-R(1),R(2),-R(3),EPSI)
     $              - TVNL( A(1), -A(2), B(3),-R(1),R(2),-R(3),EPSI)
            ELSE
               VAL = BVNMVN( A ,B, INFIN, R(1) )
            ENDIF
         ELSE IF (INFIN(2) .EQ. 0 ) THEN
            SELECT CASE (INFIN(3))
            CASE (2:)           ! % % A(2)=-INF
               VAL = TVNL( B(1), B(2), B(3),R(1),R(2),R(3),EPSI )
     $              - TVNL( A(1), B(2), B(3),R(1),R(2),R(3),EPSI )
     $              - TVNL( B(1), B(2), A(3),R(1),R(2),R(3),EPSI )
     $              + TVNL( A(1), B(2), A(3),R(1),R(2),R(3),EPSI )
            CASE (1)            !%  % A(2)=-INF B(3) = INF
               VAL = TVNL( B(1), B(2), -A(3),R(1),-R(2),-R(3),EPSI)
     $              - TVNL(A(1), B(2), -A(3),R(1),-R(2),-R(3),EPSI)
            CASE (0)            ! % % A(2)=-INF A(3) = -INF
               VAL = TVNL( B(1), B(2), B(3),R(1),R(2),R(3),EPSI )
     $              - TVNL( A(1), B(2), B(3),R(1),R(2),R(3),EPSI)
            CASE DEFAULT
               VAL = BVNMVN(A,B,INFIN,R(1))
            END SELECT
         ELSE
            VAL = BVNMVN(A(1:3:2),B(1:3:2),INFIN(1:3:2),R(2))
         ENDIF
      ELSE IF ( INFIN(1) .EQ. 1 ) THEN
         SELECT CASE (INFIN(2))
         CASE (2)
            SELECT CASE (INFIN(3))
            CASE (2)            !% B(1) = INF   %OK
               VAL =  TVNL(-A(1), B(2), B(3),-R(1),-R(2),R(3),EPSI )
     $              - TVNL(-A(1), B(2), A(3),-R(1),-R(2),R(3),EPSI )
     $              - TVNL(-A(1), A(2), B(3),-R(1),-R(2),R(3),EPSI )
     $              + TVNL(-A(1), A(2), A(3),-R(1),-R(2),R(3),EPSI )
            CASE (1)            ! % B(1) = INF   B(3) = INF %OK
               VAL = TVNL(-A(1), B(2), -A(3),-R(1),R(2),-R(3),EPSI )
     $              - TVNL(-A(1), A(2), -A(3),-R(1),R(2),-R(3),EPSI)
            CASE (0)            ! % B(1) = INF   A(3) = -INF %OK
               VAL =  TVNL(-A(1), B(2), B(3),-R(1),-R(2),R(3),EPSI )
     $              - TVNL(-A(1), A(2), B(3),-R(1),-R(2),R(3),EPSI)
            CASE (-1)
               VAL = BVNMVN(A,B,INFIN,R(1))
            END SELECT
         CASE (1)               !%B(2) = INF
            SELECT CASE (INFIN(3))
            CASE (2)            ! % B(1) = INF  B(2) = INF % OK
               VAL = TVNL( -A(1), -A(2), B(3),-R(1),R(2),-R(3),EPSI )
     &              - TVNL( -A(1), -A(2),A(3),-R(1),R(2),-R(3),EPSI )
            CASE (1)            ! % B(1:3) = INF %OK
               VAL = TVNL( -A(1), -A(2), -A(3),R(1),R(2),R(3),EPSI)
            CASE (0)            !  % B(1:2) = INF A(3) = -INF %OK
               VAL = TVNL( -A(1), -A(2), B(3),R(1),-R(2),-R(3),EPSI )
            CASE (:-1)
               VAL = BVNMVN(A,B,INFIN,R(1))
            END SELECT
         CASE (0) ! A(2) = -INF
            SELECT CASE ( INFIN(3))
            CASE (2)            ! B(1) = INF , A(2) = -INF %OK
               VAL = TVNL( -A(1), B(2), B(3),-R(1),R(2),-R(3),EPSI )
     &              - TVNL( -A(1), B(2),A(3),-R(1),R(2),-R(3),EPSI )
            CASE (1)            ! B(1) = INF , A(2) = -INF  B(3) = INF % OK
               VAL = TVNL( -A(1), B(2), -A(3),-R(1),-R(2),R(3),EPSI)
            CASE (0)            !% B(1) = INF , A(2:3) = -INF
               VAL = TVNL( -A(1), B(2), B(3),-R(1),-R(2),R(3),EPSI )
            CASE (:-1)
               VAL = BVNMVN(A,B,INFIN,R(1))
            END SELECT
         CASE DEFAULT
            VAL = BVNMVN(A(1:3:2),B(1:3:2),INFIN(1:3:2),R(2))
         END SELECT
      ELSE IF ( INFIN(1) .EQ. 0 ) THEN
         SELECT CASE (INFIN(2))
         CASE (2)
            SELECT CASE (INFIN(3))
            CASE (2:)            ! A(1) = -INF %OK
               VAL =  TVNL( B(1), B(2), B(3),R(1),R(2),R(3),EPSI )
     &              - TVNL( B(1), B(2), A(3),R(1),R(2),R(3),EPSI)
     &              - TVNL( B(1), A(2), B(3),R(1),R(2),R(3),EPSI )
     &              + TVNL( B(1), A(2), A(3),R(1),R(2),R(3),EPSI )
            CASE (1) ! % A(1) = -INF , B(3) = INF %OK
               VAL =  TVNL( B(1), B(2), -A(3),R(1),-R(2),-R(3),EPSI )
     $              - TVNL( B(1), A(2), -A(3),R(1),-R(2),-R(3),EPSI )
            CASE (0)            ! A(1) = -INF , A(3) = -INF %OK
               VAL =  TVNL( B(1), B(2), B(3),R(1),R(2),R(3),EPSI )
     &              - TVNL( B(1), A(2), B(3),R(1),R(2),R(3),EPSI )
            CASE DEFAULT
               VAL = BVNMVN(A,B,INFIN,R(1))
            END SELECT
         CASE (1)               ! B(2) = INF
            SELECT CASE (INFIN(3))
            CASE (2:)            ! A(1) = -INF B(2) = INF %OK
               VAL =  TVNL( B(1), -A(2), B(3),-R(1),R(2),-R(3),EPSI)
     $              - TVNL( B(1), -A(2), A(3),-R(1),R(2),-R(3),EPSI)
            CASE (1)            ! A(1) = -INF B(2) = INF  B(3) = INF %OK
               VAL =  TVNL( B(1), -A(2), -A(3),-R(1),-R(2),R(3),EPSI)
            CASE (0)            ! % A(1) = -INF B(2) = INF  A(3) = -INF %OK
               VAL = TVNL(B(1), -A(2), B(3),-R(1),R(2),-R(3),EPSI)
            CASE DEFAULT
               VAL = BVNMVN(A,B,INFIN,R(1))
            END SELECT
         CASE (0)               ! A(2) = -INF
            SELECT CASE (INFIN(3))
            CASE (2:)            ! %  A(1:2) = -INF
               VAL =  TVNL( B(1), B(2), B(3),R(1),R(2),R(3),EPSI)
     $              - TVNL( B(1), B(2), A(3),R(1),R(2),R(3),EPSI)
            CASE (1)            ! A(1:2) = -INF B(3) = INF
               VAL =  TVNL( B(1), B(2), -A(3),R(1),-R(2),-R(3),EPSI)
            CASE (0)            !  % A(1:3) = -INF
               VAL =  TVNL( B(1), B(2), B(3),R(1),R(2),R(3),EPSI )
            CASE DEFAULT
               VAL = BVNMVN(A,B,INFIN,R(1))
            END  SELECT
         CASE DEFAULT
            VAL = BVNMVN(A(1:3:2),B(1:3:2),INFIN(1:3:2),R(2))
         END SELECT
      ELSE
         VAL = BVNMVN(A(2:3),B(2:3),INFIN(2:3),R(3))
      END IF
      CONTAINS
      DOUBLE PRECISION FUNCTION TVNL(H1,H2,H3, R12,R13,R23, EPSI )
      !Returns Trivariate Normal CDF
      DOUBLE PRECISION, INTENT(IN) :: R12,R13,R23
      DOUBLE PRECISION, INTENT(IN) :: H1,H2,H3, EPSI
!     Locals
      INTEGER, PARAMETER :: NU = 0
      DOUBLE PRECISION,DIMENSION(3) :: H,R
      H(:) = (/ H1, H2, H3 /)
      R(:) = (/ R12, R13, R23 /)
      TVNL = TVTL(NU,H,R,EPSI)
      END FUNCTION TVNL
      END  FUNCTION TVNMVN
      FUNCTION BVNMVN( LOWER, UPPER, INFIN, CORREL ) RESULT (VAL)
      IMPLICIT NONE
*
*     A function for computing bivariate normal probabilities.
*
*  Parameters
*
*     LOWER  REAL, array of lower integration limits.
*     UPPER  REAL, array of upper integration limits.
*     INFIN  INTEGER, array of integration limits flags:
*            if INFIN(I) = 0, Ith limits are (-infinity, UPPER(I)];
*            if INFIN(I) = 1, Ith limits are [LOWER(I), infinity);
*            if INFIN(I) = 2, Ith limits are [LOWER(I), UPPER(I)].
*     CORREL REAL, correlation coefficient.
*

      DOUBLE PRECISION, DIMENSION(:), INTENT (IN) :: LOWER, UPPER
      DOUBLE PRECISION,               INTENT (IN) :: CORREL
      INTEGER,          DIMENSION(:), INTENT (IN) :: INFIN
      DOUBLE PRECISION :: VAL
      DOUBLE PRECISION :: E
      SELECT CASE (INFIN(1))
      CASE (2:)
         SELECT CASE ( INFIN(2) )
         CASE (2:)
            VAL =  BVU ( LOWER(1), LOWER(2), CORREL )
     &           - BVU ( UPPER(1), LOWER(2), CORREL )
     &           - BVU ( LOWER(1), UPPER(2), CORREL )
     &           + BVU ( UPPER(1), UPPER(2), CORREL )

         CASE (1)
            VAL =  BVU ( LOWER(1), LOWER(2), CORREL )
     &           - BVU ( UPPER(1), LOWER(2), CORREL )
         CASE (0)
            VAL =  BVU ( -UPPER(1), -UPPER(2), CORREL )
     &           - BVU ( -LOWER(1), -UPPER(2), CORREL )
         CASE DEFAULT
             CALL MVNLIMITS(LOWER(1),UPPER(1),INFIN(1),E,VAL)
         END SELECT
      CASE (1)
         SELECT CASE ( INFIN(2))
         CASE ( 2: )
            VAL =  BVU ( LOWER(1), LOWER(2), CORREL )
     &           - BVU ( LOWER(1), UPPER(2), CORREL )
         CASE (1)
            VAL =  BVU ( LOWER(1), LOWER(2), CORREL )
         CASE (0)
            VAL =  BVU ( LOWER(1), -UPPER(2), -CORREL )
         CASE DEFAULT
            CALL MVNLIMITS(LOWER(2),UPPER(2),INFIN(2),E,VAL)
         END SELECT
      CASE (0)
         SELECT CASE ( INFIN(2))
         CASE ( 2: )
            VAL =  BVU ( -UPPER(1), -UPPER(2), CORREL )
     &           - BVU ( -UPPER(1), -LOWER(2), CORREL )
         CASE ( 1 )
            VAL =  BVU ( -UPPER(1), LOWER(2), -CORREL )
         CASE (0)
            VAL =  BVU ( -UPPER(1), -UPPER(2), CORREL )
         CASE DEFAULT
            CALL MVNLIMITS(LOWER(1),UPPER(1),INFIN(1),E,VAL)
         END SELECT
      CASE DEFAULT !ELSE  !INFIN(1)<0
         CALL MVNLIMITS(LOWER(2),UPPER(2),INFIN(2),E,VAL)
      END SELECT
      END  FUNCTION BVNMVN
      FUNCTION BVU( SH, SK, R ) RESULT (VAL)
!      USE GLOBALDATA, ONLY: XMAX
      IMPLICIT NONE
*
!     A function for computing bivariate normal probabilities.
!
!       Yihong Ge
!       Department of Computer Science and Electrical Engineering
!       Washington State University
!       Pullman, WA 99164-2752
!     and
!       Alan Genz
!       Department of Mathematics
!       Washington State University
!       Pullman, WA 99164-3113
!       Email : alangenz@wsu.edu
!
!    This function is based on the method described by
!        Drezner, Z and G.O. Wesolowsky, (1989),
!        On the computation of the bivariate normal integral,
!        Journal of Statist. Comput. Simul. 35, pp. 101-107,
!    with major modifications for double precision, and for |R| close to 1.
!
! BVU - calculate the probability that X > SH and Y > SK.
!       (to accuracy of 1e-16?)
!
! Parameters
!
!   SH  REAL, lower integration limit
!   SK  REAL, lower integration limit
!   R   REAL, correlation coefficient
!
!   LG  INTEGER, number of Gauss Rule Points and Weights
!
! Revised pab added check on XMAX
      DOUBLE PRECISION, INTENT(IN) :: SH, SK, R
      DOUBLE PRECISION  :: VAL
! Local variables
      DOUBLE PRECISION :: ZERO,ONE,FOUR
      DOUBLE PRECISION :: SQTWOPI ,TWOPI1,FOURPI1
      DOUBLE PRECISION :: HALF,ONETHIRD,ONEEIGHT,ONESIXTEEN
      DOUBLE PRECISION :: TWELVE, EXPMIN, XMAX
      INTEGER :: I, LG, NG
      PARAMETER ( ZERO = 0.D0,ONE=1.0D0,HALF=0.5D0)
      PARAMETER (FOUR = 4.0D0, TWELVE = 12.0D0)
      PARAMETER (EXPMIN =  -100.0D0)
      PARAMETER (ONESIXTEEN = 0.0625D0) !1/16
      PARAMETER (ONEEIGHT = 0.125D0 )  !1/8
      PARAMETER (ONETHIRD = 0.3333333333333333333333D0)
!      PARAMETER (TWOPI   = 6.283185307179586D0 )
      PARAMETER (TWOPI1  = 0.15915494309190D0 )  !1/(2*pi)
      PARAMETER (FOURPI1 = 0.0795774715459476D0 ) !/1/(4*pi)
      PARAMETER (SQTWOPI = 2.50662827463100D0) ! SQRT(2*pi)
      PARAMETER (XMAX    = 8.3D0)
      DOUBLE PRECISION, DIMENSION(10,3) :: X, W
      DOUBLE PRECISION :: AS, A, B, C, D, RS, XS
      DOUBLE PRECISION :: SN, ASR, H, K, BS, HS, HK
!     Gauss Legendre Points and Weights, N =  6
      DATA  ( W(I,1), X(I,1), I = 1,3) /
     *  0.1713244923791705D+00,-0.9324695142031522D+00,
     *  0.3607615730481384D+00,-0.6612093864662647D+00,
     *  0.4679139345726904D+00,-0.2386191860831970D+00/
!     Gauss Legendre Points and Weights, N = 12
      DATA ( W(I,2), X(I,2), I = 1,6) /
     *  0.4717533638651177D-01,-0.9815606342467191D+00,
     *  0.1069393259953183D+00,-0.9041172563704750D+00,
     *  0.1600783285433464D+00,-0.7699026741943050D+00,
     *  0.2031674267230659D+00,-0.5873179542866171D+00,
     *  0.2334925365383547D+00,-0.3678314989981802D+00,
     *  0.2491470458134029D+00,-0.1252334085114692D+00/
!     Gauss Legendre Points and Weights, N = 20
      DATA ( W(I,3), X(I,3), I = 1,10) /
     *  0.1761400713915212D-01,-0.9931285991850949D+00,
     *  0.4060142980038694D-01,-0.9639719272779138D+00,
     *  0.6267204833410906D-01,-0.9122344282513259D+00,
     *  0.8327674157670475D-01,-0.8391169718222188D+00,
     *  0.1019301198172404D+00,-0.7463319064601508D+00,
     *  0.1181945319615184D+00,-0.6360536807265150D+00,
     *  0.1316886384491766D+00,-0.5108670019508271D+00,
     *  0.1420961093183821D+00,-0.3737060887154196D+00,
     *  0.1491729864726037D+00,-0.2277858511416451D+00,
     *  0.1527533871307259D+00,-0.7652652113349733D-01/
      SAVE W, X
      VAL = ZERO
      HK = MIN(SH,SK)
      IF ( HK  < -XMAX) THEN    ! pab 24.05.2003
         VAL = FI(-MAX(SH,SK))
         RETURN
      ELSE IF ( XMAX  < MAX(SH,SK)) THEN
         RETURN
      ENDIF
      IF ( ABS(R)  <  0.3D0 ) THEN
         NG = 1
         LG = 3
      ELSE IF ( ABS(R)  <  0.75D0 ) THEN
         NG = 2
         LG = 6
      ELSE
         NG = 3
         LG = 10
      ENDIF
      H = SH
      K = SK
      HK = H*K

      IF ( ABS(R)  <  0.925D0 ) THEN
         IF (ABS(R) .GT. ZERO ) THEN
            HS = ( H*H + K*K )*HALF
            ASR = ASIN(R)
            DO I = 1, LG
               SN  = SIN(ASR*(ONE + X(I,NG))*HALF)
               VAL = VAL + W(I,NG)*EXP( ( SN*HK - HS )/( ONE - SN*SN ) )
               SN  = SIN(ASR*(ONE - X(I,NG))*HALF)
               VAL = VAL + W(I,NG)*EXP( ( SN*HK - HS )/( ONE - SN*SN ) )
            END DO
            VAL = VAL*ASR*FOURPI1
         ENDIF
         VAL = VAL + FI(-H)*FI(-K)
      ELSE
         IF ( R  <  ZERO ) THEN
            K  = -K
            HK = -HK
         ENDIF
         IF ( ABS(R)  <  ONE ) THEN
            AS  = ( ONE - R )*( ONE + R )
            A   = SQRT(AS)
            B   = ABS( H - K ) !**2
            BS  = B * B
            C   = ( FOUR - HK ) * ONEEIGHT !/8D0
            D   = ( TWELVE - HK ) * ONESIXTEEN !/16D0
            ASR =  -(BS/AS + HK)*HALF
            IF (ASR.GT.EXPMIN) THEN
               VAL = A*EXP( ASR ) *
     &              ( ONE - C*(BS - AS)*(ONE - D*BS*0.2D0)*ONETHIRD +
     &              C*D*AS*AS*0.2D0 )
            ENDIF
            IF ( HK .GT. EXPMIN ) THEN
               VAL = VAL - EXP(-HK*HALF)*SQTWOPI*FI(-B/A)*B
     +              *( ONE - C*BS*( ONE - D*BS*0.2D0 )*ONETHIRD )
            ENDIF
            A = A * HALF
            DO I = 1, LG
               XS  = ( A * (ONE + X(I,NG)) ) !**2
               XS  = XS * XS
               RS  = SQRT( ONE - XS )
               ASR = -(BS / XS + HK) * HALF
               IF (ASR.GT.EXPMIN) THEN
                  VAL = VAL + A*W(I,NG)*EXP( ASR )
     &                 * ( EXP( - HALF*HK*( ONE - RS )/( ONE + RS ) )/RS
     $                 -( ONE + C*XS*( ONE + D*XS ) ) )
               ENDIF
               XS  = ( A * (ONE - X(I,NG)) ) !**2
               XS  = XS * XS
               RS  = SQRT( ONE - XS )
               ASR = -(BS / XS + HK) * HALF
               IF (ASR.GT.EXPMIN) THEN
                  VAL = VAL + A*W(I,NG)*EXP( ASR )
     &                 *( EXP( - HALF*HK*( ONE - RS )/( ONE + RS ) )/RS-
     $                 ( ONE + C*XS*( ONE + D*XS ) ) )
               ENDIF
            END DO
            VAL = -VAL*TWOPI1
         ENDIF
         IF ( R .GT. ZERO ) THEN
            VAL =  VAL + FI( -MAX( H, K ) )
         ELSE
            VAL = -VAL
            IF ( H  <  K )  VAL = VAL +  FI(K)-FI(H)
         ENDIF
      ENDIF
      RETURN
      END FUNCTION BVU
      DOUBLE PRECISION FUNCTION STUDNT( NU, T )
      IMPLICIT NONE
!
!     Student t Distribution Function
!
!                       T
!         STUDNT = C   I  ( 1 + y*y/NU )**( -(NU+1)/2 ) dy
!                   NU -INF
!
      INTEGER, INTENT(IN) :: NU
      DOUBLE PRECISION, INTENT(IN) :: T
!     Locals
      INTEGER :: J
      DOUBLE PRECISION :: ZRO, ONE
      PARAMETER ( ZRO = 0.0D0, ONE = 1.0D0 )
      DOUBLE PRECISION, PARAMETER :: PI = 3.14159265358979D0
      DOUBLE PRECISION :: CSSTHE, SNTHE, POLYN, TT, TS, RN

      IF ( NU  <  1 ) THEN
         STUDNT = FI( T )
      ELSE IF ( NU .EQ. 1 ) THEN
         STUDNT = ( ONE + 2.0D0*ATAN(T)/PI )*0.5D0
      ELSE IF ( NU .EQ. 2 ) THEN
         STUDNT = ( ONE + T/SQRT( 2.0D0 + T*T ))*0.5D0
      ELSE
         RN = NU ! convert to double
         TT = T * T
         CSSTHE = ONE/( ONE + TT/RN )
         POLYN = 1
         DO J = NU-2, 2, -2
            POLYN = ONE + ( J - 1 )*CSSTHE*POLYN/J
         END DO

         IF ( MOD( NU, 2 ) .EQ. 1 ) THEN
            TS = T/SQRT(RN)
            STUDNT = ( ONE + 2.0D0*( ATAN(TS) +
     &           TS*CSSTHE*POLYN )/PI )*0.5D0
         ELSE
            SNTHE = T/SQRT( RN + TT )
            STUDNT = ( ONE + SNTHE*POLYN )*0.5D0
         END IF
         STUDNT = MAX( ZRO, MIN( STUDNT, ONE ) )
      ENDIF
      END FUNCTION STUDNT
      DOUBLE PRECISION FUNCTION BVTL( NU, DH, DK, R )
      IMPLICIT NONE
!*
!*     A function for computing bivariate t probabilities.
!*
!*       Alan Genz
!*       Department of Mathematics
!*       Washington State University
!*       Pullman, WA 99164-3113
!*       Email : alangenz@wsu.edu
!*
!*    This function is based on the method described by
!*        Dunnett, C.W. and M. Sobel, (1954),
!*        A bivariate generalization of Student's t-distribution
!*        with tables for certain special cases,
!*        Biometrika 41, pp. 153-169.
!*
!* BVTL - calculate the probability that X < DH and Y < DK.
!*
!* parameters
!*
!*   NU number of degrees of freedom (NOTE: NU = 0 gives bivariate normal prb)
!*   DH 1st lower integration limit
!*   DK 2nd lower integration limit
!*   R   correlation coefficient
!*
      INTEGER, INTENT(IN) ::NU
      DOUBLE PRECISION, INTENT(IN) :: DH, DK, R
!     Locals
      INTEGER :: J, HS, KS
      DOUBLE PRECISION ::  ORS, HRK, KRH, BVT
      DOUBLE PRECISION ::  DH2, DK2, SNU ,DNU, DHDK
!, BVND, STUDNT
      DOUBLE PRECISION :: GMPH, GMPK, XNKH, XNHK, QHRK, HKN, HPK, HKRN
      DOUBLE PRECISION :: BTNCKH, BTNCHK, BTPDKH, BTPDHK
      DOUBLE PRECISION :: ZERO, ONE, EPS, PI,TPI
      PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, EPS = 1.0D-15 )
      PARAMETER (PI =  3.14159265358979D0, TPI = 6.28318530717959D0)
      IF ( NU  <  1 ) THEN
         BVTL = BVU( -DH, -DK, R )
      ELSE IF ( ONE - R .LE. EPS .OR. 1.0D+16<MAX(ABS(DH),ABS(DK))) THEN
            BVTL = STUDNT( NU, MIN( DH, DK ) )
      ELSE IF ( R + ONE  .LE. EPS ) THEN
         IF ( DH .GT. -DK )  THEN
            BVTL = STUDNT( NU, DH ) - STUDNT( NU, -DK )
         ELSE
            BVTL = ZERO
         END IF
      ELSE
         !PI = ACOS(-ONE)
         !TPI = 2*PI
         DNU = NU           ! convert to double
         SNU = SQRT(DNU)
         ORS = ONE - R * R
         HRK = DH - R * DK
         KRH = DK - R * DH
         DK2 = DK * DK
         DH2 = DH * DH
         IF ( ABS(HRK) + ORS .GT. ZERO ) THEN
            XNHK = HRK**2/( HRK**2 + ORS*( DNU + DK2 ) )
            XNKH = KRH**2/( KRH**2 + ORS*( DNU + DH2 ) )
         ELSE
            XNHK = ZERO
            XNKH = ZERO
         END IF
         HS = INT(DSIGN( ONE, HRK)) !DH - R*DK )
         KS = INT(DSIGN( ONE, KRH)) !DK - R*DH )
         IF ( MOD( NU, 2 ) .EQ. 0 ) THEN
            BVT = ATAN2( SQRT(ORS), -R )/TPI
            GMPH = DH/SQRT( 16.0D0*( DNU + DH2 ) )
            GMPK = DK/SQRT( 16.0D0*( DNU + DK2 ) )
            BTNCKH = 2*ATAN2( SQRT( XNKH ), SQRT( ONE - XNKH ) )/PI
            BTPDKH = 2*SQRT( XNKH*( ONE - XNKH ) )/PI
            BTNCHK = 2*ATAN2( SQRT( XNHK ), SQRT( ONE - XNHK ) )/PI
            BTPDHK = 2*SQRT( XNHK*( ONE - XNHK ) )/PI
            DO J = 1, NU/2
               BVT = BVT + GMPH*( ONE + KS*BTNCKH )
               BVT = BVT + GMPK*( ONE + HS*BTNCHK )
               BTNCKH = BTNCKH + BTPDKH
               BTPDKH = 2*J*BTPDKH*( ONE - XNKH )/( 2*J + 1 )
               BTNCHK = BTNCHK + BTPDHK
               BTPDHK = 2*J*BTPDHK*( ONE - XNHK )/( 2*J + 1 )
               GMPH = GMPH*( 2*J - 1 )/( 2*J*( ONE + DH2/DNU ) )
               GMPK = GMPK*( 2*J - 1 )/( 2*J*( ONE + DK2/DNU ) )
            END DO
         ELSE  ! NU is ODD
            DHDK = DH*DK
            QHRK = SQRT( DH2 + DK2 - 2.0D0*R*DHDK + DNU*ORS )
            HKRN = DHDK + R*DNU
            HKN = DHDK - DNU
            HPK = DH + DK
            BVT = ATAN2( -SNU*( HKN*QHRK + HPK*HKRN ),
     &                          HKN*HKRN-DNU*HPK*QHRK )/TPI
            IF ( BVT  <  -EPS ) BVT = BVT + ONE
            GMPH = DH/( TPI*SNU*( ONE + DH2/DNU ) )
            GMPK = DK/( TPI*SNU*( ONE + DK2/DNU ) )
            BTNCKH = SQRT( XNKH )
            BTPDKH = BTNCKH
            BTNCHK = SQRT( XNHK )
            BTPDHK = BTNCHK
            DO J = 1, ( NU - 1 )/2
               BVT = BVT + GMPH*( ONE + KS*BTNCKH )
               BVT = BVT + GMPK*( ONE + HS*BTNCHK )
               BTPDKH = ( 2*J - 1 )*BTPDKH*( ONE - XNKH )/( 2*J )
               BTNCKH = BTNCKH + BTPDKH
               BTPDHK = ( 2*J - 1 )*BTPDHK*( ONE - XNHK )/( 2*J )
               BTNCHK = BTNCHK + BTPDHK
               GMPH = 2*J*GMPH/( ( 2*J + 1 )*( ONE + DH2/DNU ) )
               GMPK = 2*J*GMPK/( ( 2*J + 1 )*( ONE + DK2/DNU ) )
            END DO
         END IF
         BVTL = BVT
      END IF
      END FUNCTION BVTL

      DOUBLE PRECISION FUNCTION TVTL( NU1, H, R, EPSI )
      USE TRIVARIATEVAR !Block to transfer variables to TVTMFN
      IMPLICIT NONE
!
!     A function for computing trivariate normal and t-probabilities.
!     This function uses algorithms developed from the ideas
!     described in the papers:
!       R.L. Plackett, Biometrika 41(1954), pp. 351-360.
!       Z. Drezner, Math. Comp. 62(1994), pp. 289-294.
!     with adaptive integration from (0,0,1) to (0,0,r23) to R.
!
!      Calculate the probability that X(I) < H(I), for I = 1,2,3
!    NU   INTEGER degrees of freedom; use NU = 0 for normal cases.
!    H    REAL array of uppoer limits for probability distribution
!    R    REAL array of three correlation coefficients, R should
!         contain the lower left portion of the correlation matrix r.
!         R should contains the values r21, r31, r23 in that order.
!   EPSI  REAL required absolute accuracy; maximum accuracy for most
!          computations is approximately 1D-14
!
!    The software is based on work described in the paper
!     "Numerical Computation of Rectangular Bivariate and Trivariate
!      Normal and t Probabilities", by the code author:
!
!       Alan Genz
!       Department of Mathematics
!       Washington State University
!       Pullman, WA 99164-3113
!       Email : alangenz@wsu.edu
!
!      EXTERNAL TVTMFN
      INTEGER,                       INTENT(IN) :: NU1
      DOUBLE PRECISION,DIMENSION(:), INTENT(IN) :: H, R
      DOUBLE PRECISION,              INTENT(IN) ::  EPSI
!Locals
      DOUBLE PRECISION :: R12, R13,  TVT
      DOUBLE PRECISION :: ONE, ZERO, EPS,  PT
!      DOUBLE PRECISION RUA, RUB, AR, RUC,
!        BVTL, PHID, ADONET
      PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
      PARAMETER ( PT = 1.57079632679489661923132169163975D0 ) !pi/2

!      COMMON /TVTMBK/ H1, H2, H3, R23, RUA, RUB, AR, RUC, NU
      EPS = MAX( 1.0D-13, EPSI )
!      PT = ASIN(ONE)
      NU = NU1
      H1 = H(1)
      H2 = H(2)
      H3 = H(3)
      R12 = R(1)
      R13 = R(2)
      R23 = R(3)
!
!     Sort R's and check for special cases
!
      IF ( ABS(R12) .GT. ABS(R13) ) THEN
         H2 = H3
         H3 = H(2)
         R12 = R13
         R13 = R(1)
      END IF
      IF ( ABS(R13) .GT. ABS(R23) ) THEN
         H1 = H2
         H2 = H(1)
         R23 = R13
         R13 = R(3)
      END IF
      TVT = 0
      IF ( ABS(H1) + ABS(H2) + ABS(H3)  <  EPS ) THEN
         TVT = (ONE + (ASIN(R12) + ASIN(R13) + ASIN(R23))/PT )*0.125D0
      ELSE IF ( NU  <  1 .AND. ABS(R12) + ABS(R13)  <  EPS ) THEN
         TVT = FI(H1)*BVTL( NU, H2, H3, R23 )
      ELSE IF ( NU  <  1 .AND. ABS(R13) + ABS(R23)  <  EPS ) THEN
         TVT = FI(H3)*BVTL( NU, H1, H2, R12 )
      ELSE IF ( NU  <  1 .AND. ABS(R12) + ABS(R23)  <  EPS ) THEN
         TVT = FI(H2)*BVTL( NU, H1, H3, R13 )
      ELSE IF ( ONE - R23  <  EPS ) THEN
         TVT = BVTL( NU, H1, MIN( H2, H3 ), R12 )
      ELSE IF ( R23 + ONE  <  EPS ) THEN
         IF  ( H2 .GT. -H3 )
     &        TVT = BVTL( NU, H1, H2, R12 ) - BVTL( NU, H1, -H3, R12 )
      ELSE
!
!        Compute singular TVT value
!
         IF ( NU  <  1 ) THEN
            TVT = BVTL( NU, H2, H3, R23 )*FI(H1)
         ELSE IF ( R23 .GE. ZERO ) THEN
            TVT = BVTL( NU, H1, MIN( H2, H3 ), ZERO )
         ELSE IF ( H2 .GT. -H3 ) THEN
            TVT = BVTL( NU, H1, H2, ZERO ) - BVTL( NU, H1, -H3, ZERO )
         END IF
!
!        Use numerical integration to compute probability
!
!
         RUA = ASIN( R12 )
         RUB = ASIN( R13 )
         AR  = ASIN( R23)
         RUC = SIGN( PT, AR ) - AR
         TVT = TVT + ADONET(TVTMFN, ZERO, ONE, EPS )/( 4.D0*PT )
      END IF
      TVTL = MAX( ZERO, MIN( TVT, ONE ) )
      END FUNCTION TVTL
!      CONTAINS
!
      DOUBLE PRECISION FUNCTION TVTMFN( X )
      USE TRIVARIATEVAR
      IMPLICIT NONE
!
!     Computes Plackett formula integrands
!
      DOUBLE PRECISION, INTENT(IN) :: X
! Locals
      DOUBLE PRECISION R12, RR2, R13, RR3, R, RR, ZRO !, PNTGND
! Parameters transfeered from TRIVARIATEVAR
!      INTEGER :: NU
!      DOUBLE PRECISION :: H1, H2, H3, R23, RUA, RUB, AR, RUC
      PARAMETER ( ZRO = 0.0D0 )
!      COMMON /TVTMBK/ H1, H2, H3, R23, RUA, RUB, AR, RUC, NU
      TVTMFN = 0.0D0
      CALL SINCS( RUA*X, R12, RR2 )
      CALL SINCS( RUB*X, R13, RR3 )
      IF ( ABS(RUA) .GT. ZRO )
     &     TVTMFN = TVTMFN + RUA*PNTGND( NU, H1,H2,H3, R13,R23,R12,RR2 )
      IF ( ABS(RUB) .GT. ZRO )
     &     TVTMFN = TVTMFN + RUB*PNTGND( NU, H1,H3,H2, R12,R23,R13,RR3 )
      IF ( NU .GT. 0 ) THEN
         CALL SINCS( AR + RUC*X, R, RR )
         TVTMFN = TVTMFN - RUC*PNTGND( NU, H2, H3, H1, ZRO, ZRO, R, RR )
      END IF
      END FUNCTION TVTMFN
!
      SUBROUTINE SINCS( X, SX, CS )
!
!     Computes SIN(X), COS(X)^2, with series approx. for |X| near PI/2
!
      DOUBLE PRECISION, INTENT(IN) :: X
      DOUBLE PRECISION, INTENT(OUT) :: SX, CS
!Locals
      DOUBLE PRECISION :: EE, PT, KS, KC, ONE, SMALL, HALF, ONETHIRD
      PARAMETER (ONE = 1.0D0, SMALL = 5.0D-5, HALF = 0.5D0 )
      PARAMETER ( PT = 1.57079632679489661923132169163975D0 )
      PARAMETER ( KS = 0.0833333333333333333333333333333D0) !1/12
      PARAMETER ( KC = 0.1333333333333333333333333333333D0) !2/15
      PARAMETER ( ONETHIRD = 0.33333333333333333333333333333333D0) !1/3
      EE = ( PT - ABS(X) )
      EE = EE * EE
      IF ( EE  <  SMALL ) THEN
         SX = SIGN( ONE - EE*( ONE - EE*KS )*HALF, X )
         CS = EE *( ONE - EE*( ONE - EE*KC )*ONETHIRD)
      ELSE
         SX = SIN(X)
         CS = ONE - SX*SX
      END IF
      END SUBROUTINE SINCS
!
      DOUBLE PRECISION FUNCTION PNTGND( NU, BA, BB, BC, RA, RB, R, RR )
      IMPLICIT NONE
!
!     Computes Plackett formula integrand
!
      INTEGER, INTENT(IN) :: NU
      DOUBLE PRECISION, INTENT(IN) :: BA, BB, BC, RA, RB, R, RR
! Locals
      DOUBLE PRECISION :: DT, FT, BT,RAB2, BARB, rNU!, PHID, STUDNT
      PNTGND = 0.0D0
      FT   = ( RA - RB )
      RAB2 = FT*FT
      DT   = RR*( RR -  RAB2 - 2.0D0*RA*RB*( 1.D0 - R ) )
      IF ( DT .GT. 0.0D0 ) THEN
         BT   = ( BC*RR + BA*( R*RB - RA ) + BB*( R*RA -RB ) )/SQRT(DT)
         BARB = ( BA - R*BB )
         FT   = ( BARB * BARB ) / RR + BB * BB
         IF ( NU  <  1 ) THEN
            IF ( BT .GT. -10.0D0 .AND. FT  <  100.0D0 ) THEN
               PNTGND = EXP( -FT * 0.5D0 ) * FI( BT )
!               PNTGND = EXP( -FT*0.5D0)
!               IF ( BT  <  10.0D0 ) PNTGND = PNTGND * FI(BT)
            END IF
         ELSE
            rNU = NU
            FT  = SQRT( 1.0D0 + FT/rNU )
            PNTGND = STUDNT( NU, BT/FT )/FT**NU
         END IF
      END IF
      END  FUNCTION PNTGND
!
      DOUBLE PRECISION FUNCTION ADONET( F, A, B, TOL )
      IMPLICIT NONE
!
!     One Dimensional Globally Adaptive Integration Function
!
!      EXTERNAL F
      DOUBLE PRECISION, INTENT(IN) :: A, B, TOL
      INTEGER :: NL, I, IM, IP
      PARAMETER ( NL = 100 )
      DOUBLE PRECISION, DIMENSION(NL) ::  EI, AI, BI, FI
      DOUBLE PRECISION :: FIN, ERR !, KRNRDT
      DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0
      DOUBLE PRECISION, PARAMETER :: HALF = 0.5D0
      DOUBLE PRECISION, PARAMETER :: ONE  = 1.0D0
      DOUBLE PRECISION, PARAMETER :: FOUR = 4.0D0
      INTERFACE
         FUNCTION F(Z) RESULT (VAL)
         DOUBLE PRECISION, INTENT(IN) :: Z
         DOUBLE PRECISION :: VAL
         END FUNCTION F
      END INTERFACE
!      COMMON /ABLK/ ERR, IM
      AI(1) = A
      BI(1) = B
      ERR = ONE
      IP = 1
      IM = 1
      DO WHILE ( (ERR .GT. TOL) .AND. (IM  <  NL) )
         IM = IM + 1
         BI(IM) = BI(IP)
         AI(IM) = ( AI(IP) + BI(IP) ) * HALF
         BI(IP) = AI(IM)
         FI(IP) = KRNRDT( AI(IP), BI(IP), F, EI(IP) )
         FI(IM) = KRNRDT( AI(IM), BI(IM), F, EI(IM) )
         ERR = ZERO
         FIN = ZERO
         DO I = 1, IM
            IF ( EI(I) .GT. EI(IP) ) IP = I
            FIN = FIN + FI(I)
            ERR = ERR + EI(I)*EI(I)
         END DO
         ERR = FOUR * SQRT( ERR )
      END DO
      ADONET = FIN
      END FUNCTION ADONET
!
      DOUBLE PRECISION FUNCTION KRNRDT( A, B, F, ERR )
!
!     Kronrod Rule
!
      DOUBLE PRECISION, intent(in) :: A, B
      DOUBLE PRECISION, intent(out) :: ERR
      DOUBLE PRECISION T, CEN, FC, WID, RESG, RESK
!
!        The abscissae and weights are given for the interval (-1,1);
!        only positive abscissae and corresponding weights are given.
!
!        XGK    - abscissae of the 2N+1-point Kronrod rule:
!                 XGK(2), XGK(4), ...  N-point Gauss rule abscissae;
!                 XGK(1), XGK(3), ...  optimally added abscissae.
!        WGK    - weights of the 2N+1-point Kronrod rule.
!        WG     - weights of the N-point Gauss rule.
!
      INTEGER :: J, N
      PARAMETER ( N = 11 )
      DOUBLE PRECISION, PARAMETER :: HALF = 0.5D0
      DOUBLE PRECISION WG(0:(N+1)/2), WGK(0:N), XGK(0:N)
      SAVE WG, WGK, XGK
      INTERFACE
          FUNCTION F(Z) RESULT (VAL)
         DOUBLE PRECISION, INTENT(IN) :: Z
         DOUBLE PRECISION :: VAL
         END FUNCTION F
      END INTERFACE
      DATA WG( 0)/ 0.2729250867779007D+00/
      DATA WG( 1)/ 0.5566856711617449D-01/
      DATA WG( 2)/ 0.1255803694649048D+00/
      DATA WG( 3)/ 0.1862902109277352D+00/
      DATA WG( 4)/ 0.2331937645919914D+00/
      DATA WG( 5)/ 0.2628045445102478D+00/
!
      DATA XGK( 0)/ 0.0000000000000000D+00/
      DATA XGK( 1)/ 0.9963696138895427D+00/
      DATA XGK( 2)/ 0.9782286581460570D+00/
      DATA XGK( 3)/ 0.9416771085780681D+00/
      DATA XGK( 4)/ 0.8870625997680953D+00/
      DATA XGK( 5)/ 0.8160574566562211D+00/
      DATA XGK( 6)/ 0.7301520055740492D+00/
      DATA XGK( 7)/ 0.6305995201619651D+00/
      DATA XGK( 8)/ 0.5190961292068118D+00/
      DATA XGK( 9)/ 0.3979441409523776D+00/
      DATA XGK(10)/ 0.2695431559523450D+00/
      DATA XGK(11)/ 0.1361130007993617D+00/
!
      DATA WGK( 0)/ 0.1365777947111183D+00/
      DATA WGK( 1)/ 0.9765441045961290D-02/
      DATA WGK( 2)/ 0.2715655468210443D-01/
      DATA WGK( 3)/ 0.4582937856442671D-01/
      DATA WGK( 4)/ 0.6309742475037484D-01/
      DATA WGK( 5)/ 0.7866457193222764D-01/
      DATA WGK( 6)/ 0.9295309859690074D-01/
      DATA WGK( 7)/ 0.1058720744813894D+00/
      DATA WGK( 8)/ 0.1167395024610472D+00/
      DATA WGK( 9)/ 0.1251587991003195D+00/
      DATA WGK(10)/ 0.1312806842298057D+00/
      DATA WGK(11)/ 0.1351935727998845D+00/
!
!           Major variables
!
!           CEN  - mid point of the interval
!           WID  - half-length of the interval
!           RESG - result of the N-point Gauss formula
!           RESK - result of the 2N+1-point Kronrod formula
!
!           Compute the 2N+1-point Kronrod approximation to
!            the integral, and estimate the absolute error.
!
      WID = ( B - A ) * HALF
      CEN = ( B + A ) * HALF
      FC  = F(CEN)
      RESG = FC * WG(0)
      RESK = FC * WGK(0)
      DO J = 1, N
         T  = WID * XGK(J)
         FC = F( CEN - T ) + F( CEN + T )
         RESK = RESK + WGK(J) * FC
         IF( MOD( J, 2 ) .EQ. 0 ) RESG = RESG + WG(J/2) * FC
      END DO
      KRNRDT = WID * RESK
      ERR = ABS( WID * ( RESK - RESG ) )
      END FUNCTION KRNRDT
!      END FUNCTION TVTL

      FUNCTION GAUSINT (X1, X2, A, B, C, D) RESULT (value)
!      USE GLOBALDATA,ONLY:  xCutOff
      IMPLICIT NONE
      DOUBLE PRECISION, INTENT(in) :: X1,X2,A,B,C,D
      DOUBLE PRECISION             :: value
!Local variables
      DOUBLE PRECISION             :: Y1,Y2,Y3
      DOUBLE PRECISION, PARAMETER :: SQTWOPI1=3.9894228040143d-1 !=1/sqrt(2*pi)
      DOUBLE PRECISION, PARAMETER :: XMAX = 37.d0
      ! Let  X  be standardized Gaussian variable,
      ! i.e., X=N(0,1). The function calculate the
      !  following integral E[I(X1<X<X2)(A+BX)(C+DX)
      ! where I(X1<X<X2) is an indicator function of
      ! the set {X1<X<X2}.
      IF (X1.GE.X2) THEN
         value = 0.d0
         RETURN
      ENDIF
      IF (ABS (X1) .GT.XMAX) THEN
         Y1 = 0.d0
      ELSE
         Y1 = (A * D+B * C + X1 * B * D) * EXP ( - 0.5d0 * X1 * X1)
      ENDIF
      IF (ABS (X2) .GT.XMAX) THEN
         Y2 = 0.d0
      ELSE
         Y2 = (A * D+B * C + X2 * B * D) * EXP ( - 0.5d0 * X2 * X2)
      ENDIF
      Y3 = (A * C + B * D) * (FI (X2) - FI (X1) )
      value = Y3 + SQTWOPI1 * (Y1 - Y2)
      RETURN
      END FUNCTION GAUSINT


      FUNCTION GAUSINT2 (X1, X2, A, B) RESULT (value)
      IMPLICIT NONE
      DOUBLE PRECISION, INTENT(in) :: X1,X2,A,B
      DOUBLE PRECISION             :: value
! Local variables
      DOUBLE PRECISION             :: X0,Y0,Y1,Y2
      DOUBLE PRECISION, PARAMETER :: SQTWOPI1=3.9894228040143d-1 !=1/sqrt(2*pi)
!     Let  X  be standardized Gaussian variable,
!     i.e., X=N(0,1). The function calculate the
!     following integral E[I(X1<X<X2)ABS(A+BX)]
!     where I(X1<X<X2) is an indicator function of
!     the set {X1<X<X2}.
      IF (X1.GE.X2) THEN
         value = 0.d0
         RETURN
      ENDIF
      IF (ABS(B).EQ.0.d0) THEN
         value = ABS(A)*(FI(X2)-FI(X1))
         RETURN
      ENDIF

      Y1 = -A*FI(X1)+SQTWOPI1*B*EXP(-0.5d0*X1*X1)
      Y2 = A*FI(X2)-SQTWOPI1*B*EXP(-0.5d0*X2*X2)
      IF ((B*X1 < -A).AND.(-A < B*X2))THEN
         X0 = -A/B
         Y0 = 2.d0*(A*FI(X0)-SQTWOPI1*B*EXP(-0.5d0*X0*X0))
         value=ABS(Y2-Y1-Y0)
      ELSE
         value=ABS(Y1+Y2)
      ENDIF
      RETURN
      END FUNCTION GAUSINT2

      SUBROUTINE EXLMS(A, X1, X2, INFIN, LOWER, UPPER, Ca,Pa)
      IMPLICIT NONE
      DOUBLE PRECISION, INTENT(in) :: A, X1, X2
      DOUBLE PRECISION, INTENT(out) :: LOWER, UPPER,Ca,Pa
      INTEGER,INTENT(in) :: INFIN
      ! Local variables
      DOUBLE PRECISION :: P1
!      DOUBLE PRECISION, PARAMETER :: AMAX = 5D15

      ! Let  X  be standardized Gaussian variable,
      ! i.e., X=N(0,1). The function calculate the
      !  following integral E[I(X1<X<X2)ABS(A+X)]/C(A)
      ! where I(X1<X<X2) is an indicator function of
      ! the set {X1<X<X2} and C(A) is a normalization factor
      ! i.e. E(I(-inf<X<inf)*ABS(A+X)) = C(A)
      !
      ! Pa = probability at the inflection point of the CDF
      ! Ca = C(A) normalization parameter
      ! A  = location parameter of the inflection point
      ! X1,X2 = integration limits < infinity
      LOWER = 0.d0
      UPPER = 1.d0

      P1 = EXFUN(-A,A)
      Ca = A+2D0*P1
      Pa = P1/Ca
      IF ( INFIN  <  0 ) RETURN
      IF ( INFIN .NE. 0 ) THEN
         IF (X1.LE.-A) THEN
            LOWER = EXFUN(X1,A)/Ca
         ELSE
            LOWER = 1D0-EXFUN(-X1,-A)/Ca
         ENDIF
      ENDIF
      IF ( INFIN .NE. 1 ) THEN
         IF (X2.LE.-A) THEN
            UPPER = EXFUN(X2,A)/Ca
         ELSE
            UPPER = 1D0-EXFUN(-X2,-A)/Ca
         ENDIF
      ENDIF

      RETURN
      CONTAINS
      FUNCTION EXFUN(X,A) RESULT (P)
      DOUBLE PRECISION, INTENT(in) :: X,A
      DOUBLE PRECISION :: P
      DOUBLE PRECISION, PARAMETER :: SQTWOPI1=3.9894228040143d-1 !=1/sqrt(2*pi)
      P = EXP(-X*X*0.5D0)*SQTWOPI1-A*FI(X)
      END FUNCTION EXFUN
      END SUBROUTINE EXLMS


      FUNCTION EXINV(P,A,Ca,Pa) RESULT (VAL)
!EXINV calculates the inverse of the CDF of abs(x+A)*exp(-x^2/2)/SQRT(2*pi)/C(A)
!      where C(A) is a normalization parameter depending on A
!
!   CALL:  val   = exinv(P,A,Ca,Pa)
!
!     val  = quantiles
!     p    = probabilites
!     A    = location parameter
!     Ca   = normalization parameter
!     Pa   = excdf(-A,A) probability at the inflection point of the CDF
      double precision, intent(in) :: P,A,Ca,Pa
      double precision :: val
! local variables
      double precision, parameter :: amax = 5.D15
      double precision, parameter :: epsl = 5.D-15
      double precision, parameter :: xmax = 8
      double precision :: P1,Xk,Ak,Zk,SGN

!      if (P<0.D0.OR.P.GT.1.D0) PRINT *,'warning P<0 or P>1'

! The inverse cdf of 0 is -inf, and the inverse cdf of 1 is inf.
      if (P.LE.EPSL.OR.P+EPSL.GE.1.D0) THEN
         VAL = SIGN(xmax,P-0.5D0)
         return
      endif
      Ak = ABS(A)
      if (EPSL < Ak .AND. Ak < amax) THEN
         IF (ABS(p-Pa).LE.EPSL) THEN
            VAL = SIGN(MIN(Ak,xmax),-A)
            RETURN
         ENDIF
         IF (Ak < 1D-2) THEN   ! starting guess always less than 0.2 from the true value
            IF (P.GE.0.5D0) THEN
               xk = SQRT(-2D0*log(2D0*(1D0-P)))
            ELSE
               xk = -SQRT(-2D0*log(2D0*P))
            ENDIF
         ELSE
            xk = FIINV(P)       ! starting guess always less than 0.8 from the true value
            ! Modify starting guess if possible in order to speed up Newtons method
            IF (1D-3.LE.P.AND. P.LE.0.99D0.AND.
     &           3.5.LE.Ak.AND.Ak.LE.1D3 ) THEN
               SGN = SIGN(1.d0,-A)
               Zk = xk*SGN
               xk = SGN*(Zk+((1D0/(64.9495D0*Ak-178.3191D0)-0.02D0/Ak)*
     &              Zk+1D0/(-0.99679234298211D0*Ak-0.07195350071872D0))/
     &            (Zk/(-1.48430620263825D0*Ak-0.33340759016175D0)+1D0))
            ELSEIF ((P < 1D-3.AND.A.LE.-3.5D0).OR.
     &              (3.5D0.LE.A.AND.P.GT.0.99D0)) THEN
               SGN = SIGN(1.d0,-A)
               Zk = xk*SGN
               P1 = -2.00126182192701D0*Ak-2.57306603933111D0
               xk = SGN*Zk*(1D0+
     &              P1/((-0.99179258785909D0*Ak-0.21359746002397D0)*
     &              (Zk+P1)))
            ENDIF
         ENDIF
         ! Check if the starting guess is on the correct side of the inflection point
         IF (xk.LE.-A .AND. P.GT.Pa) xk = 1.D-2-A
         IF (xk.GE.-A .AND. P < Pa) xk = -1.D-2-A


         IF (P < Pa) THEN
            VAL = funca(xk,A,P*Ca)
         ELSE  ! exploit the symmetry of the CDF
            VAL = -funca(-xk,-A,(1.D0-P)*Ca)
         ENDIF

      ELSEIF (ABS(A).LE.EPSL) THEN
         IF (P>=0.5D0) THEN
            VAL = SQRT(-2D0*log(2D0*(1.D0-P)))
         ELSE
            VAL = -SQRT(-2D0*log(2D0*P))
         ENDIF
      ELSE  ! ABS(A) > AMAX
         VAL = FIINV(P)
      ENDIF
      !CALL EXLMS(A,0.d0,VAL,0,ak,P1,zk,sgn)
      !If (ABS(p-P1).GT.0.0001) PRINT *,'excdf(x,a)-p',p-P1
      RETURN

      CONTAINS

      function funca(xk0,ak,CaP) RESULT (xk)
      double precision, intent(in) :: xk0,ak,CaP ! =Ca*P
      DOUBLE PRECISION :: xk
!Local variables
      INTEGER,          PARAMETER :: ixmax = 25
      double precision, parameter :: crit = 7.1D-08 ! = sqrt(1e-15)
      double precision, parameter :: SQTWOPI1 = 0.39894228040143D0 !=1/SQRT(2*pi)
      double precision, parameter :: SQTWOPI = 2.50662827463100D0 !=SQRT(2*pi)
      INTEGER :: IX
      DOUBLE PRECISION :: H,H1,tmp0,tmp1,XNEW
      ! Newton's Method or Fixed point iteration to find the inverse of the EXCDF.
      ! Assumption: xk0 < -ak and xk < -ak
      ! Permit no more than IXMAX iterations.
      IX = 0
      H  = 1.D0
      xk = xk0    ! starting guess for the iteration


!     Break out of the iteration loop for the following:
!     1) The last update is very small (compared to x).
!     2) The last update is very small (compared to sqrt(eps)=crit).
!     3) There are more than 15 iterations. This should NEVER happen.
      IF (.TRUE..OR.ABS(ak) < 1.D-2) THEN
      ! Newton's method
      !~~~~~~~~~~~~~~~~~
      DO WHILE( ABS(H).GT.MIN(crit*ABS(xk),crit).AND.IX < IXMAX)

         IX = IX+1
                                !print *,'Iteration ',IX

         tmp0  = FI(xk)
         tmp1  = EXP(-xk*xk*0.5D0)*SQTWOPI1 ! =normpdf(x)
         H1 = (tmp1-ak*tmp0-CaP)/(ABS(xk+ak)*tmp1)
         H  = DSIGN(MIN(ABS(H1),0.7D0/DBLE(IX)),H1) ! Only allow smaller and smaller steps

         xnew = xk - H
                                ! Make sure that the current guess is less than -a.
                                ! When Newton's Method suggests steps that lead to -a guesses
                                ! take a step 9/10ths of the way to -a:
         IF (xnew.GT.-ak-crit) THEN
            xnew = (xk - 9.D0*ak)*1D-1
            H    = xnew - xk
         ENDIF
         xk = xnew
      END DO
      ELSE                      ! FIXED POINT iteration
                                !~~~~~~~~~~~~~~~~~~~~~~~
         DO WHILE (ABS(H).GT.MIN(crit*ABS(xk),crit).AND.IX < IXMAX)
            IX   = IX+1
            tmp0 = SQTWOPI1*EXP(-xk*xk*0.5D0)/FI(xk)
            tmp1 = -2.D0*LOG(SQTWOPI*CaP*tmp0/(tmp0-ak))
            SGN  = sign(1.D0,tmp1)
            xnew = -SQRT(SGN*tmp1)*SGN
            ! Make sure that the current guess is less than -a.
            ! When this method suggests steps that lead to -a guesses
            ! take a step 9/10ths of the way to -a:
            IF (xnew.GT.-ak-crit) xnew = (xk - 9.D0*ak)*1.D-1

            H  = xnew - xk
            xk = xnew
         END DO
      ENDIF

      !print *,'EXINV total number of iterations ',IX
      if (IX.GE.IXMAX) THEN
!         print *, 'Warning: EXINV did not converge. Cap=',Cap
!         print *, 'The last step was:  ', h, ' value=,',xk,' ak=',ak
      endif
      return
      END FUNCTION FUNCA
      END FUNCTION EXINV
      END MODULE FIMOD