From 70eef0af8a540a8dc46985abee16f948820a093a Mon Sep 17 00:00:00 2001
From: Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>
Date: Tue, 1 Dec 2020 11:13:24 +0100
Subject: [PATCH] Trying to fix NaN for lmax > 30.

It is impossible to compute the crosssection if LMAX > 30.
This is due to a lack of precision in the djmn.f subroutine which
is written in simple precision. This commit is a first step to promote
this subroutine and all its dependensices to doucle precision.
---
 src/SConstruct                                | 11 ++++++-----
 src/msspec/spec/fortran/memalloc/allocation.f |  2 ++
 src/msspec/spec/fortran/memalloc/dim_mod.f    |  6 ++++--
 src/msspec/spec/fortran/memalloc/modules.f    | 14 ++++++++++++++
 4 files changed, 26 insertions(+), 7 deletions(-)

diff --git a/src/SConstruct b/src/SConstruct
index 45dae48..3739618 100644
--- a/src/SConstruct
+++ b/src/SConstruct
@@ -106,11 +106,12 @@ else:
     env.Append(F2PY=p.stdout.decode('utf8').strip())
 
 if GetOption('dbg'):
-    gfortran_env.Append(FORTRANFLAGS = ['-g', '-Wall', '-Wextra', '-Warray-temporaries', 
-                                        '-Wconversion', '-fbacktrace', '-ffree-line-length-0', 
-                                        '-fcheck=all', '-ffpe-trap=zero,overflow,underflow,invalid,denormal', 
-                                        '-finit-real=nan'],
-                        F2PY_OPTS = ['--noopt', '--debug-capi', '--debug'])
+    flags = ['-g', '-Wall', '-Wextra', '-Warray-temporaries', 
+             '-Wconversion', '-fbacktrace', '-ffree-line-length-0', 
+             '-fcheck=all', '-ffpe-trap=zero,overflow,underflow,invalid,denormal', 
+             '-finit-real=nan']
+    gfortran_env.Append(FORTRANFLAGS = flags,
+                        F2PY_OPTS = ['--noopt', '--debug-capi', '--debug', '--f77flags=-g -fcheck=all -fbacktrace -ffpe-trap=zero,overflow,underflow,invalid,denormal'])
 
 if GetOption('verbose'):
     env.Replace(FORTRANCOMSTR = "")
diff --git a/src/msspec/spec/fortran/memalloc/allocation.f b/src/msspec/spec/fortran/memalloc/allocation.f
index dd8a841..1ce3b3d 100644
--- a/src/msspec/spec/fortran/memalloc/allocation.f
+++ b/src/msspec/spec/fortran/memalloc/allocation.f
@@ -44,6 +44,7 @@
         USE COEFRLM_MOD
         USE EXPROT_MOD
         USE EXPFAC_MOD
+        USE DEXPFAC_MOD
         USE LOGAMAD_MOD
         USE EXPFAC2_MOD
         USE FACTSQ_MOD
@@ -155,6 +156,7 @@
         CALL ALLOC_COEFRLM()
         CALL ALLOC_EXPROT()
         CALL ALLOC_EXPFAC()
+        CALL ALLOC_DEXPFAC()
         CALL ALLOC_LOGAMAD()
         CALL ALLOC_EXPFAC2()
         CALL ALLOC_FACTSQ()
diff --git a/src/msspec/spec/fortran/memalloc/dim_mod.f b/src/msspec/spec/fortran/memalloc/dim_mod.f
index 475e2bd..e0432a4 100644
--- a/src/msspec/spec/fortran/memalloc/dim_mod.f
+++ b/src/msspec/spec/fortran/memalloc/dim_mod.f
@@ -58,8 +58,10 @@ C       ===============================================================
         
           NCG_M=4*LI_M+2
 
-          N_BESS=100*NL_M
-          N_GAUNT=5*NL_M
+C          N_BESS=100*NL_M
+C          N_GAUNT=5*NL_M
+          N_BESS=200*NL_M
+          N_GAUNT=10*NL_M
 
           NLTWO=2*NL_M
         END SUBROUTINE INIT_DIM
diff --git a/src/msspec/spec/fortran/memalloc/modules.f b/src/msspec/spec/fortran/memalloc/modules.f
index 3362336..9e8ab0d 100644
--- a/src/msspec/spec/fortran/memalloc/modules.f
+++ b/src/msspec/spec/fortran/memalloc/modules.f
@@ -764,6 +764,20 @@ C=======================================================================
         END SUBROUTINE ALLOC_EXPFAC
       END MODULE EXPFAC_MOD
 
+C=======================================================================
+      MODULE DEXPFAC_MOD
+        IMPLICIT NONE
+        REAL*8, ALLOCATABLE, DIMENSION(:,:) :: DEXPF
+      CONTAINS
+        SUBROUTINE ALLOC_DEXPFAC()
+          USE DIM_MOD
+          IF (ALLOCATED(DEXPF)) THEN
+            DEALLOCATE(DEXPF)
+          ENDIF
+          ALLOCATE(DEXPF(0:2*NL_M-2,0:2*NL_M-2))
+        END SUBROUTINE ALLOC_DEXPFAC
+      END MODULE DEXPFAC_MOD
+
 C=======================================================================
       MODULE LOGAMAD_MOD
         IMPLICIT NONE