我不知道这个问题是这里的主题或其他地方（或不上的话题在所有任意位置）。 我继承的Fortran 90代码，不会牛顿拉夫逊插值，其中温度的对数抵抗的压力对数内插。
内插的类型为
t = a ln(p) + b
和其中a，b被定义为
a = ln(tup/tdwn)/(alogpu - alogpd)
和
b = ln T - a * ln P
下面是测试程序。 它仅示出单次迭代。 但实际的程序运行了三个超过K，J和我的循环。 在现实中pthta是3D阵列（K，J，i）和thta是一维数组（k）的
program test
implicit none
integer,parameter :: dp = SELECTED_REAL_KIND(12,307)
real(kind=dp) kappa,interc,pres,dltdlp,tup,tdwn
real(kind=dp) pthta,alogp,alogpd,alogpu,thta,f,dfdp,p1
real(kind=dp) t1,resid,potdwn,potup,pdwn,pup,epsln,thta1
integer i,j,kout,n,maxit,nmax,resmax
kappa = 2./7.
epsln = 1.
potdwn = 259.39996337890625
potup = 268.41687198359159
pdwn = 100000.00000000000
pup = 92500.000000000000
alogpu = 11.43496392350051
alogpd = 11.512925464970229
thta = 260.00000000000000
alogp = 11.512925464970229
! known temperature at lower level
tdwn = potdwn * (pdwn/100000.)** kappa
! known temperature at upper level
tup = potup *(pup/100000.)** kappa
! linear change of temperature wrt lnP between different levels
dltdlp = dlog(tup/tdwn)/(alogpu-alogpd)
! ln(T) value(intercept) where Pressure is 1 Pa and assume a linear
! relationship between P and T
interc = dlog(tup) - dltdlp*alogpu
! Initial guess value for pressure
pthta = exp((dlog(thta)-interc-kappa*alogp)/(dltdlp-kappa))
n=0
1900 continue
!First guess of temperature at intermediate level
t1 = exp(dltdlp * dlog(pthta)+interc)
!Residual error when calculating Newton Raphson iteration(Pascal)
resid = pthta - 100000.*(t1/thta)**(1./kappa)
print *, dltdlp,interc,t1,resid,pthta
if (abs(resid) .gt. epsln) then
n=n+1
if (n .le. nmax) then
! First guess of potential temperature given T1 and
! pressure level guess
thta1 = t1 * (100000./pthta)**kappa
f= thta - thta1
dfdp = (kappa-dltdlp)*(100000./pthta)**kappa*exp(interc + (dltdlp -1.)*dlog(pthta))
p1 = pthta - f/dfdp
if (p1 .le. pdwn) then
if (p1 .ge. pup) then
pthta = p1
goto 1900
else
n = nmax
end if
end if
else
if (resid .gt. resmax) resmax = resid
maxit = maxit+1
goto 2100
end if
end if
2100 continue
end program test
当你从一个数据文件的真实数据运行此程序渣油的值以下
2.7648638933897018E-010
它不会对整个执行太大的差别。 大部分值的范围是
1E-10 and 1E-12
因此，考虑这些值IF条件下
IF (abs(resid) .gt. epsln)
不会被调用和牛顿Raphson迭代永远不会被执行。 所以，我看着两种方式来得到这个工作。 一是删除这两个步骤中的指数通话
pthta = exp((dlog(thta)-interc-kappa*alogp)/(dltdlp-kappa))
t1 = exp(dltdlp * dlog(pthta)+interc)
即保持对数空间中的一切，牛顿Raphson迭代完成后采取的指数。 这部分没有收敛没有问题。
我努力使这项工作的另一种方式是截断
t1 = exp(dltdlp * dlog(pthta)+interc)
当我把它截断成整数的渣油的变化值显着从1E-10到813，我不知道如何截断该函数调用导致如此大的价值变动。 截断该结果确实得到圆满完成。 所以，我不知道这是更好的方式来进一步进行。
我如何决定这将是更好的方式来处理这个？
从研究的角度来看，我会说你的第一个解决方案很可能是更合适的方法。 在物理模拟，每个人都应该与通过定义总是积极性质的对数工作。 在上面的代码中，这些将是温度和压力。 严格正定的物理变量往往导致溢出和下溢的计算，无论你用Fortran语言或任何其他的编程语言，或任何可能的变量类型的。 如果事情都可能发生，它会发生。
这是一个关于其他物理量也同样如此，例如，能量（伽马射线脉冲串的典型能量为〜10 ^ 54尔格）中，任意尺寸的对象的体积（半径的100维球的体积10米是〜10 ^ 100），或者甚至概率（在许多统计问题的似然函数可采取〜10 ^值{ - 1000}或更小）。 正定变量与工作日志变换将使你的代码来处理数字一样大〜10 ^ 10 ^ 307（一个双精度变量）。
几点注意事项还就在你的代码中使用的Fortran语言的语法：
变量RESMAX在代码中使用，无需初始化。
当为变量分配值，它指定的那种字面常量恰当，否则，程序结果可能受影响的是很重要的。 例如，这里是你的原代码与英特尔Fortran编译器2018在调试模式下进行编译的输出：
-0.152581477302743 7.31503025786548 259.608693509165 -3.152934473473579E-002 99474.1999921620
这里是相同的代码的输出，但与那种参数后缀的所有文字常数_dp （见下面的代码的修订版）：
-0.152580456940175 7.31501855886952 259.608692604963 -8.731149137020111E-011 99474.2302854451
从在此答案订正代码的输出是从原始代码中的上述问题的输出略有不同。
有没有必要使用.gt. ， .ge. ， .le. ， .lt. ，...，进行比较。 这些都是传统的FORTRAN语法，据我所知。 使用而不是更吸引人符号（ < ， > ， <= >= ， == ）用于比较。
有没有必要使用GOTO在Fortran程序语句。 这又是传统的FORTRAN。 通常情况下，简洁大方的DO-循环和if块可以代替GOTO语句，就像下面的修改后的代码。
有没有必要使用特定类型的内部函数Fortran语言了（如dexp ， dlog ，...对于双精度）。 Fortran语言的几乎所有的（也许是全部）内部函数具有通用名称（ exp ， log在目前Fortran标准，...）。
以下是该方案在这个问题上的修订，能够解决上述所有过时的语法，以及应对极端大或小正定变量的问题（我可能走得太远了对数变换的一些变量会不会导致溢出或下溢，但在这里，我的目的是为了只显示后面正定变量以及如何处理与他们的算术对数变换的逻辑，而可能造成溢/下溢/ error_in_results）。
program test
implicit none
integer,parameter :: dp = SELECTED_REAL_KIND(12,307)
real(kind=dp) kappa,interc,pres,dltdlp,tup,tdwn
real(kind=dp) pthta,alogp,alogpd,alogpu,thta,f,dfdp,p1
real(kind=dp) t1,resid,potdwn,potup,pdwn,pup,epsln,thta1
integer i,j,kout,n,maxit,nmax,resmax
real(kind=dp) :: log_resmax, log_pthta, log_t1, log_dummy, log_residAbsolute, sign_of_f
real(kind=dp) :: log_epsln, log_pdwn, log_pup, log_thta, log_thta1, log_p1, log_dfdp, log_f
logical :: residIsPositive, resmaxIsPositive, residIsBigger
log_resmax = log(log_resmax)
resmaxIsPositive = .true.
kappa = 2._dp/7._dp
epsln = 1._dp
potdwn = 259.39996337890625_dp
potup = 268.41687198359159_dp
pdwn = 100000.00000000000_dp
pup = 92500.000000000000_dp
alogpu = 11.43496392350051_dp
alogpd = 11.512925464970229_dp
thta = 260.00000000000000_dp
alogp = 11.512925464970229_dp
log_epsln = log(epsln)
log_pup = log(pup)
log_pdwn = log(pdwn)
log_thta = log(thta)
! known temperature at lower level
tdwn = potdwn * (pdwn/1.e5_dp)**kappa
! known temperature at upper level
tup = potup *(pup/1.e5_dp)** kappa
! linear change of temperature wrt lnP between different levels
dltdlp = log(tup/tdwn)/(alogpu-alogpd)
! ln(T) value(intercept) where Pressure is 1 Pa and assume a linear
! relationship between P and T
interc = log(tup) - dltdlp*alogpu
! Initial guess value for pressure
!pthta = exp( (log(thta)-interc-kappa*alogp) / (dltdlp-kappa) )
log_pthta = ( log_thta - interc - kappa*alogp ) / ( dltdlp - kappa )
n=0
MyDoLoop: do
!First guess of temperature at intermediate level
!t1 = exp(dltdlp * log(pthta)+interc)
log_t1 = dltdlp * log_pthta + interc
!Residual error when calculating Newton Raphson iteration(Pascal)
!resid = pthta - 1.e5_dp*(t1/thta)**(1._dp/kappa)
log_dummy = log(1.e5_dp) + ( log_t1 - log_thta ) / kappa
if (log_pthta>=log_dummy) then
residIsPositive = .true.
log_residAbsolute = log_pthta + log( 1._dp - exp(log_dummy-log_pthta) )
else
residIsPositive = .false.
log_residAbsolute = log_dummy + log( 1._dp - exp(log_pthta-log_dummy) )
end if
print *, "log-transformed values:"
print *, dltdlp,interc,log_t1,log_residAbsolute,log_pthta
print *, "non-log-transformed values:"
if (residIsPositive) print *, dltdlp,interc,exp(log_t1),exp(log_residAbsolute),exp(log_pthta)
if (.not.residIsPositive) print *, dltdlp,interc,exp(log_t1),-exp(log_residAbsolute),exp(log_pthta)
!if (abs(resid) > epsln) then
if ( log_residAbsolute > log_epsln ) then
n=n+1
if (n <= nmax) then
! First guess of potential temperature given T1 and
! pressure level guess
!thta1 = t1 * (1.e5_dp/pthta)**kappa
log_thta1 = log_t1 + ( log(1.e5_dp)-log_pthta ) * kappa
!f = thta - thta1
if ( log_thta>=thta1 ) then
log_f = log_thta + log( 1._dp - exp( log_thta1 - log_thta ) )
sign_of_f = 1._dp
else
log_f = log_thta + log( 1._dp - exp( log_thta - log_thta1 ) )
sign_of_f = 1._dp
end if
!dfdp = (kappa-dltdlp)*(1.e5_dp/pthta)**kappa*exp(interc + (dltdlp -1._dp)*log(pthta))
! assuming kappa-dltdlp>0 is TRUE always:
log_dfdp = log(kappa-dltdlp) + kappa*(log(1.e5_dp)-log_pthta) + interc + (dltdlp -1._dp)*log_pthta
!p1 = pthta - f/dfdp
! p1 should be, by definition, positive. Therefore:
log_dummy = log_f - log_dfdp
if (log_pthta>=log_dummy) then
log_p1 = log_pthta + log( 1._dp - sign_of_f*exp(log_dummy-log_pthta) )
else
log_p1 = log_dummy + log( 1._dp - sign_of_f*exp(log_pthta-log_dummy) )
end if
!if (p1 <= pdwn) then
if (log_p1 <= log_pdwn) then
!if (p1 >= pup) then
if (log_p1 >= log_pup) then
log_pthta = log_p1
cycle MyDoLoop
else
n = nmax
end if
end if
else
!if (resid > resmax) resmax = resid
residIsBigger = ( residIsPositive .and. resmaxIsPositive .and. log_residAbsolute>log_resmax ) .or. &
( .not.residIsPositive .and. .not.resmaxIsPositive .and. log_residAbsolute<log_resmax ) .or. &
( residIsPositive .and. .not. resmaxIsPositive )
if ( residIsBigger ) then
log_resmax = log_residAbsolute
resmaxIsPositive = residIsPositive
end if
maxit = maxit+1
end if
end if
exit MyDoLoop
end do MyDoLoop
end program test
下面是此方案，该方案与原始代码的输出非常一致的样本输出：
log-transformed values:
-0.152580456940175 7.31501855886952 5.55917546888014
-22.4565579499410 11.5076538974964
non-log-transformed values:
-0.152580456940175 7.31501855886952 259.608692604963
-1.767017293116268E-010 99474.2302854451
为了便于比较，这里是从原代码的输出：
-0.152580456940175 7.31501855886952 259.608692604963
-8.731149137020111E-011 99474.2302854451
