!**************************************************/
program f90_BTS
!**************************************************/
    implicit none
    integer,parameter::NDmax=1000
    integer,parameter::NYE=14
    integer,parameter::NBTS=10000
    real(8)::data_org(1:NDmax)
    real(8)::pyear(1:NYE)
    real(8)::qq(1:NYE)
    real(8)::qbts(1:NBTS,1:NYE)
    real(8)::dataw(1:NBTS)
    integer::nd                         !Number of data
    real(8),allocatable::datax(:)       !Observed data
    real(8),allocatable::datab(:)       !data for bootstrap
    integer,allocatable::idata(:)       !
    real(8)::p,qave,work
    integer::i,j,k,nnn,mmm,kl,ku
    integer::imethod
    character(len=50)::fnameR,fnameW,fnameH
    character(len=200)::strcom
    character(len=1000)::linebuf
    character(len=50)::dummy
    integer::io

    pyear(1)=2.0D0
    pyear(2)=5.0D0
    pyear(3)=10.0D0
    pyear(4)=20.0D0
    pyear(5)=30.0D0
    pyear(6)=50.0D0
    pyear(7)=100.0D0
    pyear(8)=200.0D0
    pyear(9)=300.0D0
    pyear(10)=500.0D0
    pyear(11)=1000.0D0
    pyear(12)=2000.0D0
    pyear(13)=3000.0D0
    pyear(14)=5000.0D0

    !Reading input file and output file from command line
    call getarg(1,dummy)  !Method of analysis
    call getarg(2,fnameR) !Input data file
    call getarg(3,fnameW) !Output data file
    call getarg(4,fnameH) !Histogram data file
    read(dummy,*) imethod !convert character 'dummy' to integer 'imethod'

    ! imethod
    ! 1  Gumbel         Gumbel distribution
    ! 2  GEV            Generalized Extreme Value distribution
    ! 3  Weibull (Goda) Weibull distribution (L-moments by GODA)
    ! 4  Weibull (LSM)  Weibull distribution (least square method)
    ! 5  Weibull (MLM)  Weibull distribution (Maximum Likelihood method)
    ! 6  SQRT-ET        SQRT exponential-type distribution of maximum
    ! 7  LP3            Log Pearson type III distribution
    ! 8  LN3 (Iwai)     Logarithmic Normal distribution with 3 variables (quantile)
    ! 9  LN3 (Moment)   Logarithmic Normal distribution with 3 variables (moment)
    ! 10  LN3 (Trial)    Logarithmic Normal distribution with 3 variables (trial)
    ! 11 Exponential    Exponential distribution
    ! 12 GPD            Generalized Pareto distribution

    !========================*/
    ! Data input */
    !========================*/
    open(11,file=fnameR,status='old')
    read(11,'(a)') strcom
    do i=1,NDmax
        read(11,*,iostat=io) data_org(i)
        if(io/=0)exit
    end do
    close(11)
    nd=i-1
    allocate(datax(1:nd))
    allocate(datab(1:nd))
    allocate(idata(1:nd))
    do i=1,nd
        datax(i)=data_org(i)
    end do
    !=============================*/
    ! Calculation on BTS method */
    !=============================*/
    ! Result: Quantile qbts(1:NBTS,1:NYE)
    do nnn=1,NBTS
        call IRAND(nd,idata)
        do i=1,nd
            datab(i)=datax(idata(i))
        end do
        call CAL_QQ(imethod,nd,datab,pyear,qq,NYE)
        do mmm=1,NYE
            qbts(nnn,mmm)=qq(mmm)
        end do
    end do

    !=============================*/
    ! Output of result */
    !=============================*/
    if(imethod==1) strcom="Gumbel bootstrap"
    if(imethod==2) strcom="GEV bootstrap"
    if(imethod==3) strcom="Weibull (Goda) bootstrap"
    if(imethod==4) strcom="Weibull (LSM) bootstrap"
    if(imethod==5) strcom="Weibull (MLM) bootstrap"
    if(imethod==6) strcom="SQRT-ET bootstrap"
    if(imethod==7) strcom="LP3 bootstrap"
    if(imethod==8) strcom="LN3 (Iwai) bootstrap"
    if(imethod==9) strcom="LN3 (Moment) bootstrap"
    if(imethod==10)strcom="LN3 (Trial) bootstrap"
    if(imethod==11)strcom="Exponential bootstrap"
    if(imethod==12)strcom="GPD bootstrap"

    ! Calculation using original data
    call CAL_QQ(imethod,nd,datax,pyear,qq,NYE)

    open(12,file=fnameW,status='replace')
    write(12,'(a)') trim(strcom)
    write(12,'(a)') 'Year,Prob.,Qorg,Qprd,Q95,Q05'
    do mmm=1,NYE
        do i=1,NBTS
            dataw(i)=qbts(i,mmm)
        end do
        call SSORT(NBTS,dataw)
        do i=1,NBTS
            qbts(i,mmm)=dataw(i)
        end do
        qave=0.0D0
        do i=1,NBTS
            qave=qave+qbts(i,mmm)
        end do
        qave=qave/dble(NBTS)
        kl=int(0.025*dble(NBTS))
        ku=NBTS-kl
        p=1.0D0-1.0D0/pyear(mmm)
        write(linebuf,'(f7.0,",",e10.5,4(",",f7.1))') &
            pyear(mmm),p,qq(mmm),qave,qbts(ku,mmm),qbts(kl,mmm)
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    close(12)

    !=============================*/
    ! Output for Histogram */
    !=============================*/
    mmm=NYE
    dummy=' (5000year)'
    if(imethod==1) strcom="Gumbel bootstrap"
    if(imethod==2) strcom="GEV bootstrap"
    if(imethod==3) strcom="Weibull (Goda) bootstrap"
    if(imethod==4) strcom="Weibull (LSM) bootstrap"
    if(imethod==5) strcom="Weibull (MLM) bootstrap"
    if(imethod==6) strcom="SQRT-ET bootstrap"
    if(imethod==7) strcom="LP3 bootstrap"
    if(imethod==8) strcom="LN3 (Iwai) bootstrap"
    if(imethod==9) strcom="LN3 (Moment) bootstrap"
    if(imethod==10)strcom="LN3 (Trial) bootstrap"
    if(imethod==11)strcom="Exponential bootstrap"
    if(imethod==12)strcom="GPD bootstrap"
    strcom=trim(adjustl(strcom))//trim(adjustl(dummy))
    open(12,file=fnameH,status='replace')
    write(12,'(a)') trim(strcom)
    do i=1,NBTS
        write(linebuf,'(f7.1)') qbts(i,mmm)
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    close(12)

