[Module] module的调用问题

program main_sieve
use design
    implicit none
    real Ci,H0,LV,C,W,Di
    real ::H0_f,LV_f
    real,intrinsic::log
    write(*,*)"请依次输入液相流量,气相流量,液相密度,气相密度,初估板间距,初估泡沫层高,粘度:"
    read(*,*) Ls,Vs,Rl,Rv,Ht,hli,st
    write(*,*)"液相流量",Ls
    write(*,*)"气相流量",Vs
    write(*,*)"液相密度",Rl
    write(*,*)"气相密度",Rv
    write(*,*)"初估板间距",Ht
    write(*,*)"初估泡沫层高",hli
    write(*,*)"粘度",st
    H0=H0_f(Ht,hli)
    LV=LV_f(Ls,Vs,Rl,Rv)
    Ci = exp(-4.531+1.6562*(H0)+5.5496*(H0)**2-6.4695*(H0)**3+&
        (-0.474675+0.079*(H0)-1.39*(H0)**2+1.3212*(H0)**3)*log(LV)&
        +(-0.07291+0.088307*(H0)-0.49123*(H0)**2+0.43196*(H0)**3)*&
        (log(LV))**2)
    print *,"标准系数Ci:",Ci
    C = Ci/(20/st)**0.2
    print *,"经验系数C:",C
    W = 0.8*C*sqrt((Rl-Rv)/Rv)
    print *,"空速W:",W
    Di = sqrt(Vs/0.785/W)
    print *,"塔径初步估计为Di:",Di
    write(*,*)"开始对塔径进行校核"
    call C2
    write(*,*)"开始计算降液管"
    call C3
    write(*,*)"开始计算孔布置"
    call C4
    write(*,*)"开始计算压降"
    call C5
    write(*,*)"开始计算漏液,雾沫夹带"
    call C6
    write(*,*)"计算完成"
end

function H0_f(Ht,hli) 
implicit none
real Ht,hli
real H0_f
H0_f = Ht-hli
end function


function LV_f(Ls,Vs,Rl,Rv) 
implicit none
real Ls,Vs,Rl,Rv
real LV_f
LV_f = (Ls/Vs)*sqrt(Rl/Rv)
end function

module design
implicit none
real Ls,Vs,Rl,Rv,st,Ht
save Ls,Vs,Rl,Rv,st,Ht
contains

subroutine C2
implicit none
    real At,Af,yin,hf,hli,t,Ht,D,s
    real Wg,ent
    write(*,*)"基于塔径的预估结果，要进行圆整"
    read *,D
    write(*,*)"圆整后的塔径大小D:",D
    write(*,*)"输入比例因子"
    read *,yin
    At = 3.1415926*0.25*D**2
    s = (asin(yin)-yin*sqrt(1-yin**2))/3.1415926
    Af = At*s
    write(*,*)"精馏塔截面积为",At,"降液管面积为",Af
    Wg = Vs/(At-Af)
    write(*,*)"空速为",Wg
    hf = 2.5*hli
    ent = 0.22*(73/st)*(Wg/(12*(Ht-hf)))**3.2
    t = Ht*Af/Ls
    write(*,*)"雾沫夹带量",ent,"停留时间为",t
    if ((ent.LT.0.1).and.(t.GT.5)) then
        print *,"塔径为:",D,"是合理的,下面将做进一步的核算"
    else
        print *,"重复塔径估算,主要是修改板间距的大小"
    end if
end subroutine

subroutine C3
implicit none
    real lw,D,how,hl,hli,h0,Ls,Ei,E,hwi,hw,yin
        lw = D * yin
        Ei = 3600*Ls/(D*yin)**2.5
        write(*,*)"计算出的Ei值为",Ei
    write(*,*)"根据计算出的Ei值找出相应的E值"
        E = 1+(0.00012-0.00386*log(yin))*Ei
    write(*,*)"相应的E(收缩系数)值为",E
        how = 0.00284*E*(3600*Ls/lw)**0.6666667
        write(*,*)"堰上清液层高度为",how
    write(*,*)"下面进行堰高和板上清液层高度的计算,首先进行预估计算"
        hwi = hli-how
        write(*,*)"预估的堰高的值为",hwi
    write(*,*)"圆整堰高的值"
        read *,hw
        write(*,*)"圆整后的值为",hw
        hl = hw+how
        h0 = hw-0.013
        write(*,*)"板上清液层高度为",hl,"降液管底隙为",h0
