大家好，负数取小数指数的问题
Program Gas exchangeC real AaO2, AbO2, XNO2, DLO2, O2MM, PIO2,PaCO,CVO2
C real AaCO,AbCO, XNCO, DLCO, COMM, PICO, PaO2,CVCO
C real AaCO2,AbCO2,XNCO2,DLCO2,CO2MM,PICO2,CVCO2
COMMON /parametersO2/ AaO2, AbO2, XNO2, DLO2, O2MM,PIO2,PaCO0,CVO2
COMMON /parametersCO/ AaCO,AbCO, XNCO, DLCO, COMM, PICO,PaO20,CVCO
COMMON /parametersCO2/ AaCO2,AbCO2,XNCO2,DLCO2,CO2MM,PICO2,CVCO2
COMMON /VAR/QV, QB, QE, QI ,Cos
DIMENSION PBloodO2(20)
DIMENSION PBloodCO(20)
DIMENSION PBloodCO2(20)
DIMENSION CIO2(20)
DIMENSION CICO(20)
DIMENSION CICO2(20)
! Read the constants
OPEN(UNIT=1, FILE='CONSTANT.dat', STATUS='OLD')
OPEN(UNIT=2, FILE='PAlveolarO2.dat', STATUS='UNKNOWN')
OPEN(UNIT=3, FILE='PAlveolarCO.dat', STATUS='UNKNOWN')
OPEN(UNIT=4, FILE='PAlveolarCO2.dat', STATUS='UNKNOWN')
OPEN(UNIT=5, FILE='COHB.dat', STATUS='UNKNOWN')
OPEN(UNIT=6, FILE='PBLOODO2.dat', STATUS='UNKNOWN')
OPEN(UNIT=7, FILE='PBLOODCO.dat', STATUS='UNKNOWN')
OPEN(UNIT=8, FILE='PBLOODCO2.dat', STATUS='UNKNOWN')
OPEN(UNIT=9, FILE='CIO2.dat', STATUS='UNKNOWN')
OPEN(UNIT=10, FILE='CICO.dat', STATUS='UNKNOWN')
OPEN(UNIT=11, FILE='CICO2.dat', STATUS='UNKNOWN')
write (2,*) "O2 partial pressure in Alveolar"
write (3,*) "CO partial pressure in Alveolar"
write (4,*) "CO2 partial pressure in Alveolar"
write (5,*) "COHB percentage in blood"
Read(1,*)AaO2,AbO2,XNO2,DLO2,O2MM
Read(1,*)AaCO,AbCO,XNCO,DLCO,COMM
Read(1,*)AaCO2, AbCO2, XNCO2, DLCO2, CO2MM
! Control parameters
Write (*,*) "O2 exposure concentration =(unit in mmHg)"
Read(*,*) PIO2
Write (*,*) "CO exposure concentration =(unit in mmHg)"
Read(*,*) PICO
Write (*,*) "CO2 exposure concentration = (unit in mmHg)"
Read(*,*) PICO2
Write (*,*) "Volume of blood in the body = (unit in ml)"
Read(*,*) QV
Write (*,*) "Blood flow rate = (unit in ml/min)"
Read(*,*) QB
Write (*,*) "Ventilation rate = (unit in ml/min)"
Read(*,*) QE
! Control parameters
Write (*,*) "Cycle number = "
Read (*,*) ICycle
Write (*,*) "Iterative number = "
Read (*,*) Iteration
Write (*,*) "Vascular Compartments = "
Read (*,*) Cos
CVO2 = 0.155742
CVCO = 0.0000000001
CVCO2 = 0.54382
DO I = 1, ICycle
DO J = 1, Cos
PaO20= 0.0000000001
PaCO0= 0.0000000001
PaCO20 = 0.0000000001
DO K = 1, Iteration
PaO2= PaO20- FunO2(PaO20)/ dFunO2(PaO20)
PaCO= PaCO0- FunCO(PaCO0)/ dFunCO(PaCO0)
PaCO2 = PaCO20 - FunCO2(PaCO20)/ dFunCO2(PaCO20)
PaO20= PaO2
PaCO0= PaCO
PaCO20 = PaCO2
ENDDO
SO2= ST(PaO20+COMM*PaCO0) / ( 1+(COMM*PaCO0)/PaO20 )
CaO2 = AbO2*PaO20 + XNO2*SO2
CIO2(J) = CaO2
CVO2=CIO2(J)
SCO= ST(PaO20+COMM*PaCO0) / ( 1+PaO20/(COMM*PaCO0) )
CaCO = AbCO*PaCO0 + XNCO*SCO
CICO(J) = CaCO
CVCO=CICO(J)
SCO2= 0.063*PaCO20**(5000/60181)
+ /(1+0.63*PaCO20**(5000/60181)+0.275*ST(PaO20+COMM*PaCO0) )
CaCO2 = AbCO2*PaCO20 + XNCO2*SCO2
CICO2(J) = CaCO2
CVCO2=CICO2(J)
PBLOODO2(J) = PaO20
PBLOODCO(J) = PaCO0
PBLOODCO2(J) = PaCO20
WRITE (6,*) J, PBLOODO2(J)
WRITE (7,*) J, PBLOODCO(J)
WRITE (8,*) J, PBLOODCO2(J)
WRITE (9,*) J, CIO2(J)
WRITE (10,*) J, CICO(J)
WRITE (11,*) J, CICO2(J)
END DO
QI = QE + QB*( (CaO2-CVO2)+(CaCO-CVCO)
+ + (CaCO2-CVCO2) ) / 0.826184
CVO2 = CaO2 - 250/QB
CVCO = CaCO + 0.007/QB
CVCO2 = CaCO2 + 250/QB
PaO2AVER=SUM(PBLOODO2)/Cos
PAlveolarO2 = ( QI*AaO2*PIO2 + DLO2*PaO2AVER) /(QE*AaO2 +DLO2)
PaCOAVER=SUM(PBLOODCO)/Cos
PAlveolarCO = ( QI*AaCO*PICO + DLCO*PaCOAVER )/(QE*AaCO +DLCO)
PaCO2AVER=SUM(PBLOODCO2)/Cos
PAlveolarCO2 = ( QI*AaCO2*PICO2 + DLCO2*PaCO2AVER )
+ / ( QE*AaCO2 + DLCO2 )
Write (2,*) I, PAlveolarO2
Write (3,*) I, PAlveolarCO
Write (4,*) I, PAlveolarCO2
Write (5,*) I, SCO
ENDDO
End Program
FUNCTION FunO2(x)
COMMON /parametersO2/ AaO2, AbO2, XNO2, DLO2, O2MM,PIO2,PaCO0,CVO2
COMMON /VAR/QV,QB, QE, QI,Cos
FunO2 = ( QB*AbO2 + DLO2/Cos*QE*AaO2 / (DLO2/Cos+QE*AaO2) )*x
+ +XNO2*x*QB / (X+O2MM*PaCO0)
+ *( (X+O2MM*PaCO0)**3 + 150*(X+O2MM*PaCO0) )
+ /( (X+O2MM*PaCO0)**3 + 150*(X+O2MM*PaCO0) + 23400 )
+ -DLO2/Cos*QI*AaO2*PIO2 / (DLO2/Cos+QE*AaO2) - QB*CVO2
END FUNCTION
FUNCTION dFunO2(x)
COMMON /parametersO2/ AaO2, AbO2, XNO2, DLO2, O2MM,PIO2,PaCO0,CVO2
COMMON /VAR/QV,QB, QE, QI,Cos
dFunO2 =( QB*AbO2 + DLO2/Cos*QE*AaO2 / (DLO2/Cos+QE*AaO2) )
+ +XNO2*QB
+ *( ( (X+O2MM*PaCO0)**2 + 150 )
+ * ( (X+O2MM*PaCO0)**3 + 150*(X+O2MM*PaCO0) + 23400 )
+ + 2*X*(X+O2MM*PaCO0)
+ * ( (X+O2MM*PaCO0)**3 + 150*(X+O2MM*PaCO0) + 23400 )
+ -X*( (X+O2MM*PaCO0)**2 + 150 )
+ *( 3*(X+O2MM*PaCO0)**2 + 150) )
+ / ( (X+O2MM*PaCO0)**3 + 150*(X+O2MM*PaCO0) + 23400 )**2
END FUNCTION
FUNCTION FunCO(x)
COMMON /parametersCO/ AaCO, AbCO, XNCO, DLCO, COMM,PICO,PaO20,CVCO
COMMON /VAR/QV, QB, QE, QI,Cos
FunCO = ( QB*AbCO + DLCO/Cos*QE*AaCO / (DLCO/Cos+QE*AaCO) )*x
+ +XNCO*COMM*QB*x / (PaO20+COMM*x)
+ *( (PaO20+COMM*x)**3 + 150*(PaO20+COMM*x) )
+ /( (PaO20+COMM*x)**3 + 150*(PaO20+COMM*x) +23400 )
+ -DLCO/Cos*QI*AaCO*PICO/ (DLCO/Cos+QE*AaCO)
+ -QB*CVCO
END FUNCTION
FUNCTION dFunCO(x)
COMMON /parametersCO/ AaCO, AbCO, XNCO, DLCO, COMM,PICO,PaO20,CVCO
COMMON /VAR/QV,QB, QE, QI,Cos
dFunCO = ( QB*AbCO + DLCO/Cos*QE*AaCO / (DLCO/Cos+QE*AaCO) )
+ + XNCO*COMM*QB
+ *( ( (PaO20+COMM*x)**2 + 2*COMM*x*(PaO20+COMM*x)+ 150)
+ *( (PaO20+COMM*x)**3 + 150*(PaO20+COMM*x) + 23400 )
+ - ( (PaO20+COMM*x)**2 + 150 )
+ *x*COMM*( 3*(PaO20+COMM*x)**2 + 150 ) )
+ /( (PaO20+COMM*x)**3 + 150*(PaO20+COMM*x) + 23400 )**2
END FUNCTION
FUNCTION FunCO2(X)
COMMON /parametersCO2/ AaCO2,AbCO2,XNCO2,DLCO2,CO2MM,PICO2,CVCO2
COMMON /parametersCO/ AaCO,AbCO, XNCO, DLCO, COMM, PICO,PaO20,CVCO
COMMON /parametersO2/ AaO2, AbO2, XNO2, DLO2, O2MM,PIO2,PaCO0,CVO2
COMMON /VAR/QV,QB, QE, QI,Cos
FunCO2 = ( QB*AbCO2 + (DLCO2/Cos*QE*AaCO2)/(DLCO2/Cos+QE*AaCO2))*X
+ + QB*XNCO2*(0.063*x**(50000/60181))
+ /( 1 + 0.63*x**(50000/60181) + 0.275*ST(PaO20+CO2MM*PaCO0) )
+ - DLCO2/Cos*QI*AaCO2*PICO2 / (DLCO2/Cos+QE*AaCO2)
+ - QB*CVCO2
END FUNCTION
FUNCTION dFunCO2(X)
COMMON /parametersCO2/ AaCO2,AbCO2,XNCO2,DLCO2,CO2MM,PICO2,CVCO2
COMMON /parametersCO/ AaCO,AbCO, XNCO, DLCO, COMM, PICO,PaO20,CVCO
COMMON /parametersO2/ AaO2, AbO2, XNO2, DLO2, O2MM,PIO2,PaCO0,CVO2
COMMON /VAR/QV,QB, QE, QI,Cos
dFunCO2 = QB*AbCO2 + (DLCO2/Cos*QE*AaCO2) / (DLCO2/Cos + QE*AaCO2)
+ + QB*XNCO2*(0.05234*( 1 + 0.063*x**(50000/60181)
+ + 0.275*ST(PaO20+CO2MM*PaCO0 ) - 0.063*x**(50000/60181) )
+ / ( x**0.16917*( 1 + 0.063*x**(50000/60181)
+ + 0.275*ST(PaO20+CO2MM*PaCO0) )**2 ) )
END FUNCTION
FUNCTION ST(x)
ST = ( x**3 + 150*x ) / ( x**3 + 150*x + 23400 )
END FUNCTION
由于是固定格式，程序在我这边编译的时候，只有部分行长度超过了允许长度，需要往左侧挪一下，挪完之后，编译再没出现任何错误，至于运行期错误，楼主没有提供数据文件，所以无法知道是什么错误。
ps 本论坛的代码可以高亮显示的，楼主可以编辑下帖子 代码在 编辑菜单的 “<> ”输入 哦，知道啦，谢谢你的指导。还有一个要读入的constant文件，屏幕输入的数据依次是0,0,149,5000,5000,5000,600,10。
我这里也能运行程序，就是只能循环两次自动跳出来，然后就出现First-chance exception in GAS.exe (KERNELBASE.DLL): 0xC0000090: Float Invalid Operation.检查了好几天了也不知道哪里出问题。希望你能帮我解答一下:-lol 9个屏幕输入，你只给了8个。 vvt 发表于 2014-12-9 11:59
9个屏幕输入，你只给了8个。
0，0,149,5000,5000,5000,600,100,10谢谢 珊瑚虫 发表于 2014-12-8 20:10
由于是固定格式，程序在我这边编译的时候，只有部分行长度超过了允许长度，需要往左侧挪一下，挪完之后，编 ...
哦，知道啦，谢谢你的指导。还有一个要读入的constant文件，屏幕输入的数据依次是0,0,149,5000,5000,5000,600,10。
我这里也能运行程序，就是只能循环两次自动跳出来，然后就出现First-chance exception in GAS.exe (KERNELBASE.DLL): 0xC0000090: Float Invalid Operation.检查了好几天了也不知道哪里出问题。希望你能帮我解答一下
你的问题在于：
FunCO2 函数中
FunCO2 = ( QB*AbCO2 + (DLCO2/Cos*QE*AaCO2)/(DLCO2/Cos+QE*AaCO2))*x
+ + QB*XNCO2*(0.063*x**0.83083)这里有 x ** 0.83083 这个语句。 Fortran 中的指数为小数时，x 不能为负数。
这个 0.83083 是什么含义？可否写成分数的形式？
另外，为了安全，我建议你把 **2.0 都改为 **2，**3.0 都改为 **3
指数如果可以是整数，就不要写成小数。
vvt 发表于 2014-12-9 23:20
你的问题在于：
FunCO2 函数中
我按照要求改了一下，还是出现一样的问题。0.83083是该函数中给定的值，没有整数形式，:-P 本帖最后由 vvt 于 2014-12-10 16:06 编辑
负数不能开偶次方，所以负数也就不能取偶数倒数的指数。
这是数学问题决定的。不是 fortran 决定的。
解决方法：
1.寻求 0.83083 的分数形式。比如 多少除以多少这种。
2.保证 x 大于 0。
没有别的办法。
如果这两个办法你都无法做到，那我就帮不了你了。 vvt 发表于 2014-12-10 16:04
负数不能开偶次方，所以负数也就不能取偶数倒数的指数。
这是数学问题决定的。不是 fortran 决定的。
非常感谢你的指导。:-P第一点我改正了。第二点由于X是通过牛顿迭代法求出来的，不知道是不是负值。可是根据实际情况不可能是负的。
页:
[1]
2