contains
!---------------------------------------------------------------------------*/
    subroutine CAL_QQ(imethod,nd,datat,pyear,qq,NYE)
        !Gumbel
        integer,intent(in)::imethod,nd,NYE
        real(8),intent(in)::pyear(1:NYE)
        real(8),intent(inout)::datat(1:nd)
        real(8),intent(out)::qq(1:NYE)
        integer::i
        real(8)::datay(1:nd),alpha,p
        real(8)::a_gum,c_gum
        real(8)::k_gev,a_gev,c_gev
        real(8)::k_wgo,a_wgo,c_wgo
        real(8)::k_wls,a_wls,c_wls,k_wml,a_wml,c_wml
        real(8)::a_sqr,b_sqr
        real(8)::a_lp3,b_lp3,c_lp3
        real(8)::a_lnq,m_lnq,s_lnq
        real(8)::a_lnm,m_lnm,s_lnm
        real(8)::a_lnt,m_lnt,s_lnt
        real(8)::a_exp,c_exp
        real(8)::k_gpd,a_gpd,c_gpd

        call SSORT(nd,datat)
        alpha=0.4D0
        do i=1,nd
            datay(i)=(dble(i)-alpha)/(dble(nd+1)-2.0D0*alpha)
        end do

        select case(imethod)
        case(1)
            call GUM_COEF(nd,datat,a_gum,c_gum)
            do i=1,NYE
                p=1.0D0-1.0D0/pyear(i)
                qq(i)=XP_GUM(p,a_gum,c_gum)
            end do
        case(2)
            call GEV_COEF(nd,datat,k_gev,a_gev,c_gev)
            do i=1,NYE
                p=1.0D0-1.0D0/pyear(i)
                qq(i)=XP_GEV(p,k_gev,a_gev,c_gev)
            end do
        case(3)
            call WGO_COEF(nd,datat,k_wgo,a_wgo,c_wgo)
            do i=1,NYE
                p=1.0D0-1.0D0/pyear(i)
                qq(i)=XP_WEI(p,k_wgo,a_wgo,c_wgo)
            end do
        case(4)
            call WML_COEF(nd,datat,datay,k_wls,a_wls,c_wls,k_wml,a_wml,c_wml,0)
            do i=1,NYE
                p=1.0D0-1.0D0/pyear(i)
                qq(i)=XP_WEI(p,k_wls,a_wls,c_wls)
            end do
        case(5)
            call WML_COEF(nd,datat,datay,k_wls,a_wls,c_wls,k_wml,a_wml,c_wml,1)
            do i=1,NYE
                p=1.0D0-1.0D0/pyear(i)
                qq(i)=XP_WEI(p,k_wml,a_wml,c_wml)
            end do
        case(6)
            call SQR_COEF(nd,datat,a_sqr,b_sqr)
            do i=1,NYE
                p=1.0D0-1.0D0/pyear(i)
                qq(i)=XP_SQR(p,a_sqr,b_sqr,datat(nd))
            end do
        case(7)
            call LP3_COEF(nd,datat,a_lp3,b_lp3,c_lp3)
            do i=1,NYE
                p=1.0D0-1.0D0/pyear(i)
                if(b_lp3<1.0D4)then
                    qq(i)=XP_LP3(p,a_lp3,b_lp3,c_lp3)
                else
                    qq(i)=XP_LP3_SP(nd,datax,p)
                end if
            end do
        case(8)
           call LNQ_COEF(nd,datat,a_lnq,m_lnq,s_lnq)
            do i=1,NYE
                p=1.0D0-1.0D0/pyear(i)
                qq(i)=XP_LN3(p,a_lnq,m_lnq,s_lnq)
            end do
        case(9)
           call LNM_COEF(nd,datat,a_lnm,m_lnm,s_lnm)
            do i=1,NYE
                p=1.0D0-1.0D0/pyear(i)
                qq(i)=XP_LN3(p,a_lnm,m_lnm,s_lnm)
            end do
        case(10)
           call LNT_COEF(nd,datat,datay,a_lnt,m_lnt,s_lnt)
            do i=1,NYE
                p=1.0D0-1.0D0/pyear(i)
                qq(i)=XP_LN3(p,a_lnt,m_lnt,s_lnt)
            end do
        case(11)
           call EXP_COEF(nd,datat,a_exp,c_exp)
            do i=1,NYE
                p=1.0D0-1.0D0/pyear(i)
                qq(i)=XP_EXP(p,a_exp,c_exp)
            end do
        case(12)
           call GPD_COEF(nd,datat,k_gpd,a_gpd,c_gpd)
            do i=1,NYE
                p=1.0D0-1.0D0/pyear(i)
                qq(i)=XP_GPD(p,k_gpd,a_gpd,c_gpd)
            end do
        end select

    end subroutine CAL_QQ
