Fortran Coder

查看: 477|回复: 4

[数值问题] fortran运行结果不合理,求大佬们的帮忙

[复制链接]

1

帖子

1

主题

0

精华

新人

F 币
12 元
贡献
5 点
发表于 2020-1-21 09:40:56 | 显示全部楼层 |阅读模式
[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次方呢,可我又不知道为什么,本人小白一枚,求各位大佬支招


回复

使用道具 举报

1543

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1454 元
贡献
956 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2020-1-22 21:16:59 | 显示全部楼层
Totalmolarvolume 函数里 xh2s,T,n 这三个变量,没有定义,也没有赋值。就直接使用,当然会出错。

建议:
1. 要写 implicit none,每个程序单元都要写。
2. 学习 debug 单步调试。

167

帖子

14

主题

0

精华

大宗师

F 币
4011 元
贡献
818 点
发表于 2020-1-24 01:56:56 | 显示全部楼层
fcode 发表于 2020-1-22 21:16
Totalmolarvolume 函数里 xh2s,T,n 这三个变量,没有定义,也没有赋值。就直接使用,当然会出错。

建议 ...

实际上,只要主程序和每个module里的开头写上就可以了吧
我实际测试的结果是,只要这样写了,之后的所有子程序和函数里都不能有未声明的变量了

1543

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1454 元
贡献
956 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2020-1-24 09:04:34 | 显示全部楼层
liudy02 发表于 2020-1-24 01:56
实际上,只要主程序和每个module里的开头写上就可以了吧
我实际测试的结果是,只要这样写了,之后的所有 ...

Implicit None 如果写在module里,那么module里contains的过程都有效。

但是楼主玩的都是外部子程序,需要每个程序单元都写

167

帖子

14

主题

0

精华

大宗师

F 币
4011 元
贡献
818 点
发表于 2020-1-25 02:14:40 | 显示全部楼层
fcode 发表于 2020-1-24 09:04
Implicit None 如果写在module里,那么module里contains的过程都有效。

但是楼主玩的都是外部子程序,需 ...

原来如此,谢谢~
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

QQ|捐赠本站|Archiver|关于我们 About Us|群聊|Fcode

GMT+8, 2020-9-19 08:10

Powered by Discuz! X3.2

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表