fortran运行结果不合理，求大佬们的帮忙

sunjiyue

[Fortran] 纯文本查看 复制代码

program enthalpy
    implicit none
!This program aim to calculate the enthalpy of h2s hydrate dissociation in the presence of water


    double precision P,T,xh2s,n
    double precision dPdT,H1,V,Hdiss,Hsolve

    print*,"please imput temperature (K)"
    read*,T
    print*,"please imput pressure (MPa)"
    read*,P
    print*,"please imput hydation number"
    read*,n
    print*,"please imput mole fraction of h2s in liquid"
    read*,xh2s

    call dissenthalpy(Hdiss,T,dPdT,V)
    call solveenthalpy(Hsolve,xh2s,n)

    H1=Hdiss+Hsolve

    print*,"H1=",H1
    pause
!    return
    end

    subroutine Derivative(dPDT,T)
    double precision dPdT,T,e(4)
    data e/-436.1,4.6,-0.018,0.0000196/
    dPdT=e(1)+e(2)*T+e(3)*T**2+e(4)*T**3    !MPa/K
    print*,"dPdT=",dPdT
    return
    end

    subroutine Totalmolarvolume(V,Vliquid,Vgas,Vhydrate)     
    double precision V,Vliquid,Vgas,Vhydrate,xh2s,n,parameter1,parameter2,density2,l(4)                 !density1,molecularweight2,V4,molecularweight,h(6) 
    !data h/44.202,2544.564,-35.957,0.191,-4.488E-4,3.965E-7/
    data l/21.90583,-0.21947,7.33695E-4,-8.18347E-7/
    !calculating the molar volume change of solution
    parameter1=0.0000045
    parameter2=0.000018
    Vliquid=xh2s*parameter1+parameter2   !m3/mol

    !calculating the molar volume change of hydrate
   ! density1=1050  !kg/m3   !this value is guess
   ! molecularweight2=0.1461 !kg/mol
    Vhydrate=1.391E-4    !m3/mol, Vhydrate=molecularweight/density1

    !calculating the molar volume change of vapor
    !density2=h(1)+h(2)*T+h(3)*T**2+h(4)*T**3+h(5)*T**4+h(6)*T**5   !kg/m3
    !molecularweight=0.034 !kg/mol
   ! Vgas=molecularweight/density2   !m3/mol
     Vgas=l(1)+l(2)*T+l(3)*T**2+l(4)*T**3

    !calculating the total molar volume
    V=(1-(n*xh2s)/(1-xh2s))*Vgas+n*Vliquid-Vhydrate   !m3/mol

    print*,"Vliquid=",Vliquid
    print*,"Vhydrate=",Vhydrate
    print*,"Vgas=",Vgas
    print*,"V=",V
    return
    end

    subroutine dissenthalpy(Hdiss,T,dPdT,V)
    double precision Hdiss,T,dPdT,V,Vliquid,Vgas,Vhydrate
    call Derivative(dPDT,T)
    call Totalmolarvolume(V,Vliquid,Vgas,Vhydrate)
    Hdiss=1000000.0*dPdT*T*V
    print*,"Hdiss=",Hdiss
    return
    end

    subroutine solveenthalpy(Hsolve,xh2s,n)
    double precision Hsolve,xh2s,n
    Hsolve=17960.0*n*xh2s/(1-xh2s)   !J/mol
    print*,"Hsolve=",Hsolve
    return
    end

按照上面的运行结果，我输入T=274.15 K, P=0.118 MPa，n=6.2, xh2s=0.004的情况下，输出结果为：
dPdT=  -124.007293931873
Vliquid= -4.165183558045402E+056
Vhydrate=  1.390999968862161E-004
Vgas=  1.013072811542869E+018
V=  3.855278546349224E+118
Hdiss= -1.310663612078491E+129
Hsolve=   447.196787148594
H1= -1.310663612078491E+129

上述高亮部分结果不对，为什么会出现10的129次方呢，可我又不知道为什么，本人小白一枚，求各位大佬支招