!---------------------------------------------------------------------------*/
    subroutine IRAND(nd,idata)
        integer,intent(in)::nd
        integer,intent(out)::idata(1:nd)
        integer::i
        real(8)::a,b,x,r

        a=1.0
        b=dble(nd)
        do i=1,nd
            call random_number(x)
            r=a+(b-a)*x
            idata(i)=int(r+0.5)
        end do
    end subroutine IRAND
!---------------------------------------------------------------------------*/
    subroutine SSORT(nd,datat)
        integer,intent(in)::nd
        real(8),intent(inout)::datat(:)
        integer::i,kh,kv
        real(8)::temp

        kv=1
        do
            kv=3*kv+1
            if(nd<kv)exit
        end do
        do
            kv=kv/3
            do i=kv+1,nd
                temp=datat(i)
                kh=i
                do while(datat(kh-kv)>temp)
                    datat(kh)=datat(kh-kv)
                    kh=kh-kv
                    if(kh<kv)exit
                end do
                datat(kh)=temp
            end do
            if(kv==1)exit
        end do
    end subroutine SSORT
!---------------------------------------------------------------------------*/
! Gumbel
!    subroutine GUM_COEF(nd,datat,a_gum,c_gum)
!    real(8) function XP_GUM(p,a_gum,c_gum)
!---------------------------------------------------------------------------*/
    subroutine GUM_COEF(nd,datat,a_gum,c_gum)
        ! Gumbel distribution  */
        integer,intent(in)::nd
        real(8),intent(in)::datat(1:nd)
        real(8),intent(out)::a_gum,c_gum
        integer::i
        real(8)::b0,b1
        real(8)::lam1,lam2

        ! datat() is sorted data */
        b0=0.0D0
        b1=0.0D0
        do i=1,nd
            b0=b0+datat(i)
            b1=b1+dble(i-1)*datat(i)
        end do
        b0=b0/dble(nd) !Average*/
        b1=b1/dble(nd*(nd-1))
        lam1=b0
        lam2=2.0D0*b1-b0

        !Gumbel F(x)=exp{-exp(-(x-c)/a)} T=1/(1-F(x)) */
        a_gum=lam2/log(2.0D0)
        c_gum=lam1-0.5772D0*a_gum

    end subroutine GUM_COEF


    real(8) function XP_GUM(p,a_gum,c_gum)
        real(8),intent(in)::p,a_gum,c_gum
        XP_GUM=c_gum-a_gum*log(-log(p))
    end function XP_GUM
!---------------------------------------------------------------------------*/
! GEV
!    subroutine GEV_COEF(nd,datat,k_gev,a_gev,c_gev)
!    real(8) function XP_GEV(p,k_gev,a_gev,c_gev)
!---------------------------------------------------------------------------*/
    subroutine GEV_COEF(nd,datat,k_gev,a_gev,c_gev)
        integer,intent(in)::nd
        real(8),intent(in)::datat(1:nd)
        real(8),intent(out)::k_gev,a_gev,c_gev
        ! Generalized Extreme Value distribution  */
        integer::i
        real(8)::b0,b1,b2,d
        real(8)::lam1,lam2,lam3
        real(8)::x,gammax
        real(8)::pi=3.14159265358979323846D0

        ! datat() is sorted data */
        b0=0.0D0
        b1=0.0D0
        b2=0.0D0
        do i=1,nd
            b0=b0+datat(i)
            b1=b1+dble(i-1)*datat(i)
            b2=b2+dble((i-1)*(i-2))*datat(i)
        end do
        b0=b0/dble(nd) !Average*/
        b1=b1/dble(nd*(nd-1))
        b2=b2/dble(nd*(nd-1)*(nd-2))
        lam1=b0
        lam2=2.0D0*b1-b0
        lam3=6.0D0*b2-6.0D0*b1+b0

        !GEV F(x)=exp{-(1-k*(x-c)/a)^(1/k)} T=1/(1-F(x)) */
        d=2.0D0*lam2/(lam3+3.0D0*lam2)-log(2.0D0)/log(3.0D0)
        k_gev=7.8590D0*d+2.9554D0*d*d
        x=1.0D0+k_gev
        gammax=sqrt(2.0D0*pi/x)*(x/exp(1.0D0)*sqrt(x*sinh(1.0D0/x)+1.0D0/810.0D0/x**6.0D0))**x
        a_gev=k_gev*lam2/(1.0D0-2.0D0**(-k_gev))/gammax
        c_gev=lam1-a_gev/k_gev*(1.0D0-gammax)

    end subroutine GEV_COEF


    real(8) function XP_GEV(p,k_gev,a_gev,c_gev)
        real(8),intent(in)::p,k_gev,a_gev,c_gev
        XP_GEV=c_gev+(a_gev/k_gev)*(1.0D0-(-log(p))**k_gev)
    end function XP_GEV
