Module EffPotFunctions

Use Control
Use Mathematics
! Use Model_Data_MSSM
Use Settings

! squared masses lesser than MIN_MASS are set to zero
Real(dp), parameter :: MIN_MASS = 1.0_dp
! for squared masses lesser than SMALL_MASS expansions are made
Real(dp), parameter :: SMALL_MASS = 1.0_dp
!Real(dp) :: hstep_pn=0.1_dp, hstep_sa=0.001_dp

 Real(dp) :: epsD = 1E-8_dp, epsM = 1.0E-8_dp, epsMnull = 1.0E-10_dp

 Interface SmallD
  Module Procedure SmallDR, SmallDC
 End Interface

 Interface SmallM
  Module Procedure SmallMR, SmallMC, SmallMR2, SmallMC2
 End Interface

  Interface VerySmallM
  Module Procedure VerySmallMR, VerySmallMC, VerySmallMR2, VerySmallMC2
 End Interface

Contains

! Function h(x,Q) for 1-loop effective potential
! following S.P. Martin, arXiv:hep-ph/0111209v2
! "Two-loop effective potential for a general renormalizable theory
! and softly broken supersymmetry"

Logical Function VerySmallMR(x,y,z,Q)
Implicit None
Real(dp), Intent(in) :: x,y,z
Real(dp), Intent(in) :: Q
If (Max(Max(Abs(x),Abs(y)),Abs(z)).le.epsMnull) Then
 VerySmallMR = .True.
Else If ((Abs(x)/Max(Q,Max(Abs(y),Abs(z)))).le.epsMnull) Then
 VerySmallMR = .True.
Else
 VerySmallMR = .False.
End if
End Function VerySmallMR

Logical Function VerySmallMC(x,y,z,Q)
Implicit None
Complex(dp), Intent(in) :: x,y,z
Real(dp), Intent(in) :: Q
If (Max(Max(Abs(x),Abs(y)),Abs(z)).le.epsMnull) Then
 VerySmallMC = .True.
Else If ((Abs(x)/Max(Q,Max(Abs(y),Abs(z)))).lt.epsMnull) Then
 VerySmallMC = .True.
Else
 VerySmallMC = .False.
End if
End Function VerySmallMC

Logical Function VerySmallMR2(x,y)
Implicit None
Real(dp), Intent(in) :: x,y
If (Max(Abs(x),Abs(y)).le.epsMnull) Then
 VerySmallMR2 = .True.
Else If (Abs(x/y).le.epsMnull) Then
 VerySmallMR2 = .True.
Else
 VerySmallMR2 = .False.
End if
End Function VerySmallMR2

Logical Function VerySmallMC2(x,y)
Implicit None
Complex(dp), Intent(in) :: x,y
If (Max(Abs(x),Abs(y)).le.epsMnull) Then
 VerySmallMC2 = .True.
Else If (Abs(x/y).le.epsMnull) Then
 VerySmallMC2 = .True.
Else
 VerySmallMC2 = .False.
End if
End Function VerySmallMC2

recursive Complex(dp) Function I_SMartin_ex(x_in,y_in,z_in,Q2) result(res)
Implicit None
! x_in,y_in,z_in are supposed to be squared masses
! symmetric under interchange of any two of these
Real(dp),Intent(in) :: x_in,y_in,z_in,Q2
Integer :: ncontrol
Complex(dp) :: x,y,z,tempI
x = x_in
y = y_in
z = z_in
res = Cmplx(0._dp,0._dp)
ncontrol = 0

If (Q2.le.0._dp) Then
  write(*,*) 'invalid renormalization scale!'
  return
end if

Call SortArguments3C(x,y,z) ! sort x<=y<=z

If (VerySmallM(x,y,z,Q2)) Then ! x=0
!  Write(*,*) x
  If (VerySmallM(y,z,z,Q2)) Then ! y=0 (z >= y holds)
    If (VerySmallM(z,z,z,Q2)) then ! z=0
      res = ZeroC ! I(0,0,0)
      ncontrol = ncontrol +1
    Else ! z!=0
         ! I(0,0,z)
         res = -0.5_dp*z*(log(z/Q2))**2 &
         & + 2*z*log(z/Q2)-2.5_dp*z-pi**2*z/6._dp
         ncontrol = ncontrol +1
!       res = 0._dp
    End If
  Else ! y != 0
    If (VerySmallM(Abs(y-z),Abs(z))) Then ! y=z
	 ! I(0,y,y)
         res = -y*(zlog(y/Q2))**2+4._dp*y*zlog(y/Q2)-5._dp*y
         ncontrol = ncontrol +1
! res = 0._dp
    Else ! z != y
          ! I(0,y,z)
	res = (z-y)*(CLi2_mod(y/z)-log(z/y)*log((z-y)/Q2) &
	& +0.5_dp*(log(z/Q2))**2-pi**2/6._dp)  &
	& -2.5_dp*(z+y)+2*y*log(y/Q2)+2*z*log(z/Q2) &
	& -z*log(z/Q2)*log(y/Q2)
	ncontrol = ncontrol +1
! res = 0._dp
    End If
  End If
Else If (SmallM(x,y,z,Q2)) then ! x small
  If (SmallM(y,z,z,Q2)) then ! y small (z >= y holds)

      ! I(eps1,eps2,z)
      tempI = I_SMartin_ex(0._dp,0._dp,Real(z,dp),Q2)
      res = tempI + x/z*(-z-tempI + J_SMartin_1x(Real(z,dp),Q2)*(3._dp-zlog(x/Q2))) &  !! SWITCH 2
      & + y/z*(-z-tempI + J_SMartin_1x(Real(z,dp),Q2)*(3._dp-log(y/Q2))) &
      & + x*y/z**2*(-2._dp*tempI+4._dp*J_SMartin_1x(Real(z,dp),Q2) - 2.*z &
      & - (z+J_SMartin_1x(Real(z,dp),Q2))*(log(x/Q2)+log(y/Q2))+z*log(x/Q2)*log(y/Q2))
      ncontrol = ncontrol +1

  else ! y large (implies z also large)
    If (SmallD(y,z)) then ! y=z (large)
      ! I(eps,y,y)
      res = 2._dp
      res = 2._dp*J_SMartin_1x(Real(y,dp),Q2) - 2._dp*y &
      & - J_SMartin_2x(Real(y,dp),Real(y,dp),Q2)/y &
      & + x*(4._dp+J_SMartin_2x(Real(y,dp),Real(y,dp),Q2)/(2._dp*y**2) &  !! SWITCH 1
      & + 3._dp/y*J_SMartin_1x(Real(y,dp),Q2) &
      & -(1._dp+J_SMartin_1x(Real(y,dp),Q2)/y)*log(x/Q2))
      ncontrol = ncontrol +1

    else ! y != z
      ! I(eps,y,z)
      tempI = I_SMartin_ex(0._dp,Real(y,dp),Real(z,dp),Q2)
      res = tempI + x*(-(y+z)*tempI-2.*J_SMartin_2x(Real(y,dp),Real(z,dp),Q2) &  !! SWITCH 3
      & +3._dp*y*J_SMartin_1x(Real(y,dp),Q2) + 3.*z*J_SMartin_1x(Real(z,dp),Q2) &
      & -z*J_SMartin_1x(Real(y,dp),Q2)-y*J_SMartin_1x(Real(z,dp),Q2) - (y+z)**2 &
      & + (y-z)*(J_SMartin_1x(Real(z,dp),Q2)-J_SMartin_1x(Real(y,dp),Q2))*log(x))/(y-z)**2
      ncontrol = ncontrol +1
    end if
  end if

else ! x large

  ! most general case
  ! I(x,y,z)
  res = Cmplx(0._dp,0._dp)
  res=res + 0.5_dp*(x-y-z)*log(y/Q2)*log(z/Q2)
  res=res + 0.5_dp*(y-x-z)*log(x/Q2)*log(z/Q2)
  res=res + 0.5_dp*(z-x-y)*log(x/Q2)*log(y/Q2)
  res=res + 2._dp*x*log(x/Q2)
  res=res + 2._dp*y*log(y/Q2)
  res=res + 2._dp*z*log(z/Q2)
  res=res - 2.5_dp*(x+y+z)
  res=res - 0.5*XiFunction(x,y,z)
  ncontrol = ncontrol +1

end if

! end of all ifs

! control variable against sloppy programming
if (ncontrol.gt.1) then
  write(*,*) 'in I_SMartin_ex more than one case applies. sloppy programming!'
end if
End Function I_SMartin_ex

Complex(dp) Function I_SMartinC(x,y,z,Q2)
Implicit None
Complex(dp), Intent(in) :: x,y,z
Real(dp), Intent(in) :: Q2
I_SMartinC = I_SMartin(Real(x,dp),Real(y,dp),Real(z,dp),Q2)
End Function

Complex(dp) Function h_SMartin(m2,Q2)
Implicit None
Complex(dp) :: z
Real(dp),Intent(in) :: m2,Q2
! input: squared mass m2, squared renormalization scale Q2.
If (Q2.le.0._dp) Then
  write(*,*) 'invalid renormalization scale!'
  return
end if
z = CMPLX(m2,0._dp)
If (abs(m2).le.MIN_MASS) Then
   h_SMartin = 0._dp
Else
   h_SMartin = 0.25_dp * (m2**2) * (log(z/Q2)-1.5_dp)
End If
End Function h_SMartin

Real(dp) Function DumbDiff2(func,x,k1,k2,h)
! finite difference method
Implicit None
Real(dp),external :: func
Integer,intent(in) :: k1,k2
Real(dp), Intent(in) :: x(2),h
Real(dp) :: res
if (k1.eq.k2) then
 DumbDiff2 = (func1D(x(k1)-h)-2._dp*func1D(x(k1))+func1D(x(k1)+h))/(h**2.)
else
 res = func((/x(1)+h,x(2)+h/))
 res = res-func((/x(1)+h,x(2)-h/))
 res = res-func((/x(1)-h,x(2)+h/))
 res = res+func((/x(1)-h,x(2)-h/))
 res = res/(4._dp*h**2)
 DumbDiff2=res
end if
Contains
Real(dp) function func1D(y)
Real(dp),Intent(in) :: y
Real(dp) :: ytemp(2)
ytemp=x
ytemp(k1)=y
func1D = func(ytemp)
end function
End function

Real(dp) function partialDiff(func,x,k,n,h)
! wrapper function to calculate the partial derivative of func w.r.t. to k-th variable
! evaluated at the n-dim. point x, with stepsize h.
Implicit None
Integer,Intent(in) :: n,k
Real(dp),External :: func
Real(dp),Intent(in) :: x(n),h
Real(dp) :: xtemp(n)
partialDiff=DumbDiff(func_k,x(k),h)
  Contains
  Real(dp) function func_k(y)
  Implicit None
  Real(dp) :: y
  xtemp = x
  xtemp(k)=y
  func_k = func(xtemp)
  End function
End function

Real(dp) Function DumbDiff(f,x,h)
implicit none
Real(dp),intent(in) :: x,h
real(dp),external :: f
DumbDiff = (-f(x+2*h)+8*f(x+h)-8*f(x-h)+f(x-2*h))/(12.*h)
end function

Real(dp) Function dfridr(func,x,h,err,iout)
! algorithm based on
! Ridders, C.J.F. 1982,
! Advances in Engineering Software
! vol. 4, no. 2, pp. 75–76.

! Returns the derivative of a function func at a point x
! by Ridders' method of polynomial extrapolation.
! The value h is input as an estimated initial stepsize; it need not be small, but rather should
! be an increment in x over which func changes substantially.
! An estimate of the error in the derivative is returned as err.
! Parameters: Stepsize is decreased by CON at each iteration. Max size of tableau is set by NTAB. Return when error is SAFE worse than the best so far.
Implicit None
INTEGER, parameter :: NTAB=40
! original: CON=1.4
Real(dp), parameter :: CON=1.4_dp,CON2=CON*CON,BIG=1.D30,safe=20.0_dp
REAL(dp),Intent(in) :: h,x
Real(dp) :: err_int
Real(dp),Intent(out) :: err
Real(dp), EXTERNAL :: func
INTEGER i,j
INTEGER,Intent(out) :: iout
REAL(dp) :: errt,fac,hh,a(NTAB,NTAB)
if (h.eq.0.) write(*,*) 'h must be nonzero in dfridr'
hh=h
a(1,1)=(func(x+hh)-func(x-hh))/(2._dp*hh)
dfridr = a(1,1) ! the zero-th order solution
err_int=BIG
iout=2
do i=2,NTAB
! Successive columns in the Neville tableau will go to smaller
! stepsizes and higher orders of extrapolation.
  hh=hh/CON
  a(1,i)=(func(x+hh)-func(x-hh))/(2._dp*hh)
  ! Try new, smaller stepsize.
  fac=CON2
  do j=2,i
    ! Compute extrapolations of various orders, requiring no new
    ! function evaluations.
    a(j,i)=(a(j-1,i)*fac-a(j-1,i-1))/(fac-1._dp)
    fac=CON2*fac
    errt=max(abs(a(j,i)-a(j-1,i)),abs(a(j,i)-a(j-1,i-1)))
    ! The error strategy is to compare each new extrapolation to one order lower, both at
    ! the present stepsize and the previous one.
    if (errt.le.err_int) then
    ! If error is decreased, save the improved answer.
       err_int=errt
      dfridr=a(j,i)
    endif
  enddo
  if(abs(a(i,i)-a(i-1,i-1)).ge.SAFE*err_int) then
  iout=2*i
  exit
  end if
  ! If higher order is worse by a signicant factor SAFE, then quit early.
