C     Version  July 2007
C
C  The MREG module provide 3 programs.
C
C    1) MREG
C    2) RIND
C    3) FI - normal CDF
C 
C MREG and RIND are explained in the following:
C
C
C CALL MREG(F,R,B,DB,AA,BB,A,DA,VDER,M,N,NIT,INFR)
C
C  F    = expectation
C  R    = Covariance R(i+(j-1)*N) = Cov( Delta(T(i)), Delta(T(j)),    length RDIM (in)          
C  B    = Covariance  B(i) = Cov(Delta(T(i)), XN),  B(N+1)=Var(XN)    length NMAX (in)
C  DB   = Covariance DB(i) = Cov(Delta(T(i)), Y0), DB(N+1)=Cov(XN,Y0) length NMAX (in)
C  AA   = Regression matrix coefficients    size MMAX x MMAX
C  BB   = Regression vector coefficients    length MMax + 1
C  A    = Slepian model coefficients,  length (MMax + 1) * NMAX                                                                                                                                                       
C  DA   = Slepian model coefficients,  length  MMax + 1 
C  VDER = variance of Y0, Var(Y0)
C  M    = Number of regressors       ( 0 < M < MMAX)
C  N    = dimension of the problem   (     N < NMAX)
C  NIT  = 0,1,2..., maximum # of iterations/integrations done by quadrature 
C         to calculate the indicator function 
C  INFR = 1 means all input are the same as in the previous call except BB, A and DA
C         0 indicate new input
C
C  The program MREG computes the following problem:
C
C  We consider a process X(I)=X(T(I)) at the grid of  N  points  T(1),...,T(N),
C
C          X(I) = -A(I) + Z*A(I+N) + Sum Xj*A(I+(j+1)*N) + Delta(I), j=1,...,M-1
C
C  where the sum disappears if M=1. We assume that Z,Xj are independent
C  standard Rayleigh, Gaussian distributed rv. and independent of the zero
C  mean Gaussian residual process, with covariance structure given in  R,
C
C          R(i+(j-1)N) = Cov (Delta(T(i)), Delta(T(j))).
C
C  Additionally we have a zero mean Gaussian variable XN, 
C  independent of Z,Xj with  covariance structure defined by 
C  B(i)= Cov (Delta(T(i)),XN), i=1,...,N,  B(N+1)=Var(XN). 
C  Furthermore XN and Z,Xj satisfies the following equation system
C
C     (BB + (XN,0,...,0)^T =  AA*(Z,X1,...,Xm-1)^T       (***)
C
C  where AA is (M,M) matrix, BB is M-vector. We rewrite this equation, by
C  introducing a variable Xm=XN/SQRT(Var(XN)) and construct   new matrix AA1
c  by adding the column (SQRT(Var(XN)),0,...,0) and the row with only zeros.
C  The equations (***) writtes
C
C     (BB,0)^T  =  AA1*(Z,X1,...,Xm-1,Xm)^T       (****)
C
C  where AA1 is (M+1,M+1) matrix, We assume that the rank of AA1 is M,
C  otherwise the density is singular and we give a output F=0.CC
C
C  Let Y0 be a zero-mean Gaussian variable independent of Z,Xj 
C  with  covariance structure defined by 
C  DB(i)= Cov (Delta(T(i)),Y0), i=1,...,N,  DB(N+1)=Cov(XN,Y0), Var(Y0)=VDER. 
C  Let Y be defined by
C
C     Y=-DA(1) + Z*DA(2) + Sum Xj*DA(2+j) +Y0, j=1,...,M-1.
C
C  The MREG program computes:
C
C  F = E[ Y^+ *1{ HH<X(I)<0 for all I, I=1,...,N}|Z,X1,...,Xm-1 solves (***)]
C      *f_{Z,X1,....,Xm-1}(***).
C
C  In the simplest case NIT=0 we define (Delta(1),...,Delta(N),XN)=0.0d0 
C
C  We renormalize  vectors AN and DA, the covariance fkn R, DB
C  and VDER. Then by renormalization we choose the Gaussian variable X such
C  that F is written in the form
C
C  F = E[(D0(1)+X*D0(2)+Y1)^+*(PC+X*PD)^+*1{HH <A0(I)+X*B0(I)+Delta1(I)<0}]
C
C  Observe, PC+X*PD>0 defines integration region for X.
C  In the simplest case NIT=0 we define (Delta(1),...,Delta(N),Y1)=0.0d0 
C  For NIT=1 only (Delta(1),...,Delta(N))=0, i.e. we have to compute
C  a one dimensional integral. Finally by conditioning on X the problem is
C  put in the format of RIND-problem.
C
C  INF indicates whether one
C  has already called the subroutine before and ONLY! inputs BB, DA or A
C  was changed.
C
C  Observe the limitations are :  N<NMAX = 201, 0<M< 5 = MMAX.
C
C
c                  3-IX-93
C
C  CALL  RIND(XIND,R,BU,DBUN,DB,SQ,VDER,NIT,N,INFR)
C
C  XIND = expectation/density (inout)
C  R    = Covariances Delta(ti), Delta(tj), size RDIM x 1 (in)
C         R(i+(j-1)*N) = Cov( Delta(T(i)), Delta(T(j)) 
C  BU   = expectation of y(t), i.e., E(y(t)), size NMAX x 1  (in)
C  DBUN = expectation of Y, i.e., E(Y) 
C  DB   = Covariances Delta(T(i)), Y,    size NMAX x 1  (in)
C         DB(i) = Cov(Delta(T(i)), Y)
C  SQ   = standard deviations of Delta(T(i)), size NMAX x 1
C         SQ(I) = SQRT (R(I+(I-1)*N)) 
C  VDER = variance of Y, Var(Y)
C  NIT  = 0,1,2..., maximum # of iterations/integrations done by quadrature 
C         to calculate the indicator function 
C  N    = dimension of the problem
C  INFR = 1 means R, DB and SQ are the same as in the previous
C         0 indicate new R, DB and SQ.  
C
C  The program RIND computes the following problem:
C
C  Let the process  y(t)=BU(t)+Delta(t), where  Delta(t)  is zero mean
C  Gaussian process and BU(t) be expectation of  y(t). Consider the process  x
C  at the grid T(1),...,T(N), N=0,...,50, (N=0 means the empty grid T).
C  Let  y(I) = BU(T(I)) + Delta(T(I)).  Observe we do not assume that the
C  points T(I) are ordered or from, e.g. T(I) are in R^K.
C
C  The covariance fkn of Delta at the points T(1),...,T(N), are given in
C  the vector R;  Cov( Delta(T(i)), Delta(T(j)) = R(i+(j-1)*N),
C  furter E[y(T(i))] = BU(i). Hence the dimension of R must be N*N=2500.
C  The vector SQ(1), ...., SQ(N) contains standard deviations of the residual
C  Delta(T(I)), e.g. SQ(I) = SQRT (R(I+(I-1)*N)). However the small values
C  of SQ could be corrupted by nummerical errors especially if the covariance
C  structure was computed using the FFT algorithm. IF R(I+(I-1)*N)<EPS
C  SQ(I)=0.    and is used as an indicator that one is not allowed to
C  condition on Delta(T(I)). Further when one have conditioned on the point
C  T(I) the variance is put to zero.
C
C  Consider in addition to y(t) a Gaussian variable  Y, E[Y]=DBUN, Var(Y)=VDER,
C  DB(I)=Cov(Delta(T(I)),Y).
C
C  *** XIND - is the result;  XIND=E[Y^+1{ HH<y(I)<0 for all I, I=1,...,N}] ***
C
C In the speccial case by choosing  DB(I)=0, VDER=0 and DBUN=1, IAC=0,1,
C (if IAC=0 VDER can take any positive value), the output XIND is equal to
C XIND=Prob( HH < y(I) < 0 for all I, I=1,...,N).
C
C
C  Some control variables:
C  INFR=1 means that both R, DB and SQ are the same as in the previous
C  call of RIND subroutine, INFR=0 indicates the new R, DB and SQ. The history
C  of the conditioning is saved in the vectors INF(5), INFO(5): INF(1), ...,
C  INF(5) are the times one has conditioned in the subroutines RINDT1,...,
C  RINDT5, respectively. After conditioning INFO(i)=INF(i). Now if INF=INFO
C  then the conditioning tree has not be changed and one not need to compute
C  the conditonal covariances. This is really time saving trick. We are assume
C  that the program saves all the time the matrices at the same fysical
C  location, e.g. the values of variables are saved during the execution.
C  This has all the time be checked when new compilator will be used.
C
C  The variable ISQ marks which type of conditioning will be used ISQ=0
C  means random time where the probability is minimum, ISQ=1 is the time
C  where the variance of the residual process is minimal(ISQ=1 is faster).
C
C  NIT, IAC are  described in CROSSPACK paper, EPS0 is the accuracy constant
C  used in choosing the number of nodes in numerical integrations
C  (XX1, H1 vectors). The nodes and weights and other parameters are
C  read in the subroutine INITINTEG from files Z.DAT, H.DAT and ACCUR.DAT.
C
C
C    NIT=0, IAC=1 then one uses RIND0 - subroutine, all other cases
C    goes through RIND1, ...,RIND5. NIT=0, here means explicite formula
C    approximation for XIND=E[Y^+1{ HH<BU(I)<0 for all I, I=1,...,N}], where
C    BU(I) is deterministic function.
C
C    NIT=1, leads tp call RIND1, IAC=0 is also explicit form approximation,
C    while IAC=1 leads to maximum one dimensional integral.
C    .......
C    NIT=5, leads tp call RIND5, IAC is maximally 4-dimensional integral,
C    while IAC=1 leads to maximum 5 dimensional integral.
C
!
! Revised pab August 2007
! - replaced a call to SVDCMP with DSVDC derived from Lapack
! Revised pab July 2007
!
! - fixed some bugs 
! - reimplemented as module mregmodule
! - moved the functions/subroutines in twog.f into rindmod and renamed it to MREG. -> mreg and rind publicly available
! - All commonblocks are replaced with a corresponding module

! References
! Rychlik, I and Lindgren, G (1993)
! "CROSSREG - A Technique for First Passage and Wave Density Analysis"
! Probability in the Engineering and Informational Sciences, Vol 7, pp 125--148
!
! Lindgren, G and Rychlik, I (1991)
! "Slepian Models and Regression Approximations in Crossing  and xtreme value Theory",
! International Statistical Review, Vol 59, 2, pp 195--225


      MODULE SIZEMOD
      IMPLICIT NONE
      INTEGER, PARAMETER :: MMAX = 6, NMAX = 201
      INTEGER, PARAMETER :: RDIM = NMAX*NMAX 
      END MODULE SIZEMOD
      
      MODULE EPSMOD
      IMPLICIT NONE             
                      ! Constants determining accuracy of integration
                      !-----------------------------------------------
                      !if the conditional variance are less than: 
C      DOUBLE PRECISION :: EPS2=1.d-4    !- EPS2, the variable is 
                                        !  considered deterministic 
      DOUBLE PRECISION :: EPS = 1.d-2   ! SQRT(EPS2)
C      DOUBLE PRECISION :: XCEPS2=1.d-16 ! if Var(Xc) is less return NaN
      DOUBLE PRECISION :: EPSS = 5.d-5  ! accuracy of Indicator 
C      DOUBLE PRECISION :: CEPSS=0.99995d0 ! accuracy of Indicator 
      DOUBLE PRECISION :: EPS0 = 5.d-5 ! used in GAUSSLE1 to implicitly 
                                       ! determ. # nodes  
   
C      DOUBLE PRECISION :: fxcEpss=1.d-20 ! if less do not compute E(...|Xc)
C      DOUBLE PRECISION :: xCutOff=5.d0  ! upper/lower truncation limit of the 
                                       ! normal CDF 
C      DOUBLE PRECISION :: FxCutOff  = 0.99999942669686d0 
C      DOUBLE PRECISION :: CFxCutOff = 5.733031438470704d-7       ! 1-FxCutOff, 
      
      END MODULE EPSMOD

      MODULE RINTMOD
      DOUBLE PRECISION, save :: C = 4.5d0
      DOUBLE PRECISION, save :: FC = 0.999993204653751d0 
C     COMMON /RINT/   C,FC
      END MODULE RINTMOD

      MODULE TBRMOD
      USE SIZEMOD
      IMPLICIT NONE  
      DOUBLE PRECISION, DIMENSION(NMAX) :: HH
      END MODULE TBRMOD

      MODULE EXPACCMOD
      DOUBLE PRECISION,PARAMETER:: PMAX = 40.0d0
C     COMMON /EXPACC/ PMAX
      END MODULE EXPACCMOD

      MODULE INFCMOD
      IMPLICIT NONE  
      INTEGER, save :: ISQ = 0, IAC=1
      INTEGER, DIMENSION(10) :: INF,INFO
C     DOUBLE PRECISION, DIMENSION(10):: 
C      COMMON /INFC/   ISQ,INF,INFO
      END MODULE INFCMOD
      MODULE CHECKMOD
      IMPLICIT NONE
C      III01,III11,... - variables,counts how many times one calls
C      subroutine RIND0,RIND1,..., III*1 are also modified in the
C      subroutines RIND*. This gives us statistics over the complexity of
C      numerical calculations.
      INTEGER :: III01,III11,III21,III31,III41,III51
      INTEGER :: III61,III71,III81,III91,III101
      INTEGER :: III0
      END MODULE CHECKMOD


      MODULE QUADRMOD
      IMPLICIT NONE         ! Quadratures available: Legendre
      INTEGER            :: I

C      BLOCK DATA inithermite

      INTEGER, PARAMETER :: NNW = 13
      INTEGER, DIMENSION(25) :: NN
      REAL*8  Z(126),H(126)
      
      
C      COMMON /QUADR/  Z,H,NN,NNW
c      COMMON /EXPACC/ PMAX
C      COMMON /RINT/   C,FC
      
C      DATA NNW /13/
      DATA (NN(I),I=1,NNW)/2,3,4,5,6,7,8,9,10,12,16,20,24/
C      DATA PMAX/40./
C      DATA C/4.5/
      DATA (H(I),I=1,61)/1.0d0,1.0d0,0.555555555555556d0,
     * 0.888888888888889d0,
     * 0.555555555555556d0,0.347854845137454d0,0.652145154862546d0,
     * 0.652145154862546d0,0.347854845137454d0,0.236926885056189d0,
     * 0.478628670499366d0,0.568888888888889d0,0.478628670499366d0,
     * 0.236926885056189d0,0.171324492379170d0,0.360761573048139d0,
     * 0.467913934572691d0,0.467913934572691d0,0.360761573048139d0,
     * 0.171324492379170d0,0.129484966168870d0,0.279705391489277d0,
     * 0.381830050505119d0,0.417959183673469d0,0.381830050505119d0,
     * 0.279705391489277d0,0.129484966168870d0,0.101228536290376d0,
     * 0.222381034453374d0,0.313706645877887d0,0.362683783378362d0,
     * 0.362683783378362d0,0.313706645877887d0,0.222381034453374d0,
     * 0.101228536290376d0,0.081274388361574d0,0.180648160694857d0,
     * 0.260610696402935d0,0.312347077040003d0,0.330239355001260d0,
     * 0.312347077040003d0,0.260610696402935d0,0.180648160694857d0,
     * 0.081274388361574d0,0.066671344308688d0,0.149451349150581d0,
     * 0.219086362515982d0,0.269266719309996d0,0.295524224714753d0,
     * 0.295524224714753d0,0.269266719309996d0,0.219086362515982d0,
     * 0.149451349150581d0,0.066671344308688d0,0.047175336386512d0,
     * 0.106939325995318d0,0.160078328543346d0,0.203167426723066d0,
     * 0.233492536538355d0,0.249147048513403d0,0.249147048513403d0/
      DATA (H(I),I=62,101)/0.233492536538355d0,0.203167426723066d0,
     * 0.160078328543346d0,0.106939325995318d0,
     * 0.047175336386512d0,0.027152459411754094852d0,
     * 0.062253523938647892863d0,0.095158511682492784810d0,
     * 0.124628971255533872052d0,0.149595988816576732081d0,
     * 0.169156519395002538189d0,0.182603415044923588867d0, 
     * 0.189450610455068496285d0,0.189450610455068496285d0,
     * 0.182603415044923588867d0,0.169156519395002538189d0,
     * 0.149595988816576732081d0,0.124628971255533872052d0,
     * 0.095158511682492784810d0,0.062253523938647892863d0,
     * 0.027152459411754094852d0,0.017614007139152118312d0,
     * 0.040601429800386941331d0,0.062672048334109063570d0,
     * 0.083276741576704748725d0,0.101930119817240435037d0,
     * 0.118194531961518417312d0,0.131688638449176626898d0,
     * 0.142096109318382051329d0,0.149172986472603746788d0,
     * 0.152753387130725850698d0,0.152753387130725850698d0,
     * 0.149172986472603746788d0,0.142096109318382051329d0,
     * 0.131688638449176626898d0,0.118194531961518417312d0,
     * 0.101930119817240435037d0,0.083276741576704748725d0,
     * 0.062672048334109063570d0,0.040601429800386941331d0/
      DATA (H(I),I=102,126)/0.017614007139152118312d0,
     * 0.012341229799987199547d0, 0.028531388628933663181d0,
     * 0.044277438817419806169d0, 0.059298584915436780746d0,
     * 0.073346481411080305734d0, 0.086190161531953275917d0,
     * 0.097618652104113888270d0, 0.107444270115965634783d0,
     * 0.115505668053725601353d0, 0.121670472927803391204d0,
     * 0.125837456346828296121d0, 0.127938195346752156974d0,
     * 0.127938195346752156974d0, 0.125837456346828296121d0,
     * 0.121670472927803391204d0, 0.115505668053725601353d0,
     * 0.107444270115965634783d0, 0.097618652104113888270d0,
     * 0.086190161531953275917d0, 0.073346481411080305734d0,
     * 0.059298584915436780746d0, 0.044277438817419806169d0,
     * 0.028531388628933663181d0, 0.012341229799987199547d0/ 
        
      DATA (Z(I),I=1,58)/-0.577350269189626d0,0.577350269189626d0,
     * -0.774596669241483d0,0.0d0,
     *  0.774596669241483d0, -0.861136311594053d0, -0.339981043584856d0,
     *  0.339981043584856d0,  0.861136311594053d0, -0.906179845938664d0,
     * -0.538469310105683d0,0.0d0,
     *  0.538469310105683d0,  0.906179845938664d0, -0.932469514203152d0,
     * -0.661209386466265d0, -0.238619186083197d0,  0.238619186083197d0,
     *  0.661209386466265d0,  0.932469514203152d0, -0.949107912342759d0,
     * -0.741531185599394d0, -0.405845151377397d0, 0.0d0,
     *  0.405845151377397d0,  0.741531185599394d0,  0.949107912342759d0,
     * -0.960289856497536d0, -0.796666477413627d0, -0.525532409916329d0,
     * -0.183434642495650d0,  0.183434642495650d0,  0.525532409916329d0,
     *  0.796666477413627d0,  0.960289856497536d0, -0.968160239507626d0,
     * -0.836031107326636d0, -0.613371432700590d0, -0.324253423403809d0,
     *  0.0d0, 
     *  0.324253423403809d0,  0.613371432700590d0,  0.836031107326636d0,
     *  0.968160239507626d0, -0.973906528517172d0, -0.865063366688985d0,
     * -0.679409568299024d0, -0.433395394129247d0, -0.148874338981631d0,
     *  0.148874338981631d0,  0.433395394129247d0,  0.679409568299024d0,
     *  0.865063366688985d0,  0.973906528517172d0, -0.981560634246719d0,
     * -0.904117256370475d0, -0.769902674194305d0, -0.587317954286617d0/
      DATA (Z(I),I=59,99)/-0.367831498198180d0, -0.125233408511469d0,
     *      0.125233408511469d0,  0.367831498198180d0,  
     *      0.587317954286617d0,  0.769902674194305d0,
     *      0.904117256370475d0,  0.981560634246719d0,
     *     -0.989400934991649932596d0,
     *     -0.944575023073232576078d0, -0.865631202387831743880d0,
     *     -0.755404408355003033895d0, -0.617876244402643748447d0,
     *     -0.458016777657227386342d0, -0.281603550779258913230d0,
     *     -0.095012509837637440185d0,  0.095012509837637440185d0,
     *      0.281603550779258913230d0,  0.458016777657227386342d0,
     *      0.617876244402643748447d0,  0.755404408355003033895d0,
     *      0.865631202387831743880d0,  0.944575023073232576078d0,
     *      0.989400934991649932596d0, -0.993128599185094924786d0,
     *     -0.963971927277913791268d0, -0.912234428251325905868d0,
     *     -0.839116971822218823395d0, -0.746331906460150792614d0,
     *     -0.636053680726515025453d0, -0.510867001950827098004d0,
     *     -0.373706088715419560673d0, -0.227785851141645078080d0,
     *     -0.076526521133497333755d0,  0.076526521133497333755d0,
     *      0.227785851141645078080d0,  0.373706088715419560673d0,
     *      0.510867001950827098004d0,  0.636053680726515025453d0,
     *      0.746331906460150792614d0,  0.839116971822218823395d0/
      DATA (Z(I),I=100,126)/0.912234428251325905868d0,
     *      0.963971927277913791268d0,  0.993128599185094924786d0,
     *     -0.995187219997021360180d0, -0.974728555971309498198d0,
     *     -0.938274552002732758524d0, -0.886415527004401034213d0,
     *     -0.820001985973902921954d0, -0.740124191578554364244d0,
     *     -0.648093651936975569252d0, -0.545421471388839535658d0,
     *     -0.433793507626045138487d0, -0.315042679696163374387d0,
     *     -0.191118867473616309159d0, -0.064056892862605626085d0,
     *      0.064056892862605626085d0,  0.191118867473616309159d0,
     *      0.315042679696163374387d0,  0.433793507626045138487d0,
     *      0.545421471388839535658d0,  0.648093651936975569252d0,
     *      0.740124191578554364244d0,  0.820001985973902921954d0,
     *      0.886415527004401034213d0,  0.938274552002732758524d0,
     *      0.974728555971309498198d0,  0.995187219997021360180d0/
      END MODULE QUADRMOD
      

C
      MODULE MREGMOD
      IMPLICIT NONE
      PRIVATE
      PUBLIC :: RIND, MREG, FI

      INTERFACE RIND
      MODULE PROCEDURE RIND
      END INTERFACE

      INTERFACE MREG
      MODULE PROCEDURE MREG
      END INTERFACE

      INTERFACE FI
      MODULE PROCEDURE FI
      END INTERFACE

      INTERFACE C1_C2
      MODULE PROCEDURE C1_C2
      END INTERFACE

      INTERFACE GAUSS1
      MODULE PROCEDURE GAUSS1
      END INTERFACE

      INTERFACE GAUSINT
      MODULE PROCEDURE GAUSINT
      END INTERFACE

      INTERFACE PYTHAG
      MODULE PROCEDURE PYTHAG
      END INTERFACE


      CONTAINS


      SUBROUTINE RIND(XIND,R,BU,DBUN,DB,SQ,VDER,NIT,N,INFR)
      USE TBRMOD
      USE INFCMOD
      USE CHECKMOD
      USE EPSMOD
      USE SIZEMOD
      IMPLICIT NONE
      REAL*8, intent(inout) :: XIND,DBUN,VDER
      REAL*8, DIMENSION(RDIM), intent(inout) :: R
      REAL*8, DIMENSION(NMAX), intent(inout) :: BU,DB, SQ
      INTEGER, intent(in) :: NIT,N,INFR
      REAL*8 SDER
      INTEGER, save :: NNIT
      INTEGER I,III
C     DIMENSION R(1),BU(1),SQ(1),DB(1)
C     DIMENSION INF(10),INFO(10),HH(101)
C     COMMON /TBR/ HH
C     COMMON /INFC/   ISQ,INF,INFO
C     COMMON /CHECK1/ III01,III11,III21,III31,III41,III51
C    *,III61,III71,III81,III91,III101
C     COMMON /EPS/    EPS,EPSS,CEPSS
C
C      III01,III11,... - variables,counts how many times one calls
C      subroutine RIND0,RIND1,..., III*1 are also modified in the
C      subroutines RIND*. This gives us statistics over the complexity of
C      numerical calculations.
C
      XIND=0.0d0
      IF (N.lt.1) go to 99
      
      IF (INFR.EQ.0) THEN
         NNIT=MIN(NIT,N)
         if (NNIT.gt.10) NNIT=10
         DO I=1,10
            INF(I)=0
            INFO(I)=0
         enddo
         III=0
         DO I=1,N
          IF (SQ(I).GT.EPS) then
            III=1
            else
             IF(BU(I).GT.0.0d0) THEN
             RETURN
             END IF
             IF(BU(I).LT.HH(I)) THEN
             RETURN
             END IF
          END IF
         enddo
        END IF
        IF (III.eq.0) go to 99

!      GO TO (10,20,30,40,50,60,70,80,90,100) NNIT
      SELECT CASE (NNIT)
      CASE (1)
       CALL RIND1(XIND,R,BU,DBUN,DB,SQ,VDER,N)
       iii11=iii11+1
      CASE(2)
       CALL RIND2(XIND,R,BU,DBUN,DB,SQ,VDER,N)
       iii21=iii21+1
      CASE(3)
       CALL RIND3(XIND,R,BU,DBUN,DB,SQ,VDER,N)
       iii31=iii31+1
      CASE(4)
       CALL RIND4(XIND,R,BU,DBUN,DB,SQ,VDER,N)
       iii41=iii41+1
      CASE(5)
       CALL RIND5(XIND,R,BU,DBUN,DB,SQ,VDER,N)
       iii51=iii51+1
      CASE(6)
       CALL RIND6(XIND,R,BU,DBUN,DB,SQ,VDER,N)
       iii61=iii61+1
      CASE(7)
       CALL RIND7(XIND,R,BU,DBUN,DB,SQ,VDER,N)
       iii71=iii71+1
      CASE(8)
       CALL RIND8(XIND,R,BU,DBUN,DB,SQ,VDER,N)
       iii81=iii81+1
      CASE (9)
       CALL RIND9(XIND,R,BU,DBUN,DB,SQ,VDER,N)
      iii91=iii91+1
      CASE (10)
       CALL RIND10(XIND,R,BU,DBUN,DB,SQ,VDER,N)
       iii101=iii101+1
      CASE DEFAULT
       CALL RIND0(XIND,BU,DBUN,VDER,N)
       iii01=iii01+1
      END SELECT
      RETURN

 99   continue   
      SDER=0.0d0
      IF(VDER.GT.EPS) SDER=SQRT(VDER)
      XIND=PMEAN(DBUN,SDER)
      return
      END SUBROUTINE RIND

      SUBROUTINE RIND0(XIND,BU,DBUN,VDER,N)
      USE TBRMOD
      USE EPSMOD
      USE SIZEMOD
      IMPLICIT NONE
      INTEGER, intent(in) :: N
      REAL*8, intent(inout) ::  XIND,DBUN,VDER
      REAL*8, DIMENSION(NMAX), intent(inout) :: BU
      REAL*8  SDER
      INTEGER I
!      DIMENSION BU(NMAX)
C      DIMENSION HH(101)
C      COMMON /EPS/ EPS,EPSS,CEPSS
C      COMMON /TBR/ HH
      
      IF (N.LT.1) GO TO 20
      XIND=0.0d0
      IF(DBUN.LT.0.0d0) THEN
      RETURN
      END IF
      DO I=1,N
      IF(BU(I).GT.0.0d0) THEN
      RETURN
      END IF
      IF(BU(I).LT.HH(I)) THEN
      RETURN
      END IF
      enddo
20    CONTINUE
      SDER=0.0d0
      IF(VDER.GT.EPS) SDER=SQRT(VDER)
      XIND=PMEAN(DBUN,SDER)
      RETURN
      END SUBROUTINE RIND0

      SUBROUTINE RIND1(XIND,R,BU,DBUN,DB,SQ,VDER,N)
      USE SIZEMOD
      USE TBRMOD
      USE INFCMOD
      USE CHECKMOD
      USE EPSMOD
      USE RINTMOD
      IMPLICIT NONE
      REAL*8, intent(inout) :: XIND,DBUN,VDER
      REAL*8, DIMENSION(RDIM), intent(inout) :: R
      REAL*8, DIMENSION(NMAX), intent(inout) :: BU,DB,SQ
      INTEGER, intent(in) :: N
      REAL*8, DIMENSION(NMAX), save :: B1,SQ1
      REAL*8, DIMENSION(24) :: XX1, H1
      REAL*8 XMI,XMA,DER,SDER
      REAL*8, save :: DB1N,SDER1,VDER1, SS0
      REAL*8 XFF,XF,X,XH, SQ0, HHB 
      INTEGER I,III,J,II0,N1
!     INTEGER IAC,N
!      real*8 XIND,R,BU,DBUN,DB,SQ,VDER
!     REAL*8 XX1,H1, B1,SQ1,XMI,XMA, SDER,DB1N , DER, SDER1
!     REAL*8 XFF,X, XH, SQ0, HHB, SS0, VDER1, XF
!     INTEGER I,J,III,II0,N1
C      DIMENSION R(1),BU(1),SQ(1),DB(1),B1(NMAX),SQ1(NMAX)
!     DIMENSION R(RDIM),BU(NMAX),SQ(NMAX),DB(NMAX) ,B1(NMAX),SQ1(NMAX)
!      DIMENSION XX1(24),H1(24)
C      DIMENSION HH(101),INF(10),INFO(10)
      
C      COMMON /EPS/    EPS,EPSS,CEPSS
C      COMMON /RINT/   C,FC
C      COMMON /TBR/    HH
C      COMMON /INFC/   ISQ,INF,INFO
C      COMMON /CHECK1/ III01,III11,III21,III31,III41,III51
C     *,III61,III71,III81,III91,III101
      
c      print *,'Topp of R1:',sq(1),sq(2),sq(3)
      XIND=0.0d0

C Choice of the time for conditioning, two methods
C
C ISQ=1; INF(1)=II0 is the point where the  SQ(I) obtaines its maximum, SQ0
C is the maximal st. deviation of the residual.
C
C ISQ=0; INF(1) is the time point when the probability P(hh<BU(I)+Delta(I)<0)
C obtains its minimum which is denoted by XFF.
C
      XFF=1.0d0
      IF (N.LT.1) GO TO 11
C
C  If N<1 means empty grid we can not condition on any point and hence GO TO 11
C  then XFF=1.0d0 and the XIND will be approximated by E[Y^+]=PMEAN(DBUN,SDER).
C  Obs. E[Y]=DBUN and Var(Y)=VDER.
C
      SQ0=0.0d0
      DO I=1,N
      SQ1(i)=0.0d0
      IF (SQ(I).LE.eps) GO TO 1
      HHB=HH(I)-BU(I)
C
C Obs. SQ(I)<=EPS idicates that the point is not good for conditioning.
C There can be two reasons for it: 1 Variance of the residual is too small
C or the point was already used before.
C
      if (ISQ.GT.1) then
      SS0  =R(I+(I-1)*N)
      DB1N =DB(I)
      VDER1=DB1N*DB1N/SS0
      IF (VDER1.gt.SQ0) Then
      SQ0=VDER1
      II0=I
      END IF
      ELSE
      IF (SQ(I).GT.SQ0) THEN
      SQ0=SQ(I)
      II0=I
      END IF
      END IF

      X=-BU(I)/SQ(I)
      XH=HHB/SQ(I)
      XF=FI(X)-FI(XH)
      IF(XF.LT.XFF) THEN
      INF(1)=I
      XFF=XF
      END IF
1     CONTINUE
      enddo
11    CONTINUE
C
C   If the minimum probability XFF is close to 0 !!! then the indicator 1{}
C   can be bounded by EPSS leading to the approximation of XIND=0 and RETURN.
C
      IF(XFF.LT.EPSS) RETURN
C
C  We are stoping because we assume that for all sample pathes I(0,t)=1
C
C   If the minimum probability XFF is close to one then the indicator 1{}
C   can be bounded by 1 leading to the approximation of XIND by E[Y^+]=
C   PMEAN(DBUN,SDER), if IAC=1 or with E[Y]^+=MAX(0,DBUN).
C   This reduces the order of integrations.
C
      IF (XFF.GT.0.999d0*FC) THEN
      SDER=0.
      IF(VDER.GT.EPS) SDER=SQRT(VDER)
      XIND=MAX(DBUN,0.0d0)
      IF(IAC.LT.1) RETURN
      XIND=PMEAN(DBUN,SDER)
      RETURN
      END IF
C
C  We are conditioning on the point T(INF(1)). If ISQ=1 INF(1)=ii0.
C  Obviously, X(INF(1))=BU(INF(1))+Delta(INF(1)), where SQ0 is a standard
C  deviation of Delta(IN(1)), hence  if X(INF(1))>XMA or X(INF(1))<XM1
C  1{}=0. Hence the values of X(INF(1)) are truncated to the interval [xmi,xma].

      IF(ISQ.EQ.1) INF(1)=II0
      SQ0=SQ(INF(1))
      XMA=-BU(INF(1))/SQ0
      XMI=XMA+HH(INF(1))/SQ0
      XMI=MAX(-C,XMI)
      XMA=MIN(C,XMA)
      IF (XMI.GT.XMA) XMA=-C
C
C   Now we are checking whether INF(I)=INFO(I), I=1,..,5, what indicates
C   that all conditional covariances and variancec are unchanged since the
C   last call of this subroutine (R,SQ,DB) as well as
C
      III=0
      DO I=1,10
      III=III + ABS(INF(I)-INFO(I))
      enddo
      IF (III.EQ.0) GO TO 99
      DO I=1,N
       B1(I) = R(I+(INF(1)-1)*N)
      enddo
      SS0 = B1(INF(1))
      DB1N = DB(INF(1))
      INFO(1) = INF(1)
      VDER1 = VDER-DB1N*DB1N/SS0
      SDER1 = 0.0d0
      IF(VDER1.GT.EPS) SDER1=SQRT(VDER1)
      DB1N = DB1N/SQ0
      DO I=1,N
       B1(I) = B1(I)/SQ0
      enddo
99    CONTINUE
C
C  Here conditioning is done
C
      CALL C1_C2(XMI,XMA,BU,B1,DBUN,DB1N,0.0d0,SQ1,N)
      IF(FI(XMA)-FI(XMI).LT.EPSS) RETURN
C *******************************************************
C
C  In this special case RIND1 if IAC=0 one can explicitly compute XIND
C  and stop.
      IF (IAC.LT.1) THEN
      XIND=GAUSINT(XMI,XMA,DBUN,DB1N,1.0d0,0.0d0)
      RETURN
      END IF
      CALL GAUSS1(N1,H1,XX1,XMI,XMA,EPS0)
c      print *,XMI,XMA,EPS0,N1
c      write(11,*) XMI,XMA,EPS0,N1
      XIND=0.0d0
      DO J=1,N1
      DER=DBUN+XX1(J)*DB1N
      XIND=XIND+PMEAN(DER,SDER1)*H1(J)
      III01=III01+1
c      IF (N.eq.15) then
c      print *,'der,dbun,db1n,sder1',der,dbun,db1n,sder1
c      write(11,*) der,dbun,db1n,sder1
c      end if
10    CONTINUE
      enddo
c      IF (N.eq.15) then
c      do 999 iii=1,N
c      print *,iii,sq(Iii)
c999   continue
c      write(11,*) XIND,INF(1),INF(2),INF(3),inf(4)
c      pause
c      end if
      RETURN
      END SUBROUTINE RIND1

      SUBROUTINE RIND2(XIND,R,BU,DBUN,DB,SQ,VDER,N)
      USE SIZEMOD
      USE TBRMOD
      USE INFCMOD
      USE CHECKMOD
      USE EPSMOD
      USE RINTMOD
      IMPLICIT NONE
      REAL*8, intent(inout) :: XIND,DBUN,VDER
      REAL*8, DIMENSION(RDIM), intent(inout) :: R
      REAL*8, DIMENSION(NMAX), intent(inout) :: BU,DB,SQ
      INTEGER, intent(in) :: N
      REAL*8, DIMENSION(RDIM), save :: R1
      REAL*8, DIMENSION(NMAX), save :: BU1,DB1,B1,SQ1
      REAL*8, DIMENSION(24) :: XX1, H1
      REAL*8, save :: DB1N,SDER1,VDER1, SS0
      REAL*8 XMI,XMA,DER,SDER
      REAL*8 XIND1,XFF,XF,X,XH, XR1, SQ0 
      INTEGER I,III,J,II0,N1

!     INTEGER IAC,N
!     INTEGER I,J, II0,III, N1
!     REAL*8 XIND,R,BU,DBUN,DB,SQ,VDER
!      real*8 XX1,H1,R1,B1,DB1,BU1,SQ1
!     REAL*8 SS0, DB1N, VDER1, XR1
!     REAL*8 XFF,X, XH, SQ0, SDER, XF, XMI, XMA
!     REAL*8 SDER1, DER,XIND1
C      DIMENSION R(1),BU(1),SQ(1),DB(1) 
!     DIMENSION R(RDIM),BU(NMAX),SQ(NMAX),DB(NMAX)
!      DIMENSION R1(RDIM),B1(NMAX),DB1(NMAX),BU1(NMAX),SQ1(NMAX)
!      DIMENSION XX1(24),H1(24)
C      DIMENSION HH(101),INF(10),INFO(10)
      
C      COMMON /EPS/    EPS,EPSS,CEPSS
C      COMMON /RINT/   C,FC
C      COMMON /TBR/    HH
C      COMMON /INFC/   ISQ,INF,INFO
c      COMMON/CHECK/III0,III1,III2,III3,III4
C      COMMON /CHECK1/ III01,III11,III21,III31,III41,III51
C     *,III61,III71,III81,III91,III101
      
c      PRINT *,'Topp of 2:',sq(1),sq(2),sq(3)
      XIND=0.0d0
      XFF=1.0d0
      IF (N.LT.1) GO TO 11
      SQ0=0.0d0
      DO I=1,N
         IF (SQ(I).LE.EPS) THEN
            SQ1(I)=SQ(I)
            GO TO 1
         END IF
c      CSQ=C*SQ(I)
c      IF (BU(I).LT.-CSQ.AND.BU(I).GT.HH(I)+CSQ) GO TO 1
c      IF (BU(I).LT.-CSQ.AND.BU(I).GT.HH(I)+CSQ) SQ1(I)=EPS1
C      IF (BU(I).GT. CSQ.OR.BU(I).LT.HH(I)-CSQ) RETURN
         IF (SQ(I).GT.SQ0) THEN
            SQ0=SQ(I)
            II0=I
         END IF
         X=-BU(I)/SQ(I)
         XH=X+HH(I)/SQ(I)
         XF=FI(X)-FI(XH)
c      IF(XF.GT.CEPSS) SQ1(I)=EPS1
         IF (XF.LT.XFF) THEN
            INF(2)=I
            XFF=XF
         END IF
 1    CONTINUE
      ENDDO
11    CONTINUE
      IF(XFF.LT.EPSS) RETURN
      IF (XFF.GT.0.9999d0*FC) THEN
      SDER=0.0d0
      IF (VDER.GT.EPS) SDER=SQRT(VDER)
      XIND=MAX(DBUN,0.0d0)
      IF(IAC.LT.1) RETURN
      XIND=PMEAN(DBUN,SDER)
      RETURN
      END IF
      IF(ISQ.EQ.1) INF(2)=II0
      SQ0=SQ(INF(2))
      XMA=-BU(INF(2))/SQ0
      XMI=XMA+HH(INF(2))/SQ0
      XMI=MAX(-C,XMI)
      XMA=MIN(C,XMA)
      IF (XMI.GT.XMA) XMA=-C
C *************************************************
C
C  We are conditioning on the point T(INF(2)) and write new model
C  BU(I)+X*B1(I)+Delta1(I), I=1,...,N  (Obs. we do not use I=1,N)
C  SQ1(I) is standard deviation of Delta1 DBUN=BU'(N), DB1N=B1'(N) and X is
C  N(0,1) independent of Delta1, SDER1 is standard deviation of Delta1'(N).
C
      III=0
      DO I=2,10
      III=III + ABS(INF(I)-INFO(I))
      ENDDO
      IF (III.EQ.0) GO TO 99
      CALL M_COND(R1,B1,DB1,R,DB,INF(2),N)
C      III1=III1+1
      SS0=B1(INF(2))
      INFO(2)=INF(2)
      DB1N=DB(INF(2))
      VDER1=VDER-DB1N*DB1N/SS0
      SDER1=0.0d0
      IF (VDER1.GT.EPS) SDER1=SQRT(VDER1)
C      SQ0=SQRT(SS0)
      DB1N=DB1N/SQ0
      SQ1(INF(2))=0.0d0
      DO I=1,N
C
C IF (SQ1(I).EQ.EPS1) GO TO 3 - the .EQ. can not be raplaced with .LT. without
C some general changes in the strategy of SQ and SQ1 values. More exactly
C SQ can not be changed in this subroutine when for some  I  we would like to
C put SQ1(I)=EPS1 in the first loop. This SQ1 should not be changed here and
C thus we have GO TO 3 statement. Observe that the other SQ1 values are
C
      IF (SQ(I).LE.EPS) GO TO 3
C      IF (SQ1(I).EQ.EPS1.OR.SQ(I).LE.EPS1) GO TO 3
      XR1=R1(I+(I-1)*N)
c      IF(XR1.LT.0.0d0) CALL ERROR(I,N,-1)
      SQ1(I)=0.0d0
      IF (XR1.GT.EPS) SQ1(I)=SQRT(XR1)
3     CONTINUE
      ENDDO
      DO I=1,N
5     B1(I)=B1(I)/SQ0
      ENDDO
99    CONTINUE
C
C ***********************************************************
C
C We shall condition on the values of X, XMI<X<XMA, but for some
C X values XIND will be zero leading to reduced accuracy. Hence we try
C to exclude them and narrow the interval [XMI,XMA]
C
c      PRINT *,'2:**  ',XMI,XMA
      CALL C1_C2(XMI,XMA,BU,B1,DBUN,DB1N,SDER1,SQ1,N)
c      PRINT *,'2:****',XMI,XMA
      IF(FI(XMA)-FI(XMI).LT.EPSS)  THEN
c         print *, 'Leaving R2: Exit 4', XIND,XIND1,VDER1
         RETURN
      ENDIF
      CALL GAUSS1(N1,H1,XX1,XMI,XMA,EPS0)
      DO J=1,N1
        DO  I=1,N
20       BU1(I)=BU(I)+XX1(J)*B1(I)
        ENDDO
      DER=DBUN+XX1(J)*DB1N
c      print *,'R2: before calling R1: (SQ1):',sq1(1),sq1(2),sq1(3)
      CALL RIND1(XIND1,R1,BU1,DER,DB1,SQ1,VDER1,N)
      III11=III11+1
10    XIND=XIND+XIND1*H1(J)
      ENDDO
c      print *, 'Leaving R2: Exit 5', XIND,XIND1,VDER1
      RETURN
      END SUBROUTINE RIND2

      SUBROUTINE RIND3(XIND,R,BU,DBUN,DB,SQ,VDER,N)
      USE SIZEMOD
      USE TBRMOD
      USE INFCMOD
      USE CHECKMOD
      USE EPSMOD
      USE RINTMOD
      IMPLICIT NONE
      REAL*8, intent(inout) :: XIND,DBUN,VDER
      REAL*8, DIMENSION(RDIM), intent(inout) :: R
      REAL*8, DIMENSION(NMAX), intent(inout) :: BU,DB,SQ
      INTEGER, intent(in) :: N
      REAL*8, DIMENSION(RDIM), save :: R1
      REAL*8, DIMENSION(NMAX), save :: BU1,DB1,B1,SQ1
      REAL*8, DIMENSION(24) :: XX1, H1
      REAL*8, save :: DB1N,SDER1,VDER1, SS0
      REAL*8 XMI,XMA,DER,SDER
      REAL*8 XIND1,XFF,XF,X,XH, XR1, SQ0 
      INTEGER I,III,J,II0,N1

!     INTEGER IAC,N
!     INTEGER I,J, II0,III, N1
!     REAL*8 XIND,R,BU,DBUN,DB,SQ,VDER
!      real*8 XX1,H1
!     REAL*8 R1,B1,DB1,BU1,SQ1
!     REAL*8 XFF, SQ0, X, XH,XF, SDER, XMI,XMA
!     REAL*8 SS0, DB1N, VDER1, XR1
!     REAL*8 SDER1, DER,XIND1
C      DIMENSION R(1),BU(1),SQ(1),DB(1)
!     DIMENSION R(RDIM),BU(NMAX),SQ(NMAX),DB(NMAX)
!      DIMENSION R1(RDIM),B1(NMAX),DB1(NMAX),BU1(NMAX),SQ1(NMAX)
!      DIMENSION XX1(24),H1(24)
C     DIMENSION HH(101),INF(10),INFO(10)
      
C      COMMON /EPS/    EPS,EPSS,CEPSS
C      COMMON /RINT/   C,FC
C      COMMON /TBR/    HH
C      COMMON /INFC/   ISQ,INF,INFO
C      COMMON/CHECK/III0,III1,III2,III3,III4
C      COMMON /CHECK1/ III01,III11,III21,III31,III41,III51
C     *,III61,III71,III81,III91,III101
      
c      PRINT *,'Topp of 3:',sq(1),sq(2),sq(3)
      XIND=0.0d0
      XFF=1.0d0
      IF (N.LT.1) GO TO 11
      SQ0=0.0d0
      DO I=1,N
         IF (SQ(I).LE.EPS) THEN
            SQ1(I)=SQ(I)
            GO TO 1
         END IF
         IF (SQ(I).GT.SQ0) THEN
            SQ0=SQ(I)
            II0=I
         END IF
         X=-BU(I)/SQ(I)
         XH=X+HH(I)/SQ(I)
         XF=FI(X)-FI(XH)
         IF(XF.LT.XFF) THEN
            INF(3)=I
            XFF=XF
         END IF
 1    CONTINUE
      ENDDO
 11   CONTINUE
      IF (XFF.LT.EPSS) RETURN
      IF (XFF.GT.0.9999d0*FC) THEN
         SDER=0.0d0
         IF(VDER.GT.EPS) SDER=SQRT(VDER)
         XIND=MAX(DBUN,0.0d0)
         IF (IAC.LT.1) RETURN
         XIND=PMEAN(DBUN,SDER)
         RETURN
      END IF
      IF(ISQ.EQ.1) INF(3)=II0
      SQ0=SQ(INF(3))
      XMA=-BU(INF(3))/SQ0
      XMI=XMA+HH(INF(3))/SQ0
      XMI=MAX(-C,XMI)
      XMA=MIN(C,XMA)
      IF (XMI.GT.XMA) XMA=-C
      III=0
      DO I=3,10
         III=III + ABS(INF(I)-INFO(I))
      ENDDO
      IF (III.EQ.0) GO TO 99
      CALL M_COND(R1,B1,DB1,R,DB,INF(3),N)
      SS0=B1(INF(3))
      DB1N=DB(INF(3))
      VDER1=VDER-DB1N*DB1N/SS0
      SDER1=0.0d0
      IF (VDER1.GT.EPS) SDER1=SQRT(VDER1)
      INFO(3)=INF(3)
      DB1N=DB1N/SQ0
      SQ1(INF(3))=0.0d0
      DO I=1,N
         IF (SQ(I).LE.EPS) GO TO 3
         XR1=R1(I+(I-1)*N)
         SQ1(I)=0.0d0
         IF (XR1.GT.EPS) SQ1(I)=SQRT(XR1)
3     CONTINUE
      ENDDO
      DO I=1,N
 5    B1(I)=B1(I)/SQ0
      ENDDO
99    CONTINUE
c      PRINT *,'3:**  ',XMI,XMA
      CALL C1_C2(XMI,XMA,BU,B1,DBUN,DB1N,SDER1,SQ1,N)
c      PRINT *,'3:****',XMI,XMA,EPSS
      IF (FI(XMA)-FI(XMI).LT.EPSS) RETURN
      CALL GAUSS1(N1,H1,XX1,XMI,XMA,EPS0)
      DO J=1,N1
         DO I=1,N
20        BU1(I)=BU(I)+XX1(J)*B1(I)
         ENDDO
      DER=DBUN+XX1(J)*DB1N
c      print *,'R3: before calling R2: (SQ1):',sq1(1),sq1(2),sq1(3)
      CALL RIND2(XIND1,R1,BU1,DER,DB1,SQ1,VDER1,N)
      III21=III21+1
10    XIND=XIND+XIND1*H1(J)
      ENDDO
      RETURN
      END SUBROUTINE RIND3

      SUBROUTINE RIND4(XIND,R,BU,DBUN,DB,SQ,VDER,N)
      USE SIZEMOD
      USE TBRMOD
      USE INFCMOD
      USE CHECKMOD
      USE EPSMOD
      USE RINTMOD
      IMPLICIT NONE
      REAL*8, intent(inout) :: XIND,DBUN,VDER
      REAL*8, DIMENSION(RDIM), intent(inout) :: R
      REAL*8, DIMENSION(NMAX), intent(inout) :: BU,DB,SQ
      INTEGER, intent(in) :: N
      REAL*8, DIMENSION(RDIM), save :: R1
      REAL*8, DIMENSION(NMAX), save :: BU1,DB1,B1,SQ1
      REAL*8, DIMENSION(24) :: XX1, H1
      REAL*8, save :: DB1N,SDER1,VDER1, SS0
      REAL*8 XMI,XMA,DER,SDER
      REAL*8 XIND1,XFF,XF,X,XH, XR1, SQ0 
      INTEGER I,III,J,II0,N1

!     INTEGER IAC,N
!     INTEGER I,J,II0,III, N1
!     REAL*8 XIND,R,BU,DBUN,DB,SQ,VDER
!      real*8 XX1,H1
!     REAL*8 R1,B1,DB1,BU1,SQ1
!     REAL*8 XFF, SQ0, X, XH,XF, SDER, XMI,XMA
!     REAL*8 SS0, DB1N, VDER1, XR1
!     REAL*8 SDER1, DER,XIND1
C      DIMENSION R(1),BU(1),SQ(1),DB(1)
!     DIMENSION R(RDIM),BU(NMAX),SQ(NMAX),DB(NMAX)
!      DIMENSION R1(RDIM),B1(NMAX),DB1(NMAX),BU1(NMAX),SQ1(NMAX)
!      DIMENSION XX1(24),H1(24)
C      DIMENSION HH(101),INF(10),INFO(10)
      
C      COMMON /EPS/    EPS,EPSS,CEPSS
C      COMMON /RINT/   C,FC
C      COMMON /TBR/    HH
C      COMMON /INFC/   ISQ,INF,INFO
C      COMMON/CHECK/III0,III1,III2,III3,III4
C      COMMON /CHECK1/ III01,III11,III21,III31,III41,III51
C     *,III61,III71,III81,III91,III101
      
c      PRINT *,'Topp of 4:',SQ(1),SQ(2),SQ(3)
      XIND=0.0d0
      XFF=1.0d0
      IF (N.LT.1) GO TO 11
      SQ0=0.0d0
      DO I=1,N
         IF (SQ(I).LE.EPS) THEN
            SQ1(I)=SQ(I)
            GO TO 1
         END IF
         IF (SQ(I).GT.SQ0) THEN
            SQ0=SQ(I)
            II0=I
         END IF
         X=-BU(I)/SQ(I)
         XH=X+HH(I)/SQ(I)
         XF=FI(X)-FI(XH)
         IF (XF.LT.XFF) THEN
            INF(4)=I
            XFF=XF
         END IF
 1    CONTINUE
      ENDDO
 11   CONTINUE
      IF (XFF.LT.EPSS) RETURN
      IF (XFF.GT.0.9999d0*FC) THEN
         SDER=0.0d0
         IF(VDER.GT.EPS) SDER=SQRT(VDER)
         XIND=MAX(DBUN,0.0d0)
         IF (IAC.LT.1) RETURN
         XIND=PMEAN(DBUN,SDER)
         RETURN
      END IF
      IF(ISQ.EQ.1) INF(4)=II0
      SQ0=SQ(INF(4))
      XMA=-BU(INF(4))/SQ0
      XMI=XMA+HH(INF(4))/SQ0
      XMI=MAX(-C,XMI)
      XMA=MIN(C,XMA)
      IF (XMI.GT.XMA) XMA=-C
      III=0
      DO I=4,10
         III=III + ABS(INF(I)-INFO(I))
      ENDDO
      IF (III.EQ.0) GO TO 99
      CALL M_COND(R1,B1,DB1,R,DB,INF(4),N)
      SS0=B1(INF(4))
      DB1N=DB(INF(4))
      VDER1=VDER-DB1N*DB1N/SS0
      SDER1=0.0d0
      IF (VDER1.GT.EPS) SDER1=SQRT(VDER1)
      INFO(4)=INF(4)
      DB1N=DB1N/SQ0
      SQ1(INF(4))=0.0d0
      DO I=1,N
         IF (SQ(I).LE.EPS) GO TO 3
         XR1=R1(I+(I-1)*N)
         SQ1(I)=0.0d0
         IF (XR1.GT.EPS) SQ1(I)=SQRT(XR1)
3     CONTINUE
      ENDDO
      DO I=1,N
5     B1(I)=B1(I)/SQ0
      ENDDO
99    CONTINUE
C      PRINT *,'**',XMI,XMA
      CALL C1_C2(XMI,XMA,BU,B1,DBUN,DB1N,SDER1,SQ1,N)
C      PRINT *,INF(4),XMI,XMA
      IF(FI(XMA)-FI(XMI).LT.EPSS) RETURN
      CALL GAUSS1(N1,H1,XX1,XMI,XMA,EPS0)
      DO  J=1,N1
      DO  I=1,N
20    BU1(I)=BU(I)+XX1(J)*B1(I)
      ENDDO
      DER=DBUN+XX1(J)*DB1N
      CALL RIND3(XIND1,R1,BU1,DER,DB1,SQ1,VDER1,N)
      III31=III31+1
10    XIND=XIND+XIND1*H1(J)
      ENDDO
      RETURN
      END SUBROUTINE RIND4

      SUBROUTINE RIND5(XIND,R,BU,DBUN,DB,SQ,VDER,N)
      USE SIZEMOD
      USE TBRMOD
      USE INFCMOD
      USE CHECKMOD
      USE EPSMOD
      USE RINTMOD
      IMPLICIT NONE
      REAL*8, intent(inout) :: XIND,DBUN,VDER
      REAL*8, DIMENSION(RDIM), intent(inout) :: R
      REAL*8, DIMENSION(NMAX), intent(inout) :: BU,DB,SQ
      INTEGER, intent(in) :: N
      REAL*8, DIMENSION(RDIM), save :: R1
      REAL*8, DIMENSION(NMAX), save :: BU1,DB1,B1,SQ1
      REAL*8, DIMENSION(24) :: XX1, H1
      REAL*8, save :: DB1N,SDER1,VDER1, SS0
      REAL*8 XMI,XMA,DER,SDER
      REAL*8 XIND1,XFF,XF,X,XH, XR1, SQ0 
      INTEGER I,III,J,II0,N1

!     INTEGER IAC,N
!     INTEGER I,J,III,II0,N1
!     REAL*8 XIND,R,BU,DBUN,DB,SQ,VDER
!      real*8 XX1,H1
!     REAL*8 R1,B1,DB1,BU1,SQ1
!     REAL*8 XFF, SQ0, X, XH,XF, SDER, XMI,XMA
!     REAL*8 SS0, DB1N, VDER1, XR1
!     REAL*8 SDER1, DER,XIND1
C     DIMENSION R(1),BU(1),SQ(1),DB(1)
!     DIMENSION R(RDIM),BU(NMAX),SQ(NMAX),DB(NMAX)
!      DIMENSION R1(RDIM),B1(NMAX),DB1(NMAX),BU1(NMAX),SQ1(NMAX)
      
!      DIMENSION XX1(24),H1(24)
C      DIMENSION INF(10),INFO(10),HH(101)

C      COMMON /EPS/    EPS,EPSS,CEPSS
C      COMMON /RINT/   C,FC
C      COMMON /TBR/    HH
C      COMMON /INFC/   ISQ,INF,INFO
C      COMMON/CHECK/III0,III1,III2,III3,III4
C      COMMON /CHECK1/ III01,III11,III21,III31,III41,III51
C     *,III61,III71,III81,III91,III101
      
      XIND=0.0d0
      XFF=1.0d0
      IF (N.LT.1) GO TO 11
      SQ0=0.0d0
      DO I=1,N
         IF (SQ(I).LE.EPS) THEN
            SQ1(I)=SQ(I)
            GO TO 1
         END IF
         IF (SQ(I).GT.SQ0) THEN
            SQ0=SQ(I)
            II0=I
         END IF
         X=-BU(I)/SQ(I)
         XH=X+HH(I)/SQ(I)
         XF=FI(X)-FI(XH)
         IF (XF.LT.XFF) THEN
            INF(5)=I
            XFF=XF
         END IF
 1    CONTINUE
      ENDDO
 11   CONTINUE
      IF (XFF.LT.EPSS) RETURN
      IF (XFF.GT.0.9999d0*FC) THEN
         SDER=0.d0
         IF(VDER.GT.EPS) SDER=SQRT(VDER)
         XIND=MAX(DBUN,0.0d0)
         IF (IAC.LT.1) RETURN
         XIND=PMEAN(DBUN,SDER)
         RETURN
      END IF
      IF(ISQ.EQ.1) INF(5)=II0
      SQ0=SQ(INF(5))
      XMA=-BU(INF(5))/SQ0
      XMI=XMA+HH(INF(5))/SQ0
      XMI=MAX(-C,XMI)
      XMA=MIN(C,XMA)
      IF (XMI.GT.XMA) XMA=-C
      III=0
      DO I=5,10
         III=III + ABS(INF(I)-INFO(I))
      ENDDO
      IF (III.EQ.0) GO TO 99
      CALL M_COND(R1,B1,DB1,R,DB,INF(5),N)
      SS0=B1(INF(5))
      DB1N=DB(INF(5))
      VDER1=VDER-DB1N*DB1N/SS0
      SDER1=0.0d0
      IF (VDER1.GT.EPS) SDER1=SQRT(VDER1)
      INFO(5)=INF(5)
      DB1N=DB1N/SQ0
      SQ1(INF(5))=0.0d0
      DO I=1,N
         IF (SQ(I).LE.EPS) GO TO 3
         XR1=R1(I+(I-1)*N)
         SQ1(I)=0.0d0
         IF (XR1.GT.EPS) SQ1(I)=SQRT(XR1)
3     CONTINUE
      ENDDO
      DO I=1,N
5     B1(I)=B1(I)/SQ0
      ENDDO
99    CONTINUE
      CALL C1_C2(XMI,XMA,BU,B1,DBUN,DB1N,SDER1,SQ1,N)
      IF(FI(XMA)-FI(XMI).LT.EPSS) RETURN
      CALL GAUSS1(N1,H1,XX1,XMI,XMA,EPS0)
      DO J=1,N1
       DO I=1,N
        BU1(I)=BU(I)+XX1(J)*B1(I)
       ENDDO
       DER=DBUN+XX1(J)*DB1N
       CALL RIND4(XIND1,R1,BU1,DER,DB1,SQ1,VDER1,N)
       III41=III41+1
       XIND=XIND+XIND1*H1(J)
      ENDDO
      RETURN
      END SUBROUTINE RIND5
C
      SUBROUTINE RIND6(XIND,R,BU,DBUN,DB,SQ,VDER,N)
      USE SIZEMOD
      USE TBRMOD
      USE INFCMOD
      USE CHECKMOD
      USE EPSMOD
      USE RINTMOD
      IMPLICIT NONE
      REAL*8, intent(inout) :: XIND,DBUN,VDER
      REAL*8, DIMENSION(RDIM), intent(inout) :: R
      REAL*8, DIMENSION(NMAX), intent(inout) :: BU,DB,SQ
      INTEGER, intent(in) :: N
      REAL*8, DIMENSION(RDIM), save :: R1
      REAL*8, DIMENSION(NMAX), save :: BU1,DB1,B1,SQ1
      REAL*8, DIMENSION(24) :: XX1, H1
      REAL*8, save :: DB1N,SDER1,VDER1, SS0
      REAL*8 XMI,XMA,DER,SDER
      REAL*8 XIND1,XFF,XF,X,XH, XR1, SQ0 
      INTEGER I,III,J,II0,N1

!     INTEGER IAC,N
!     INTEGER I,J,III,II0,N1
!     REAL*8 XIND,R,BU,DBUN,DB,SQ,VDER
!      real*8 XX1,H1
!     REAL*8 R1,B1,DB1,BU1,SQ1
!     REAL*8 XFF, SQ0, X, XH,XF, SDER, XMI,XMA
!     REAL*8 SS0, DB1N, VDER1, XR1
!     REAL*8 SDER1, DER,XIND1
C      DIMENSION R(1),BU(1),SQ(1),DB(1)
!     DIMENSION R(RDIM),BU(NMAX),SQ(NMAX),DB(NMAX)
!      DIMENSION R1(RDIM),B1(NMAX),DB1(NMAX),BU1(NMAX),SQ1(NMAX)
!      DIMENSION XX1(24),H1(24)
C      DIMENSION HH(101),INF(10),INFO(10)
      
C      COMMON /EPS/    EPS,EPSS,CEPSS
C      COMMON /RINT/   C,FC
C      COMMON /TBR/    HH
C      COMMON /INFC/   ISQ,INF,INFO
C      COMMON /CHECK1/ III01,III11,III21,III31,III41,III51
C     *,III61,III71,III81,III91,III101
      
      XIND=0.0d0
      XFF=1.0d0
      IF (N.LT.1) GO TO 11
      SQ0=0.0d0
      DO I=1,N
         IF (SQ(I).LE.EPS) THEN
            SQ1(I)=SQ(I)
            GO TO 1
         END IF
c      CSQ=C*SQ(I)
c      IF (BU(I).LT.-CSQ.AND.BU(I).GT.HH(I)+CSQ) GO TO 1
c      IF (BU(I).LT.-CSQ.AND.BU(I).GT.HH(I)+CSQ) SQ1(I)=EPS1
C      IF (BU(I).GT. CSQ.OR.BU(I).LT.HH(I)-CSQ) RETURN
         IF (SQ(I).GT.SQ0) THEN
            SQ0=SQ(I)
            II0=I
         END IF
         X=-BU(I)/SQ(I)
         XH=X+HH(I)/SQ(I)
         XF=FI(X)-FI(XH)
c      IF(XF.GT.CEPSS) SQ1(I)=EPS1
         IF (XF.LT.XFF) THEN
            INF(6)=I
            XFF=XF
         END IF
 1    CONTINUE
      ENDDO
11    CONTINUE
      IF(XFF.LT.EPSS) RETURN
      IF (XFF.GT.0.9999d0*FC) THEN
      SDER=0.0d0
      IF (VDER.GT.EPS) SDER=SQRT(VDER)
      XIND=MAX(DBUN,0.0d0)
      IF(IAC.LT.1) RETURN
      XIND=PMEAN(DBUN,SDER)
      RETURN
      END IF
      IF(ISQ.EQ.1) INF(6)=II0
      SQ0=SQ(INF(6))
      XMA=-BU(INF(6))/SQ0
      XMI=XMA+HH(INF(6))/SQ0
      XMI=MAX(-C,XMI)
      XMA=MIN(C,XMA)
      IF (XMI.GT.XMA) XMA=-C
C *************************************************
C
C  We are conditioning on the point T(INF(2)) and write new model
C  BU(I)+X*B1(I)+Delta1(I), I=1,...,N  (Obs. we do not use I=1,N)
C  SQ1(I) is standard deviation of Delta1 DBUN=BU'(N), DB1N=B1'(N) and X is
C  N(0,1) independent of Delta1, SDER1 is standard deviation of Delta1'(N).
C
      III=0
      DO I=6,10
      III=III + ABS(INF(I)-INFO(I))
      ENDDO
      IF (III.EQ.0) GO TO 99
      CALL M_COND(R1,B1,DB1,R,DB,INF(6),N)
C      III1=III1+1
      SS0=B1(INF(6))
      INFO(6)=INF(6)
      DB1N=DB(INF(6))
      VDER1=VDER-DB1N*DB1N/SS0
      SDER1=0.0d0
      IF (VDER1.GT.EPS) SDER1=SQRT(VDER1)
C      SQ0=SQRT(SS0)
      DB1N=DB1N/SQ0
      SQ1(INF(6))=0.0d0
      DO I=1,N
C
C IF (SQ1(I).EQ.EPS1) GO TO 3 - the .EQ. can not be raplaced with .LT. without
C some general changes in the strategy of SQ and SQ1 values. More exactly
C SQ can not be changed in this subroutine when for some  I  we would like to
C put SQ1(I)=EPS1 in the first loop. This SQ1 should not be changed here and
C thus we have GO TO 3 statement. Observe that the other SQ1 values are
C
      IF (SQ(I).LE.EPS) GO TO 3
C      IF (SQ1(I).EQ.EPS1.OR.SQ(I).LE.EPS1) GO TO 3
      XR1=R1(I+(I-1)*N)
c      IF(XR1.LT.0.0d0) CALL ERROR(I,N,-1)
      SQ1(I)=0.0d0
      IF (XR1.GT.EPS) SQ1(I)=SQRT(XR1)
3     CONTINUE
      ENDDO
      DO I=1,N
5     B1(I)=B1(I)/SQ0
      ENDDO
99    CONTINUE
C
C ***********************************************************
C
C We shall condition on the values of X, XMI<X<XMA, but for some
C X values XIND will be zero leading to reduced accuracy. Hence we try
C to exclude them and narrow the interval [XMI,XMA]
C
      CALL C1_C2(XMI,XMA,BU,B1,DBUN,DB1N,SDER1,SQ1,N)
      IF(FI(XMA)-FI(XMI).LT.EPSS)  THEN
         RETURN
      ENDIF
      CALL GAUSS1(N1,H1,XX1,XMI,XMA,EPS0)
      DO J=1,N1
       DO I=1,N
20       BU1(I)=BU(I)+XX1(J)*B1(I)
       ENDDO
      DER=DBUN+XX1(J)*DB1N
      CALL RIND5(XIND1,R1,BU1,DER,DB1,SQ1,VDER1,N)
      III51=III51+1
10    XIND=XIND+XIND1*H1(J)
      ENDDO
      RETURN
      END SUBROUTINE RIND6

      SUBROUTINE RIND7(XIND,R,BU,DBUN,DB,SQ,VDER,N)
      USE SIZEMOD
      USE TBRMOD
      USE INFCMOD
      USE CHECKMOD
      USE EPSMOD
      USE RINTMOD
      IMPLICIT NONE
      REAL*8, intent(inout) :: XIND,DBUN,VDER
      REAL*8, DIMENSION(RDIM), intent(inout) :: R
      REAL*8, DIMENSION(NMAX), intent(inout) :: BU,DB,SQ
      INTEGER, intent(in) :: N
      REAL*8, DIMENSION(RDIM), save :: R1
      REAL*8, DIMENSION(NMAX), save :: BU1,DB1,B1,SQ1
      REAL*8, DIMENSION(24) :: XX1, H1
      REAL*8, save :: DB1N,SDER1,VDER1, SS0
      REAL*8 XMI,XMA,DER,SDER
      REAL*8 XIND1,XFF,XF,X,XH, XR1, SQ0 
      INTEGER I,III,J,II0,N1

!     INTEGER IAC,N
!     INTEGER I,J,III,II0,N1
!     REAL*8 XIND,R,BU,DBUN,DB,SQ,VDER
!      real*8 XX1,H1
!     REAL*8 R1,B1,DB1,BU1,SQ1
!     REAL*8 XFF, SQ0, X, XH,XF, SDER, XMI,XMA
!     REAL*8 SS0, DB1N, VDER1, XR1
!     REAL*8 SDER1, DER,XIND1
C      DIMENSION R(1),BU(1),SQ(1),DB(1)
!     DIMENSION R(RDIM),BU(NMAX),SQ(NMAX),DB(NMAX)
!      DIMENSION R1(RDIM),B1(NMAX),DB1(NMAX),BU1(NMAX),SQ1(NMAX)
!      DIMENSION XX1(24),H1(24)
C      DIMENSION HH(101),INF(10),INFO(10)
      
C      COMMON /EPS/    EPS,EPSS,CEPSS
C      COMMON /RINT/   C,FC
C      COMMON /TBR/    HH
C      COMMON /INFC/   ISQ,INF,INFO
C      COMMON /CHECK1/ III01,III11,III21,III31,III41,III51
C     *,III61,III71,III81,III91,III101
      
      XIND=0.0d0
      XFF=1.0d0
      IF (N.LT.1) GO TO 11
      SQ0=0.0d0
      DO I=1,N
         IF (SQ(I).LE.EPS) THEN
            SQ1(I)=SQ(I)
            GO TO 1
         END IF
c      CSQ=C*SQ(I)
c      IF (BU(I).LT.-CSQ.AND.BU(I).GT.HH(I)+CSQ) GO TO 1
c      IF (BU(I).LT.-CSQ.AND.BU(I).GT.HH(I)+CSQ) SQ1(I)=EPS1
C      IF (BU(I).GT. CSQ.OR.BU(I).LT.HH(I)-CSQ) RETURN
         IF (SQ(I).GT.SQ0) THEN
            SQ0=SQ(I)
            II0=I
         END IF
         X=-BU(I)/SQ(I)
         XH=X+HH(I)/SQ(I)
         XF=FI(X)-FI(XH)
c      IF(XF.GT.CEPSS) SQ1(I)=EPS1
         IF (XF.LT.XFF) THEN
            INF(7)=I
            XFF=XF
         END IF
 1    CONTINUE
      ENDDO
11    CONTINUE
      IF(XFF.LT.EPSS) RETURN
      IF (XFF.GT.0.9999d0*FC) THEN
      SDER=0.0d0
      IF (VDER.GT.EPS) SDER=SQRT(VDER)
      XIND=MAX(DBUN,0.0d0)
      IF(IAC.LT.1) RETURN
      XIND=PMEAN(DBUN,SDER)
      RETURN
      END IF
      IF(ISQ.EQ.1) INF(7)=II0
      SQ0=SQ(INF(7))
      XMA=-BU(INF(7))/SQ0
      XMI=XMA+HH(INF(7))/SQ0
      XMI=MAX(-C,XMI)
      XMA=MIN(C,XMA)
      IF (XMI.GT.XMA) XMA=-C
C *************************************************
C
C  We are conditioning on the point T(INF(2)) and write new model
C  BU(I)+X*B1(I)+Delta1(I), I=1,...,N  (Obs. we do not use I=1,N)
C  SQ1(I) is standard deviation of Delta1 DBUN=BU'(N), DB1N=B1'(N) and X is
C  N(0,1) independent of Delta1, SDER1 is standard deviation of Delta1'(N).
C
      III=0
      DO I=7,10
      III=III + ABS(INF(I)-INFO(I))
      ENDDO
      IF (III.EQ.0) GO TO 99
      CALL M_COND(R1,B1,DB1,R,DB,INF(7),N)
C      III1=III1+1
      SS0=B1(INF(7))
      INFO(7)=INF(7)
      DB1N=DB(INF(7))
      VDER1=VDER-DB1N*DB1N/SS0
      SDER1=0.0d0
      IF (VDER1.GT.EPS) SDER1=SQRT(VDER1)
C      SQ0=SQRT(SS0)
      DB1N=DB1N/SQ0
      SQ1(INF(7))=0.0d0
      DO I=1,N
C
C IF (SQ1(I).EQ.EPS1) GO TO 3 - the .EQ. can not be raplaced with .LT. without
C some general changes in the strategy of SQ and SQ1 values. More exactly
C SQ can not be changed in this subroutine when for some  I  we would like to
C put SQ1(I)=EPS1 in the first loop. This SQ1 should not be changed here and
C thus we have GO TO 3 statement. Observe that the other SQ1 values are
C
      IF (SQ(I).LE.EPS) GO TO 3
C      IF (SQ1(I).EQ.EPS1.OR.SQ(I).LE.EPS1) GO TO 3
      XR1=R1(I+(I-1)*N)
c      IF(XR1.LT.0.0d0) CALL ERROR(I,N,-1)
      SQ1(I)=0.0d0
      IF (XR1.GT.EPS) SQ1(I)=SQRT(XR1)
3     CONTINUE
      ENDDO
      DO I=1,N
       B1(I)=B1(I)/SQ0
      ENDDO
99    CONTINUE
C
C ***********************************************************
C
C We shall condition on the values of X, XMI<X<XMA, but for some
C X values XIND will be zero leading to reduced accuracy. Hence we try
C to exclude them and narrow the interval [XMI,XMA]
C
      CALL C1_C2(XMI,XMA,BU,B1,DBUN,DB1N,SDER1,SQ1,N)
      IF(FI(XMA)-FI(XMI).LT.EPSS)  THEN
         RETURN
      ENDIF
      CALL GAUSS1(N1,H1,XX1,XMI,XMA,EPS0)
      DO J=1,N1
       DO I=1,N
        BU1(I)=BU(I)+XX1(J)*B1(I)
       ENDDO
       DER=DBUN+XX1(J)*DB1N
       CALL RIND6(XIND1,R1,BU1,DER,DB1,SQ1,VDER1,N)
       III61=III61+1
       XIND=XIND+XIND1*H1(J)
      ENDDO
      RETURN
      END SUBROUTINE RIND7

      SUBROUTINE RIND8(XIND,R,BU,DBUN,DB,SQ,VDER,N)
      USE SIZEMOD
      USE TBRMOD
      USE INFCMOD
      USE CHECKMOD
      USE EPSMOD
      USE RINTMOD
          IMPLICIT NONE
      REAL*8, intent(inout) :: XIND,DBUN,VDER
      REAL*8, DIMENSION(RDIM), intent(inout) :: R
      REAL*8, DIMENSION(NMAX), intent(inout) :: BU,DB,SQ
      INTEGER, intent(in) :: N
      REAL*8, DIMENSION(RDIM), save :: R1
      REAL*8, DIMENSION(NMAX), save :: BU1,DB1,B1,SQ1
      REAL*8, DIMENSION(24) :: XX1, H1
      REAL*8, save :: DB1N,SDER1,VDER1, SS0
      REAL*8 XMI,XMA,DER,SDER
      REAL*8 XIND1,XFF,XF,X,XH, XR1, SQ0 
      INTEGER I,III,J,II0,N1

!     INTEGER IAC,N
!     INTEGER I,J,III,II0,N1
!     REAL*8 XIND,R,BU,DBUN,DB,SQ,VDER
!      real*8 XX1,H1
!     REAL*8 R1,B1,DB1,BU1,SQ1
!     REAL*8 XFF, SQ0, X, XH,XF, SDER, XMI,XMA
!     REAL*8 SS0, DB1N, VDER1, XR1
!     REAL*8 SDER1, DER,XIND1
C      DIMENSION R(1),BU(1),SQ(1),DB(1)
!     DIMENSION R(RDIM),BU(NMAX),SQ(NMAX),DB(NMAX)
!      DIMENSION R1(RDIM),B1(NMAX),DB1(NMAX),BU1(NMAX),SQ1(NMAX)
!      DIMENSION XX1(24),H1(24)
C      DIMENSION HH(101),INF(10),INFO(10)
      
C      COMMON /EPS/    EPS,EPSS,CEPSS
C      COMMON /RINT/   C,FC
C      COMMON /TBR/    HH
C      COMMON /INFC/   ISQ,INF,INFO
C      COMMON /CHECK1/ III01,III11,III21,III31,III41,III51
C     *,III61,III71,III81,III91,III101
      
      XIND=0.0d0
      XFF=1.0d0
      IF (N.LT.1) GO TO 11
      SQ0=0.0d0
      DO I=1,N
         IF (SQ(I).LE.EPS) THEN
            SQ1(I)=SQ(I)
            GO TO 1
         END IF
c      CSQ=C*SQ(I)
c      IF (BU(I).LT.-CSQ.AND.BU(I).GT.HH(I)+CSQ) GO TO 1
c      IF (BU(I).LT.-CSQ.AND.BU(I).GT.HH(I)+CSQ) SQ1(I)=EPS1
C      IF (BU(I).GT. CSQ.OR.BU(I).LT.HH(I)-CSQ) RETURN
         IF (SQ(I).GT.SQ0) THEN
            SQ0=SQ(I)
            II0=I
         END IF
         X=-BU(I)/SQ(I)
         XH=X+HH(I)/SQ(I)
         XF=FI(X)-FI(XH)
c      IF(XF.GT.CEPSS) SQ1(I)=EPS1
         IF (XF.LT.XFF) THEN
            INF(8)=I
            XFF=XF
         END IF
 1    CONTINUE
      ENDDO
11     CONTINUE
      IF(XFF.LT.EPSS) RETURN
      IF (XFF.GT.0.9999d0*FC) THEN
      SDER=0.0d0
      IF (VDER.GT.EPS) SDER=SQRT(VDER)
      XIND=MAX(DBUN,0.0d0)
      IF(IAC.LT.1) RETURN
      XIND=PMEAN(DBUN,SDER)
      RETURN
      END IF
      IF(ISQ.EQ.1) INF(8)=II0
      SQ0=SQ(INF(8))
      XMA=-BU(INF(8))/SQ0
      XMI=XMA+HH(INF(8))/SQ0
      XMI=MAX(-C,XMI)
      XMA=MIN(C,XMA)
      IF (XMI.GT.XMA) XMA=-C
C *************************************************
C
C  We are conditioning on the point T(INF(2)) and write new model
C  BU(I)+X*B1(I)+Delta1(I), I=1,...,N  (Obs. we do not use I=1,N)
C  SQ1(I) is standard deviation of Delta1 DBUN=BU'(N), DB1N=B1'(N) and X is
C  N(0,1) independent of Delta1, SDER1 is standard deviation of Delta1'(N).
C
      III=0
      DO I=8,10
      III=III + ABS(INF(I)-INFO(I))
      ENDDO
      IF (III.EQ.0) GO TO 99
      CALL M_COND(R1,B1,DB1,R,DB,INF(8),N)
C      III1=III1+1
      SS0=B1(INF(8))
      INFO(8)=INF(8)
      DB1N=DB(INF(8))
      VDER1=VDER-DB1N*DB1N/SS0
      SDER1=0.0d0
      IF (VDER1.GT.EPS) SDER1=SQRT(VDER1)
C      SQ0=SQRT(SS0)
      DB1N=DB1N/SQ0
      SQ1(INF(8))=0.0d0
      DO I=1,N
C
C IF (SQ1(I).EQ.EPS1) GO TO 3 - the .EQ. can not be raplaced with .LT. without
C some general changes in the strategy of SQ and SQ1 values. More exactly
C SQ can not be changed in this subroutine when for some  I  we would like to
C put SQ1(I)=EPS1 in the first loop. This SQ1 should not be changed here and
C thus we have GO TO 3 statement. Observe that the other SQ1 values are
C
      IF (SQ(I).LE.EPS) GO TO 3
C      IF (SQ1(I).EQ.EPS1.OR.SQ(I).LE.EPS1) GO TO 3
      XR1=R1(I+(I-1)*N)
c      IF(XR1.LT.0.0d0) CALL ERROR(I,N,-1)
      SQ1(I)=0.0d0
      IF (XR1.GT.EPS) SQ1(I)=SQRT(XR1)
3     CONTINUE
      ENDDO
      DO I=1,N
       B1(I)=B1(I)/SQ0
      ENDDO
99    CONTINUE
C
C ***********************************************************
C
C We shall condition on the values of X, XMI<X<XMA, but for some
C X values XIND will be zero leading to reduced accuracy. Hence we try
C to exclude them and narrow the interval [XMI,XMA]
C
      CALL C1_C2(XMI,XMA,BU,B1,DBUN,DB1N,SDER1,SQ1,N)
      IF(FI(XMA)-FI(XMI).LT.EPSS)  THEN
         RETURN
      ENDIF
      CALL GAUSS1(N1,H1,XX1,XMI,XMA,EPS0)
      DO J=1,N1
       DO I=1,N
        BU1(I)=BU(I)+XX1(J)*B1(I)
       ENDDO
       DER=DBUN+XX1(J)*DB1N
       CALL RIND7(XIND1,R1,BU1,DER,DB1,SQ1,VDER1,N)
       III71=III71+1
       XIND=XIND+XIND1*H1(J)
      ENDDO
      RETURN
      END SUBROUTINE RIND8

      SUBROUTINE RIND9(XIND,R,BU,DBUN,DB,SQ,VDER,N)
      USE SIZEMOD
      USE TBRMOD
      USE INFCMOD
      USE CHECKMOD
      USE EPSMOD
      USE RINTMOD
      IMPLICIT NONE
      REAL*8, intent(inout) :: XIND,DBUN,VDER
      REAL*8, DIMENSION(RDIM), intent(inout) :: R
      REAL*8, DIMENSION(NMAX), intent(inout) :: BU,DB,SQ
      INTEGER, intent(in) :: N
      REAL*8, DIMENSION(RDIM), save :: R1
      REAL*8, DIMENSION(NMAX), save :: BU1,DB1,B1,SQ1
      REAL*8, DIMENSION(24) :: XX1, H1
      REAL*8, save :: DB1N,SDER1,VDER1, SS0
      REAL*8 XMI,XMA,DER,SDER
      REAL*8 XIND1,XFF,XF,X,XH, XR1, SQ0 
      INTEGER I,III,J,II0,N1

!     INTEGER IAC,N
!     INTEGER I,J,III,II0,N1
!     REAL*8 XIND,R,BU,DBUN,DB,SQ,VDER
!      real*8 XX1,H1
!     REAL*8 R1,B1,DB1,BU1,SQ1
!     REAL*8 XFF, SQ0, X, XH,XF, SDER, XMI,XMA
!     REAL*8 SS0, DB1N, VDER1, XR1
!     REAL*8 SDER1, DER,XIND1
C      DIMENSION R(1),BU(1),SQ(1),DB(1)
!     DIMENSION R(RDIM),BU(NMAX),SQ(NMAX),DB(NMAX)
!      DIMENSION R1(RDIM),B1(NMAX),DB1(NMAX),BU1(NMAX),SQ1(NMAX)
!      DIMENSION XX1(24),H1(24)
C      DIMENSION HH(101),INF(10),INFO(10)
      
C      COMMON /EPS/    EPS,EPSS,CEPSS
C      COMMON /RINT/   C,FC
C      COMMON /TBR/    HH
C      COMMON /INFC/   ISQ,INF,INFO
C      COMMON /CHECK1/ III01,III11,III21,III31,III41,III51
C     *,III61,III71,III81,III91,III101
      
      XIND=0.0d0
      XFF=1.0d0
      IF (N.LT.1) GO TO 11
      SQ0=0.0d0
      DO I=1,N
         IF (SQ(I).LE.EPS) THEN
            SQ1(I)=SQ(I)
            GO TO 1
         END IF
c      CSQ=C*SQ(I)
c      IF (BU(I).LT.-CSQ.AND.BU(I).GT.HH(I)+CSQ) GO TO 1
c      IF (BU(I).LT.-CSQ.AND.BU(I).GT.HH(I)+CSQ) SQ1(I)=EPS1
C      IF (BU(I).GT. CSQ.OR.BU(I).LT.HH(I)-CSQ) RETURN
         IF (SQ(I).GT.SQ0) THEN
            SQ0=SQ(I)
            II0=I
         END IF
         X=-BU(I)/SQ(I)
         XH=X+HH(I)/SQ(I)
         XF=FI(X)-FI(XH)
c      IF(XF.GT.CEPSS) SQ1(I)=EPS1
         IF (XF.LT.XFF) THEN
            INF(9)=I
            XFF=XF
         END IF
1        CONTINUE 
      ENDDO
11    CONTINUE
      IF(XFF.LT.EPSS) RETURN
      IF (XFF.GT.0.9999d0*FC) THEN
      SDER=0.0d0
      IF (VDER.GT.EPS) SDER=SQRT(VDER)
      XIND=MAX(DBUN,0.0d0)
      IF(IAC.LT.1) RETURN
      XIND=PMEAN(DBUN,SDER)
      RETURN
      END IF
      IF(ISQ.EQ.1) INF(9)=II0
      SQ0=SQ(INF(9))
      XMA=-BU(INF(9))/SQ0
      XMI=XMA+HH(INF(9))/SQ0
      XMI=MAX(-C,XMI)
      XMA=MIN(C,XMA)
      IF (XMI.GT.XMA) XMA=-C
C *************************************************
C
C  We are conditioning on the point T(INF(2)) and write new model
C  BU(I)+X*B1(I)+Delta1(I), I=1,...,N  (Obs. we do not use I=1,N)
C  SQ1(I) is standard deviation of Delta1 DBUN=BU'(N), DB1N=B1'(N) and X is
C  N(0,1) independent of Delta1, SDER1 is standard deviation of Delta1'(N).
C
      III=0
      DO I=9,10
       III=III + ABS(INF(I)-INFO(I))
      ENDDO
      IF (III.EQ.0) GO TO 99
      CALL M_COND(R1,B1,DB1,R,DB,INF(9),N)
C      III1=III1+1
      SS0=B1(INF(9))
      INFO(9)=INF(9)
      DB1N=DB(INF(9))
      VDER1=VDER-DB1N*DB1N/SS0
      SDER1=0.0d0
      IF (VDER1.GT.EPS) SDER1=SQRT(VDER1)
C      SQ0=SQRT(SS0)
      DB1N=DB1N/SQ0
      SQ1(INF(9))=0.0d0
      DO I=1,N
C
C IF (SQ1(I).EQ.EPS1) GO TO 3 - the .EQ. can not be raplaced with .LT. without
C some general changes in the strategy of SQ and SQ1 values. More exactly
C SQ can not be changed in this subroutine when for some  I  we would like to
C put SQ1(I)=EPS1 in the first loop. This SQ1 should not be changed here and
C thus we have GO TO 3 statement. Observe that the other SQ1 values are
C
       IF (SQ(I)>EPS) THEN
        XR1=R1(I+(I-1)*N)
c      IF(XR1.LT.0.0d0) CALL ERROR(I,N,-1)
        SQ1(I)=0.0d0
        IF (XR1.GT.EPS) SQ1(I)=SQRT(XR1)
       ENDIF
      ENDDO
      DO I=1,N
       B1(I)=B1(I)/SQ0
      ENDDO
99    CONTINUE
C
C ***********************************************************
C
C We shall condition on the values of X, XMI<X<XMA, but for some
C X values XIND will be zero leading to reduced accuracy. Hence we try
C to exclude them and narrow the interval [XMI,XMA]
C
      CALL C1_C2(XMI,XMA,BU,B1,DBUN,DB1N,SDER1,SQ1,N)
      IF(FI(XMA)-FI(XMI).LT.EPSS)  THEN
         RETURN
      ENDIF
      CALL GAUSS1(N1,H1,XX1,XMI,XMA,EPS0)
      DO J=1,N1
       DO I=1,N
        BU1(I)=BU(I)+XX1(J)*B1(I)
       ENDDO
       DER=DBUN+XX1(J)*DB1N
       CALL RIND8(XIND1,R1,BU1,DER,DB1,SQ1,VDER1,N)
       III81=III81+1
       XIND=XIND+XIND1*H1(J)
      ENDDO
      RETURN
      END SUBROUTINE RIND9

      SUBROUTINE RIND10(XIND,R,BU,DBUN,DB,SQ,VDER,N)
      USE SIZEMOD
      USE TBRMOD
      USE INFCMOD
      USE CHECKMOD
      USE EPSMOD
      USE RINTMOD
      IMPLICIT NONE
      REAL*8, intent(inout) :: XIND,DBUN,VDER
      REAL*8, DIMENSION(RDIM), intent(inout) :: R
      REAL*8, DIMENSION(NMAX), intent(inout) :: BU,DB,SQ
      INTEGER, intent(in) :: N
      REAL*8, DIMENSION(RDIM), save :: R1
      REAL*8, DIMENSION(NMAX), save :: BU1,DB1,B1,SQ1
      REAL*8, DIMENSION(24) :: XX1, H1
      REAL*8, save :: DB1N,SDER1,VDER1, SS0
      REAL*8 XMI,XMA,DER,SDER
      REAL*8 XIND1,XFF,XF,X,XH, XR1, SQ0 
      INTEGER I,III,J,II0,N1

!     INTEGER IAC,N
!     INTEGER I,J,III,II0,N1
!     REAL*8 XIND,R,BU,DBUN,DB,SQ,VDER
!      real*8 XX1,H1
!     REAL*8 R1,B1,DB1,BU1,SQ1
!     REAL*8 XFF, SQ0, X, XH,XF, SDER, XMI,XMA
!     REAL*8 SS0, DB1N, VDER1, XR1
!     REAL*8 SDER1, DER,XIND1
C      DIMENSION R(1),BU(1),SQ(1),DB(1)
!     DIMENSION R(RDIM),BU(NMAX),SQ(NMAX),DB(NMAX)
!      DIMENSION R1(RDIM),B1(NMAX),DB1(NMAX),BU1(NMAX),SQ1(NMAX)
!      DIMENSION XX1(24),H1(24)
C      DIMENSION HH(101),INF(10),INFO(10)
     
C      COMMON /EPS/    EPS,EPSS,CEPSS
C      COMMON /RINT/   C,FC
C      COMMON /TBR/    HH
C      COMMON /INFC/   ISQ,INF,INFO
C      COMMON /CHECK1/ III01,III11,III21,III31,III41,III51
C     *,III61,III71,III81,III91,III101

      XIND=0.0d0
      XFF=1.0d0
      IF (N.LT.1) GO TO 11
      SQ0=0.0d0
      DO I=1,N
         IF (SQ(I).LE.EPS) THEN
            SQ1(I)=SQ(I)
            GO TO 1
         END IF
c      CSQ=C*SQ(I)
c      IF (BU(I).LT.-CSQ.AND.BU(I).GT.HH(I)+CSQ) GO TO 1
c      IF (BU(I).LT.-CSQ.AND.BU(I).GT.HH(I)+CSQ) SQ1(I)=EPS1
C      IF (BU(I).GT. CSQ.OR.BU(I).LT.HH(I)-CSQ) RETURN
         IF (SQ(I).GT.SQ0) THEN
            SQ0=SQ(I)
            II0=I
         END IF
         X=-BU(I)/SQ(I)
         XH=X+HH(I)/SQ(I)
         XF=FI(X)-FI(XH)
c      IF(XF.GT.CEPSS) SQ1(I)=EPS1
         IF (XF.LT.XFF) THEN
            INF(10)=I
            XFF=XF
         END IF
 1    CONTINUE
      ENDDO
11    CONTINUE
      IF(XFF.LT.EPSS) RETURN
      IF (XFF.GT.0.9999d0*FC) THEN
      SDER=0.0d0
      IF (VDER.GT.EPS) SDER=SQRT(VDER)
      XIND=MAX(DBUN,0.0d0)
      IF(IAC.LT.1) RETURN
      XIND=PMEAN(DBUN,SDER)
      RETURN
      END IF
      IF(ISQ.EQ.1) INF(10)=II0
      SQ0=SQ(INF(10))
      XMA=-BU(INF(10))/SQ0
      XMI=XMA+HH(INF(10))/SQ0
      XMI=MAX(-C,XMI)
      XMA=MIN(C,XMA)
      IF (XMI.GT.XMA) XMA=-C
C *************************************************
C
C  We are conditioning on the point T(INF(2)) and write new model
C  BU(I)+X*B1(I)+Delta1(I), I=1,...,N  (Obs. we do not use I=1,N)
C  SQ1(I) is standard deviation of Delta1 DBUN=BU'(N), DB1N=B1'(N) and X is
C  N(0,1) independent of Delta1, SDER1 is standard deviation of Delta1'(N).
C
      III=0
      DO I=10,10
      III=III + ABS(INF(I)-INFO(I))
      ENDDO
      IF (III.EQ.0) GO TO 99
      CALL M_COND(R1,B1,DB1,R,DB,INF(10),N)
C      III1=III1+1
      SS0=B1(INF(10))
      INFO(10)=INF(10)
      DB1N=DB(INF(10))
      VDER1=VDER-DB1N*DB1N/SS0
      SDER1=0.0d0
      IF (VDER1.GT.EPS) SDER1=SQRT(VDER1)
C      SQ0=SQRT(SS0)
      DB1N=DB1N/SQ0
      SQ1(INF(10))=0.0d0
      DO I=1,N
C
C IF (SQ1(I).EQ.EPS1) GO TO 3 - the .EQ. can not be raplaced with .LT. without
C some general changes in the strategy of SQ and SQ1 values. More exactly
C SQ can not be changed in this subroutine when for some  I  we would like to
C put SQ1(I)=EPS1 in the first loop. This SQ1 should not be changed here and
C thus we have GO TO 3 statement. Observe that the other SQ1 values are
cc
       IF (SQ(I)>EPS) THEN

        XR1=R1(I+(I-1)*N)
c      IF(XR1.LT.0.0d0) CALL ERROR(I,N,-1)
        SQ1(I)=0.0d0
        IF (XR1.GT.EPS) SQ1(I)=SQRT(XR1)
       ENDIF
      ENDDO
      DO I=1,N
       B1(I)=B1(I)/SQ0
      ENDDO
99    CONTINUE
C
C ***********************************************************
C
C We shall condition on the values of X, XMI<X<XMA, but for some
C X values XIND will be zero leading to reduced accuracy. Hence we try
C to exclude them and narrow the interval [XMI,XMA]
C
      CALL C1_C2(XMI,XMA,BU,B1,DBUN,DB1N,SDER1,SQ1,N)
      IF(FI(XMA)-FI(XMI).LT.EPSS)  THEN
        RETURN
      ENDIF
      CALL GAUSS1(N1,H1,XX1,XMI,XMA,EPS0)
      DO J=1,N1
       DO I=1,N
        BU1(I)=BU(I)+XX1(J)*B1(I)
       ENDDO
       DER=DBUN+XX1(J)*DB1N
       CALL RIND9(XIND1,R1,BU1,DER,DB1,SQ1,VDER1,N)
       III91=III91+1
       XIND=XIND+XIND1*H1(J)
      ENDDO
      RETURN
      END SUBROUTINE RIND10
      
      SUBROUTINE C1_C2(C1,C2,BU,B1,DBUN,DB1N,SDER,SQ,N)
C
C  We assume that the process  y  is of form y(I)=BU(I)+X*B1(I)+Delta(I),
C  I=1,...,N, SQ(I) is standard deviation of Delta(I), where X is  N(0,1)
C  independent of Delta.  Let Y = DBUN + DB1N*X + Z, where Z is zero-mean
C  Gaussian with standart independent of X (it can depend on Delta(I)) with
C  standart deviation SDER. Since we are truncating all Gaussian  variables to
C  the interval [-C,C], then if for any I
C
C  a) BU(I)+x*B1(I)-C*SQ(I)>0  or
C
C  b) BU(I)+x*B1(I)+C*SQ(I)<HH  then
C
C  XIND|X=x = E[Y^+1{ HH<y(I)<0 for all I, I=1,...,N}|X=x] = 0 !!!!!!!!!
C
C  Further, see discussion in comments to the subroutine PMEAN, by first upper-
C  bounding the indicator 1{} in XIND by 1, XIND|X=x = 0 if
C
C  c) DBUN+x*DB1N+4.5*SDER<0.0d0
C
C  Consequently, for increasing the accuracy (by excluding possible discon-
C  tinuouities) we shall exclude such X=x values for which XIND|X=x = 0.0d0
C  XIND=E([XIND|X]). Hence we assume that if C1<X<C2 any of the previous
C  conditions are satisfied.
C
C  OBSERVE!!, C1, C2 has to be set to upper bounds of possible values, e.g.
C  C1=-C, C2=C before calling C1_C2 subroutine.
C

C NOTE: Check that integration limits contained in TBRMOD is used correctly
C       If order of variables are changed then also the integration limits in HH should reflect this
      USE EPSMOD
      USE RINTMOD
      USE TBRMOD
      USE SIZEMOD
      IMPLICIT NONE
      REAL*8, DIMENSION(NMAX), intent(inout) :: BU,SQ,B1
      REAL*8, intent(inout) :: C1,C2
      REAL*8, intent(in) :: DBUN,DB1N,SDER
      INTEGER, intent(in) :: N
      REAL*8 CSQ, HHB, CC1, CC2, X
      INTEGER I   
C      DIMENSION BU(1),B1(1),SQ(1)
!     DIMENSION BU(NMAX),SQ(NMAX),B1(NMAX)
C      DIMENSION HH(101)
C      COMMON /EPS/EPS,EPSS,CEPSS
C      COMMON/RINT/C,FC
C      COMMON/TBR/HH
      DO I=1,N
      CSQ=C*SQ(I)
      HHB=HH(I)-BU(I)
C
C  If ABS(B1(I)) < EPS we can have overflow and hence we consider two cases
C  1) BU(I) is so large or small so we can surely assume that the probability
C     of staying between the barriers is 0, consequently C1=C2=0
C  2) we do not change the original limits.
C
      IF (ABS(B1(I)).LT. EPS) THEN
         IF (BU(I).GT.CSQ.OR.BU(I).LT.HH(I)-CSQ) THEN
            C1=0.0d0
            C2=0.0d0
            RETURN
         END IF