end subroutine

subroutine C4
implicit none
    real d0,dd,Ws,Wc,Wd,A0,Aa,x,r,D,n,ns,yin,th,bl
    write(*,*)"请依次输入孔间距,孔径的大小"
    read *,dd,d0
    write(*,*)"孔间距为",dd
    write(*,*)"孔径为",d0
    write(*,*)"请依次输入降液管出口进口安定区宽度和精馏塔边缘安定区宽度"
    read *,Ws,Wc
        write(*,*)"降液管出口进口安定区宽度为",Ws
        write(*,*)"精馏塔边缘安定区宽度为",Wc
        bl = 0.907/(dd/d0)**2
        print *,"开孔区开孔率为：",bl
        Wd = 0.5*D*(1-sqrt(1-yin*yin))
        print *,"降液管宽度为:",Wd
        x = 0.5*D - (Ws+Wd)
        r = 0.5*D - Wc
        Aa = 2*((x*sqrt(r*r-x*x))+(r*r*asin(x/r)))
        print *,"鼓泡区面积为:",Aa
        A0 = Aa*bl
        n = 250000*A0/3.1415926
        print *,"孔数为",n
end subroutine

subroutine C5
implicit none
    real hc,hp,hs,hll,th,d0,Vs,A0,W0,W0m,Rv,Rl,C0,st,K,E0,hl
    write(*,*)"计算孔流系数,先输入板厚度"
    read *,th
        C0 = 0.8806-0.0677*(d0/th)+0.00732*(d0/th)**2-0.000338*(d0/th)**3
        print *,"孔流系数C0为",C0
    write(*,*)"计算干板压降"
        W0 = Vs/A0
        if (d0.LT.12) then
            hc = 0.0512*(Rv/Rl)*(W0/C0)**2
        else
            hc = 0.0512*(Rv/Rl)*(W0/(1.15*C0))**2
        end if
        print *,"干板压降的值为",hc
    write(*,*)"计算稳定性"
    write(*,*)"首先计算表面张力压头"
        hs = 4*st/(9.81*Rl*d0)
        print *,"表面张力压头为",hs
    write(*,*)"计算漏液点下的孔速"
        if (d0.GT.12) then
            W0m = 4.4*1.15*C0*sqrt((0.001+0.13*hl-hs)*Rl/Rv)
        else if (d0.LT.3) then
            W0m = 4.4*C0*sqrt((0.0051+0.05*hl)*Rl/Rv)
        else
            W0m = 4.4*C0*sqrt((0.0056+0.13*hl-hs)*Rl/Rv)
        end if
        print *,"漏液点下的孔速为",W0m
            K = W0/W0m
        print *,"稳定系数为",K
    write(*,*)"开始计算塔板压降"
    write(*,*)"首先计算液层有效阻力系数"    
        E0 = 0.652029*(W0*sqrt(Rv))**(-0.185964)
        hll = hl*E0
        print *,"液层有效阻力为",hll
        hp = hll + hs + hc
        print *,"塔板压降为",hp
end subroutine

subroutine C6
implicit none
    real Hd,hl,hdi,hp,Ls,D,lw,h0,A0w,hw,Ht,Af
    integer types
    write(*,*)"首先计算hd"
    write(*,*)"是否设置内堰是一个关键因素,是代表1,不是代表0,用types变量声明"
        lw = D*0.66
        read *,types
        select case(types)
            case(1)
            read *,A0w
            hdi = 0.1*(Ls/A0w)**2
            case(0)
            hdi = 0.153*(Ls/(lw*h0))**2
            case default
            write(*,*)"输入错误,请重新输入"
        end select
        print *,"局部阻力为",hdi
        Hd = hl+hdi+hp
        if ((2*Hd-hw).LT.0.3) then
            print *,"该工况下不会发生液泛情况"
        else
            print *,"液泛了,亲,你要重新设计了"
        end if
    write(*,*)"继续判断工艺是否合理"
        if (((Ls/Af).LT.0.1).AND.((Ht*Af/Ls).GT.5))  then
            print *,"恭喜你,你快设计完成了"
        else
            print *,"靠,挂了,需要调整一些参数了"
        end if
end subroutine
end module