enddo
err=err_int
END function

! Subroutine dfridrMulDimDummy(dfunc,d,x,h,err,der)
! EXTERNAL :: dfunc
! Integer, Intent(in) :: d
! Real(dp), Intent(out) :: der(d), err
! REAL(dp),Intent(in) :: h,x
! INTEGER, parameter :: NTAB=40
! Real(dp), parameter :: CON=1.4_dp,CON2=CON*CON,BIG=1.D30,SAFE=20.0_dp
! Real(dp) :: err_int, fm(d),  fp(d)
! INTEGER :: i,j
! Real(dp), Allocatable :: a(:,:,:)
! REAL(dp) :: errt,fac,hh
!
! Allocate( a(d,NTAB,NTAB) )
!
!
! Write(*,*) "a"
! der = 0._dp
! End subroutine dfridrMulDimDummy

Subroutine dfridrMulDim(dfunc,d,x,h,err,der)
! algorithm based on
! Ridders, C.J.F. 1982,
! Advances in Engineering Software
! vol. 4, no. 2, pp. 75–76.

! Returns the derivative of a function func at a point x
! by Ridders' method of polynomial extrapolation.
! The value h is input as an estimated initial stepsize; it need not be small, but rather should
! be an increment in x over which func changes substantially.
! An estimate of the error in the derivative is returned as err.
! Parameters: Stepsize is decreased by CON at each iteration. Max size of tableau is set by NTAB. Return when error is SAFE worse than the best so far.
Implicit None
EXTERNAL :: dfunc
INTEGER, parameter :: NTAB=40
! original: CON=1.4
Integer, Intent(in) :: d
Real(dp), Intent(out) :: der(d)
Real(dp), parameter :: CON=1.4_dp,CON2=CON*CON,BIG=1.D30,safe=20.0_dp
REAL(dp),Intent(in) :: h,x
Real(dp) :: err_int, fm(d),  fp(d)
Real(dp),Intent(out) :: err
INTEGER :: i,j
Real(dp), Allocatable :: a(:,:,:)
REAL(dp) :: errt,fac,hh

Allocate( a(d,NTAB,NTAB) )

if (h.eq.0.) write(*,*) 'h must be nonzero in dfridr'
hh=h
Call dfunc(x+hh,fp)
Call dfunc(x-hh,fm)
a(:,1,1)=(fp-fm)/(2._dp*hh)
der = a(:,1,1) ! the zero-th order solution
err_int=BIG
do i=2,NTAB
! Successive columns in the Neville tableau will go to smaller
! stepsizes and higher orders of extrapolation.
  hh=hh/CON
Call dfunc(x+hh,fp)
Call dfunc(x-hh,fm)
  a(:,1,i)=(fp-fm)/(2._dp*hh)
  ! Try new, smaller stepsize.
  fac=CON2
  do j=2,i
    ! Compute extrapolations of various orders, requiring no new
    ! function evaluations.
    a(:,j,i)=(a(:,j-1,i)*fac-a(:,j-1,i-1))/(fac-1._dp)
    fac=CON2*fac
    errt=max(maxval(abs(a(:,j,i)-a(:,j-1,i))),maxval(abs(a(:,j,i)-a(:,j-1,i-1))))
    ! The error strategy is to compare each new extrapolation to one order lower, both at
    ! the present stepsize and the previous one.
    if (errt.le.err_int) then
    ! If error is decreased, save the improved answer.
       err_int=errt
      der=a(:,j,i)
    endif
  enddo
  if(maxval(abs(a(:,i,i)-a(:,i-1,i-1))).ge.SAFE*err_int) then
  exit
  end if
  ! If higher order is worse by a signicant factor SAFE, then quit early.
enddo
err=err_int
Deallocate(a)
END Subroutine

Real(dp) FUNCTION D2FRIDRMIXED(FUNC,X,Y,HX,HY,ERR)
IMPLICIT NONE
INTEGER NTAB
Real(dp) :: ERR,HX,HY,X,Y,FUNC,CON,CON2,BIG,SAFE
PARAMETER(CON=1.4_dp,CON2=CON*CON,BIG=1.0E+30_dp,NTAB=40,safe=20.0_dp)
EXTERNAL FUNC
INTEGER I,J
Real(dp) :: ERRT,FAC,HHX,HHY,A(NTAB,NTAB)
! IF(HX.EQ.0.0_dp) PAUSE 'HX MUST BE NONZERO IN DFRIDR'
! IF(HY.EQ.0.0_dp) PAUSE 'HY MUST BE NONZERO IN DFRIDR'
HHX=HX
HHY=HY
A(1,1)=(FUNC(X+HHX,Y+HHY)-FUNC(X+HHX,Y-HHY) &
& -FUNC(X-HHX,Y+HHY)+FUNC(X-HHX,Y-HHY))/(4*HHX*HHY)
ERR=BIG
DO I=2,NTAB
HHX=HHX/CON
HHY=HHY/CON
A(1,I)=(FUNC(X+HHX,Y+HHY)-FUNC(X+HHX,Y-HHY) &
& -FUNC(X-HHX,Y+HHY)+FUNC(X-HHX,Y-HHY))/(4.0_dp*HHX*HHY)
FAC=CON2
DO J=2,I
A(J,I)=(A(J-1,I)*FAC-A(J-1,I-1))/(FAC-1.0_dp)
FAC=CON2*FAC
ERRT=MAX1(ABS(A(J,I)-A(J-1,I)),ABS(A(J,I)-A(J-1,I-1)))
IF (ERRT.LE.ERR) THEN
ERR=ERRT
D2FRIDRMIXED=A(J,I)
ENDIF
ENDDO
IF(ABS(A(I,I)-A(I-1,I-1)).GE.SAFE*ERR) RETURN
ENDDO
RETURN
END function

Real(dp) FUNCTION D2FRIDR(FUNC,X,H,ERR)
IMPLICIT NONE
INTEGER NTAB
Real(dp) :: ERR,H,X,FUNC,CON,CON2,BIG,SAFE
PARAMETER(CON=1.4_dp,CON2=CON*CON,BIG=1.0E+30_dp,NTAB=40,safe=20.0_dp)
EXTERNAL FUNC
INTEGER I,J
Real(dp) :: ERRT,FAC,HH,A(NTAB,NTAB)

! Write(*,*) "X", X, H

! IF(H.EQ.0.0_dp) PAUSE 'H MUST BE NONZERO IN DFRIDR'
HH=H
A(1,1)=(FUNC(X+HH)-2.0_dp*FUNC(X)+FUNC(X-HH))/(HH**2.0_dp)
ERR=BIG
DO I=2,NTAB
HH=HH/CON
! Write(*,*) "X", X, HH
A(1,I)=(FUNC(X+HH)-2.0_dp*FUNC(X)+FUNC(X-HH))/(HH**2.0_dp)
! Write(*,*) "end"
FAC=CON2
DO J=2,I
A(J,I)=(A(J-1,I)*FAC-A(J-1,I-1))/(FAC-1.0_dp)
FAC=CON2*FAC
ERRT=MAX1(ABS(A(J,I)-A(J-1,I)),ABS(A(J,I)-A(J-1,I-1)))
IF (ERRT.LE.ERR) THEN
ERR=ERRT
D2FRIDR=A(J,I)
ENDIF
ENDDO
IF(ABS(A(I,I)-A(I-1,I-1)).GE.SAFE*ERR) RETURN
ENDDO
RETURN
END function

Subroutine D2FRIDRMIXEDmuldim(dFUNC,D,X,Y,HX,HY,ERR,der)
IMPLICIT NONE
EXTERNAL dFUNC
Integer, Intent(in) :: d
Real(dp), Intent(out) :: der(d)
Real(dp) :: ERR,HX,HY,X,Y,FUNC,CON,CON2,BIG,SAFE
Real(dp) :: fxpyp(d), fxpym(d), fxmyp(d), fxmym(d)
INTEGER, Parameter :: NTAB=40
PARAMETER(CON=1.4_dp,CON2=CON*CON,BIG=1.0E+30_dp,safe=20.0_dp)
Integer :: I,J
Real(dp) :: ERRT,FAC,HHX,HHY
Real(dp), Allocatable :: a(:,:,:)

Allocate( a(d,NTAB,NTAB) )

! IF(HX.EQ.0.0_dp) PAUSE 'HX MUST BE NONZERO IN DFRIDR'
! IF(HY.EQ.0.0_dp) PAUSE 'HY MUST BE NONZERO IN DFRIDR'
HHX=HX
HHY=HY

Call dFUNC(X+HHX,Y+HHY,fxpyp)
Call dFUNC(X+HHX,Y-HHY,fxpym)
Call dFUNC(X-HHX,Y+HHY,fxmyp)
Call dFUNC(X-HHX,Y-HHY,fxmym)

A(:,1,1)=(fxpyp-fxpym -fxmyp+fxmym)/(4*HHX*HHY)
ERR=BIG
DO I=2,NTAB
HHX=HHX/CON
HHY=HHY/CON

Call dFUNC(X+HHX,Y+HHY,fxpyp)
Call dFUNC(X+HHX,Y-HHY,fxpym)
Call dFUNC(X-HHX,Y+HHY,fxmyp)
Call dFUNC(X-HHX,Y-HHY,fxmym)

A(:,1,I)=(fxpyp-fxpym -fxmyp+fxmym)/(4.0_dp*HHX*HHY)
FAC=CON2
DO J=2,I
A(:,J,I)=(A(:,J-1,I)*FAC-A(:,J-1,I-1))/(FAC-1.0_dp)
FAC=CON2*FAC
ERRT=MAX1(MaxVal(ABS(A(:,J,I)-A(:,J-1,I))),MaxVal(ABS(A(:,J,I)-A(:,J-1,I-1))))
IF (ERRT.LE.ERR) THEN
ERR=ERRT
der=A(:,J,I)
ENDIF
ENDDO
IF(MaxVal(ABS(A(:,I,I)-A(:,I-1,I-1))).GE.SAFE*ERR) Exit
ENDDO

Deallocate(a)

END Subroutine

Subroutine D2FRIDRmuldim(dFUNC,d,X,H,ERR,der)
IMPLICIT NONE
EXTERNAL dFUNC
INTEGER, Parameter :: NTAB=40
Integer, Intent(in) :: d
Real(dp), Intent(out) :: der(d)
Real(dp) :: ERR,H,X,FUNC,CON,CON2,BIG,SAFE
PARAMETER(CON=1.4_dp,CON2=CON*CON,BIG=1.0E+30_dp,safe=20.0_dp)
INTEGER I,J
Real(dp) :: ERRT,FAC,HH, fn(d), f(d),fp(d)
Real(dp), Allocatable :: a(:,:,:)

Allocate( a(d,NTAB,NTAB) )
! Write(*,*) "X", X, H

! IF(H.EQ.0.0_dp) PAUSE 'H MUST BE NONZERO IN DFRIDR'
HH=H

Call dFunc(x+hh,fp)
Call dFunc(x,f)
Call dFunc(x-hh,fn)

A(:,1,1)=(fp-2.0_dp*f+fn)/(HH**2.0_dp)
ERR=BIG
DO I=2,NTAB
HH=HH/CON

Call dFunc(x+hh,fp)
Call dFunc(x,f)
Call dFunc(x-hh,fn)

A(:,1,I)=(fp-2.0_dp*f+fn)/(HH**2.0_dp)
! Write(*,*) "end"
FAC=CON2
DO J=2,I
A(:,J,I)=(A(:,J-1,I)*FAC-A(:,J-1,I-1))/(FAC-1.0_dp)
FAC=CON2*FAC
ERRT=MAX1(MaxVal(ABS(A(:,J,I)-A(:,J-1,I))),MaxVal(ABS(A(:,J,I)-A(:,J-1,I-1))))
IF (ERRT.LE.ERR) THEN
ERR=ERRT
der=A(:,J,I)
ENDIF
ENDDO
IF(MaxVal(ABS(A(:,I,I)-A(:,I-1,I-1))).GE.SAFE*ERR) RETURN
ENDDO
Deallocate(a)
END Subroutine

Real(dp) function partialDiff_Ridders(func,x,k,n,err,h,iout)
! wrapper function to calculate the partial derivative of func w.r.t. to k-th variable
! evaluated at the n-dim. point x, starting with stepsize h.
! returns error estimate err
Implicit None
Integer,Intent(in) :: n,k
Integer :: iout
Real(dp),External :: func
Real(dp),Intent(in) :: x(n), h
Real(dp),Intent(out) :: err
Real(dp) :: xtemp(n)

partialDiff_Ridders=dfridr(func_k,x(k),Min(1._dp,hstep_pn*x(k)),err,iout)
! partialDiff_Ridders=dfridr(func_k,x(k),hstep_pn*x(k),err,iout)
  Contains
  Real(dp) function func_k(y)
  ! projects out other components
  Implicit None
  Real(dp) :: y
  xtemp = x
  xtemp(k)=y
  func_k = func(xtemp)
  End function
End function

function partialDiff_RiddersMulDim(func,d,x,k,n,err)
! wrapper function to calculate the partial derivative of func w.r.t. to k-th variable
! evaluated at the n-dim. point x, starting with stepsize h.
! returns error estimate err
Implicit None
External :: func
Integer,Intent(in) :: n,k,d
Integer :: iout
Real(dp) :: partialDiff_RiddersMulDim(d)
Real(dp) :: test(d)
Real(dp),Intent(in) :: x(n)
Real(dp),Intent(out) :: err
Real(dp) :: xtemp(n),hrel(n)

! hrel(1) = hstep_sa_vd
! hrel(2) = hstep_sa_vu

! Write(*,*) x
Call dfridrMulDim(func_k,d,x(k),Min(1._dp,hstep_sa*x(k)),err,test)
! Call dfridrMulDim(func_k,d,x(k),hstep_sa*x(k),err,test)
partialDiff_RiddersMulDim = test
! dfridrMulDim(dfunc,d,x,h,err,der)

  Contains

  subroutine func_k(y,out)
  ! projects out other components
  Implicit None
  Real(dp),Intent(in) :: y
  Real(dp),Intent(out) :: out(d)
  xtemp = x
  xtemp(k)=y
  Call func(xtemp,out)
!   func_k1 = func(xtemp)
  End subroutine

End function

function partialDiffXY_RiddersMulDim(func,d,x,k1,k2,n,err)
! wrapper function to calculate second partial derivative w.r.t. k1-th variable, then k2-th
! evaluated at n-dim. point x
! applies partialDiff_Ridders twice
Implicit None
External :: func
Integer,Intent(in) :: n,k1,k2,d
Real(dp) :: partialDiffXY_RiddersMulDim(d)
Integer :: itemp=0,i1
Real(dp),Intent(in) :: x(n)
Real(dp),Intent(out) :: err
Real(dp) :: xtemp(n), hrel(n)

! hrel(1) = hstep_sa_vd
! hrel(2) = hstep_sa_vu

If (k1.eq.k2) Then
!  Call D2FRIDRmuldim(func_k1,d,x(k1),hstep_sa*x(k1),err,partialDiffXY_RiddersMulDim)
!  Call D2FRIDRmuldim(func_k1,d,x(k1),hstep_sa,err,partialDiffXY_RiddersMulDim)
 Call D2FRIDRmuldim(func_k1,d,x(k1),Min(1._dp,x(k1)*hstep_sa),err,partialDiffXY_RiddersMulDim)
Else
!  Call D2FRIDRMIXEDmuldim(func_k1k2,d,x(k2),x(k1),hstep_sa*x(k2), hstep_sa*x(k1),err,partialDiffXY_RiddersMulDim)
!  Call D2FRIDRMIXEDmuldim(func_k1k2,d,x(k2),x(k1),hstep_sa, hstep_sa,err,partialDiffXY_RiddersMulDim)
Call D2FRIDRMIXEDmuldim(func_k1k2,d,x(k2),x(k1),Min(1._dp,hstep_sa*x(k2)), Min(1._dp, hstep_sa*x(k1)),err,partialDiffXY_RiddersMulDim)
end if

Contains

  subroutine func_k1(y,out)
  ! projects out other components
  Implicit None
  Real(dp),Intent(in) :: y
  Real(dp),Intent(out) :: out(d)
  xtemp = x
  xtemp(k1)=y
  Call func(xtemp,out)
!   func_k1 = func(xtemp)
  End subroutine

  subroutine func_k1k2(y1,y2,out)
  ! projects out other components
  Implicit None
  Real(dp),Intent(in) :: y1,y2
  Real(dp),Intent(out) :: out(d)
  xtemp = x
  xtemp(k2)=y1
  xtemp(k1)=y2
  Call func(xtemp,out)
!   func_k1k2 = func(xtemp)
  End subroutine

End function

Subroutine ScanEffPot(func,func2,x0,y0,hx,hy,N)
Implicit none
Real(dp),External :: func,func2
Real(dp),Intent(in) :: x0,y0,hx,hy
Integer,intent(in) :: N
Real(dp) :: x,y,dx,dy
Integer :: i,j
open(394729,file="scaneffpot.out",status="replace")
open(394728,file="scaneffpot2L.out",status="replace")
write(394729,*) x0,y0,func((/x0,y0/))
write(394728,*) x0,y0,func2((/x0,y0/))
dx=(hx)/(N-1)
dy=(hy)/(N-1)
do i=-N+2,N
 do j=-N+2,N
  x=x0+(i-1)*dx
  y=y0+(j-1)*dy
  write(394729,*) x,y,func((/x,y/))
  write(394728,*) x,y,func2((/x,y/))
 end do
end do
close(394729)
close(394728)
End subroutine

Real(dp) function partialDiffXY_Ridders(func,x,k1,k2,n,err,h,iout)
! wrapper function to calculate second partial derivative w.r.t. k1-th variable, then k2-th
! evaluated at n-dim. point x
! applies partialDiff_Ridders twice
Implicit None
Integer,Intent(in) :: n,k1,k2
Integer :: itemp=0,i1,iout
Real(dp),External :: func
Real(dp),Intent(in) :: x(n),h
Real(dp),Intent(out) :: err
Real(dp) :: xtemp(n),hrel(n)
iout=0

! hrel(1) = hstep_sa_vd
! hrel(2) = hstep_sa_vu

! ! partialDiffXY_Ridders = partialDiff_Ridders(funcD1p,x,k2,n,err,h,i1)
If (k1.eq.k2) Then
!  write(*,*) k1,k2
! partialDiffXY_Ridders = D2FRIDR(func_k1,x(k1),hstep_pn*x(k1),err)
! partialDiffXY_Ridders = D2FRIDR(func_k1,x(k1),hstep_pn,err)
partialDiffXY_Ridders = D2FRIDR(func_k1,x(k1),Min(1._dp,hstep_pn*x(k1)),err)
Else
partialDiffXY_Ridders = D2FRIDRMIXED(func_k1k2,x(k2),x(k1),Min(1._dp,hstep_pn*x(k2)), Min(1._dp,hstep_pn*x(k1)),err)
! partialDiffXY_Ridders = D2FRIDRMIXED(func_k1k2,x(k2),x(k1),hstep_pn, hstep_pn,err)
! partialDiffXY_Ridders = D2FRIDRMIXED(func_k1k2,x(k2),x(k1),hstep_pn*x(k2), hstep_pn*x(k1),err)
end if

Contains
  Real(dp) function funcD1p(x)
  Implicit None
  Real(dp) :: x(n)
  funcD1p = partialDiff_Ridders(func,x,k1,n,err,h,itemp)
  iout = iout + itemp
  End function

  Real(dp) function func_k1(y)
  ! projects out other components
  Implicit None
  Real(dp) :: y
  xtemp = x
  xtemp(k1)=y
  func_k1 = func(xtemp)
  End function

  Real(dp) function func_k1k2(y1,y2)
  ! projects out other components
  Implicit None
  Real(dp) :: y1,y2
  xtemp = x
  xtemp(k2)=y1
  xtemp(k1)=y2
  func_k1k2 = func(xtemp)
  End function

End function

Real(dp) function partialDiffXY_Ridders2(func,x,k1,k2,n,err,h,iout)
! wrapper function to calculate second partial derivative w.r.t. k1-th variable, then k2-th
! evaluated at n-dim. point x
! applies partialDiff_Ridders twice
Implicit None
Integer,Intent(in) :: n,k1,k2
Integer :: itemp=0,i1,iout
Real(dp),External :: func
Real(dp),Intent(in) :: x(n),h(2)
Real(dp),Intent(out) :: err
iout=0
partialDiffXY_Ridders2 = partialDiff_Ridders(funcD1p,x,k2,n,err,h(k2),i1)
Contains
  Real(dp) function funcD1p(x)
  Implicit None
  Real(dp) :: x(n)
  funcD1p = partialDiff_Ridders(func,x,k1,n,err,h(k1),itemp)
  iout = iout + itemp
  End function
End function

!! These two-loop integral functions J(x),J(x,y),I(x,y,z) are taken from
!! Stephen P. Martin, arXiv:hep-ph/111209v2, 29 Mar 2002
!! FERMILAB-Pub-01/348-T

Real(dp) Function J_SMartin_1x(x,Q2)
Implicit None
! x is supposed to be a squared mass
Real(dp),Intent(in) :: x,Q2
If (Q2.le.0._dp) Then
  write(*,*) 'invalid renormalization scale!'
  return
end if
If (VerySmallM(x,Q2)) Then
   J_SMartin_1x = 0._dp
Else
   J_SMartin_1x = (x)*(zlog(Cmplx(x/Q2,0._dp,dp))-1._dp)
End If

! J_SMartin_1x = 0._dp

End Function J_SMartin_1x

Real(dp) Function J_SMartin_2x(x,y,Q2)
Implicit None
! x,y are supposed to be squared masses
! symmetric x <-> y
Real(dp),Intent(in) :: x,y,Q2
If (Q2.le.0._dp) Then
  write(*,*) 'invalid renormalization scale!'
  return
end if
If (VerySmallM(x,y,y,Q2).or.(VerySmallM(y,x,x,Q2))) Then
 J_SMartin_2x = 0
Else
 J_SMartin_2x = J_SMartin_1x(x,Q2)*J_SMartin_1x(y,Q2)
End if

! J_SMartin_2x = 0._dp

End Function J_SMartin_2x

Complex(dp) Function RFunction(x,y,z)
Implicit None
Complex(dp), Intent(in) :: x,y,z
Complex(dp) :: temp
temp = x**2+y**2+z**2-2._dp*x*y-2._dp*x*z-2._dp*y*z
RFunction = sqrt(temp)
! If (temp.lt.0._dp) Then
!   RFunction = Cmplx(0.,+sqrt(-temp),dp)
! Else
!   RFunction = Cmplx(sqrt(temp),0.,dp)
! End If
End Function RFunction

Complex(dp) Function I_SMartin(x_in,y_in,z_in,Q2)
Implicit None
! x_in,y_in,z_in are supposed to be squared masses
! symmetric under interchange of any two of these
Real(dp),Intent(in) :: x_in,y_in,z_in,Q2
Complex(dp) :: res,x,y,z
x = x_in
y = y_in
z = z_in
If (Q2.le.0._dp) Then
  write(*,*) 'invalid renormalization scale!'
  return
end if
Call SortArguments3C(x,y,z) ! sort x<=y<=z

If (SmallM(x,y,z,Q2)) Then
 I_SMartin = I_SMartin_ex(x_in,y_in,z_in,Q2)

! !  ! I(0,y,z)
! !  If (SmallM(y,x,z,Q2)) Then
! !   If (SmallM(z,y,x,Q2)) Then
! !    I_SMartin = ZeroC ! I(0,0,0)=0
! !   Else
! !    I_SMartin = -0.5_dp*z*(log(z/Q2))**2 &
! !       & + 2*z*log(z/Q2)-2.5_dp*z-pi**2*z/6._dp
! !    ! I(0,0,z)
! !   End If
! !  Else If (SmallD(z,y)) Then
! !   ! I(0,y,y)
! ! !   Write(*,*) Real(x,dp),Real(y,dp),Real(z,dp)
! !   I_SMartin = -y*(zlog(y/Q2))**2+4._dp*y*zlog(y/Q2)-5._dp*y
! ! ! I_SMartin = 0.
! !  Else
! ! !   Write(*,*) Real(x,dp),Real(y,dp),Real(z,dp)
! !   ! I(0,y,z), z>y
! !   I_SMartin = (z-y)*(CLi2_mod(y/z)-log(z/y)*log((z-y)/Q2) &
! !   & +0.5_dp*(log(z/Q2))**2-pi**2/6._dp)  &
! !   & -2.5_dp*(z+y)+2*y*log(y/Q2)+2*z*log(z/Q2) &
! !   & -z*log(z/Q2)*log(y/Q2)
! ! ! I_SMartin =0
! !  End If

! I_SMartin = 0._dp

Else ! I(x,y,z) general case

res = Cmplx(0._dp,0._dp)
res=res + 0.5_dp*(x-y-z)*zlog(y/Q2)*zlog(z/Q2)
res=res + 0.5_dp*(y-x-z)*zlog(x/Q2)*zlog(z/Q2)
res=res + 0.5_dp*(z-x-y)*zlog(x/Q2)*zlog(y/Q2)
res=res + 2._dp*x*zlog(x/Q2)
res=res + 2._dp*y*zlog(y/Q2)
res=res + 2._dp*z*zlog(z/Q2)
res=res - 2.5_dp*(x+y+z)
res=res - 0.5_dp*XiFunction(x,y,z)

! res = 0._dp

I_SMartin = res

End If

! I_SMartin = 0._dp

End Function I_SMartin

Complex(dp) Function XiFunction(x0,y0,z0)
Implicit None
Complex(dp), Intent(in) :: x0,y0,z0
Complex(dp) :: x,y,z
Complex(dp) :: res,R
x=x0
y=y0
z=z0
Call SortArguments3C(x,y,z)
R=Rfunction(x,y,z)
res = Cmplx(0.,0.,dp)
res = res + 2._dp*log((z+x-y-R)/(2._dp*z))*log((z+y-x-R)/(2._dp*z))
res = res - log(x/z)*log(y/z)
res = res - 2._dp*CLi2_mod((z+x-y-R)/(2._dp*z))
res = res - 2._dp*CLi2_mod((z+y-x-R)/(2._dp*z))
res = res + (pi**2)/3._dp
res = R*res
XiFunction = res
End Function XiFunction

!! sorts three real arguments a,b,c from small to large
!! and returns them
Subroutine SortArguments3(a,b,c)
Implicit None
Real(dp), Intent(inout) :: a,b,c
Real(dp) :: x,y,z
x=a ! assume a is smallest
y=b
z=c
If (b.le.a) Then
 x=b
 y=a
 z=c