C
C  In other cases this part follows from the description of the problem.
C
       ELSE
         IF (B1(I).GT.EPS) THEN
            CC1=(HHB-CSQ)/B1(I)
            CC2=(-BU(I)+CSQ)/B1(I)
            IF (C1.LT.CC1) C1=CC1
            IF (C2.GT.CC2) C2=CC2
          ELSE
            CC2=(HHB-CSQ)/B1(I)
            CC1=(-BU(I)+CSQ)/B1(I)
            IF (C1.LT.CC1) C1=CC1
            IF (C2.GT.CC2) C2=CC2
          END IF
      END IF
      ENDDO
      X=-DBUN-4.5d0*SDER
      IF(DB1N.GT.EPS.AND.C1.LT.X/DB1N) C1=X/DB1N
      IF(DB1N.LT.-EPS.AND.C2.GT.X/DB1N) C2=X/DB1N
      if(abs(db1n).lt.eps.and.x.gt.0.0d0) then
            C1=0.0d0
            C2=0.0d0
            RETURN
         END IF

c
c  In the following three rows we are cutting C1, C2 to the interval [-C,C].
c  Obs. all tree lines are neccessary.
c
      C1=MAX(-C,C1)
      C2=MIN( C,C2)
      IF (C1.GT.C2) C2=-C