!---------------------------------------------------------------------------*/
! Weibull (GODA, MOM, MLM)
!    subroutine WGO_COEF(nd,datat,k_wgo,a_wgo,c_wgo)
!    subroutine WML_COEF(nd,datat,datap,k_wls,a_wls,c_wls,k_wml,a_wml,c_wml,newton)
!    subroutine LSM(nd,vecx,vecy,aa,bb,rr)
!    real(8) function XP_WEI(p,k_wei,a_wei,c_wei)
!---------------------------------------------------------------------------*/
    subroutine WGO_COEF(nd,datat,k_wgo,a_wgo,c_wgo)
        integer,intent(in)::nd
        real(8),intent(in)::datat(1:nd)
        real(8),intent(out)::k_wgo,a_wgo,c_wgo
        integer::i
        real(8)::b0,b1,b2
        real(8)::lam1,lam2,lam3,t
        real(8)::x,gammax
        real(8)::pi=3.14159265358979323846D0

        ! datat() is sorted data */
        b0=0.0D0
        b1=0.0D0
        b2=0.0D0
        do i=1,nd
            b0=b0+datat(i)
            b1=b1+dble(i-1)*datat(i)
            b2=b2+dble((i-1)*(i-2))*datat(i)
        end do
        b0=b0/dble(nd) !Average*/
        b1=b1/dble(nd*(nd-1))
        b2=b2/dble(nd*(nd-1)*(nd-2))
        lam1=b0
        lam2=2.0D0*b1-b0
        lam3=6.0D0*b2-6.0D0*b1+b0

        t=lam3/lam2
        k_wgo=285.3D0*t**6.0D0-658.6D0*t**5.0D0+622.8D0*t**4.0D0-317.2D0*t**3.0D0+98.52D0*t**2.0D0-21.256D0*t+3.5160D0
        x=1.0D0+1.0D0/k_wgo
        gammax=sqrt(2.0D0*pi/x)*(x/exp(1.0D0)*sqrt(x*sinh(1.0D0/x)+1.0D0/810.0D0/x**6.0D0))**x
        a_wgo=lam2/(1.0D0-2.0**(-1.0D0/k_wgo))/gammax
        c_wgo=lam1-a_wgo*gammax

    end subroutine WGO_COEF


    subroutine WML_COEF(nd,datat,datap,k_wls,a_wls,c_wls,k_wml,a_wml,c_wml,newton)
        integer,intent(in)::nd,newton
        real(8),intent(in)::datat(1:nd),datap(1:nd)
        real(8),intent(out)::k_wls,a_wls,c_wls
        real(8),intent(out)::k_wml,a_wml,c_wml
        integer::i,j,itmax=1000
        real(8)::c,dc
        real(8)::vecx(1:nd),vecy(1:nd)
        real(8)::aa,bb,rr1,rr2
        real(8)::kk,ff,df,tt,t0,t1,t2,t3

        ! Culculation of coefficient (c) */
        dc=0.1D0
        c=datat(1)
        rr1=0.0
        do j=1,itmax
            c=c-dc
            do i=1,nd
                tt=datat(i)-c
                if(tt<=0.0D0)tt=0.1D0*dc
                vecx(i)=log(tt)
                vecy(i)=log(-log(1.0D0-datap(i)))
            end do
            call LSM(nd,vecx,vecy,aa,bb,rr2)
            if(abs(rr2)<abs(rr1))exit
            rr1=rr2
        end do
        k_wls=aa
        a_wls=exp(-bb/aa)
        c_wls=c

        k_wml=0.0D0
        a_wml=0.0D0
        c_wml=0.0D0

        if(newton/=0)then
            ! Newton-Raphson method */
            dc=0.1D0
            kk=k_wls
            do while(1.0D-5<abs(ff/df))
                t0=0.0
                t1=0.0
                t2=0.0
                t3=0.0
                do i=1,nd
                    tt=datax(i)-c_wls
                    if(tt<=0.0D0)tt=dc
                    t0=t0+log(tt)
                    t1=t1+tt**kk
                    t2=t2+log(tt)*tt**kk
                    t3=t3+log(tt)*log(tt)*tt**kk
                end do
                ff=1.0/kk+t0/dble(nd)-t2/t1
                df=-1.0/kk/kk-(t3*t1-t2*t2)/t1/t1
                kk=kk-ff/df
            end do
            k_wml=kk
            a_wml=(t1/dble(nd))**(1.0D0/kk)
            c_wml=c_wls
        end if

    end subroutine WML_COEF


    subroutine LSM(nd,vecx,vecy,aa,bb,rr)
    !--------------------------*/
    ! Regression (y=aa*x+bb)   */
    !--------------------------*/
        integer,intent(in)::nd
        real(8),intent(in)::vecx(1:nd),vecy(1:nd)
        real(8),intent(out)::aa,bb,rr
        integer::i
        real(8)::x1,y1,x2,y2,xy,xm,ym,c1,c2,c3
    
        x1=0.0D0
        y1=0.0D0
        x2=0.0D0
        xy=0.0D0
        do i=1,nd
            x1=x1+vecx(i)
            y1=y1+vecy(i)
            x2=x2+vecx(i)*vecx(i)
            xy=xy+vecx(i)*vecy(i)
        end do
        xm=x1/dble(nd)
        ym=y1/dble(nd)
        aa=(dble(nd)*xy-x1*y1)/(dble(nd)*x2-x1*x1)
        bb=(x2*y1-x1*xy)/(dble(nd)*x2-x1*x1)
        c1=0.0D0
        c2=0.0D0
        c3=0.0D0
        do i=1,nd
            c1=c1+(vecx(i)-xm)*(vecy(i)-ym)
            c2=c2+(vecx(i)-xm)*(vecx(i)-xm)
            c3=c3+(vecy(i)-ym)*(vecy(i)-ym)
        end do
        rr=c1/sqrt(c2*c3)
    end subroutine LSM


    real(8) function XP_WEI(p,k_wei,a_wei,c_wei)
        real(8),intent(in)::p,k_wei,a_wei,c_wei
        XP_wei=c_wei+a_wei*(-log(1.0D0-p))**(1.0D0/k_wei)
    end function XP_wei