End If
If (c.le.x) Then
 x=c
 y=a
 z=b
End If
! compare last two elements
If (z.le.y) Then
 a=x ! output
 b=z
 c=y
Else
 a=x !output
 b=y
 c=z
End If
End Subroutine SortArguments3

!! sorts three complex arguments a,b,c from small to large
!! and returns them
Subroutine SortArguments3C(a,b,c)
Implicit None
Complex(dp), Intent(inout) :: a,b,c
Complex(dp) :: x,y,z
x=a ! assume a is smallest
y=b
z=c
If (Abs(b).le.Abs(a)) Then
 x=b
 y=a
 z=c
End If
If (Abs(c).le.Abs(x)) Then
 x=c
 y=a
 z=b
End If
! compare last two elements
If (Abs(z).le.Abs(y)) Then
 a=x ! output
 b=z
 c=y
Else
 a=x !output
 b=y
 c=z
End If
End Subroutine SortArguments3C

 Complex(DP) Function CLI2_mod(Z)
 !--------------------------------------------------------------
 !     routine for the complex dilogarithm function
 !--------------------------------------------------------------
 Implicit None
  !-------
  ! input
  !-------
  Complex(DP), Intent(in) :: Z

  !-----------------
  ! local variables
  !-----------------
  Complex(DP) :: sumC, cy, cz, cz2, CLI2
  Real(dp) :: rz, iz, az
  Integer :: i1, jsgn, ipi12
  !---------------------------------------------
  ! Bernulli numbers / (2 n + 1)!
  !---------------------------------------------
  Real(dp), Parameter :: bf(20) = (/ - 1._dp/4._dp, + 1._dp/36._dp        &
       & , - 1._dp/36.e2_dp, + 1._dp/21168.e1_dp,  - 1._dp/108864.e2_dp   &
       & , + 1._dp/52690176.e1_dp,  - 691._dp/16999766784.e3_dp           &
       & , + 1._dp/1120863744.e3_dp,  - 3617._dp/18140058832896.e4_dp     &
       & , + 43867._dp/97072790126247936.e3_dp                            &
       & , - 174611._dp/168600109166641152.e5_dp                          &
       & , + 77683._dp/32432530090601152512.e4_dp                         &
       & , - 236364091._dp/4234560341829359173632.e7_dp                   &
       & , + 657931._dp/5025632054039239458816.e6_dp                      &
       & , - 3392780147._dp/109890470493622010006470656.e7_dp             &
       & , +172.3168255201_dp/2355349904102724211909.3102313472e6_dp      &
       & , -770.9321041217_dp/4428491985594062112714.2791446528e8_dp      &
       & ,  4.1576356446138997196178996207752266734882541595116e-29_dp    &
       & , -9.9621484882846221031940067024558388498548600173945e-31_dp    &
       & ,  2.3940344248961653005211679878937495629342791569329e-32_dp    &
       & /)

  rz = Real(z,dp)
  iz = Aimag(z)
  az = Sqrt(rz**2 + iz**2)
  !------------------
  ! special cases
  !------------------
  If (iz.Eq.0._dp) Then
    If (rz.gt.1._dp) Then
    CLi2_mod = Cmplx( Li2(rz), -pi*log(rz),dp)
    Else
    CLi2_mod = Cmplx( Li2(rz), 0.,dp)
    End If
   Return
  Else If (az.Lt.Epsilon(1._dp)) Then
   CLi2_mod = z
   Return
  End If
  !-----------------------------------------------------
  ! transformation to |z|<1, Re(z)<=0.5
  !-----------------------------------------------------
  If (rz.Le.0.5_dp) Then
   If (az.Gt.1._dp) Then
    cy = -0.5_dp * Log(-z)**2
    cz = - Log(1._dp - 1._dp / z)
    jsgn = -1
    ipi12 = -2
   Else !  (az.le.1._dp)
    cy = 0
    cz = - Log(1._dp - z)
    jsgn = 1
    ipi12 = 0
   End If
  Else ! rz.gt.0.5_dp
   If (az.Le.Sqrt(2*rz)) Then
    cz = - Log(z)
    cy = cz * Log(1._dp - z)
    jsgn = -1
    ipi12 = 2
   Else ! (az.gt.Sqrt(2*rz))
    cy = -0.5_dp * Log(-z)**2
    cz = - Log(1._dp - 1._dp / z)
    jsgn = -1
    ipi12 = -2
   End If
  End If
  !--------------------------------------------
  ! the dilogarithm
  !--------------------------------------------
  cz2 = cz**2
  sumC = cz2 * bf(20)
  Do i1=3,19
   sumC = cz2 * ( sumC + bf(22-i1) )
  End Do
  ! watch the powers of z
  sumC = cz + cz2 * (bf(1) + cz * (bf(2) + sumC))
  If (jsgn.Eq.1) Then
   CLi2 = sumC + cy
  Else
   CLi2 = - sumC + cy
  End If
  If (ipi12.Ne.0) Cli2 = Cli2 + ipi12 * pi2o12

  CLI2_mod = CLI2
 End Function CLI2_mod

!--------------------------------------------
! 2-Loop functions for different contributions
!
!--------------------------------------------

  !--------------------------------------------
  ! three scalars, two cubic couplings
  !--------------------------------------------
Complex(dp) Function Fep_SSS(x,y,z,Q2)
Implicit None
Real(dp), Intent(in) :: x,y,z,Q2
 Fep_SSS = -I_SMartin(x,y,z,Q2)
End Function Fep_SSS

  !--------------------------------------------
  ! two scalars, one quartic couplings
  !--------------------------------------------
Complex(dp) Function Fep_SS(x,y,Q2)
Implicit None
Real(dp), Intent(in) :: x,y,Q2
 Fep_SS = J_SMartin_2x(x,y,Q2)
End Function Fep_SS

  !--------------------------------------------
  ! one scalar, two fermions, two cubic couplings
  !--------------------------------------------
Complex(dp) Function Fep_FFS(x,y,z,Q2)
Implicit None
Real(dp), Intent(in) :: x,y,z,Q2
Complex(dp) :: res
res=Cmplx(0.,0.,dp)
res=res+J_SMartin_2x(x,y,Q2)-J_SMartin_2x(x,z,Q2)-J_SMartin_2x(y,z,Q2)
res=res+(x+y-z)*I_SMartin(x,y,z,Q2)
Fep_FFS=res
End Function Fep_FFS

  !--------------------------------------------
  ! like FFS, with mass insertion
  !--------------------------------------------
Complex(dp) Function Fep_FFSbar(x,y,z,Q2)
Implicit None
Real(dp), Intent(in) :: x,y,z,Q2
Fep_FFSbar=2*I_SMartin(x,y,z,Q2)
End Function Fep_FFSbar

  !--------------------------------------------
  ! two scalars, one vector, two cubic couplings
  !--------------------------------------------
Complex(dp) Function Fep_SSV(x,y,z,Q2)
Implicit None
Real(dp), Intent(in) :: x,y,z,Q2
Complex(dp) :: res
res=Cmplx(0.,0.,dp)

If (SmallM(z,x,y,Q2)) Then
  ! f_SSV(x,y,0) special case
 res=res+(x+y)**2+3._dp*(x+y)*I_SMartin(x,y,0._dp,Q2)
 res=res+3._dp*J_SMartin_2x(x,y,Q2)
 res=res-2._dp*x*J_SMartin_1x(x,Q2)
 res=res-2._dp*y*J_SMartin_1x(y,Q2)
Else
  ! f_SSV(x,y,z) general case
res = (-x**2-y**2-z**2+2._dp*(x*y+x*z+y*z))*I_SMartin(x,y,z,Q2)
res = res+I_SMartin(0._dp,x,y,Q2)*((x-y)**2)
res = res+(y-x-z)*J_SMartin_2x(x,z,Q2)+(x-y-z)*J_SMartin_2x(y,z,Q2)
res = res + z*J_SMartin_2x(x,y,Q2)
res = res/z
res = res+2._dp*(x+y-z/3._dp)*J_SMartin_1x(z,Q2)
End If

Fep_SSV = res
End Function Fep_SSV

  !--------------------------------------------
  ! one scalar, one vector, one quartic coupling
  !--------------------------------------------
Complex(dp) Function Fep_VS(x,y,Q2)
Implicit None
Real(dp), Intent(in) :: x,y,Q2
 Fep_VS = 3._dp*J_SMartin_2x(x,y,Q2)
End Function Fep_VS

  !--------------------------------------------
  ! one scalar, two vectors, two cubic couplings
  !--------------------------------------------
Complex(dp) Function Fep_VVS(x,y,z,Q2)
Implicit None
Real(dp), Intent(in) :: x,y,z,Q2
Real(dp) :: xt
Complex(dp) :: res
res=Cmplx(0.,0.,dp)
If (SmallM(x,y,z,Q2)) Then
 If (SmallM(y,x,z,Q2)) Then
  ! f_VSS(0,0,z)
  res=-3._dp*I_SMartin(0._dp,0._dp,z,Q2)
  res=res+1.5_dp*J_SMartin_1x(z,Q2)-0.25_dp*z
  Fep_VVS=res
  return
 Else
  xt=y
 End If
Else ! x>0
 If (SmallM(y,x,z,Q2)) Then
  xt=x
 Else ! x>0 and y>0 general case
   ! f_VSS(x,y,z)
 res=(-x**2-y**2-z**2)
 res=res-10._dp*x*y+2._dp*x*z+2._dp*y*z
 res=res*I_SMartin(x,y,z,Q2)
 res=res+(x-z)**2*I_SMartin(0._dp,x,z,Q2)
 res=res+(y-z)**2*I_SMartin( 0._dp,y,z,Q2)
 res=res-z**2*I_SMartin(0._dp,0._dp,z,Q2)
 res=res+(z-x-y)*J_SMartin_2x(x,y,Q2)
 res=res+y*J_SMartin_2x(x,z,Q2)
 res=res+z*J_SMartin_2x(y,z,Q2)
 res=res/(4._dp*x*y)
 res=res+0.5_dp*J_SMartin_1x(x,Q2)
 res=res+0.5_dp*J_SMartin_1x(y,Q2)
 Fep_VVS=res
 Return
 End If
End If
! case f_VSS(xt,0,z) with xt=x or y
res=xt/4._dp+0.5_dp*z
res=res-0.75_dp*z/xt*I_SMartin(0._dp,0._dp,z,Q2)
res=res+(0.75_dp*z/xt-9._dp/4._dp)*I_SMartin(0._dp,xt,z,Q2)
res=res+0.75_dp/xt*J_SMartin_2x(xt,z,Q2)
Fep_VVS=res
End Function Fep_VVS

  !--------------------------------------------
  ! two fermions, one vector, two cubic couplings
  !--------------------------------------------
Complex(dp) Function Fep_FFV(x,y,z,Q2)
Implicit None
Real(dp), Intent(in) :: x,y,z,Q2
Complex(dp) :: res
res=Cmplx(0.,0.,dp)
If (SmallM(z,x,y,Q2)) Then
 ! f_FFV(x,y,0)
 res=-(x+y)**2
 res=res+2._dp*x*J_SMartin_1x(x,Q2)
 res=res+2._dp*y*J_SMartin_1x(y,Q2)
 Fep_FFV=res
 Return
Else
 ! f_FFV(x,y,z) general case
 res=x**2+y**2-2._dp*z**2
 res=res-2._dp*x*y+x*z+y*z
 res=res*I_SMartin(x,y,z,Q2)
 res=res-(x-y)**2*I_SMartin(0._dp,x,y,Q2)
 res=res+(x-y-2._dp*z)*J_SMartin_2x(x,z,Q2)
 res=res+(y-x-2._dp*z)*J_SMartin_2x(y,z,Q2)
 res=res+2._dp*z*J_SMartin_2x(x,y,Q2)
 res=res/z
 res=res+2._dp*(-x-y+z/3._dp)*J_SMartin_1x(z,Q2)
 Fep_FFV=res
 Return
End If
End Function Fep_FFV

  !--------------------------------------------
  ! like FFV with mass insertion
  !--------------------------------------------
Complex(dp) Function Fep_FFVbar(x,y,z,Q2)
Implicit None
Real(dp), Intent(in) :: x,y,z,Q2
 Fep_FFVbar=6._dp*I_SMartin(x,y,z,Q2)
End Function Fep_FFVbar

  !--------------------------------------------
  ! sum of three diagrams:
  ! two vectors, one quartic coupling
  ! + three vectors, two cubic couplings
  ! + one vector, two ghosts, two cubic couplings
  !--------------------------------------------
Complex(dp) Function Fep_gauge(a,b,c,Q2)
Implicit None
Real(dp), Intent(in) :: a,b,c,Q2
Real(dp) :: x,y,z
 Complex(dp) :: res
 res=Cmplx(0.,0.,dp)
