diff options
| -rw-r--r-- | src/IDAxiBrillBH.F | 1 | ||||
| -rw-r--r-- | src/bhbrill.x | 10 |
2 files changed, 9 insertions, 2 deletions
diff --git a/src/IDAxiBrillBH.F b/src/IDAxiBrillBH.F index d064efc..01000ef 100644 --- a/src/IDAxiBrillBH.F +++ b/src/IDAxiBrillBH.F @@ -62,6 +62,7 @@ c Perhaps this and others should go into cctk.h real*8 o70,o71,o72,o73,o74,o75,o76,o77,o78,o79 real*8 o80,o81,o82,o83,o84,o85,o86,o87,o88,o89 real*8 o90,o91,o92,o93,o94,o95,o96,o97,o98,o99 + integer i22 real*8 pi real*8 adm integer :: nx,ny,nz diff --git a/src/bhbrill.x b/src/bhbrill.x index b7721c9..270519b 100644 --- a/src/bhbrill.x +++ b/src/bhbrill.x @@ -19,8 +19,14 @@ o19 = -(o15*o18) o20 = o16 + o17 + o19 o21 = exp(o20) - o22 = -2.00000000000000d0 + n - o23 = o8**o22 + i22 = -2 + n +c The next lines causes a floating point exception on Alphas, +c if o22 is 0.0 but o8 is negative. +c This is because of y**x = exp(x*ln(y)). +c Changing the type of the exponent to integer fixes the problem. +c o22 = -2.00000000000000d0 + n +c o23 = o8**o22 + o23 = o8**i22 o24 = n**2 o25 = 2.00000000000000d0*etagrd(i)*eta0*o15 o26 = o16 + o19 + o25 |