C      PRINT *,2,C1,C2
      RETURN
      END SUBROUTINE C1_C2

      REAL*8 FUNCTION GAUSINT(X1,X2,A,B,C,D)
C
C Let  X  be stardized Gaussian variable, i.e. X=N(0,1).
C The function calculate the followin integral E[I(X1<X<X2)(A+BX)(C+DX)],
C where I(X1<X<X2) is an indicator function of the set {X1<X<X2}.
C
      IMPLICIT NONE
      REAL*8, intent(in) :: X1,X2,A,B,C,D
      REAL*8 Y1,Y2,Y3
      REAL*8, PARAMETER:: SP = 0.398942280401433d0 
      IF(X1.GE.X2) THEN
      GAUSINT=0.0d0
      RETURN
      END IF
      Y1=(A*D+B*C+X1*B*D)*EXP(-0.5d0*X1*X1)
      Y2=(A*D+B*C+X2*B*D)*EXP(-0.5d0*X2*X2)
      Y3=(A*C+B*D)*(FI(X2)-FI(X1))
      GAUSINT=Y3+SP*(Y1-Y2)
      RETURN
      END FUNCTION GAUSINT


      SUBROUTINE GAUSS1(N,H1,XX1,XMI,XMA,EPS0) 
      USE CHECKMOD
      USE QUADRMOD
      IMPLICIT NONE
      INTEGER, intent(out) :: N
      REAL*8, DIMENSION(24), intent(out) :: H1,XX1
      REAL*8, intent(in) :: XMI,XMA,EPS0
      REAL*8, DIMENSION(24) ::  Z1
      REAL*8 SDOT, SDOT1, DIFF1
      INTEGER NNN, J, I1