!---------------------------------------------------------------------------*/
! SQRT exponential-type distribution
!    subroutine SQR_COEF(nd,datat,a_sqr,b_sqr)
!    real(8) function FSQR(nd,datat,bb,a1,a2)
!    real(8) function XP_SQR(p,a_sqr,b_sqr,xmax)
!---------------------------------------------------------------------------*/
    subroutine SQR_COEF(nd,datat,a_sqr,b_sqr)
        !Required Function: FSQR
        integer,intent(in)::nd
        real(8),intent(in)::datat(1:nd)
        real(8),intent(out)::a_sqr,b_sqr
        integer::i
        real(8)::a1,a2,b1,b2,bb,f1,f2,ff

        ! Condition for a1>0 */
        bb=0.0D0
        do i=1,nd
            bb=bb+sqrt(datat(i))
        end do
        bb=(2.0D0*dble(nd)/bb)*(2.0D0*dble(nd)/bb)+0.001D0

        ! Bisection method */
        b1=bb
        b2=b1+0.5D0
        bb=0.5D0*(b1+b2)
        f1=FSQR(nd,datat,b1,a1,a2)
        f2=FSQR(nd,datat,b2,a1,a2)
        ff=FSQR(nd,datat,bb,a1,a2)
        do while(0.001<abs(a1-a2))
            if(f1*ff<0.0D0)b2=bb
            if(ff*f2<0.0D0)b1=bb
            if(ff==0.0D0)exit
            if(0.0D0<f1*ff.and.0.0D0<ff*f2)then
                b1=b2
                b2=b1+0.5D0
            end if
            bb=0.5D0*(b1+b2)
            f1=FSQR(nd,datat,b1,a1,a2)
            f2=FSQR(nd,datat,b2,a1,a2)
            ff=FSQR(nd,datat,bb,a1,a2)
        end do
        a_sqr=0.5D0*(a1+a2)
        b_sqr=bb

    end subroutine SQR_COEF


    real(8) function FSQR(nd,datat,bb,a1,a2)
        integer,intent(in)::nd
        real(8),intent(in)::datat(1:nd),bb
        real(8),intent(out)::a1,a2
        integer::i
        real(8)::bxr,bxe,bx1e,bxre

        bxr=0.0D0
        bxe=0.0D0
        bx1e=0.0D0
        bxre=0.0D0
        do i=1,nd
            bxr=bxr+sqrt(bb*datat(i))
            bxe=bxe+exp(-sqrt(bb*datat(i)))
            bx1e=bx1e+bb*datat(i)*exp(-sqrt(bb*datat(i)))
            bxre=bxre+sqrt(bb*datat(i))*exp(-sqrt(bb*datat(i)))
        end do
        a1=(bxr-2.0D0*dble(nd))/bx1e
        a2=dble(nd)/(bxe+bxre)
        FSQR=a1-a2

    end function FSQR


    real(8) function XP_SQR(p,a_sqr,b_sqr,xmax)
        real(8),intent(in)::p,a_sqr,b_sqr,xmax
        real(8)::x1,x2,ff,df,tp

        tp=sqrt(b_sqr*xmax) ! Initial value of tp */
        ff=log(1.0D0+tp)-tp-log(-1.0D0/a_sqr*log(p))
        df=1.0D0/(1.0D0+tp)-1.0D0
        tp=tp-ff/df
        x2=tp*tp/b_sqr ! Iniial value of x2 */
        do while(0.001<abs(x2-x1))
            x1=x2
            ff=log(1.0D0+tp)-tp-log(-1.0D0/a_sqr*log(p))
            df=1.0D0/(1.0D0+tp)-1.0D0
            tp=tp-ff/df
            x2=tp*tp/b_sqr
            if(abs(x2-x1)<0.001)exit
        end do
        XP_SQR=x2

    end function XP_SQR
