program f90_ams_org
    !*********************************************
    ! Buckling pressure of steel pipe under external ground water pressure
    ! (Amstutz original)
    !*********************************************
    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,theta !Coefficients for Amstutz's formula
    real(8)::offset  !Offset at welded joint
    real(8)::DeltaD  !Initial deflection
    real(8)::m       !Parameter for offset
    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
 
            !DeltaD=0.0D0
            !offset=0.0D0

!            DeltaD=DD/300.0D0
!            offset=0.05D0*t0
             DeltaD=0.0D0; offset=0.0D0

            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

            !Initial deflection
            rm=rm*(1.0D0+1.52240775D0*DeltaD/DD)
            r0d=r0d*(1.0D0+1.52240775D0*DeltaD/DD)

            !Offset
            m=1.0D0+3.0D0*offset/t

            !Bisection method
            x1=Ess*(3.0D0*3.0D0-1.0D0)/12.0D0*t*t/rm/rm
            x2=SigF
            do
                xi=0.5D0*(x1+x2)
                SigN=x1
                call CALC(Ess,rm,t,SigN,Phi,Psi,Omega,theta)
                f1=FUNC(Ess,rm,t,k0,m,SigF,SigN,Phi,Psi)
                SigN=x2
                call CALC(Ess,rm,t,SigN,Phi,Psi,Omega,theta)
                f2=FUNC(Ess,rm,t,k0,m,SigF,SigN,Phi,Psi)
                SigN=xi
                call CALC(Ess,rm,t,SigN,Phi,Psi,Omega,theta)
                fi=FUNC(Ess,rm,t,k0,m,SigF,SigN,Phi,Psi)
                if(f1*fi<0.0D0)x2=xi
                if(fi*f2<0.0D0)x1=xi
                if(f1*f2>0.0D0)then
                    write(6,*) 'Imposible',x1,x2,f1,f2
                    stop
                end if
                if(abs(x2-x1)<1.0D-6)exit
            end do
            Pk=SigN/(rm/t)/(1.0D0+2.0D0*Omega*rm/t*(SigF-m*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

    subroutine CALC(Ess,rm,t,SigN,Phi,Psi,Omega,theta)
        real(8),parameter::pi=3.14159265358979323846D0
        real(8),intent(in)::Ess,rm,t,SigN
        real(8),intent(out)::Phi,Psi,Omega,theta
        real(8)::ang,cc,ee,alpha,delta,beta,cota,gamma

        ! ang=90->0 degree 
        ! ee*ang=270->180 degree 

        ee=sqrt(1.0D0+12.0D0*rm*rm/t/t*SigN/Ess)
        ang=270.0D0/ee
        do
            ang=ang-0.0001D0
            alpha=ang*pi/180.0D0
            cc=ee*tan(alpha)-tan(ee*alpha)
            if(0.0D0<=cc)exit
        end do
        delta=(ee*ee-1.0D0)*(1.0D0-cos(ee*alpha))
        beta=(ee-1.0D0/ee)*(ee*alpha*cos(ee*alpha)-sin(ee*alpha))
        if(ang==90.0D0)then
            cota=0.0D0
        else
            cota=1.0D0/tan(alpha)
        end if
        gamma=-ee*(ee*alpha-sin(ee*alpha)*cos(ee*alpha)&
              +ee*alpha*sin(ee*alpha)*sin(ee*alpha)/sin(alpha)/sin(alpha)&
              -ee*sin(ee*alpha)*sin(ee*alpha)*cota)

        Phi=ee*ee*ee*beta/pi/delta
        Psi=-gamma/4.0D0/beta/delta
        Omega=-cos(ee*alpha)/(1.0D0-cos(ee*alpha))
        theta=alpha/pi*180.0D0
    end subroutine CALC


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

        a=(k0/rm+SigN/Ess)*(1.0+12.0*rm*rm/t/t*SigN/Ess)**1.5
        b=2.0*Phi*rm/t*(SigF-m*SigN)/Ess*(1.0-2.0*Psi*rm/t*(SigF-m*SigN)/Ess)
        FUNC=a-b
    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_org
