program f90_ams_jpn
    !*********************************************
    ! Buckling pressure of steel pipe under external ground water pressure
    ! (Technical Standards for gates and penstocks)
    !*********************************************
    implicit none
    integer::nnn,iii
    real(8)::Es    !Elastic modulus of steel pipe
    real(8)::pos   !Poisson's ratio
    real(8)::As    !Thermal expansion coefficient
    real(8)::DD    !Internal diameter of steel pipe
    real(8)::t0    !Design thickness of steel pipe
    real(8)::eps   !Margin thickness of steel pipe
    real(8)::Sigy  !Yield stress of steel pipe
    real(8)::Siga  !Allowable stress of steel pipe
    real(8)::Jc    !Joint efficiency
    real(8)::Temp  !Temperature change
    real(8)::Bg    !Plastic deformation modulus of rock
    real(8)::Ess   !Modified elastic modulus
    real(8)::ny    !Coefficient of support effect
    real(8)::SigF  !Modified yield stress
    real(8)::t     !pipe thickness without margin
    real(8)::rm    !rm=(DD+t0)/2
    real(8)::r0d   !r0d=(DD+2*t0)/2
    real(8)::k0    !Gap size between outer surface of steel pipe and inner surface of backfilled concrete
    real(8)::SigN  !Axial stress of pipe at the moment of buckling
    real(8)::Pk    !Buckling pressure of steel pipe
    real(8)::Phi,Psi,Omega !Coefficients for Amstutz's formula
    real(8)::f1,f2,fi,x1,x2,xi !Variables for bisection method
    real(8)::pkr(1:5) !Memory of buckling pressure of steel pipe
    character::linebuf*1000,fmt1*200

    fmt1="(f10.3,',',f5.0,',',f10.3,5(',',f10.3))"

    Es=206000.0D0
    pos=0.3D0
    As=1.2D-5
    t0=30.0D0
    eps=1.5D0
    Jc=1.00D0
    Temp=20.0D0
    Bg=1.0D0

    write(6,*) '#D0,t0,rm/t0,Pk1,Pk2,Pk3,Pk4,Pk5'
    do nnn=35,140
        DD=2.0D0*dble(nnn)*t0
        do iii=1,5
            select case(iii)
                case(1) !HT100
                    Sigy=885.0D0
                    Siga=400.0D0
                case(2) !SHY685
                    Sigy=685.0D0
                    Siga=330.0D0
                case(3) !SM570
                    Sigy=450.0D0
                    Siga=240.0D0
                case(4) !SM490
                    Sigy=315.0D0
                    Siga=175.0D0
                case(5) !SM400
                    Sigy=235.0D0
                    Siga=130.0D0
            end select
            Ess=Es/(1.0D0-pos*pos)
            ny=1.5D0-0.5D0*1.0D0/(1.0D0+0.002D0*Es/Sigy)/(1.0D0+0.002D0*Es/Sigy)
            SigF=ny*Sigy/sqrt(1.0D0-pos+pos*pos)
            t=t0-eps
            rm=(DD+t0)/2.0D0
            r0d=(DD+2.0D0*t0)/2.0D0
            k0=(As*Temp+Bg*Siga*Jc/Es)*r0d/(1.0D0+Bg)
            !k0=0.4e-3*rm
            !Amstutz
            Phi=1.68D0
            Psi=0.25D0
            Omega=0.175D0
            x1=0.0D0
            x2=SigF
            !Bisection method
            do
                xi=0.5D0*(x1+x2)
                f1=func(k0,rm,x1,Ess,t,SigF,Phi,Psi)
                f2=func(k0,rm,x2,Ess,t,SigF,Phi,Psi)
                fi=func(k0,rm,xi,Ess,t,SigF,Phi,Psi)
                if(f1*fi<0.0)x2=xi
                if(fi*f2<0.0)x1=xi
                SigN=xi
                if(abs(x2-x1)<1.0D-6)exit
            end do
            Pk=SigN/(rm/t)/(1.0D0+Omega*2.0D0*rm/t*(SigF-SigN)/Ess)
            pkr(iii)=Pk
        end do
        write(linebuf,FMT=fmt1) DD,t0,0.5*DD/t0,(pkr(iii),iii=1,5)
        call del_spaces(linebuf)
        write (6,'(a)') trim(linebuf)
    end do

stop

contains

    real(8) function FUNC(k0,rm,SigN,Ess,t,SigF,Phi,Psi)
        real(8),intent(in)::k0,rm,SigN,Ess,t,SigF,Phi,Psi
        real(8)::a,b,fn

        a=(k0/rm+SigN/Ess)*(1.0D0+12.0D0*rm*rm/t/t*SigN/Ess)**1.5D0
        b=Phi*2.0D0*rm/t*(SigF-SigN)/Ess*(1.0D0-Psi*2.0D0*rm/t*(SigF-SigN)/Ess)
        fn=a-b
        FUNC=fn
    end function FUNC


    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_ams_jpn