!     DIMENSION Z1(24),XX1(1),H1(1)
C      DIMENSION Z(126),H(126)
C      DIMENSION NN(25)
C      COMMON/QUADR/ Z,H,NN,NNW
C      COMMON/CHECKQ/ III0
      REAL*8, parameter::  SP= 0.398942280401433d0
      IF (XMA.LT.XMI) THEN
      PRINT *,'Error XMIN>XMAX in GAUSS1 - stop!'
C     STOP
      END IF
      NNN=0
      DO I=1,NNW
      N=NN(I)
      DO J=1,N
       XX1(J)=0.5d0*(Z(NNN+J)*(XMA-XMI)+XMA+XMI)
       Z1(J)=XX1(J)*XX1(J)
       H1(J)=0.5d0*SP*(XMA-XMI)*H(NNN+J)*EXP(-0.5d0*Z1(J))
      ENDDO
      NNN=NNN+N
      SDOT=GAUSINT(XMI,XMA,0.0d0,1.0d0,0.0d0,1.0d0)
      SDOT1=0.d0
      DO I1=1,N
        SDOT1=SDOT1+Z1(I1)*H1(I1)
      ENDDO
      DIFF1=ABS(SDOT-SDOT1)
      IF(EPS0.LT.DIFF1) GO TO 10
      III0=III0+N