!---------------------------------------------------------------------------*/
! LP3
!    subroutine LP3_COEF(nd,datat,a_lp3,b_lp3,c_lp3)
!    real(8) function XP_LP3_SP(nd,datat,p)
!    real(8) function GDISTPP(eta,p)
!    real(8) function LINCGF(y1,z1)
!    real(8) function LGAMMA(x)
!    real(8) function TODA(rx)
!---------------------------------------------------------------------------*/
    subroutine LP3_COEF(nd,datat,a_lp3,b_lp3,c_lp3)
        ! Log Pearson III distribution (Gamma distribution) */
        integer,intent(in)::nd
        real(8),intent(in)::datat(1:nd)
        real(8),intent(out)::a_lp3,b_lp3,c_lp3
        integer::i
        real(8)::wmy,wsy,wry

        wmy=0.0D0
        do i=1,nd
            wmy=wmy+log(datat(i))
        end do
        wmy=wmy/dble(nd) !Average
        wsy=0.0D0
        do i=1,nd
            wsy=wsy+(log(datat(i))-wmy)**2.0D0
        end do
        wsy=sqrt(wsy/dble(nd-1)) !Standard deviation
        wry=0.0D0
        do i=1,nd
            wry=wry+(log(datat(i))-wmy)**3.0D0
        end do
        wry=wry*dble(nd)/dble(nd-1)/dble(nd-2)/wsy**3.0D0 !Skewness

        b_lp3=4.0D0/(wry*wry)
        a_lp3=wsy/sqrt(b_lp3)
        if(wry<0.0D0)a_lp3=-a_lp3
        c_lp3=wmy-a_lp3*b_lp3

    end subroutine LP3_COEF


    real(8) function XP_LP3(p,a_lp3,b_lp3,c_lp3)
        ! Log Pearson III distribution (Gamma distribution) */
        ! Required Function: GDISTPP
        real(8),intent(in)::p,a_lp3,b_lp3,c_lp3
        real(8)::q
        q=p
        if(a_lp3<0.0D0)q=1.0D0-p
        XP_LP3=exp(c_lp3+a_lp3*GDISTPP(b_lp3,q))
    end function XP_LP3


    real(8) function XP_LP3_SP(nd,datat,p)
        ! Log Pearson III distribution (Gamma distribution) */
        ! Required Function: TODA
        integer,intent(in)::nd
        real(8),intent(in)::datat(1:nd),p
        integer::i
        real(8)::wmy,wsy,wry,ws,wc,w1,w2
        real(8)::q,kp,ry

        wmy=0.0D0
        do i=1,nd
            wmy=wmy+log(datat(i))
        end do
        wmy=wmy/dble(nd) !Average
        wsy=0.0D0
        do i=1,nd
            wsy=wsy+(log(datat(i))-wmy)**2.0D0
        end do
        wsy=sqrt(wsy/dble(nd-1)) !Standard deviation
        wry=0.0D0
        do i=1,nd
            wry=wry+(log(datat(i))-wmy)**3.0D0
        end do
        wry=wry*dble(nd)/dble(nd-1)/dble(nd-2)/wsy**3.0D0 !Skewness
        q=p
        ry=wry
        kp=2.0D0/ry*(1.0D0+ry*TODA(q)/6.0D0-ry*ry/36.0D0)**3.0D0-2.0D0/ry
        XP_LP3_SP=exp(wmy+wsy*kp)

    end function XP_LP3_SP


    real(8) function GDISTPP(eta,p)
        ! Gamma function %-point */
        ! Required Function: LGAMMA, LINCGF
        ! eta=b_lp3,p=probability */
        real(8),intent(in)::eta,p
        integer::k
        real(8)::etad,bd,et,g1,x,dk,xx,a,b,expa,q,qd,er

        etad=eta
        bd=p
        et=log(etad)
        g1=log(bd)+LGAMMA(etad)
        x=(et+g1)*1.0D0/etad
        k=0
        do while(LINCGF(xx,et)<g1)
            k=k+1
            dk=dble(k)
            xx=x+log(dk)
            if(g1<=LINCGF(xx,et))exit
        end do
        a=0.0D0
        if(2<=k)a=x+log(dk-1.0D0)
        b=x+log(dk)
        expa=exp(a)
        if(k<=1)expa=0.0D0
        do k=1,100000
            q=0.5D0*(expa+exp(b))
            qd=log(q)
            er=LINCGF(qd,et)-g1
            if(er<0.0D0)then
                a=qd
                expa=exp(a)
            end if
            if(abs(er)<1.0D-10)exit
            if(0.0D0<er)b=qd
        end do
        GDISTPP=q

    end function GDISTPP


    real(8) function LGAMMA(x)
        real(8),intent(in)::x
        real(8)::a,b,p,q,z
        real(8),parameter::pi=3.14159265358979323846D0

        if(x<=20.0D0)then
            b=x
            a=1.0D0/(x*(x+1.0D0))
            b=b-1.0D0
            if(0.0D0<=b)then
                do while(0.0D0<=b)
                    a=a*(b+2.0D0)
                    b=b-1.0D0
                end do
            end if
            b=x-int(x)
            p=(((-45.30104765151624D0*b-261.0284950616458D0)*b-1817.799191156434D0)*b-4954.768972677267D0)*b-13572.65011103686D0
            q=((((((1.0D0*b-16.16160629773462D0)*b+93.7769173224393D0)*b-124.7947256227383D0)*b-931.1785045593485D0)*b &
               +3440.699241346454D0)*b+783.5348800559552D0)*b-13572.65011103691D0
            LGAMMA=log(a*p/q)
        else
            z=log(2.0D0*pi)-log(x)+x*(2.0D0*log(x)+log(x*sinh(1.0D0/x)+1.0D0/810.0D0/x**6.0D0)-2.0D0)
            LGAMMA=0.5D0*z
        end if
    end function LGAMMA


    real(8) function LINCGF(y1,z1)
        real(8),intent(in)::y1,z1
        integer::j
        real(8)::z,y,f,d1,d2

        z=exp(z1)
        y=exp(y1)
        f=1.0D0/z
        d1=1.0D0/z
        j=0
        do j=1,100000
            d2=d1*y/(z+dble(j))
            f=f+d2
            if(y/(z+dble(j)+1.0D0-y)*d2<f*1.0D-10.and.y/(z+dble(j)+1.0D0)<1.0D0)exit
            d1=d2
        end do
        LINCGF=z*y1-y+log(f)

    end function LINCGF


    real(8) function TODA(rx)
        ! %-point of Non-exceeding probability */
        real(8),intent(in)::rx
        integer::i
        real(8)::ry,ay,uay
        real(8)::b(0:10),sum

        if(rx==0.5D0)then
            uay=0.0D0
        else
            if(rx<0.5D0)then
                ry=rx
            else
                ry=1.0D0-rx
            end if
            ay=-log(4.0D0*ry*(1.0D0-ry))
            b(0) = 1.570796288
            b(1) = 0.03706987906
            b(2) =-0.0008364353589
            b(3) =-0.0002250947176
            b(4) = 0.000006841218299
            b(5) = 0.000005824238515
            b(6) =-0.00000104527497
            b(7) = 0.00000008360937017
            b(8) =-0.000000003231081277
            b(9) = 0.00000000003657763036
            b(10)= 0.0000000000006936233982
            sum=b(0)
            do i=1,10
                sum=sum+b(i)*ay**dble(i)
            end do
            uay=sqrt(ay*sum)
            if(rx<0.5D0)uay=-uay
        end if
        TODA=uay

    end function TODA
