From 71ea6d6a38a23051d18a4b6f8aee0ff2ab48bc32 Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Thu, 26 Oct 2023 15:02:16 +0200
Subject: [PATCH] Enable debug logging in sph

---
 src/sphere/sph.f | 33 +++++++++++++++++++++++++++------
 1 file changed, 27 insertions(+), 6 deletions(-)

diff --git a/src/sphere/sph.f b/src/sphere/sph.f
index 5eb62c2f..69d8e3b3 100644
--- a/src/sphere/sph.f
+++ b/src/sphere/sph.f
@@ -83,9 +83,9 @@ CCC   DIMENSION DC0M(NSPH,NSHL-NTL),XIV(NXI)
  6991 FORMAT('========== JXI =',I3,' ====================')
  6992 FORMAT('********** JTH =',I3,', JPH =',I3,
      1', JTHS =',I3,', JPHS =',I3,' ********************')
-      IR=5
-      IW=6
-      IT=7
+      IR=25
+      IW=26
+      IT=27
       ITIN=17
 CCC
 CCC   OTHER COMMENTS IN FILE GLOBALCOMS
@@ -270,6 +270,7 @@ C
 
       CS0=0.25D0*VK*VK*VK/PIGH
       CALL SSCR0(TFSAS,NSPH,LM,VK,EXRI)
+      PRINT *,"DEBUG: TFSAS =", TFSAS
       SQK=VK*VK*EXDC
       CALL APS(ZPV,LM,NSPH,IOG,RMI,REI,SQK,GAPS)
       CALL RABAS(INPOL,LM,NSPH,IOG,RMI,REI,TQSE,TQSPE,TQSS,TQSPS)
@@ -664,6 +665,11 @@ CCC   DIMENSION ZPV(LM,3,2,2)
       ZP=-1.0D0/DSQRT(XD)
       ZPV(L,2,1,2)=ZP
       ZPV(L,2,2,1)=ZP
+      IF(ZP.EQ.C0) GOTO 15
+C      PRINT 9010,L-1,ZP
+C      PRINT 9011,L-1,ZP
+ 9010 FORMAT("DEBUG: zpv[",I1,"][1][0][1] = ",E15.3)
+ 9011 FORMAT("DEBUG: zpv[",I1,"][1][1][0] = ",E15.3)
    15 CONTINUE
       IF(LM.EQ.1)GO TO 30
       DO 20 L=2,LM
@@ -672,6 +678,11 @@ CCC   DIMENSION ZPV(LM,3,2,2)
       ZP=DSQRT(XN/XD)
       ZPV(L,1,1,1)=ZP
       ZPV(L,1,2,2)=ZP
+      IF(ZP.EQ.C0) GOTO 20
+C      PRINT 9012,L-1,ZP
+C      PRINT 9013,L-1,ZP
+ 9012 FORMAT("DEBUG: zpv[",I1,"][0][0][0] = ",E15.3)
+ 9013 FORMAT("DEBUG: zpv[",I1,"][0][1][1] = ",E15.3)
    20 CONTINUE
       LMMO=LM-1
       DO 25 L=1,LMMO
@@ -680,6 +691,11 @@ CCC   DIMENSION ZPV(LM,3,2,2)
       ZP=-DSQRT(XN/XD)
       ZPV(L,3,1,1)=ZP
       ZPV(L,3,2,2)=ZP
+      IF(ZP.EQ.C0) GOTO 25
+C      PRINT 9014,L-1,ZP
+C      PRINT 9015,L-1,ZP
+ 9014 FORMAT("DEBUG: zpv[",I1,"][2][0][0] = ",E15.3)
+ 9015 FORMAT("DEBUG: zpv[",I1,"][2][1][1] = ",E15.3)
    25 CONTINUE
    30 CONTINUE
       RETURN
@@ -1026,7 +1042,7 @@ CCC  2RMF(LI),DRMF(LI),REF(LI),DREF(LI)
       RETURN
    42 DO 43 J=1,LIPT
       FB(J)=RFJ(J)
-   43 FN(J)=RFN(J)
+ 43   FN(J)=RFN(J)
       IF(NSH.GT.1)GO TO 65
       CRI=DC0(1)/EXDC
       DO 60 L=1,LI
@@ -1041,7 +1057,9 @@ CCC  2RMF(LI),DRMF(LI),REF(LI),DREF(LI)
       CCNC=FBI(LPO)*DFB
       CCND=FB(LPO)*DFBI
       RMI(L,I)=1.0D0+UIM*(CCNA-CCNB)/(CCNC-CCND)
-   60 REI(L,I)=1.0D0+UIM*(CRI*CCNA-CCNB)/(CRI*CCNC-CCND)
+C     PRINT *,"DEBUG: gone 60, RMI(",L,",",I,") =",RMI(L,I)
+ 60   REI(L,I)=1.0D0+UIM*(CRI*CCNA-CCNB)/(CRI*CCNC-CCND)
+C 60  PRINT *,"DEBUG: gone 60, REI(",L,",",I,") =",REI(L,I)
       RETURN
    65 DO 80 L=1,LI
       LPO=L+1
@@ -1076,6 +1094,7 @@ cccccccccccccccccccccc
    80 DREF(L)=DY2*SZ+Y2
       CRI=(1.0D0,0.0D0)
       IF(MOD(NSH,2).NE.0)CRI=DC0(IC)/EXDC
+C     PRINT *,"DEBUG: going 90, AREX =",AREX
       DO 90 L=1,LI
       LPO=L+1
       LTPO=LPO+L
@@ -1087,11 +1106,13 @@ cccccccccccccccccccccc
       CCNC=RMF(L)*DFB
       CCND=DRMF(L)*FB(LPO)*SZ*LTPO
       RMI(L,I)=1.0D0+UIM*(CCNA-CCNB)/(CCNC-CCND)
+C      PRINT *,"DEBUG: gone 90, RMI(",L,",",I,") =",RMI(L,I)
       CCNA=REF(L)*DFN
       CCNB=DREF(L)*FN(LPO)*SZ*LTPO
       CCNC=REF(L)*DFB
       CCND=DREF(L)*FB(LPO)*SZ*LTPO
-   90 REI(L,I)=1.0D0+UIM*(CRI*CCNA-CCNB)/(CRI*CCNC-CCND)
+ 90   REI(L,I)=1.0D0+UIM*(CRI*CCNA-CCNB)/(CRI*CCNC-CCND)
+C 90  PRINT *,"DEBUG: gone 90, REI(",L,",",I,") =",REI(L,I)
       RETURN
       END
       SUBROUTINE RKT(NPNTMO,STEP,X,LPO,Y1,Y2,DY1,DY2)
-- 
GitLab