x=a
y=b
z=c
Call SortArguments3(x,y,z)
If (SmallM(x,y,z,Q2)) Then
 If (SmallM(y,x,z,Q2)) Then
  If (SmallM(z,y,x,Q2)) Then
   res=0._dp
  Else
  ! f_gauge(0,0,z)
  res=13*z*I_SMartin(0._dp,0._dp,z,Q2)
  res=res-71._dp/6._dp*z*J_SMartin_1x(z,Q2)
  res=res+19._dp/4._dp*z**2
  End If
 Else
 ! f_gauge(0,y,z)
 res=43._dp*z**2*y+43._dp*z*y**2-7._dp*z**3-7._dp*y**3
 res=res*I_SMartin(0._dp,z,y,Q2)
 res=res+(2._dp*y+7._dp*z)*z**2*I_SMartin(0._dp,0._dp,z,Q2)
 res=res+(2._dp*z+7._dp*y)*y**2*I_SMartin(0._dp,0._dp,y,Q2)
 res=res+(34*z*y-7*z**2-7*y**2)*J_SMartin_2x(z,y,Q2)
 res=res/(4._dp*z*y)
 res=res+3*z**2+3*y**2+5.5_dp*z*y
 res=res-25._dp/3._dp*z*J_SMartin_1x(z,Q2)
 res=res-25._dp/3._dp*y*J_SMartin_1x(y,Q2)
 res=res+5*y*J_SMartin_1x(z,Q2)
 res=res+5*z*J_SMartin_1x(y,Q2)
 End If
Else
! f_gauge(x,y,z) general case
 res=fgauge_aux(x,y,z,Q2)
 res=res+fgauge_aux(y,x,z,Q2)
 res=res+fgauge_aux(z,y,x,Q2)
End If

Fep_gauge=res
End Function Fep_gauge

Complex(dp) Function fgauge_aux(x,y,z,Q2)
 Real(dp), Intent(in) :: x,y,z,Q2
 Complex(dp) :: res
 res=Cmplx(0.,0.,dp)
 res=-x**4-8*x**3*(y+z)+32*x**2*y*z+18*y**2*z**2
 res=res*I_SMartin(x,y,z,Q2)
 res=res+(y-z)**2*(y**2+10*y*z+z**2)*I_SMartin(0._dp,y,z,Q2)
 res=res+x**2*(2*y*z-x**2)*I_SMartin(0._dp,0._dp,x,Q2)
 res=res+(x**2-9*(y**2+z**2-x*y-x*z)+14*y*z)*x*J_SMartin_2x(y,z,Q2)
 res=res+(22*(y+z)-40*x/3._dp)*x*y*z*J_SMartin_1x(x,Q2)
 res=res/(4._dp*x*y*z)
 fgauge_aux=res
End Function fgauge_aux

Logical Function SmallDC(x,y)
Implicit None
Complex(dp), Intent(in) :: x,y
If (Max(Abs(x),Abs(y)).lt.epsD) Then
 SmallDC = .True.
Else If ((Abs(x-y)/Max(Abs(x),Abs(y))).le.epsD) Then
 SmallDC = .True.
Else
 SmallDC = .False.
End if
End Function SmallDC

Logical Function SmallDR(x,y)
Implicit None
Real(dp), Intent(in) :: x,y
If (Max(Abs(x),Abs(y)).lt.epsD) Then
 SmallDR = .True.
Else If ((Abs(x-y)/Max(Abs(x),Abs(y))).le.epsD) Then
 SmallDR = .True.
Else
 SmallDR = .False.
End if
End Function SmallDR

Logical Function SmallMR(x,y,z,Q)
Implicit None
Real(dp), Intent(in) :: x,y,z
Real(dp), Intent(in) :: Q
If (Max(Max(Abs(x),Abs(y)),Abs(z)).le.epsM) Then
 SmallMR = .True.
Else If ((Abs(x)/Max(Q,Max(Abs(y),Abs(z)))).le.epsM) Then
 SmallMR = .True.
Else
 SmallMR = .False.
End if
End Function SmallMR

Logical Function SmallMC(x,y,z,Q)
Implicit None
Complex(dp), Intent(in) :: x,y,z
Real(dp), Intent(in) :: Q
If (Max(Max(Abs(x),Abs(y)),Abs(z)).le.epsM) Then
 SmallMC = .True.
Else If ((Abs(x)/Max(Q,Max(Abs(y),Abs(z)))).lt.epsM) Then
 SmallMC = .True.
Else
 SmallMC = .False.
End if
End Function SmallMC

Logical Function SmallMR2(x,y)
Implicit None
Real(dp), Intent(in) :: x,y
If (Max(Abs(x),Abs(y)).le.epsM) Then
 SmallMR2 = .True.
Else If (Abs(x/y).le.epsM) Then
 SmallMR2 = .True.
Else
 SmallMR2 = .False.
End if
End Function SmallMR2

Logical Function SmallMC2(x,y)
Implicit None
Complex(dp), Intent(in) :: x,y
If (Max(Abs(x),Abs(y)).le.epsM) Then
 SmallMC2 = .True.
Else If (Abs(x/y).le.epsM) Then
 SmallMC2 = .True.
Else
 SmallMC2 = .False.
End if
End Function SmallMC2

End Module EffPotFunctions


text:
subroutine Poisson_2D_double_precision(Lx,Ly,Nnx,Nny,lambda,f,chii)

! Include modules defined by mkl_poisson.f90 and mkl_dfti.f90 header files
! use mkl_poisson

implicit none
integer, intent(in) :: Nnx,Nny
real(8), intent(in) :: Lx,Ly,lambda
real(8), dimension(Nnx,Nny), intent(inout)  :: f
real(8), dimension(Nnx,Nny), intent(out) :: chii

double precision pi
parameter(pi=3.14159265358979324D0)

integer nx, ny, stat
integer AA
DOUBLE PRECISION ptrb
double precision ax, bx, ay, by, q
double precision bd_ax(Nny), bd_bx(Nny), bd_ay(Nnx), bd_by(Nnx)
!double precision dpar(13*(Nnx-1)/2+7)
DOUBLE PRECISION, allocatable,dimension (:) :: workarea
!type(DFTI_DESCRIPTOR), pointer :: xhandle
character(4) BCtype

nx=Nnx-1
ny=Nny-1
AA = (1./log(2.))*log(1.0 * Nnx)
ALLOCATE(workarea(4*(Nnx) + (13 + nx)*(Nny)))
! Defining the rectangular domain 0<x<Lx, 0<y<Ly for 2D Poisson Solver
ax=0.0D0
bx=Lx
ay=0.0D0
by=Ly

!*******************************************************************************
! Setting the coefficient q to 0.
! Note that this is the way to use Helmholtz Solver to solve Poisson problem!
!*******************************************************************************
q=-(lambda*lambda)

!******************************************************************************
! NOTE THE SIGN CONVENTION USED FOR THE POISSON EQUATION!!!!!!!!!!!!!!!!!!!!!!!
!******************************************************************************

!******************************************************************************
! BOUNDARY CONDITIONS
!******************************************************************************
BCtype = 'DDDD'
bd_ax(:) = 0.0D0
bd_bx(:) = 0.0D0
bd_ay(:) = 0.0D0
bd_by(:) = 0.0D0

f(1,:) = -chii(1,:)
f(Nnx,:) = -chii(Nnx,:)
f(:,1) = -chii(:,1)
f(:,Nny) = -chii(:,Nny)
! NOTE:
!print *, ' The SOLUTION depends NOT on boundary values of the rhs for Dirichlet !!!!!!!!!!!!!!'
!f(1,:) = 1.d10
!write(*,*) f(1,:)
!print *, ' The SOLUTION depends on boundary values of the rhs for Neumann !!!!!!!!!!!!!!!!!!!!'
!f(:,1) = 1.d10
!write(*,*) f(:,1)
f = -1.0d0 * f

ptrb = 0.0d0
CALL HWSCRT(&
&    ax, bx, nx, 1, bd_ax, bd_bx, &
&    ay, by, ny, 1, bd_ay, bd_by, &
&    q, &
&    f, nx+1,&
&    ptrb, stat, workarea &
&)

!print *, f
DEALLOCATE(workarea)
! Now we can use xhandle to solve another 2D Poisson problem
chii = f

! Free MKL memory if any was allocated
!call mkl_free_buffers

end

subroutine mkl_set_dynamic(input)
implicit none
integer, intent(in) :: input

end subroutine

subroutine mkl_set_num_threads(input)
implicit none
integer, intent(in) :: input

end subroutine

max_stars_repo_name: ipelupessy/omuse

max_issues_repo_name: ipelupessy/omuse

max_forks_repo_name: ipelupessy/omuse


module exe_frontend

use, intrinsic :: iso_fortran_env, only : compiler_version, stderr=>error_unit, compiler_options
use phys_consts, only : wp
use config, only : get_compiler_vendor, gemini_cfg, read_configfile
use hwloc_ifc, only : get_cpu_count
use filesystem, only : parent, file_name, assert_is_dir, expanduser, suffix
use timeutils, only : date_filename,dateinc

implicit none (type, external)
private
public :: clean_output, cli_parser, get_Ncpu, help_gemini_bin, help_gemini_run, help_magcalc_bin, help_magcalc_run

contains

subroutine cli_parser(plan, Ncpu, path, exe, mpiexec, extra)

logical, intent(out) :: plan
integer, intent(out) :: Ncpu
character(:), allocatable, intent(out) :: path, exe, mpiexec, extra

character(1000) :: buf
integer :: argc, i, j, ierr, L

argc = command_argument_count()

call get_command_argument(1, buf, status=i)
if (i/=0) call help_run()

if(buf(1:1) == '-') then
!! not running sim, checking parameters
  select case(buf)
  case ('-h', '-help')
    call help_run()
  case ('-compiler')
    print '(A)', get_compiler_vendor()
  case ('-compiler_version')
    print '(A)', compiler_version()
  case ('-compiler_options')
    print '(A)', compiler_options()
  case ('-git')
    print '(A)', "@git_rev@"
  case ('-features')
    print '(A)', "@gemini_features@"
  case default
    write(stderr,*) "Gemini3D: unknown option: ", trim(buf)
    call help_run()
  end select

  stop
endif

!> simulation data directory
path = trim(buf)
call assert_is_dir(path)

plan = .false.
extra = ""
Ncpu = 0

do i = 2, argc
  call get_command_argument(i, buf)
  if(buf(1:1) /= "-") cycle  !< assume previous option was a -flag

  select case (buf)

  case ('-n')
    call get_command_argument(i+1, buf, length=L, status=ierr)
    if(ierr /= 0 .or. L==0 .or. buf(1:1) == "-") error stop trim(buf) // " -n missing parameter"
    read(buf, '(I6)') Ncpu
  case ('-exe')
    call get_command_argument(i+1, buf, length=L, status=ierr)
    if(ierr /= 0 .or. L==0 .or. buf(1:1) == "-") error stop trim(buf) // " -exe missing parameter"
    exe = find_exe(trim(buf))
  case ('-mpiexec')
    call get_command_argument(i+1, buf, length=L, status=ierr)
    if(ierr /= 0 .or. L==0 .or. buf(1:1) == "-") error stop "-mpiexec was specified without an executable path"
    mpiexec = find_mpiexec(trim(buf))
  case ('-plan')
    plan = .true.
  !> options passed to child executable
  case ('-dryrun', '-debug', '-nooutput')
    !! flags with no parameters
    extra = extra // ' ' // trim(buf)
  case ('-out_format')
    !! flags with one parameter
    extra = extra // ' ' // trim(buf)
    call get_command_argument(i+1, buf, length=L, status=ierr)
    if(ierr /= 0 .or. L==0 .or. buf(1:1) == "-") error stop trim(buf) // " -out_format missing parameter"
    extra = extra // ' ' // trim(buf)
  case ('-manual_grid')
    !! flags with two parameters
    extra = extra // ' ' // trim(buf)
    do j = 1,2
      call get_command_argument(i+j, buf, length=L, status=ierr)
      if(ierr /= 0 .or. L==0 .or. buf(1:1) == "-") error stop trim(buf) // " -manual_grid expected two parameters"
      extra = extra // ' ' // trim(buf)
    enddo
  case ('-start_time', '-end_time')
    !! flags with four parameters
    extra = extra // ' ' // trim(buf)
    do j = 1,4
      call get_command_argument(i+j, buf, length=L, status=ierr)
      if(ierr /= 0 .or. L==0 .or. buf(1:1) == "-") error stop trim(buf) // " -start_time expected four parameters"
      extra = extra // ' ' // trim(buf)
    enddo
  case default
    error stop "Gemini3D: unknown option: " // trim(buf)
  end select
end do

if(.not.allocated(exe)) exe = find_exe("")
if(.not.allocated(mpiexec)) mpiexec = find_mpiexec("")

end subroutine cli_parser

integer function get_Ncpu(N) result(Ncpu)

integer, intent(in) :: N
integer :: i
character(1000) :: buf

if (N > 0) then
  Ncpu = N
  print '(A,I0)', "run: CLI specified CPU count: ", Ncpu
  return
endif

call get_environment_variable("GEMINI_CPU", buf, status=i)
if (i==0) then
  read(buf,'(I6)', iostat=i) Ncpu
  if (i==0) print '(A,I0)', "run: GEMINI_CPU CPU job count: ", Ncpu
endif
if(i/=0) then
  call get_environment_variable("NSLOTS", buf, status=i)
  if (i==0) then
    read(buf,'(I6)', iostat=i) Ncpu
    if (i==0) print '(A,I0)', "run: SGE CPU job count: ", Ncpu
  endif
endif
if(i/=0) then
  call get_environment_variable("SLURM_NTASKS", buf, status=i)
  if(i==0) then
    read(buf, '(I6)', iostat=i) Ncpu
    if(i==0) print '(A,I0)', "run: SLURM CPU job count: ", Ncpu
  endif