!---------------------------------------------------------------------------*/
! LN3 (LNQ[Iwai], LNM[Moment], LNT[Trial])
!    subroutine LNQ_COEF(nd,datat,a_lnq,m_lnq,s_lnq)
!    subroutine LNM_COEF(nd,datat,a_lnm,m_lnm,s_lnm)
!    subroutine LNT_COEF(nd,datat,datap,a_lnt,m_lnt,s_lnt)
!    real(8) function XP_LN3(p,a_ln3,m_ln3,s_ln3)
!---------------------------------------------------------------------------*/
    subroutine LNQ_COEF(nd,datat,a_lnq,m_lnq,s_lnq)
        ! Log-normal distribution 3-parameter quantile-method */
        integer,intent(in)::nd
        real(8),intent(in)::datat(1:nd)
        real(8),intent(out)::a_lnq,m_lnq,s_lnq
        integer::i
        real(8)::xm,w

        xm=0.5D0*(datat(int(nd/2)+1)+datat(nd-int(nd/2)))
        a_lnq=(datat(1)*datat(nd)-xm*xm)/(datat(1)+datat(nd)-2.0D0*xm)
        w=0.0D0
        do i=1,nd
            w=w+log(datat(i)-a_lnq)
        end do
        m_lnq=w/dble(nd)
        w=0.0D0
        do i=1,nd
            w=w+(log(datat(i)-a_lnq)-m_lnq)*(log(datat(i)-a_lnq)-m_lnq)
        end do
        s_lnq=sqrt(w/dble(nd))

    end subroutine LNQ_COEF


    subroutine LNM_COEF(nd,datat,a_lnm,m_lnm,s_lnm)
        ! Log-normal distribution 3-parameter moment-method */
        integer,intent(in)::nd
        real(8),intent(in)::datat(1:nd)
        real(8),intent(out)::a_lnm,m_lnm,s_lnm
        integer::i
        real(8)::wmx,wsx,wrx
        real(8)::b,xx

        wmx=0.0D0
        do i=1,nd
            wmx=wmx+(datat(i))
        end do
        wmx=wmx/dble(nd) !Average
        wsx=0.0D0
        do i=1,nd
            wsx=wsx+((datat(i))-wmx)**2.0D0
        end do
        wsx=sqrt(wsx/dble(nd-1)) !Standard deviation
        wrx=0.0D0
        do i=1,nd
            wrx=wrx+(datat(i)-wmx)**3.0D0
        end do
        wrx=wrx*dble(nd)/dble(nd-1)/dble(nd-2)/wsx**3.0D0 !Skewness
        ! Solution of equation
        b=1.0D0+wrx*wrx/2.0D0
        xx=(b+sqrt(b*b-1.0D0))**(1.0D0/3.0D0)+(b-sqrt(b*b-1.0D0))**(1.0D0/3.0D0)-1.0D0

        s_lnm=sqrt(log(xx))
        m_lnm=log(wsx/sqrt(xx*(xx-1)))
        a_lnm=wmx-exp(m_lnm)*exp(0.5*s_lnm*s_lnm)
    end subroutine LNM_COEF


    subroutine LNT_COEF(nd,datat,datap,a_lnt,m_lnt,s_lnt)
        ! Log-normal distribution 3-parameter quantile-method */
        ! Trial method  */
        integer,intent(in)::nd
        real(8),intent(in)::datat(1:nd),datap(1:nd)
        real(8),intent(out)::a_lnt,m_lnt,s_lnt
        integer::i,j,itmax=1000
        real(8)::c,dc,tt,aa,bb,rr1,rr2
        real(8)::vecx(1:nd),vecy(1:nd)

        ! Culculation of coefficient (c) */
        dc=1.0D0
        c=datat(1)
        rr1=0.0
        do j=1,itmax
            c=c-dc
            do i=1,nd
                tt=datat(i)-c
                if(tt<=0.0D0)tt=0.1D0*dc
                vecx(i)=log(tt)
                vecy(i)=TODA(datap(i))
            end do
            call LSM(nd,vecx,vecy,aa,bb,rr2)
            if(abs(rr2)<abs(rr1))exit
            rr1=rr2
        end do
        a_lnt=c
        m_lnt=-bb/aa
        s_lnt=1.0D0/aa

    end subroutine LNT_COEF


    real(8) function XP_LN3(p,a_ln3,m_ln3,s_ln3)
        real(8),intent(in)::p,a_ln3,m_ln3,s_ln3
        XP_LN3=a_ln3+exp(m_ln3+s_ln3*TODA(p))
    end function XP_LN3