C      PRINT *,'N. of nodes',III0
      RETURN
10    CONTINUE
      ENDDO
      END SUBROUTINE GAUSS1

      
      SUBROUTINE M_COND(Syy_cd,Syyii,Syx_cd,Syy,Syx,ii,N)
C
C               INPUT:
C
C   ii     IS THE INDEX OF THE TIME ON WHICH WE ARE CONDITIONING.
C   N      number of variables in covariance matrix  Syy
C
C   Covariance matrix  Syy(I+(J-1)*N)=Cov(Yi,Yj) (is unchanged)
C   Covariance vector  Syx(I)=Cov(Yi,X) (is unchanged)
C
C              OUTPUT:
C
C   Covariance matrix Syy_cd(I+(J-1)*N)=Cov(Xi,Xj|Xii)
C   Covariance vector Syyii(I)=Cov(Xi,Xii)
C   Covariance vector Syx_cd(I)=Cov(Xi,Y|Xii)
C   Variance          Q1=Var(Xii)=Syyii(ii)
c   Obs. If   Q1<EPS there is no conditioning
C
C M_COND(R1,B1,DB1,R,DB,INF(10),N)
      USE EPSMOD
      USE SIZEMOD
      IMPLICIT NONE
      INTEGER, intent(in) :: II,N
      REAL*8, DIMENSION(RDIM), intent(inout) :: Syy_cd,Syy
      REAL*8, DIMENSION(NMAX), intent(inout)  :: Syyii,Syx_cd,Syx
      REAL*8 Q1
      INTEGER I,J
