From a research perspective, I'd say your first solution is likely the more appropriate approach. In a physical simulation, one should always work with the logarithm of the properties that are by-definition always positive. In the above code, these would be temperature and pressure. Strictly positive-definite physical variables often result in overflow and underflow in computation, whether you use Fortran or any other programming language, or any possible variable kind. If something can happen, it will happen.
This is true about other physical quantities as well, for example, energy (the typical energy of a Gamma-Ray-Burst is ~10^54 ergs), volume of objects in arbitrary dimensions (the volume of a 100-dimensional sphere of radius 10meters is ~ 10^100), or even probability (the likelihood function in many statistical problems can take values of ~10^{-1000} or less). Working with log-transform of positive-definite variables would enable your code to handle numbers as big as ~10^10^307 (for a double precision variable).
A few notes also regarding the Fortran syntax used in your code:
The variable RESMAX is used in your code without initialization.
When assigning values to variables, it is important to specify the kind of the literal constants appropriately, otherwise, the program results might be affected. For example, here is the output of your original code compiled with Intel Fortran Compiler 2018 in debug mode:
-0.152581477302743 7.31503025786548 259.608693509165
-3.152934473473579E-002 99474.1999921620
And here is the same code's output, but with all literal constants suffixed with the kind parameter _dp (see the revised version of your code below):
-0.152580456940175 7.31501855886952 259.608692604963
-8.731149137020111E-011 99474.2302854451
The output from the revised code in this answer is slightly different from the output of the original code in the above question.
There is no need to use .gt., .ge., .le., .lt., ..., for comparison. These are legacy FORTRAN syntax, as far as I am aware. Use instead the more attractive symbols ( < , > , <= , >= , == ) for comparison.
There is no necessity to use a GOTO statement in a Fortran program. This is again legacy FORTRAN. Frequently, simple elegant do-loops and if-blocks can replace GOTO statements, just as in the revised code below.
There is no need to use kind-specific intrinsic functions in Fortran anymore (such as dexp, dlog, ... for double precision). Almost all (and perhaps all) of Fortran intrinsic functions have generic names (exp, log, ...) in the current Fortran standard.
The following is a revision of the program in this question, that resolves all of the above obsolete syntax, as well as the problem of dealing with extremely large or small positive-definite variables (I probably went too far in log-transforming some variables that would never cause overflow or underflow, but my purpose here was to just show the logic behind log-transformation of positive-definite variables and how to deal with their arithmetics without potentially causing overflow/underflow/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
Here is a sample output of this program, which agrees well with the output of the original code:
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
For comparison, here is the output from the original code:
-0.152580456940175 7.31501855886952 259.608692604963
-8.731149137020111E-011 99474.2302854451