!---------------------------------------------------------------------------*/
! Exponential
!    subroutine EXP_COEF(nd,datat,a_exp,c_exp)
!    real(8) function XP_EXP(p,a_exp,c_exp)
!---------------------------------------------------------------------------*/
    subroutine EXP_COEF(nd,datat,a_exp,c_exp)
        ! Gumbel distribution  */
        integer,intent(in)::nd
        real(8),intent(in)::datat(1:nd)
        real(8),intent(out)::a_exp,c_exp
        integer::i
        real(8)::b0,b1
        real(8)::lam1,lam2

        ! datat() is sorted data */
        b0=0.0D0
        b1=0.0D0
        do i=1,nd
            b0=b0+datat(i)
            b1=b1+dble(i-1)*datat(i)
        end do
        b0=b0/dble(nd) !Average*/
        b1=b1/dble(nd*(nd-1))
        lam1=b0
        lam2=2.0D0*b1-b0

        !Exponential F(x)=1-exp{-(x-c)/a)} */
        a_exp=2.0D0*lam2
        c_exp=lam1-a_exp

    end subroutine EXP_COEF


    real(8) function XP_EXP(p,a_exp,c_exp)
        real(8),intent(in)::p,a_exp,c_exp
        XP_EXP=c_exp-a_exp*log(1.0D0-p)
    end function XP_EXP
!---------------------------------------------------------------------------*/
! GPD
!    subroutine GPD_COEF(nd,datat,k_gpd,a_gpd,c_gpd)
!    real(8) function XP_GPD(p,k_gpd,a_gpd,c_gpd)
!---------------------------------------------------------------------------*/
    subroutine GPD_COEF(nd,datat,k_gpd,a_gpd,c_gpd)
        integer,intent(in)::nd
        real(8),intent(in)::datat(1:nd)
        real(8),intent(out)::k_gpd,a_gpd,c_gpd
        ! Generalized Pareto distribution  */
        integer::i
        real(8)::b0,b1,b2
        real(8)::lam1,lam2,lam3

        ! datat() is sorted data */
        b0=0.0D0
        b1=0.0D0
        b2=0.0D0
        do i=1,nd
            b0=b0+datat(i)
            b1=b1+dble(i-1)*datat(i)
            b2=b2+dble((i-1)*(i-2))*datat(i)
        end do
        b0=b0/dble(nd) !Average*/
        b1=b1/dble(nd*(nd-1))
        b2=b2/dble(nd*(nd-1)*(nd-2))
        lam1=b0
        lam2=2.0D0*b1-b0
        lam3=6.0D0*b2-6.0D0*b1+b0

        !GPD F(x)=1-{1-k(x-c)/a}^(1/k) */
        k_gpd=(lam2-3.0D0*lam3)/(lam2+lam3)
        a_gpd=(1.0D0+k_gpd)*(2.0D0+k_gpd)*lam2
        c_gpd=lam1-(2.0D0+k_gpd)*lam2

    end subroutine GPD_COEF


    real(8) function XP_GPD(p,k_gpd,a_gpd,c_gpd)
        real(8),intent(in)::p,k_gpd,a_gpd,c_gpd
        XP_GPD=c_gpd+a_gpd/k_gpd*(1.0D0-(1.0D0-p)**k_gpd)
    end function XP_GPD

!---------------------------------------------------------------------------*/
    subroutine del_spaces(s)
        character (*), intent (inout) :: s
        character (len=len(s)) tmp
        integer i, j
        j = 1
        do i = 1, len(s)
            if (s(i:i)==' ') cycle
            tmp(j:j) = s(i:i)
            j = j + 1
        end do
        s = tmp(1:j-1)
    end subroutine del_spaces
!-----------------------------------------------------------------------*/
end program f90_BTS