!     DIMENSION Syy_cd(RDIM),Syyii(NMAX),Syx_cd(NMAX),Syy(RDIM),Syx(NMAX)
C      DIMENSION Syy_cd(1),Syyii(1),Syx_cd(1),Syy(1),Syx(1)
C      COMMON /EPS/    EPS,EPSS,CEPSS
      IF (II.LE.0.OR.II.GT.N) THEN
       PRINT *,'The conditioning time in M_COND is out of range, stop!'
       STOP
      END IF
C
C   Q1=Var(Xii)=Syyii(ii)
C
      Q1=Syy(II+(II-1)*N)
      IF(Q1.LE.eps) then
      DO I=1,N
       Syyii(I)=0.0d0
      ENDDO
      Q1=1.0d0
      else
      DO I=1,N
       Syyii(I)=Syy(I+(II-1)*N)
      ENDDO
      end if
      DO I=1,N
       DO J=1,N
        Syy_cd(I+(J-1)*N)=Syy(I+(J-1)*N)-Syy(II+(J-1)*N)*Syyii(I)/Q1
       ENDDO
       Syx_cd(I)=Syx(I)-Syx(II)*Syyii(I)/Q1
      ENDDO
      RETURN
      END SUBROUTINE M_COND
      


      REAL*8 FUNCTION PMEAN(XX,SS)