endif
if (i/=0) then
  Ncpu = get_cpu_count()

  print '(A,I0)', 'gemini3d.run: detected CPU count: ', Ncpu
endif

end function get_Ncpu

function find_exe(name) result(exe)

character(*), intent(in) :: name

character(:), allocatable :: exe, work
logical :: exists
character(1000) :: buf

if(len_trim(name) > 0) then
  inquire(file=name, exist=exists)
  if(exists) then
    exe = name
    return
  endif
  inquire(file=name // '.exe', exist=exists)
  if(exists) then
    exe = name // '.exe'
    return
  endif
endif

call get_command_argument(0, buf)

if(len_trim(name) == 0) then
  if (index(buf, "gemini3d.run.debug") > 0) then
    work = "gemini.bin.debug"
  elseif (index(buf, "gemini3d.run") > 0) then
    work = "gemini.bin"
  elseif (index(buf, "magcalc.run.debug") > 0) then
    work = "magcalc.bin.debug"
  elseif (index(buf, "magcalc.run") > 0) then
    work = "magcalc.bin"
  else
    error stop "frontend: find_exe: please specify the program name you seek"
  endif

  work = parent(buf) // '/' // work

  inquire(file=work, exist=exists)
  if (exists) then
    exe = work
    return
  else
    inquire(file=work // '.exe', exist=exists)
    if(exists) then
      exe = work // '.exe'
      return
    endif
  endif
endif

error stop "gemini3d.run: did not find " // exe // " from " // name // &
  " : please specify path to MPI runnable executable with option 'gemini3d.run -exe path/to/my.bin'"

end function find_exe

function find_mpiexec(exe) result(mpiexec)

character(*), intent(in) :: exe
character(:), allocatable :: mpiexec

character(1000) :: buf
integer :: i, L

if(len_trim(exe) > 0) then
  if(check_mpiexec(expanduser(exe))) mpiexec = expanduser(exe)
else
  call get_environment_variable("MPI_ROOT", buf, length=L, status=i)
  if (i==0 .and. L>0) then
    if(check_mpiexec(expanduser(buf) // "/bin/mpiexec")) mpiexec = expanduser(buf) // "/bin/mpiexec"
  endif
endif

if(.not.allocated(mpiexec)) mpiexec = "mpiexec"

end function find_mpiexec

logical function check_mpiexec(exe) result(ok)

character(*), intent(in) :: exe

if(exe(1:1) == "-") error stop "gemini3d.run: -mpiexec was missing an executable, got: " // exe

inquire(file=exe, exist=ok)

if(ok) return

write(stderr,"(A,/,A)") "gemini3d.run: MPIexec file not found " // exe, &
"If simulation hangs or operates incorrectly, specify -mpiexec option or set MPI_ROOT environment variable."

end function check_mpiexec

subroutine help_gemini_bin() bind(C)

print '(/,A,/)', 'GEMINI-3D: gemini.bin ' // "@git_rev@"
print '(A)', 'by Matthew Zettergren'
print '(A)', 'GLOW and auroral interfaces by Guy Grubbs'
print '(A)', 'build system and software engineering by Michael Hirsch'
print '(A)', 'Compiler vendor: '// get_compiler_vendor()
print '(A)', 'Compiler version: ' // compiler_version()
print '(/,A,/)', 'the first and only positional argument is simulation output directory.'
print '(A)', 'Optional arguments:'
print '(a,t25,a)', '-dryrun', 'allows quick check of first time step'
print '(a,t25,a)', '-manual_grid lx2 lx3', 'defines the number of MPI processes along x2 and x3.'
print '(t25,a)', '  If -manual_grid is not specified, the MPI processes are auto-assigned along x2 and x3.'
stop 'EOF: gemini.bin'

end subroutine help_gemini_bin

subroutine help_run()

character(1000) :: buf

call get_command_argument(0, buf)
if(index(buf, "gemini3d.run") > 0) call help_gemini_run()
if(index(buf, "magcalc.run") > 0) call help_magcalc_run()

error stop "help_run: unknown runner"

end subroutine help_run

subroutine help_gemini_run()

print '(/,a,/)', 'GEMINI-3D: gemini3d.run ' // "@git_rev@"
print '(a)', 'Compiler vendor: '// get_compiler_vendor()
print '(a)', 'Compiler version: ' // compiler_version()
print '(/,a,/)', 'the first and only positional argument is simulation output directory.'
print '(a)', 'Optional arguments:'
print '(a,t20,a)', '-plan', 'print MPI partition x2,x3 for given CPU count'
print '(a,t20,a)', '-dryrun', 'allows quick check of first time step'
print '(a,t20,a)', '-n', 'manually specify number of MPI images (default: auto-calculate MPI grid partitioning)'
print '(a,t20,a)', '-exe', 'specify path to gemini.bin (default: same directory as gemini3d.run)'
print '(a,t20,a)', '-mpiexec', 'specify path to mpiexec'
print '(a,t20,a)', '-compiler', 'tell the compiler family e.g. GNU, Intel, IntelLLVM.  This allows knowing what compiler was used.'
print '(a,t20,a)', '-compiler_version', 'like -compiler, and also tell the Fortran compiler version.'
print '(a,t20,a)', '-git', 'print git revision it was built from. This is not perfect, to be sure use fresh build directory.'
print '(a,t20,a)', '-compiler_options', 'print compiler flags used to build the executable.'
print '(a,t20,a)', '-features', 'print Gemini3D external features enabled as a string, machine/human-readable.'
stop 'EOF: gemini3d.run'

end subroutine help_gemini_run

subroutine help_magcalc_bin()

print '(/,A,/)', 'GEMINI-3D: magcalc.bin ' // "@git_rev@"
print '(A)', 'Compiler vendor: '// get_compiler_vendor()
print '(A)', 'Compiler version: ' // compiler_version()
print '(/,A)', 'First two positional arguments must be input directory and fieldpoint file. Example:'
print '(/,A,/)', 'mpiexec -n 4 build/magcalc.bin test2d_fang test2d_fang/fieldpoint'
print '(A)', 'Optional arguments:'
print '(A)', '-dryrun option allows quick check of first time step'
stop 'EOF: magcalc.bin'

end subroutine help_magcalc_bin

subroutine help_magcalc_run()

print '(/,A,/)', 'GEMINI-3D: magcalc.run ' // "@git_rev@"
print '(A)', 'Compiler vendor: '// get_compiler_vendor()
print '(A)', 'Compiler version: ' // compiler_version()
print '(/,A,/)', 'the first and only positional argument is simulation output directory.'
print '(A)', 'Optional arguments:'
print '(A)', '-plan  print MPI partition x2,x3 for given CPU count'
print '(A)', '-dryrun    allows quick check of first time step'
print '(A)', '-n   manually specify number of MPI images (default auto-calculate)'
print '(A)', '-exe   path to magcalc.bin'
print '(A)', '-mpiexec   path to mpiexec'
stop 'EOF: magcalc.run'

end subroutine help_magcalc_run

subroutine clean_output(path)

character(*), intent(in) :: path

type(gemini_cfg) :: cfg
integer, dimension(3) :: ymd
real(wp) :: UTsec
character(:), allocatable :: fn
logical :: exists

cfg%outdir = path
cfg%infile = path // '/inputs/config.nml'
inquire(file=cfg%infile, exist=exists)
if(.not.exists) error stop 'gemini3d.run: not a file: ' // cfg%infile

call read_configfile(cfg)

ymd = cfg%ymd0
UTsec = cfg%UTsec0

fn = date_filename(cfg%outdir, ymd, UTsec) // suffix(cfg%indatsize)

do
  !! new filename, add the 1 if it is the first
  fn = date_filename(cfg%outdir, ymd, UTsec) // suffix(cfg%indatsize)

  inquire(file=fn, exist=exists)
  if ( .not. exists ) exit
  !! last output file
  print *, 'delete: ', fn
  call unlink(fn)

  !! next time
  call dateinc(cfg%dtout, ymd,UTsec)
end do

end subroutine clean_output

subroutine unlink(path)
character(*), intent(in) :: path
integer :: i
logical :: e

inquire(file=path, exist=e)
if (.not.e) return

open(newunit=i, file=path, status='old')
close(i, status='delete')
end subroutine unlink

end module exe_frontend


text: -f uvm/uvm_express/cov.f

ext: f

max_stars_repo_name: umarcor/svunit

max_issues_repo_name: umarcor/svunit

max_forks_repo_name: umarcor/svunit


text:
*DECK CS1S2
      SUBROUTINE CS1S2 (ZR, S1, S2, NZ, ASCLE, ALIM, IUF)
C***BEGIN PROLOGUE  CS1S2
C***SUBSIDIARY
C***PURPOSE  Subsidiary to CAIRY and CBESK
C***LIBRARY   SLATEC
C***TYPE      ALL (CS1S2-A, ZS1S2-A)
C***AUTHOR  Amos, D. E., (SNL)
C***DESCRIPTION
C
C     CS1S2 TESTS FOR A POSSIBLE UNDERFLOW RESULTING FROM THE
C     ADDITION OF THE I AND K FUNCTIONS IN THE ANALYTIC CON-
C     TINUATION FORMULA WHERE S1=K FUNCTION AND S2=I FUNCTION.
C     ON KODE=1 THE I AND K FUNCTIONS ARE DIFFERENT ORDERS OF
C     MAGNITUDE, BUT FOR KODE=2 THEY CAN BE OF THE SAME ORDER
C     OF MAGNITUDE AND THE MAXIMUM MUST BE AT LEAST ONE
C     PRECISION ABOVE THE UNDERFLOW LIMIT.
C
C***SEE ALSO  CAIRY, CBESK
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   830501  DATE WRITTEN
C   910415  Prologue converted to Version 4.0 format.  (BAB)
C***END PROLOGUE  CS1S2
      COMPLEX CZERO, C1, S1, S1D, S2, ZR
      REAL AA, ALIM, ALN, ASCLE, AS1, AS2, XX
      INTEGER IUF, NZ
      DATA CZERO / (0.0E0,0.0E0) /
C***FIRST EXECUTABLE STATEMENT  CS1S2
      NZ = 0
      AS1 = ABS(S1)
      AS2 = ABS(S2)
      AA = REAL(S1)
      ALN = AIMAG(S1)
      IF (AA.EQ.0.0E0 .AND. ALN.EQ.0.0E0) GO TO 10
      IF (AS1.EQ.0.0E0) GO TO 10
      XX = REAL(ZR)
      ALN = -XX - XX + ALOG(AS1)
      S1D = S1
      S1 = CZERO
      AS1 = 0.0E0
      IF (ALN.LT.(-ALIM)) GO TO 10
      C1 = CLOG(S1D) - ZR - ZR
      S1 = CEXP(C1)
      AS1 = ABS(S1)
      IUF = IUF + 1
   10 CONTINUE
      AA = MAX(AS1,AS2)
      IF (AA.GT.ASCLE) RETURN
      S1 = CZERO
      S2 = CZERO
      NZ = 1
      IUF = 0
      RETURN
      END


text:
!##############################################################################
!# ****************************************************************************
!# <name> elementbase </name>
!# ****************************************************************************
!#
!# <purpose>
!# This module contains a set of constant definitions and structures which
!# are used by the different element modules.
!#
!# </purpose>
!##############################################################################

module elementbase

!$ use omp_lib
  use basicgeometry
  use derivatives
  use fsystem
  use perfconfig
  use transformation
  use triangulation

  implicit none

  private

!<constants>

!<constantblock description="Element evaluation tags. Defines the basis information that is \
!  necessary to evaluate an element. All tags define bits of a bitfield and can be \
!  combined by an OR command.">

  ! Calculate the coordinates of teh points that form the element.
  ! (Usually the corners of the element.)
  integer(I32), parameter, public :: EL_EVLTAG_COORDS       = 2_I32**0

  ! Prepare the coordinates on the reference element of the points where to evaluate.
  ! A set of coordinates on one element must be specified in the preparation routine.
  ! These coordinates are transferred to all elements where to evaluate.
  ! For the preparation routine on one element: Transfer the source coordinates
  ! on the reference element to the evaluation structure.
  integer(I32), parameter, public :: EL_EVLTAG_REFPOINTS    = 2_I32**1

  ! Calculate the real coordinates of the points where to evaluate.
  ! The real coordinates are calculated based on the corresponding mapping from the
  ! reference to the real element.
  integer(I32), parameter, public :: EL_EVLTAG_REALPOINTS   = 2_I32**2

  ! Calculate the Jacobian matrix of the points where to evaluate
  integer(I32), parameter, public :: EL_EVLTAG_JAC          = 2_I32**3

  ! Calculate the Jacobian determinant in the points where to evaluate
  integer(I32), parameter, public :: EL_EVLTAG_DETJ         = 2_I32**4

  ! Calculate the twist indices
  integer(I32), parameter, public :: EL_EVLTAG_TWISTIDX     = 2_I32**5

!</constantblock>

!</constants>

!<types>

!<typeblock>

  ! This structure collects information that is necessary by an element
  ! to be evaluated on a single point on a single cell.
  type t_evalElement

    ! Array with coordinates of the corners that form the real element.
    ! DIMENSION(#space dimensions)
    !  Dcoords(1)=x-coordinate,
    !  Dcoords(2)=y-coordinate.
    !  Dcoords(3)=z-coordinate (only in 3D).
    real(DP), dimension(NDIM3D,TRIA_MAXNVE) :: Dcoords

    ! Values of the Jacobian matrix that defines the mapping between the
    ! reference element and the real elements. For every point i (e.g. in 2D):
    !  Djac(1) = J_i(1,1)
    !  Djac(2) = J_i(2,1)
    !  Djac(3) = J_i(1,2)
    !  Djac(4) = J_i(2,2)
    ! ...
    ! Remark: Only used for calculating derivatives; can be set to 0.0
    ! when derivatives are not used.
    real(DP), dimension(NDIM3D*NDIM3D) :: Djac

    ! Determinant of the mapping from the reference element to the real
    ! elements in the point to be evaluated.
    ! Remark: Only used for calculating derivatives; can be set to 1.0
    ! when derivatives are not needed. Must not be set to 0.0!
    real(DP) :: ddetj

    ! Array with coordinates of the points where to evaluate.
    ! The coordinates are expected on the reference element in the corresponding
    ! coordinate system. E.g. for real coordinates, this means:
    !  Dpoints(1)=x-coordinates,
    !  Dpoints(2)=y-coordinates,
    !  Dpoints(3)=z-coordinates (in 3D).
    ! but the entries here can also be e.g. in barycentric coordinates,
    ! if the coordinate system on the reference element prescribes that.
    real(DP), dimension(TRAFO_MAXDIMREFCOORD) :: DpointRef

    ! Array with coordinates of the points where to evaluate.
    ! The coordinates are expected on the real element.
    ! It is assumed that:
    !  Dpoints(1)=x-coordinates,
    !  Dpoints(2)=y-coordinates.
    !  Dpoints(3)=z-coordinates (only 3d).
    real(DP), dimension(NDIM3D) :: DpointReal

    ! Twist index entry of the element.
    integer(I32) :: itwistIndex

  end type

  public :: t_evalElement
!</typeblock>

!<typeblock>

  ! This structure collects information that is necessary by an element
  ! to be evaluated on a set of points on a set of element cells.
  type t_evalElementSet

    ! Number of points on every element where to evalate the basis functions.
    integer :: npointsPerElement = 0

    ! Number of elements, this element set consists of
    integer  :: nelements = 0

    ! Array with coordinates of the corners that form the real element.
    ! DIMENSION(#space dimensions,NVE,nelements)
    !  Dcoords(1,.,.)=x-coordinates,
    !  Dcoords(2,.,.)=y-coordinates.
    !  Dcoords(3,.,.)=z-coordinates (only in 3D).
    ! furthermore:
    !  Dcoords(:,i,.) = Coordinates of vertex i
    ! furthermore:
    !  Dcoords(:,:,j) = Coordinates of all corner vertices of element j
    real(DP), dimension(:,:,:), pointer :: p_Dcoords => null()

    ! Values of the Jacobian matrix that defines the mapping between the
    ! reference element and the real elements. For every point i (e.g. in 2D):
    !  Djac(1,i,.) = J_i(1,1,.)
    !  Djac(2,i,.) = J_i(2,1,.)
    !  Djac(3,i,.) = J_i(1,2,.)
    !  Djac(4,i,.) = J_i(2,2,.)
    ! Remark: Only used for calculating derivatives; can be set to 0.0
    ! when derivatives are not used.
    !  Djac(:,:,j) refers to the determinants of the points of element j.
    real(DP), dimension(:,:,:), pointer :: p_Djac => null()

    ! Determinant of the mapping from the reference element to the real
    ! elements for every of the npointsPerElement points on all the elements.
    !  Ddetj(i,.) = Determinant of point i
    !  Ddetj(:,j) = determinants of all points on element j
    ! Remark: Only used for calculating derivatives; can be set to 1.0
    ! when derivatives are not needed. Must not be set to 0.0!
    real(DP), dimension(:,:), pointer :: p_Ddetj => null()

    ! Array with coordinates of the points where to evaluate.
    ! The coordinates are expected on the reference element.
    ! It is assumed that:
    !  Dpoints(1,.)=x-coordinates,
    !  Dpoints(2,.)=y-coordinates.
    !  Dpoints(3,.)=y-coordinates.
    ! furthermore:
    !  Dpoints(:,i,.) = Coordinates of point i
    ! furthermore:
    !  Dpoints(:,:,j) = Coordinates of all points on element j
    real(DP), dimension(:,:,:), pointer :: p_DpointsRef => null()

    ! Array with coordinates of the points where to evaluate.
    ! The coordinates are expected on the real element.
    ! It is assumed that:
    !  Dpoints(1,.)=x-coordinates,
    !  Dpoints(2,.)=y-coordinates.
    !  Dpoints(3,.)=z-coordinates (only 3d).
    ! furthermore:
    !  Dpoints(:,i,.) = Coordinates of point i
    ! furthermore:
    !  Dpoints(:,:,j) = Coordinates of all points on element j
    real(DP), dimension(:,:,:), pointer :: p_DpointsReal => null()

    ! Twist index array of the elements.
    ! Array with DIMENSION(nelements)
    integer(I32), dimension(:), pointer :: p_ItwistIndex => null()

    ! .TRUE., if the array with coordinates on the reference element
    ! is maintained by the caller; prevents release of memory in the
    ! cleanup routines.
    logical :: bforeignPointsRef = .false.

    ! .TRUE., if the array with coordinates on the real element
    ! is maintained by the caller; prevents release of memory in
    ! cleanup routines.
    logical :: bforeignPointsReal = .false.

    ! .TRUE., if the array p_Dcoords is maintained by the caller; prevents
    ! release of memory in cleanup routines.
    logical :: bforeignCoords = .false.

    ! Pointer to a performance pointer
    type(t_perfconfig), pointer :: p_rperfconfig => null()

  end type

  public :: t_evalElementSet

!</typeblock>

!</types>

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
      program main
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C                                                                      %
C Copyright (C) 1996, The Board of Trustees of the Leland Stanford     %
C Junior University.  All rights reserved.                             %
C                                                                      %
C The programs in GSLIB are distributed in the hope that they will be  %
C useful, but WITHOUT ANY WARRANTY.  No author or distributor accepts  %
C responsibility to anyone for the consequences of using them or for   %
C whether they serve any particular purpose or work at all, unless he  %
C says so in writing.  Everyone is granted permission to copy, modify  %
C and redistribute the programs in GSLIB, but only under the condition %
C that this notice and the above copyright notice remain intact.       %
C                                                                      %
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c-----------------------------------------------------------------------
c
c                      Univariate Transformation
c                      *************************
c
c Reads in a reference distribution and a number of other distributions
c and then transforms the values in each of the second distributions
c such that their histograms match that of the reference distribution.
c
c
c
c INPUT/OUTPUT Parameters:
c
c   ivtype      variable type (1=continuous, 0=categorical)
c   distin      file with reference distribution
c   ivr,iwt     columns for variable and weight(0=none)
c   datafl      file with uncorrected distributions
c   ivr,iwt     columns for variable and weight(0=none)
c   tmin,tmax   trimming limits
c   outfl       file for revised distributions
c   nsim        size to transform, number of realizations
c   nx, ny, nz  size of categorical variable realizations to transform
c   wx, wy, wz  window size for breaking ties
c   nxyz        size to of continuous variable data set to transform
c   zmin,zmax   minimum and maximum data values
c   ltail,ltpar lower tail: option, parameter
c   utail,utpar upper tail: option, parameter
c   ldata       honor local data (1=yes, 0=no)
c   localfl     file with estimation variance
c   ikv         column number
c   wtfac       control parameter
c
c
c
c AUTHOR: Clayton V. Deutsch                             DATE: 1989-1999
c-----------------------------------------------------------------------
      use       msflib
      parameter(MAXCAT=24,MV=500,EPSLON=1.0e-12,VERSION=2.905)

      character distin*512,datafl*512,outfl*512,localfl*512,str*512
      real      catcdf(MAXCAT),var(MV),ltpar,utpar
      integer   ltail,utail,nx,ny,nz,wx,wy,wz,category(MAXCAT),test
      logical   testfl
c
c Declare dynamic arrays:
c
      real, allocatable :: dcdf(:),dvr(:),indx(:),fuzzcat(:),
     +                     rcdf(:),rvr(:)
c
c For the random number generator:
c
      parameter(KORDEI=12,MAXOP1=KORDEI+1,MAXINT=2**30)
      real*8  acorni
      common /iaco/ ixv(MAXOP1)
      data   ixv/MAXOP1*0.0/
      data   lin/1/,lout/2/,lkv/3/
c
c Note VERSION number before anything else:
c
      write(*,9999) VERSION
 9999 format(/' TRANS Version: ',f5.3/)
c
c Get the name of the parameter file - try the default name if no input:
c
      do i=1,512
            str(i:i) = ' '
      end do
      call getarg(1,str)
      if(str(1:1).eq.' ')then
            write(*,*) 'Which parameter file do you want to use?'
            read (*,'(a)') str
      end if
      if(str(1:1).eq.' ') str(1:20) = 'trans.par           '
      inquire(file=str,exist=testfl)
      if(.not.testfl) then
            write(*,*) 'ERROR - the parameter file does not exist,'
            write(*,*) '        check for the file and try again  '
            write(*,*)
            if(str(1:20).eq.'trans.par           ') then
                  write(*,*) '        creating a blank parameter file'
                  call makepar
                  write(*,*)
            end if
            stop
      endif
      open(lin,file=str,status='OLD')
c
c Find Start of Parameters:
c
 1    read(lin,'(a4)',end=97) str(1:4)
      if(str(1:4).ne.'STAR') go to 1
c
c Read Input Parameters:
c

      read(lin,*,err=97) ivtype
      write(*,*) ' variable type = ',ivtype

      read(lin,'(a512)',err=97) distin
      call chknam(distin,512)
      write(*,*) ' reference distribution = ',distin(1:40)

      read(lin,*,err=97) ivr,iwt
      write(*,*) ' columns = ',ivr,iwt

      read(lin,'(a512)',err=97) datafl
      call chknam(datafl,512)
      write(*,*) ' data file = ',datafl(1:40)

      read(lin,*,err=97) ivrd,iwtd
      write(*,*) ' columns = ',ivrd,iwtd

      read(lin,*,err=97) tmin,tmax
      write(*,*) ' trimming limits = ',tmin,tmax

      read(lin,'(a512)',err=97) outfl
      call chknam(outfl,512)
      write(*,*) ' output file = ',outfl(1:40)

      read(lin,*,err=97) nsim
      write(*,*) ' number of sets to transform = ',nsim

      read(lin,*,err=97) nx,ny,nz
      write(*,*) ' size of model to transform = ',nx,ny,nz

      read(lin,*,err=97) wx,wy,wz
      write(*,*) ' window size of model to transform = ',wx,wy,wz

      read(lin,*,err=97) nxyz
      write(*,*) ' size of continuous variable data set = ',nxyz

      read(lin,*,err=97) zmin,zmax
      write(*,*) ' data limits = ',zmin,zmax

      read(lin,*,err=97) ltail,ltpar
      write(*,*) ' lower tail option = ',ltail,ltpar

      read(lin,*,err=97) utail,utpar
      write(*,*) ' upper tail option = ',utail,utpar

      read(lin,*,err=97) ldata
      write(*,*) ' account for local data (1=yes) = ',ldata

      if(ldata.eq.1) then

            read(lin,'(a512)',err=97) localfl
            call chknam(localfl,512)
            write(*,*) ' local file = ',localfl(1:40)

            read(lin,*,err=97) icoll
            write(*,*) ' column for kriging variance = ',icoll

            read(lin,*,err=97) wtfac
            write(*,*) ' scaling factor = ',wtfac
c
c Scale from 0 to 1  --> 0.33 to 3.0
c
            wtfac = 0.33 + wtfac*(3.0-0.33)

            read(lin,*,err=97) ixv(1)
            write(*,*) ' random number seed = ',ixv(1)
            do i=1,10000
                  rn = real(acorni(idum))
            end do

      end if

      close(lin)
c
c Check for error situation:
c
      if(ltail.ne.1.and.ltail.ne.2) then
            write(*,*) 'ERROR invalid lower tail option ',ltail
            write(*,*) '      only allow 1 or 2 - see manual '
            stop
      endif
      if(utail.ne.1.and.utail.ne.2.and.utail.ne.4) then
            write(*,*) 'ERROR invalid upper tail option ',ltail
            write(*,*) '      only allow 1,2 or 4 - see manual '
            stop
      endif
      if(utail.eq.4.and.utpar.lt.1.0) then
            write(*,*) 'ERROR invalid power for hyperbolic tail',utpar
            write(*,*) '      must be greater than 1.0!'
            stop
      endif
c
c Read in the reference distribution:
c
      inquire(file=distin,exist=testfl)
      if(.not.testfl) then
            write(*,*) 'ERROR: No reference distribution file'
            stop
      endif
      open(lin,file=distin,status='UNKNOWN')
c
c Proceed with reading in distribution:
c
      read(lin,*,err=98)
      read(lin,*,err=98)     nvari
      do i=1,nvari
            read(lin,*)
      end do
      maxref = 0
 22   read(lin,*,end =44,err=99)(var(j),j=1,nvari)
      if(var(ivr).lt.tmin.or.var(ivr).ge.tmax)go to 22
      maxref = maxref + 1
      go to 22
 44   continue
c
c Allocate the needed memory:
c
      allocate (rcdf(maxref),stat = test)
            if (test.ne.0) then
                  write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
                  stop
            end if
c
      allocate (rvr(maxref),stat = test)
            if (test.ne.0) then
                  write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
                  stop
            end if
c
      rewind(lin)
      read(lin,'(a)',err=98) str
      read(lin,*,err=98)     nvari
      do i=1,nvari
            read(lin,'()',err=98)
      end do
c
c Read as much data for reference distribution as possible:
c
      ncut = 0
      tcdf = 0
 2    read(lin,*,end=3,err=98) (var(j),j=1,nvari)
      if(var(ivr).lt.tmin.or.var(ivr).ge.tmax) go to 2

      vr = var(ivr)
      wt = 1.0
      if(iwt.ge.1) wt = var(iwt)

      if(ivtype.eq.1) then
            ncut = ncut + 1
            rvr(ncut)  = vr
            rcdf(ncut) = wt
      else
            do i=1,ncut
                  if(int(vr+0.5).eq.int(rvr(i))) then
                        icut = i
                        go to 4
                  end if
            end do
            ncut = ncut + 1
            if(ncut.gt.MAXCAT) then
                  write(*,*) 'ERROR: exceeded available storage for'
                  write(*,*) '       categories, available: ',MAXCAT
                  stop
            endif
            rvr(ncut)  = vr
            rcdf(ncut) = 0.0
            icut       = ncut
 4          continue
            rcdf(icut) = rcdf(icut) + wt
      end if
      tcdf = tcdf + wt

c
c Go back for another data?
c
      go to 2
 3    close(lin)
c
c Sort the Reference Distribution and Check for error situation:
c
      call sortem(1,ncut,rvr,1,rcdf,c,d,e,f,g,h)
      if(ncut.le.1.or.tcdf.le.EPSLON) then
            write(*,*) 'ERROR: too few data or too low weight'
            stop
      endif
      if(ivtype.eq.1.and.utail.eq.4.and.rvr(ncut).le.0.0) then
            write(*,*) 'ERROR can not use hyperbolic tail with '
            write(*,*) '      negative values! - see manual '
            stop
      endif
c
c Turn the (possibly weighted) distribution into a cdf that is useful:
c
      tcdf  = 1.0 / tcdf
      if(ivtype.eq.1) then
            oldcp = 0.0
            cp    = 0.0
            do i=1,ncut
                  cp     = cp + rcdf(i) * tcdf
                  rcdf(i) =(cp + oldcp) * 0.5
                  oldcp  = cp
            end do
      else
            do i=1,ncut
                  rcdf(i) = rcdf(i) * tcdf
            end do
      end if
c
c Write Some of the Statistics to the screen:
c
      if(ivtype.eq.1) then
            call locate(rcdf,ncut,1,ncut,0.5,j)
            gmedian = powint(rcdf(j),rcdf(j+1),rvr(j),rvr(j+1),0.5,1.0)
            write(*,900) ncut,gmedian
 900        format(/' There are ',i8,' data in reference dist,',/,
     +              '   median value        = ',f12.5)
      else
            write(*,*)
            do i=1,ncut
                  catcdf(i) = 0.0
                  write(*,901) i,rvr(i),rcdf(i)
 901              format(' Category number ',i2,' code ',f6.0,
     +                   ' proportion ',f6.4)
                  if(i.gt.1) rcdf(i) = rcdf(i) + rcdf(i-1)
            end do
      end if
c
c Get the output and distribution files ready:
c
      inquire(file=datafl,exist=testfl)
      if(.not.testfl) then
            write(*,*) 'ERROR ',datafl,' does not exist!'
            write(*,*) '  you need some distributions to work with'
            stop
      endif
      open(lin,file=datafl,status='OLD')
      open(lout,file=outfl,status='UNKNOWN')
      read(lin,*,err=96)
      read(lin,*,err=96) nvari
      do i=1,nvari
            read(lin,*)
      end do
      maxdat = 0
 20   read(lin,*,end =40,err=99)(var(j),j=1,nvari)
      if(var(ivr).lt.tmin.or.var(ivr).ge.tmax)go to 20
      maxdat = maxdat + 1
      go to 20
 40   continue
c
c Allocate the needed memory:
c
      allocate (dcdf(maxdat),stat = test)
            if (test.ne.0) then
                  write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
                  stop
            end if
c
      allocate (dvr(maxdat),stat = test)
            if (test.ne.0) then
                  write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
                  stop
            end if
c
      allocate (indx(maxdat),stat = test)
            if (test.ne.0) then
                  write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
                  stop
            end if
c
      allocate (fuzzcat(maxdat),stat = test)
            if (test.ne.0) then
                  write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
                  stop
            end if
c
      rewind(lin)
      read(lin,'(a40)',err=96) str
      read(lin,*,err=96) nvari
      do i=1,nvari
            read(lin,'(a40)',err=96) str
      end do
      write(lout,110)
 110  format('Trans output',/,'1',/,'transformed variable')
c
c MAIN Loop over all the increments to transform:
c
      if(ivtype.eq.0) then
            nxyz = nx*ny*nz
            nxy  = nx*ny
      end if
      do isim=1,nsim
c
c Read in the data values:
c
            tcdf = 0.0
            num  = 0
            do i=1,nxyz
                  read(lin,*,end=5,err=98) (var(j),j=1,nvari)
                  num = num + 1
                  dvr(num)  = var(ivrd)
                  indx(num) = real(num)
                  wt        = 0.0
                  if(dvr(num).ge.tmin.and.dvr(num).lt.tmax) then
                        wt = 1.0
                        if(iwtd.ge.1) wt = var(iwtd)
                        dcdf(num) = wt
                        tcdf      = tcdf + wt
                  endif
c
c Keep track of the proportions if working with a categorical variable:
c
                  if(dvr(num).ge.tmin.and.dvr(num).lt.tmax.and.
     +               ivtype.eq.0) then
                        do j=1,ncut
                           if(int(dvr(num)+0.5).eq.int(rvr(j))) then
                                 icut = j
                                 go to 6
                           end if
                        end do
                        write(*,*) 'Found a code not found in reference'
                        write(*,*) dvr(num)
                        stop
 6                      continue
                        catcdf(icut)  = catcdf(icut) + wt
                  end if
            end do
 5          continue
            if(tcdf.le.EPSLON) then
                  write(*,*) 'ERROR: no data'
                  stop
            endif
            if(ivtype.eq.0.and.num.ne.nxyz) then
                  write(*,*) 'ERROR: you must have ',nxyz
                  write(*,*) '       for a categorical transformation'
                  stop
            endif
c
c For categorical transformation we need a fuzzy category which is the
c local average category - this is for breaking ties in the
c transformation procedure.
c
            if(ivtype.eq.0) then
                  do i=1,nxyz
                        fuzzcat(i) = 0.0
                        nloc       = 0
                        iz = int((i-1)/nxy) + 1
                        iy = int((i-(iz-1)*nxy-1)/nx) + 1
                        ix = i- (iz-1)*nxy - (iy-1)*nx
                        do iix=-wx,wx
                        do iiy=-wy,wy
                        do iiz=-wz,wz
                              jx = ix + iix
                              jy = iy + iiy
                              jz = iz + iiz
                              if(jx.ge.1.and.jx.le.nx.and.
     +                           jy.ge.1.and.jy.le.ny.and.
     +                           jz.ge.1.and.jz.le.nz) then
                                 j = jx + (jy-1)*nx + (jz-1)*nxy
                                 if(dvr(j).ge.tmin.and.
     +                              dvr(j).lt.tmax) then
                                       nloc = nloc + 1
                                       fuzzcat(i) = fuzzcat(i) + dvr(j)
                                 end if
                              end if
                        end do
                        end do
                        end do
                        if(nloc.gt.0) then
                              fuzzcat(i) = fuzzcat(i) / real(nloc)
                        else
                              fuzzcat(i) = 2*tmax
                        end if
                  end do
            end if
c
c Turn the (possibly weighted) data distribution into a useful cdf:
c
            if(ivtype.eq.1) then
                  call sortem(1,num,dvr,2,dcdf,indx,d,e,f,g,h)
            else
                  call sortem(1,num,fuzzcat,3,dcdf,dvr,indx,e,f,g,h)
            end if
            oldcp = 0.0
            cp    = 0.0
            tcdf  = 1.0 / tcdf
            do i=1,num
                  if(dvr(i).ge.tmin.and.dvr(i).lt.tmax) then
                        cp     = cp + dcdf(i) * tcdf
                        dcdf(i) =(cp + oldcp) * 0.5
                        oldcp  = cp
                  endif
            end do
c
c Now, get the right order back:
c
            if(ivtype.eq.1) then
                  call sortem(1,num,indx,2,dcdf,dvr,d,e,f,g,h)
            else
                  call sortem(1,num,indx,3,dcdf,dvr,fuzzcat,e,f,g,h)
            end if
c
c Read in the kriging variance to array "indx" if we have to honor
c local data:
c
            if(ldata.eq.1) then
                  open(lkv,file=localfl,err=95,status='OLD')
                  read(lkv,'()',err=95)
                  read(lkv,*,   err=95) nvarik
                  do i=1,nvarik
                        read(lkv,'()',err=95)
                  end do
                  evmax = -1.0e21
                  do i=1,num
                        read(lkv,*,err=95) (var(j),j=1,nvarik)
                        indx(i) = var(icoll)
                        if(indx(i).ge.0.0.and.indx(i).lt.tmax) then
                              indx(i) = sqrt(max(indx(i),0.0))
                              if(indx(i).gt.evmax) evmax = indx(i)
                        end if
                  end do
                  close(lkv)
            end if
c
c Go through all the data back transforming them to the reference CDF:
c
            do i=1,num
                  if(dvr(i).ge.tmin.and.dvr(i).lt.tmax) then
                        if(ivtype.eq.1) then
                              zval = getz(dcdf(i),ncut,rvr,rcdf,zmin,
     +                                    zmax,ltail,ltpar,utail,utpar)
                        else
                              zval = rvr(1)
                              do j=2,ncut
                                    if(dcdf(i).gt.rcdf(j-1))
     +                              zval = rvr(j)
                              end do
                        end if
c
c Now, do we have to honor local data?
c
                        if(ldata.eq.1) then
                              if(indx(i).ge.0.and.indx(i).lt.tmax) then
                                    wt = (indx(i)/evmax)**wtfac
                                    if(ivtype.eq.1) then
                                          zval = dvr(i)+wt*(zval-dvr(i))
                                    else
                                          rn = real(acorni(idum))
                                          if(rn.gt.wt) zval = dvr(i)
                                    end if
                              end if
                        end if
                  else
                        zval = -999.0
                  endif
                  if(ivtype.eq.1) then
                        call numtext(zval,str(1:12))
                        write(lout,'(a12)') str(1:12)
                  else
                        write(lout,'(f4.0)') zval
                  end if
            end do
c
c END Main loop:
c
      end do
c
c Finished:
c
      write(*,9998) VERSION
 9998 format(/' TRANS Version: ',f5.3, ' Finished'/)
      stop
 95   stop 'ERROR in kriging variance file!'
 96   stop 'ERROR in distribution file!'
 97   stop 'ERROR in parameter file!'
 98   stop 'ERROR in global data file!'
 99   stop 'ERROR in data file!'
      end

      subroutine makepar
c-----------------------------------------------------------------------
c
c                      Write a Parameter File
c                      **********************
c
c
c
c-----------------------------------------------------------------------
      lun = 99
      open(lun,file='trans.par',status='UNKNOWN')
      write(lun,10)
 10   format('                  Parameters for TRANS',/,
     +       '                  ********************',/,/,
     +       'START OF PARAMETERS:')

      write(lun,11)
 11   format('1                         ',
     +       '-1=continuous, 0=categorical')
      write(lun,12)
 12   format('../data/true.dat          ',
     +       '-file with reference distribution')
      write(lun,13)
 13   format('1   0                     ',
     +       '-   columns for variable and weight(0=none)')
      write(lun,14)
 14   format('../data/cluster.dat       ',
     +       '-file with original distributions')
      write(lun,15)
 15   format('3   0                     ',
     +       '-   columns for variable and weight(0=none)')
      write(lun,16)
 16   format('-1.0e21  1.0e21           ',
     +       '-trimming limits')
      write(lun,17)
 17   format('trans.out                 ',
     +       '-file for transformed distributions')
      write(lun,18)
 18   format('1                         ',
     +       '-number of realizations or "sets" to trans')
      write(lun,19)
 19   format('50  50  1                 ',
     +       '-categorical: nx, ny, nz: size of 3-D model')
      write(lun,20)
 20   format(' 2   2  0                 ',
     +       '-   wx, wy, wz: window size for tie-breaking')
      write(lun,21)
 21   format('1000                      ',
     +       '-continuous: number to transform per "set"')
      write(lun,22)
 22   format('0.0   75.0                ',
     +       '-   minimum and maximum values')
      write(lun,23)
 23   format('1      1.0                ',
     +       '-   lower tail: option, parameter')
      write(lun,24)
 24   format('1     75.0                ',
     +       '-   upper tail: option, parameter')
      write(lun,25)
 25   format('0                         ',
     +       '-honor local data? (1=yes, 0=no)')
      write(lun,26)
 26   format('kt3d.out                  ',
     +       '-   file with estimation variance')
      write(lun,27)
 27   format('2                         ',
     +       '-   column number')
      write(lun,28)
 28   format('0.5                       ',
     +       '-   control parameter ( 0.33 < w < 3.0 )')
      write(lun,29)
 29   format('69069                     ',
     +       '-   random number seed (conditioning cat.)')

      close(lun)
      return
      end

ext: for