!**************************************************/
program f90_QQ
!**************************************************/
    implicit none
    integer,parameter::NDmax=1000
    real(8)::data_org(1:NDmax)
    integer::i,nd,io,imethod
    real(8),allocatable::datax(:),datay(:)
    real(8)::alpha
    real(8)::p,yy
    real(8)::rr
    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
    real(8)::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
    character(len=50)::fnameR,fnameW,fnameT
    character(len=100)::dat
    character(len=50)::dummy

    !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
    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

    open(11,file=fnameR,status='old')
    read(11,'()')
    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),datay(1:nd))
    do i=1,nd
        datax(i)=data_org(i)
    end do
    !========================*/
    ! Sorting in small order */
    !========================*/
    call SSORT(nd,datax)
    !------------------------------------------*/
    ! Non-exceeding probability                */
    ! Plotting and positioned formula          */
    !    Weibull : alpha=0 (Mean Rank)         */
    !   (Median Rank : alpha=0.3 )             */
    !    Cunnane : alpha=0.4                   */
    !    Hazen   : alpha=0.5 (Symmetrical CDF) */
    !------------------------------------------*/
    alpha=0.4D0
    do i=1,nd
        datay(i)=(dble(i)-alpha)/(dble(nd+1)-2.0D0*alpha)
    end do
    !=============================*/
    ! Calculation of coefficients */
    !=============================*/
    select case(imethod)
    case(1)
        call GUM_COEF(nd,datax,a_gum,c_gum)
    case(2)
        call GEV_COEF(nd,datax,k_gev,a_gev,c_gev)
    case(3)
        call WGO_COEF(nd,datax,k_wgo,a_wgo,c_wgo)
    case(4)
        call WML_COEF(nd,datax,datay,k_wls,a_wls,c_wls,k_wml,a_wml,c_wml,0)
    case(5)
        call WML_COEF(nd,datax,datay,k_wls,a_wls,c_wls,k_wml,a_wml,c_wml,1)
    case(6)
        call SQR_COEF(nd,datax,a_sqr,b_sqr)
    case(7)
        call LP3_COEF(nd,datax,a_lp3,b_lp3,c_lp3)
    case(8)
        call LNQ_COEF(nd,datax,a_lnq,m_lnq,s_lnq)
    case(9)
        call LNM_COEF(nd,datax,a_lnm,m_lnm,s_lnm)
    case(10)
        call LNT_COEF(nd,datax,datay,a_lnt,m_lnt,s_lnt)
    case(11)
        call EXP_COEF(nd,datax,a_exp,c_exp)
    case(12)
        call GPD_COEF(nd,datax,k_gpd,a_gpd,c_gpd)
    end select
    !===============================*/
    ! Calculation of quantile value */
    !===============================*/
    fnameT='_temp.txt'
    open(12,file=fnameT,status='replace')
    do i=1,nd
        p=datay(i)
        select case(imethod)
        case(1)
            yy=XP_GUM(p,a_gum,c_gum)
        case(2)
            yy=XP_GEV(p,k_gev,a_gev,c_gev)
        case(3)
            yy=XP_WEI(p,k_wgo,a_wgo,c_wgo)
        case(4)
            yy=XP_WEI(p,k_wls,a_wls,c_wls)
        case(5)
            yy=XP_WEI(p,k_wml,a_wml,c_wml)
        case(6)
            yy=XP_SQR(p,a_sqr,b_sqr,datax(nd))
        case(7)
            if(b_lp3<1.0D4)then
                yy=XP_LP3(p,a_lp3,b_lp3,c_lp3)
            else
                yy=XP_LP3_SP(nd,datax,p)
            end if
        case(8)
            yy=XP_LN3(p,a_lnq,m_lnq,s_lnq)
        case(9)
            yy=XP_LN3(p,a_lnm,m_lnm,s_lnm)
        case(10)
            yy=XP_LN3(p,a_lnt,m_lnt,s_lnt)
        case(11)
            yy=XP_EXP(p,a_exp,c_exp)
        case(12)
            yy=XP_GPD(p,k_gpd,a_gpd,c_gpd)
        end select
        write(12,'(f10.4,2(",",f10.2))') p,datax(i),yy
    end do
    close(12)
    !============================*/
    ! Calculation of correlation */
    !============================*/
    open(11,file=fnameT,status='old')
    do i=1,nd
        read(11,*) p,datax(i),datay(i)
    end do
    close(11)
    rr=COC(nd,datax,datay)
 
    open(12,file=fnameW,status='replace')
    write(12,'(a)') 'Input file : '//trim(adjustl(fnameR))
    write(12,'("nd=,",i6)') nd
    select case(imethod)
    case(1)
        write(12,'("GUM,(*_a_b)_r=",4(",",f10.4))'),0.0D0,a_gum,c_gum,rr
    case(2)
        write(12,'("GEV,(k_a_b)_r=",4(",",f10.4))'),k_gev,a_gev,c_gev,rr
    case(3)
        write(12,'("WGO,(k_a_c)_r=",4(",",f10.4))'),k_wgo,a_wgo,c_wgo,rr
    case(4)
        write(12,'("WLS,(k_a_c)_r=",4(",",f10.4))'),k_wls,a_wls,c_wls,rr
    case(5)
        write(12,'("WML,(k_a_c)_r=",4(",",f10.4))'),k_wml,a_wml,c_wml,rr
    case(6)
        write(12,'("SQR,(a_b_*)_r=",4(",",f10.4))'),a_sqr,b_sqr,0.0D0,rr
    case(7)
        write(12,'("LP3,(a_b_c)_r=",4(",",f10.4))'),a_lp3,b_lp3,c_lp3,rr
    case(8)
        write(12,'("LNQ,(a_m_s)_r=",4(",",f10.4))'),a_lnq,m_lnq,s_lnq,rr
    case(9)
        write(12,'("LNM,(a_m_s)_r=",4(",",f10.4))'),a_lnm,m_lnm,s_lnm,rr
    case(10)
        write(12,'("LNT,(a_m_s)_r=",4(",",f10.4))'),a_lnt,m_lnt,s_lnt,rr
    case(11)
        write(12,'("EXP,(*_a_c)_r=",4(",",f10.4))'),0.0D0,a_exp,c_exp,rr
    case(12)
        write(12,'("GPD,(k_a_c)_r=",4(",",f10.4))'),k_gpd,a_gpd,c_gpd,rr
    end select
    write(12,'(a)') 'Prob.,Obs,Calc'
    open(11,file=fnameT,status='old')
    do i=1,nd
        read(11,'(a)') dat
        write(12,'(a)') trim(dat)
    end do
    close(11)
    close(12)

stop

contains
!---------------------------------------------------------------------------*/
    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
!---------------------------------------------------------------------------*/
    real(8) function COC(nd,datat,datap)
        !Coefficient of correlation*/
        integer,intent(in)::nd
        real(8),intent(in)::datat(1:nd),datap(1:nd)
        integer::i
        real(8)::x1,y1,x2,y2,xy,xm,ym,sx,sy,cv,rr

        x1=0.0D0
        y1=0.0D0
        x2=0.0D0
        y2=0.0D0
        xy=0.0D0
        do i=1,nd
            x1=x1+datat(i)
            y1=y1+datap(i)
            x2=x2+datat(i)*datat(i)
            y2=y2+datap(i)*datap(i)
            xy=xy+datat(i)*datap(i)
        end do
        xm=x1/dble(nd)
        ym=y1/dble(nd)
        sx=sqrt(x2/dble(nd)-xm*xm)
        sy=sqrt(y2/dble(nd)-ym*ym)
        cv=xy/dble(nd)-xm*ym

        COC=cv/sx/sy;

    end function COC
!---------------------------------------------------------------------------*/
! 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
        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

!---------------------------------------------------------------------------*/
end program f90_QQ
!---------------------------------------------------------------------------*/