C
C   PMEAN is the positive mean of a Gaussian variable with mean  XX  and
C   standart deviation  SS,  i.e.  PMEAN=SS*FIFUNK(XX/SS), where
C   FIFUNK(x)=f(x)+x*FI(x), f  and  FI  are density and distribution
C   functions of  N(0,1)  variable, respectively. We have modified the
C   algorithm 209 from CACAM for evaluation of  FIFUNK, to avoid the operations
C   of type  SS*XX/SS, which can give numerical errors when  SS  and  XX  are
C   both small.
C
C                     **   NUMERICAL ACCURACY  **
C
C
C  Obs. that, our general assumption is that the process is normalized, i.e.
C  Var (X(t)) = Var (X'(t)) = 1.0d0  Consequently all conditional variances
C  are less than  1, and usualy close to zero, i.e.  SS<1.0d0. Now since  SS<1.0d0
C  and  FIFUNK(x)<0.0000001  for  x<-4.5  we have defined  FIFUNK(x)=0.0d0
C  if  x<-4.5  and  FIFUNK(x)=x  if  x>
C  4.5.  Under we have a table with
C  exact values of  FIFUNK.
C
C    x       FIFUNK(x)
C
C   -5.0    0.00000005233
C   -4.5    0.00000069515
C   -4.0    0.00000711075
C    4.0    4.00000700000
C    4.5    4.50000100000
C
C Obviously the tresholds  -4.5  and  4.5  can be increased.
C
C
      IMPLICIT NONE
      REAL*8, intent(in) :: XX, SS
      REAL*8 X,Y,W,Z
      REAL*8, parameter :: SP = 0.398942280401433d0
      IF(XX.LT.4.5d0*SS) GO TO 1
      PMEAN=XX
      RETURN
1     IF(XX.GT.-4.5d0*SS) GO TO 3
      PMEAN=0.0d0
      RETURN
3     continue
      if (SS .LT. 0.0000001d0) then
      PMEAN=0.0d0
      RETURN
      end if

      X=XX/SS
      
      IF(X==0) goto 8
      Y=0.5d0*ABS(X)
      IF(Y<1.0d0) then
       W=Y*Y
       Z=((((((((0.000124818987d0*W-0.001075204047d0)*W
     1   +0.005198775019d0)*W-0.019198292d0)*W+0.05905403564d0)*W 
     2   -0.15196875136d0)*W+0.3191529327d0)*W-0.5319230073d0)*W
     3   +0.7978845606d0)*Y*2.0d0
      else
       Y=Y-2.0d0
       Z=(((((((((((((-0.000045255659d0*Y+0.00015252929d0)*Y
     *  -0.000019538132d0)*Y-0.000676904986d0)*Y
     1 +0.001390604284d0)*Y-0.000794620820d0)*Y
     2 -0.002034254874d0)*Y+0.006549791214d0)*Y-0.010557625006d0)*Y+
     3 0.011630447319d0)*Y-0.009279453341d0)*Y+0.005353579108d0)*Y-
     4 0.002141268741d0)*Y+0.000535310849d0)*Y+0.9999366575d0
      endif
      IF(X.GT.0.0d0) PMEAN=SS*SP*EXP(-0.5d0*X*X)+XX*0.5d0*(Z+1.0d0)
      IF(X.LT.0.0d0) PMEAN=SS*SP*EXP(-0.5d0*X*X)+XX*0.5d0*(1.0d0-Z)
      RETURN
8     PMEAN=SS*SP
      RETURN
      END FUNCTION PMEAN


      REAL*8 FUNCTION FI(XX)
C
C   Algorithm 209 from CACAM.
C   FI(xx)  is a distribution functions of  N(0,1)  variable.
C
      IMPLICIT NONE
      REAL*8, intent(in) :: XX
      REAL*8 X, Y,Z, W
      X=XX
      IF(X==0) then
        FI=0.5d0
        RETURN
      endif
      Y=0.5d0*ABS(X)
      IF(Y>3.0d0) then
       IF(X.GT.0.0d0) FI=1.0d0
       IF(X.LT.0.0d0) FI=0.0d0
       RETURN
      endif
      IF (Y<1.0d0) then
       W=Y*Y
       Z=((((((((0.000124818987d0*W-0.001075204047d0)*W
     1   +0.005198775019d0)*W-0.019198292d0)*W+0.05905403564d0)*W 
     2   -0.15196875136d0)*W+0.3191529327d0)*W-0.5319230073d0)*W
     3   +0.7978845606d0)*Y*2.0d0
      ELSE
       Y=Y-2.0d0
       Z=(((((((((((((-0.000045255659d0*Y+0.00015252929d0)*Y
     1  -0.000019538132d0)*Y-0.000676904986d0)*Y+0.001390604284d0)*Y
     2  -0.000794620820d0)*Y-0.002034254874d0)*Y+0.006549791214d0)*Y
     3  -0.010557625006d0)*Y+0.011630447319d0)*Y-0.009279453341d0)*Y
     4  +0.005353579108d0)*Y-0.002141268741d0)*Y+0.000535310849d0)*Y
     5  +0.9999366575d0
      endif
100   IF(X.GT.0.0d0) FI=0.5d0*(Z+1.0d0)
      IF(X.LT.0.0d0) FI=0.5d0*(1.0d0-Z)
      RETURN
      END FUNCTION FI

C     Version  1991-XII-14

C  The MREG program.
C
C
C  We consider a process X(I)=X(T(I)) at the grid of  N  points  T(1),...,T(N),
C
C          X(I) = -A(I) + Z*A(I+N) + Sum Xj*A(I+(j+1)*N) + Delta(I),
C
C  the sum disapears if M=1, j=1,...,M-1. We assume that Z,Xj are independend
C  standart Rayleigh, Gaussian distributed rv. and independent of the zero
C  mean Gaussian residual process,  with covariance structure given in  R,
C
C          R(i+(j-1)N) = Cov (Delta(T(i)), Delta(T(j))).
C
C  Additionally we have a zero mean Gaussian variable XN, 
C  independent of Z,Xj with  covariance structure defined by 
C  B(i)= Cov (Delta(T(i)),XN), i=1,...,N,  B(N+1)=Var(XN). 
C  Furthermore XN and Z,Xj satisfies the following equation system
C
C     (BB + (XN,0,...,0)^T =  AA*(Z,X1,...,Xm-1)^T       (***)
C
C  where AA is (M,M) matrix, BB is M-vector. We rewrite this equation, by
C  introducing a variable X_M=XN/SQRT(XN) and construct   new matrix AA1
c  by adding the column (SQRT(Var(XN)),0,...,0) and the row with only zeros.
C  The equations (***) writtes
C
C     (BB,0)^T  =  AA1*(Z,X1,...,Xm-1,Xm)^T       (****)
C
C  where AA1 is (M+1,M+1) matrix, We assume that the rank of AA1 is M,
C  otherwise the density is singular and we give a output F=0.CC
C
C  Let Y0 be a zero-mean Gaussian variable independent of Z,Xj 
C  with  covariance structure defined by 
C  DB(i)= Cov (Delta(T(i)),Y0), i=1,...,N,  DB(N+1)=Cov(XN,Y0), Var(Y0)=VDER. 
C  Let Y be defined by
C
C     Y=-DA(1) + Z*DA(2) + Sum Xj*DA(2+j) +Y0,
C
C  j=1,...,M-1. The program computes:
C
C  F = E[ Y^+ *1{ HH<X(I)<0 for all I, I=1,...,N}|Z,X1,...,X_M-1 solves (***)]
C      *f_{Z,X1,....,XM-1}(***).
C
C  In the simplest case NIT=0 we define (Delta(1),...,Delta(N),XN)=0.0d0 
C
C  We renormalize  vectors AN and DA, the covariance fkn R, DB
C  and VDER. Then by renormalization we choose the Gaussian variable X such
C  that F is written in the form
C
C  F = E[(D0(1)+X*D0(2)+Y1)^+*(PC+X*PD)^+*1{HH <A0(I)+X*B0(I)+Delta1(I)<0}]
C
C  Observe, PC+X*PD>0 defines integration region for X.
C  In the simplest case NIT=0 we define (Delta(1),...,Delta(N),Y1)=0.0d0 
C  For NIT=1 only (Delta(1),...,Delta(N))=0, i.e. we have to compute
C  a one dimensional integral. Finally by conditioning on X the problem is
C  put in the format of RIND-problem.
C
C   INF indicates whether one
C  has already called the subroutine before and ONLY! inputs BB, DA or A
C  was changed.
C
C  Observe the limitations are :  N<=100, 0<M <= 5 = MMAX.
C
C     revised pab 2007
C     - replaced all common blocks with modules
C     - fixed some bugs 

      SUBROUTINE MREG(F,R,B,DB,AA,BB,A,DA,VDER,M,N,NIT,INFR)
      USE SIZEMOD
      USE EPSMOD
      USE RINTMOD 
      USE INFCMOD
      USE SVD
      IMPLICIT NONE
!     INTEGER, PARAMETER :: MMAX = 5, NMAX = 101, RDIM = 10201
      INTEGER, PARAMETER :: Nw = MMAX+1
      INTEGER, intent(in) :: M,N,NIT,INFR
      REAL*8, intent(in) :: VDER
      REAL*8, intent(inout) :: F
      REAL*8, intent(inout) :: R(RDIM),B(NMAX),DB(NMAX),BB(Nw)
      REAL*8, intent(inout) :: A(Nw*NMAX),DA(Nw),AA(MMAX-2,MMAX-2)
      REAL*8, DIMENSION(Nw),    save :: AO, DA0, W1, E1
      REAL*8, DIMENSION(Nw,Nw), save :: AA1, U1,V1 
      REAL*8, DIMENSION(NMAX),  save :: DB1, SQ
      REAL*8, DIMENSION(RDIM),  save :: R1
      REAL*8, DIMENSION(Nw*NMAX), save :: A1
      REAL*8, DIMENSION(NMAX) :: A0,B0,B1
      REAL*8, DIMENSION(24) :: XX1,H1
      REAL*8, DIMENSION(2)  :: D0
      REAL*8, save :: QD, DB1N,VDER1,SQD, SDER1,DET1
      REAL*8  XIND, DB0N,FR1, XR1,XMI,XMA
      REAL*8  CC, PC, PD, P, X
      INTEGER INF1, I, I1, J, N1,I2,infoID
      INTEGER, save :: IDET, NNIT         
C     DIMENSION AA1(6,6),V1(6,6)
C     DIMENSION W1(6),AO(6),DA0(6)
C      DIMENSION A0(2*Nmax),B0(2*Nmax),B1(2*Nmax),DB1(2*Nmax)
!      DIMENSION D0(2)
!      DIMENSION XX1(24),H1(24)
C      COMMON /EPS/EPS,EPSS,CEPSS
C      COMMON/RINT/ C,FC
C
C  If INFR=0 we have to initiate conditioning and renormalization transf.
C
      INF1=0
      F=0.0d0
      IF (N.LT.0) RETURN
      IF (N.EQ.0) GO TO 2
      DO I=1,N
       DO I1=1,M+1
        A1(I+(I1-1)*N)=A(I+(I1-1)*N)
       ENDDO
      ENDDO
2     CONTINUE
      DO I=1,M+1
      DA0(I)=DA(I)
      ENDDO

C
C  Renormalization
C
      IF (INFR.EQ.1) GO TO 105
      NNIT=MIN(NIT,N)

      DO i=1,n
         SQ(i)=0.0d0
      ENDDO
      
      DO I=1,M
      DO J=1,M
         AA1(I,J)=AA(I,J)
      ENDDO
      ENDDO
      NNIT=MIN(NIT,N)
      
      QD   =B(N+1)
      DB1N =DB(N+1)
      VDER1=VDER

      IF(QD.le.eps) then

         DB1N = 0.0d0
         SQD  = 0.0d0
         NNIT = 0
         DO  I=1,N
            DB1(I) = DB(I)
            A1(I+(M+1)*N)=0.0d0
            SQ(I) = 0.0d0
         ENDDO

      else
           SQD   = SQRT(QD)            
           VDER1 = VDER1-DB1N*DB1N/QD
           DB1N  = DB1N/SQD
           DO I=1,N
              DB1(I) = DB(I)-DB(N+1)*(B(I)/QD)
              A1(I+(M+1)*N) = B(I)/SQD
          ENDDO
      end if

      SDER1 = 0.0d0
      IF (VDER1.GT.EPS) SDER1=SQRT(VDER1)
C      print *,'sqd,SDER1',sqd,SDER1     
C PAB: BUG DA0 M+2 can be larger than NW = MMAX+1
      DA0(M+2) = DB1N
      BB(M+1) = 0.0d0

      
      AA1(1,M+1)=SQD
      DO I=1,M
       AA1(I+1,M+1)=0.0d0
       AA1(M+1,I)=0.0d0
      ENDDO
! New call to avoid calling Numerical recipess SVDCMP
      CALL dsvdc(AA1,M+1,M+1,W1, E1, U1, V1, 11, infoID)
!      CALL SVDCMP(AA1,M+1,M+1,NW,NW,W1,V1)


      DET1 = 1.0d0
      idet = 0
      DO I=1,M+1
       IF ( W1(I).LT.0.00001d0 ) THEN
          idet  = idet+1
          W1(I) = 0.0d0
         DO J=1,M+1
            AO(J)=V1(J,I)
         ENDDO
        GO TO 35
        END IF
        DET1=DET1*W1(I)
35    CONTINUE
      ENDDO
C      print *,'det1',det1

      IF(DET1.LT.0.001d0) NNIT=0
c
c    Obs. QD can not be small since NNIT>0
c      
      IF (NNIT.GT.1) THEN
      DO I=1,N
         XR1=R(I+(I-1)*N)-B(I)*(B(I)/QD) 
         IF(XR1.GT.EPS)  THEN
              SQ(I)=SQRT(XR1)
          ENDIF
      ENDDO
           
      DO I=1,N
       DO J=1,N
         R1(J+(I-1)*N)=R(J+(I-1)*N)-B(I)*(B(J)/QD)
       ENDDO
      ENDDO
      
      
      END IF


105   CONTINUE
      if (idet.gt.1) return
C
C  Renormalization is done
C
      CALL R_ORT(CC,PC,PD,U1,V1,W1,AO,BB,A1,A0,B0,DA0,D0,DET1,M+1,N)
      IF(CC.LT.0.0d0) RETURN
      XMI=-C
      XMA= C
      IF(ABS(PD).LE.EPS.AND.PC.LT.0.0d0) RETURN
      IF(ABS(PD).LE.EPS) GO TO 102
      X=-PC/PD
      IF(PD.GT.0.0d0.AND.XMI.LT.X) XMI=X
      IF(PD.LT.0.0d0.AND.XMA.GT.X) XMA=X
102   CONTINUE
c      PRINT *,'XMI,XMA',XMI,XMA
      IF(NNIT.eq.1.AND.IAC.LT.1.OR.NNIT.eq.0.OR.XMI.GE.XMA) THEN
      CALL C1_C2(XMI,XMA,A0,B0,D0(1),D0(2),0.0d0,SQ,N)
c         PRINT *,'XMI,XMA',XMI,XMA
         F=GAUSINT(XMI,XMA,D0(1),D0(2),PC,PD)*CC
c         print *,'return',f,cc
         RETURN
      END IF
C
C ***********************************************************
C
C We shall condition on the values of X, XMI<X<XMA, but for some
C X values XIND will be zero leading to reduced accuracy. Hence we try
C to exclude them and narrow the interval [XMI,XMA]
c
c      PRINT *,XMI,XMA

      CALL C1_C2(XMI,XMA,A0,B0,D0(1),D0(2),SDER1,SQ,N)
c      PRINT *,XMI,XMA
      IF(FI(XMA)-FI(XMI).LT.EPSS) RETURN
      CALL GAUSS1(N1,H1,XX1,XMI,XMA,EPS0)
      DO I2=1,N1
       FR1=CC*H1(I2)*(PC+XX1(I2)*PD)
        DO I=1,N
         B1(I)=A0(I)+XX1(I2)*B0(I)
        ENDDO
       DB0N=D0(1)+XX1(I2)*D0(2)

C
C  INF1=1 means that both R1 and SQ are the same as in the previous
C  call of TWOREG subroutine, INF=0 indicates the new R and SQ.
c      
c      print *,'go in rind'
       CALL RIND(XIND,R1,B1,DB0N,DB1,SQ,VDER1,NNIT-1,N,INF1)
c      if (n.gt.29) print *,XIND,DB0N,VDER1,b1(N),b1(n-1)
       INF1=1
       F=F+FR1*XIND
      ENDDO
      RETURN
      END SUBROUTINE MREG



      SUBROUTINE R_ORT(C,PC,PD,U1,V1,W1,AO,BB,A0,A,B,DA,D0,DET,M,N)
      USE EXPACCMOD
      USE SIZEMOD
      IMPLICIT NONE
      INTEGER, PARAMETER :: Nw = MMAX+1
      INTEGER, intent(in) :: M,N
      REAL*8, DIMENSION(Nw,Nw), intent(inout)  :: U1, V1
      REAL*8, DIMENSION(Nw   ), intent(inout)  :: W1,AO,BB, DA
      REAL*8, DIMENSION(Nw*NMAX), intent(inout):: A0
      REAL*8, DIMENSION(NMAX), intent(inout) :: A, B
      REAL*8, DIMENSION(2),    intent(inout) :: D0
      REAL*8, intent(inout):: C,PC,PD,DET
      REAL*8, DIMENSION(Nw   ) :: XO
      REAL*8  DER0,DER1,P
      INTEGER I,J
      
      
C      COMMON/EXPACC/PMAX
      REAL*8, parameter :: SP = 0.398942280401433d0
      C=-999.
      CALL SVBKSB(U1,W1,V1,M,M,Nw,Nw,BB,XO)
      P    = 0.0d0
      DER0 = -DA(1)
      DER1 = 0.0d0
      DO I=1,M
      P    = P + XO(I) * XO(I)
      DER0 = DER0 + XO(I) * DA(I+1)
      DER1 = DER1 + AO(I) * DA(I+1)
      ENDDO
      IF (P.GT.PMAX) RETURN
      C=(SP**(M-2))*EXP(-0.5d0*P)/ABS(DET)
c      print *,'XO',XO(1),XO(2),XO(3),XO(4)
c      print *,'AO',AO(1),AO(2),AO(3),AO(4)
      if(N.lt.1) go to 100
      DO I=1,N
       A(I) = -A0(I)
       B(I) = 0.0d0
       DO J=1,M
        B(I) = B(I)+AO(J)*A0(I+J*N)
        A(I) = A(I)+XO(J)*A0(I+J*N)
       ENDDO
      ENDDO
100   continue
      D0(1)=DER0
      D0(2)=DER1
      PC=XO(1)
      PD=AO(1)
      RETURN
      END SUBROUTINE R_ORT


      REAL*8 FUNCTION pythag(a,b)
      IMPLICIT NONE
      REAL*8, intent(in) :: a,b
      REAL*8 absa,absb
      absa=abs(a)
      absb=abs(b)
      IF (absa.GT.absb) THEN
         pythag=absa*SQRT(1.0d0+(absb/absa)**2)
      ELSE
         IF (absb.EQ.0.0d0) THEN
            pythag=0.0d0
         ELSE
            pythag=absb*SQRT(1.0d0+(absa/absb)**2)
         ENDIF
      ENDIF
      RETURN
      END FUNCTION PYTHAG


      SUBROUTINE SVBKSB(U,W,V,M,N,MP,NP,B,X)
C
C   Solves  AX=B  for a vector  X, where  A  is specified by the arrays
C   U, W, V  as returned by SVDCMP.  M  and  N  are the logical
C   dimensions of  A, and will be equal for a square matrices.  MP  and  NP
C   are the physical dimensions of  A.  B  is the input right-hand side.
C   X  is the output solution vector. No input quantities are destroyed,
C   so the routine may be called sequentialy with different  B's.
C
      USE SIZEMOD
      IMPLICIT NONE
C      INTEGER, PARAMETER :: NMAX=100
C   Maximum anticipated value of N
      INTEGER, intent(in) :: M,N,MP,NP
      INTEGER ::  J,I,JJ
      REAL*8, intent(inout) :: U,W,V,B,X
      REAL*8 TMP, S 
      DIMENSION U(MP,NP),W(NP),V(NP,NP),B(MP),X(NP),TMP(NMAX)
      DO J=1,N
C   Cumulate U^T*B
         S=0.0d0
         IF (W(J).NE.0.0d0) THEN
C   Nonzero rezult only if  wj  is nonzero
           DO I=1,M
             S=S+U(I,J)*B(I)
           ENDDO 
           S=S/W(J)
C   This is the divide by  wj
         ENDIF
         TMP(J)=S
      ENDDO
      DO J=1,N
         S=0.0d0
         DO JJ=1,N
           S=S+V(J,JJ)*TMP(JJ)
         ENDDO
         X(J)=S
      ENDDO
      RETURN
      END SUBROUTINE SVBKSB



      SUBROUTINE SVDCMP(A,M,N,MP,NP,W,V)
C
C  Given a matrix  A, with logical dimensions  M  by  N  and physical
C  dimensions  MP  by  NP, this routine computes its singular value
C  decomposition,  A=U.W.V^T, see Numerical Recipes, by Press W.,H.
C  Flannery, B. P., Teukolsky S.A. and Vetterling W., T. Cambrige
C  University Press 1986, Chapter 2.9. The matrix  U  replaces A  on
C  output. The diagonal matrix of singular values  W  is ouyput as a vector
C  W. The matrix  V (not the transpose  V^T) is output as  V.  M  must be
C  greater or equal to  N; if it is smaller, then  A  should be filled up
C  to square with zero rows.
C
      USE SIZEMOD
      IMPLICIT NONE
C      PARAMETER (NMAX=100)
      INTEGER, intent(in) :: M,N,MP,NP
      INTEGER :: I,L,K,J, ITS, NM
      REAL*8, intent(inout) :: A,W,V
      REAL*8 RV1,G, S,SCALE,ANORM, F,H, C, Y,Z,X
C  Maximum anticipated values of  N
      DIMENSION A(MP,NP),W(NP),V(NP,NP),RV1(NMAX)
      IF(M.LT.N) THEN
          print  *, 'You must augment  A  with extra zero rows. stop'
          stop
      ENDIF
C  Householder reduction to bidiagonal form
      G=0.0d0
      SCALE=0.0d0
      ANORM=0.0
      DO I=1,N
        L=I+1
        RV1(I)=SCALE*G
        G=0.0d0
        S=0.0d0
        SCALE=0.0d0
        IF (I.LE.M) THEN
          DO K=I,M
            SCALE=SCALE+ABS(A(K,I))
          ENDDO
          IF (SCALE.NE.0.0d0) THEN
            DO K=I,M
              A(K,I)=A(K,I)/SCALE
              S=S+A(K,I)*A(K,I)
            ENDDO
            F=A(I,I)
            G=-SIGN(SQRT(S),F)
            H=F*G-S
            A(I,I)=F-G
            IF (I.NE.N) THEN
              DO J=L,N
                S=0.0d0
                DO K=I,M
                  S=S+A(K,I)*A(K,J)
                ENDDO
                F=S/H
                DO K=I,M
                  A(K,J)=A(K,J)+F*A(K,I)
                ENDDO
              ENDDO
            ENDIF
            DO K=I,M
              A(K,I)=SCALE*A(K,I)
            ENDDO
          ENDIF
        ENDIF
        W(I)=SCALE*G
        G=0.0d0
        S=0.0d0
        SCALE=0.0d0
        IF ((I.LE.M).AND.(I.NE.N)) THEN
          DO K=L,N
            SCALE=SCALE+ABS(A(I,K))
          ENDDO
          IF (SCALE.NE.0.0d0) THEN
            DO K=L,N
              A(I,K)=A(I,K)/SCALE
              S=S+A(I,K)*A(I,K)
            ENDDO
            F=A(I,L)
            G=-SIGN(SQRT(S),F)
            H=F*G-S
            A(I,L)=F-G
            DO K=L,N
              RV1(K)=A(I,K)/H
            ENDDO
            IF (I.NE.M) THEN
              DO J=L,M
                S=0.0d0
                DO K=L,N
                  S=S+A(J,K)*A(I,K)
                ENDDO
                DO K=L,N
                  A(J,K)=A(J,K)+S*RV1(K)
                ENDDO
              ENDDO
            ENDIF
            DO K=L,N
              A(I,K)=SCALE*A(I,K)
            ENDDO
          ENDIF
        ENDIF
        ANORM=MAX(ANORM,(ABS(W(I))+ABS(RV1(I))))
      ENDDO
c        print *,'25'
C   Accumulation of right-hand transformations.
      DO I=N,1,-1
        IF (I.LT.N) THEN
          IF (G.NE.0.0d0) THEN
            DO J=L,N
              V(J,I)=(A(I,J)/A(I,L))/G
C   Double division to avoid posible underflow.
            ENDDO
            DO J=L,N
              S=0.0d0
              DO K=L,N
                S=S+A(I,K)*V(K,J)
              ENDDO
              DO K=L,N
                V(K,J)=V(K,J)+S*V(K,I)
              ENDDO
            ENDDO
          ENDIF
          DO J=L,N
            V(I,J)=0.0d0
            V(J,I)=0.0d0
         ENDDO
        ENDIF
        V(I,I)=1.0d0
        G=RV1(I)
        L=I
32    ENDDO
c        print *,'32'

C  Accumulation of the left-hang transformation
      DO I=N,1,-1
        L=I+1
        G=W(I)
        IF (I.LT.N) THEN
          DO J=L,N
            A(I,J)=0.0d0
          ENDDO
        ENDIF
        IF (G.NE.0.0d0) THEN
          G=1.0d0/G
          IF (I.NE.N) THEN
            DO J=L,N
              S=0.0d0
              DO K=L,M
                S=S+A(K,I)*A(K,J)
              ENDDO
              F=(S/A(I,I))*G
              DO K=I,M
                A(K,J)=A(K,J)+F*A(K,I)
              ENDDO
            ENDDO
          ENDIF
          DO J=I,M
            A(J,I)=A(J,I)*G
37        ENDDO
        ELSE
          DO J=I,M
            A(J,I)=0.0d0
38        ENDDO
        ENDIF
        A(I,I)=A(I,I)+1.0d0
39    ENDDO
c        print *,'39'

C      Diagonalization of the bidiagonal form
C      Loop over singular values
      DO K=N,1,-1
C   Loop allowed iterations
        DO ITS=1,30
C   Test for spliting
          DO L=K,1,-1
            NM=L-1
C   Note that RV1(1) is always zero
            IF((ABS(RV1(L))+ANORM).EQ.ANORM) GO TO 2
            IF((ABS(W(NM))+ANORM).EQ.ANORM) GO TO 1
          ENDDO
c          print *,'41'
1         C=0.0d0
          S=1.0d0
          DO I=L,K
            F=S*RV1(I)
            IF ((ABS(F)+ANORM).NE.ANORM) THEN
              G=W(I)
              H=SQRT(F*F+G*G)
              W(I)=H
              H=1.0d0/H
              C= (G*H)
              S=-(F*H)
              DO J=1,M
                Y=A(J,NM)
                Z=A(J,I)
                A(J,NM)=(Y*C)+(Z*S)
                A(J,I)=-(Y*S)+(Z*C)
              ENDDO
            ENDIF
          ENDDO
c          print *,'43'
2         Z=W(K)
          IF (L.EQ.K) THEN
C     Convergence
            IF (Z.LT.0.0d0) THEN
C   Singular values are made nonnegative
              W(K)=-Z
              DO J=1,N
                V(J,K)=-V(J,K)
              ENDDO
            ENDIF
            GO TO 3
          ENDIF
          IF (ITS.EQ.30) THEN
              print *, 'No convergence in 30 iterations. stop.'
              stop
          ENDIF
          X=W(L)
          NM=K-1
          Y=W(NM)
          G=RV1(NM)
          H=RV1(K)
          F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.0*H*Y)
          G=SQRT(F*F+1.0d0)
          F=((X-Z)*(X+Z)+H*((Y/(F+SIGN(G,F)))-H))/X
C   Next  QR  transformation
          C=1.0d0
          S=1.0d0
          DO J=L,NM
            I=J+1
            G=RV1(I)
            Y=W(I)
            H=S*G
            G=C*G
            Z=SQRT(F*F+H*H)
            RV1(J)=Z
            C=F/Z
            S=H/Z
            F= (X*C)+(G*S)
            G=-(X*S)+(G*C)
            H=Y*S
            Y=Y*C
            DO NM=1,N
              X=V(NM,J)
              Z=V(NM,I)
              V(NM,J)= (X*C)+(Z*S)
              V(NM,I)=-(X*S)+(Z*C)
            ENDDO
c            print *,'45',F,H
            Z=pythag(F,H)
            W(J)=Z
C   Rotation can be arbitrary if  Z=0.
            IF (Z.NE.0.0d0) THEN
c            print *,1/Z
              Z=1.0d0/Z
c              print *,'*'
              C=F*Z
              S=H*Z
            ENDIF
            F= (C*G)+(S*Y)
            X=-(S*G)+(C*Y)
            DO NM=1,M
              Y=A(NM,J)
              Z=A(NM,I)
              A(NM,J)= (Y*C)+(Z*S)
              A(NM,I)=-(Y*S)+(Z*C)
            ENDDO
c          print *,'46'

          ENDDO
c          print *,'47'
          RV1(L)=0.0d0
          RV1(K)=F
          W(K)=X
        ENDDO
3     CONTINUE
      ENDDO
c        print *,'49'
       
      RETURN
      END SUBROUTINE SVDCMP

      END MODULE MREGMOD