Commit 93634994 authored by Markus Pahlow's avatar Markus Pahlow

replace cut-off with approximation for large arguments of W0

The argument of the exp function cannot be much greater than about 700, so we
had introduced a limit in the subroutine tchss.  The exp function is here part
of the argument of the lambertw function, and a good approximation for
lambertw(x) exists for large x, which has the advantage that we need not call
the exp function anymore.
parent 0ff28e31
......@@ -336,14 +336,22 @@ CONTAINS
CLASS(ocm), INTENT(INOUT) :: phy
REAL(dp), INTENT(INOUT) :: RC
REAL(dp), INTENT(IN) :: phin, zN
REAL(dp) :: paravg, PARb, aim, depfrac
REAL(dp) :: paravg, PARb, aim, depfrac, l1, l2
IF ((phy%PARs.GT.phy%I0).AND.(phy%mu0.GT.1E-30_dp)) THEN ! (A4)
PARb = MAX(phy%I0, phy%PARb)
depfrac = LOG(phy%PARs/PARb)
paravg = (phy%PARs - PARb)/depfrac
aim = MIN(phy%alpha*paravg/phy%mu0, 7.e2*phy%zC)
phy%tch = 1.E0_dp/phy%zC + (1E0_dp - lambertw((1.E0_dp + RC/(phy%daylen*phy%mu0))&
*EXP(1E0_dp + aim/phy%zC), 0, 0))/aim
aim = phy%alpha*paravg/phy%mu0
IF (aim.GT.7.e2_dp*phy%zC) THEN
! approximation of lambertw(x) for large x: W0(x) = l1 - l2 + l2/l1
! with l1 = ln(x), l2 = ln(l1)
l1 = LOG(1.E0_dp + RC/(phy%daylen*phy%mu0)) + 1E0_dp + aim/phy%zC
l2 = LOG(l1)
phy%tch = 1.E0_dp/phy%zC + (1E0_dp - l1 + l2 - l2/l1)/aim
ELSE
phy%tch = 1.E0_dp/phy%zC + (1E0_dp - lambertw((1.E0_dp + RC/(phy%daylen*phy%mu0))&
*EXP(1E0_dp + aim/phy%zC), 0, 0))/aim
END IF
phy%SI = physi(phy, phin, phy%PARs, PARb, depfrac)
phy%theta = phy%tch*phy%fC
RC = RC*depfrac/phin ! for calculation of phy%A in SR grow
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment