!
!-----------------------------------------------------------------------
!
!Purpose:
! Calculation of the collision probabilities for the multicell
! surfacic approximation.
!
!Copyright:
! Copyright (C) 2025 Ecole Polytechnique de Montreal
! This library is free software; you can redistribute it and/or
! modify it under the terms of the GNU Lesser General Public
! License as published by the Free Software Foundation; either
! version 2.1 of the License, or (at your option) any later version
!
!Author(s): A. Hebert
!
!Parameters: input
! IPTRK   pointer to the tracking (L_TRACK signature).
! IFTRAK  tracking file unit.
! IMPX    print flag (equal to zero for no print).
! NREGIO  total number of merged blocks for which specific values
!         of the neutron flux and reactions rates are required.
! NBMIX   number of mixtures (NBMIX=max(MAT(i))).
! MAT     index-number of the mixture type assigned to each volume.
! VOL     volumes.
! SIGT0   total macroscopic cross sections ordered by mixture.
! SIGW0   P0 within-group scattering macroscopic cross sections
!         ordered by mixture.
! NELPIJ  number of elements in pij matrix.
! ILK     leakage flag (=.true. if neutron leakage through external
!         boundary is present).
! NBATCH  number of tracks dispached in eack OpenMP core.
! TITREC  title.
! NALBP   number of multigroup physical albedos.
! ALBP    multigroup physical albedos.
!
!Parameters: output
! PIJ     reduced and symmetrized collision probabilities.
!
!-----------------------------------------------------------------------
!
SUBROUTINE MUSP(IPTRK,IFTRAK,IMPX,NREGIO,NBMIX,MAT,VOL,SIGT0,SIGW0,NELPIJ, &
  & ILK,NBATCH,TITREC,NALBP,ALBP,PIJ)
  USE GANLIB
  !----
  !  SUBROUTINE ARGUMENTS
  !----
  LOGICAL ILK
  TYPE(C_PTR) IPTRK
  INTEGER IFTRAK,IMPX,NREGIO,NBMIX,MAT(NREGIO),NELPIJ,NBATCH,NALBP
  REAL VOL(NREGIO),SIGT0(0:NBMIX),SIGW0(0:NBMIX),PIJ(NELPIJ),ALBP(NALBP)
  CHARACTER TITREC*72
  !----
  !  LOCAL VARIABLES
  !----
  PARAMETER (EPS1=1.0E-4,NSTATE=40)
  TYPE(C_PTR) JPTRK,KPTRK
  INTEGER ISTATT(NSTATE),NNPSYS(1)
  CHARACTER TITRE2*72
  logical LSKIP
  !----
  !  ALLOCATABLE ARRAYS
  !----
  INTEGER, ALLOCATABLE, DIMENSION(:) :: MATALB,NMC_NODE,NMC_SURF,MAT2,IGEN,INUM, &
    & IFR,MIX,IMAC
  REAL, ALLOCATABLE, DIMENSION(:) :: SIGT2,SIGW2,PIJW,PISW,PSJW,PSSW,WORK,ALB,DVX
  REAL, ALLOCATABLE, DIMENSION(:,:) :: VOLSUR,PP,PSSB
  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: DPROB,DPROBX
  !
  IND(I,J) = MAX(I+J3+1,J+J3+1)*(MAX(I+J3+1,J+J3+1)-1)/2 &
  & + MIN(I+J3+1,J+J3+1)
  !
  WPR(I,J)= REAL(DPROB( IND(I,J),1 ) / DPROB( IND(I,0),1 ))
  !----
  !  BICKLEY FLAG
  !----
  SAVE IBICKL
  DATA IBICKL/0/
  !----
  !  RECOVER BICKLEY TABLES
  !----
  IF(IBICKL.EQ.0) THEN
    CALL XDRTA2
    IBICKL=1
  ENDIF
  !----
  !  RECOVER SALT SPECIFIC PARAMETERS
  !----
  REWIND IFTRAK
  CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATT)
  IF(NREGIO.NE.ISTATT(1)) THEN
     CALL XABORT('MUSP: STATE VECTOR HAS INVALID # OF ZONES.')
  ENDIF
  NMACRO=ISTATT(24) ! NGEN
  NMCEL=NMACRO
  NMERGE=NMACRO
  NGEN=NMACRO
  ALLOCATE(IGEN(NMERGE),INUM(NMCEL))
  DO IK=1,NMERGE
    IGEN(IK)=IK
  ENDDO
  DO IK=1,NMCEL
    INUM(IK)=IK
  ENDDO
  IF(NMACRO.EQ.0) CALL XABORT('MUSP: MUST MODULE TRACKING IS MANDATORY.')
  ALLOCATE(NMC_NODE(NMACRO+1),NMC_SURF(NMACRO+1))
  JPTRK=LCMGID(IPTRK,'MACRO-TRACK')
  CALL LCMGET(IPTRK,'NMC_NODE',NMC_NODE)
  CALL LCMGET(IPTRK,'NMC_SURF',NMC_SURF)
  NMIX=NMC_SURF(NMACRO+1)
  NIFR=NMC_SURF(NMACRO+1)
  !----
  ! LOOP OVER MACRO GEOMETRIES AND COMPUTE PIJ MATRICES USING EXCELP
  !----
  J1=0
  NMIX=0
  NPIJ=0
  NPIS=0
  NPSS=0
  DO IMACRO=1,NMACRO
    J2=NMC_NODE(IMACRO+1)-NMC_NODE(IMACRO)
    J3=NMC_SURF(IMACRO+1)-NMC_SURF(IMACRO)
    J1=J1+J2
    NMIX=NMIX+J3
    NPIJ=NPIJ+J2*J2
    NPIS=NPIS+J2*J3
    NPSS=NPSS+J3*J3
  ENDDO
  IF(J1.NE.NREGIO) CALL XABORT('MUSP: INVALID NREGIO.')
  IF(NMIX.NE.NMC_SURF(NMACRO+1)) CALL XABORT('MUSP: INVALID NMIX.')
  ALLOCATE(PIJW(NPIJ),PISW(NPIS),PSJW(NPIS),PSSW(NPSS))
  J1=0
  IPIJ=0
  IPIS=0
  IPSS=0
  DO IMACRO=1,NMACRO
    J2=NMC_NODE(IMACRO+1)-NMC_NODE(IMACRO)
    J3=NMC_SURF(IMACRO+1)-NMC_SURF(IMACRO)
    N2PRO=(J2+J3+1)**2
    WRITE(TITRE2,'(A,9H -- MACRO,I5.5)') TRIM(TITREC),IMACRO
    KPTRK=LCMGIL(JPTRK,IMACRO)
    KNORM=4 ! use HELIOS-type normalization
    NNPSYS(1)=1
    ALLOCATE(MAT2(J2),SIGT2(J2),SIGW2(J2))
    ALLOCATE(MATALB(-J3:J2),VOLSUR(-J3:J2,1),DPROB(N2PRO,1),DPROBX(N2PRO,1))
    CALL LCMGET(KPTRK,'MATCOD',MAT2)
    CALL EXCELP(KPTRK,IFTRAK,IMPX,J3,J2,NBMIX,MAT2,KNORM,SIGT0,1,N2PRO, &
    & 1,NNPSYS(1),NBATCH,TITRE2,NALBP,ALBP,MATALB,VOLSUR,DPROB,DPROBX)
    !----
    !  CHECK IF SCATTERING REDUCTION IS REQUIRED 
    !----
    DO I=1,J2
      SIGT2(I)=SIGT0(MAT2(I)) ! sigt by node
      SIGW2(I)=SIGW0(MAT2(I)) ! sigw by node
    ENDDO
    LSKIP=.TRUE.
    DO I=1,J2
      LSKIP=LSKIP.AND.(SIGW2(I).EQ.0.0)
    ENDDO
    !----
    !  SCATTERING REDUCTION IF LSKIP=.FALSE.
    !----
    IF(LSKIP) THEN
      ! DO NOT PERFORM SCATTERING REDUCTION.
      DO I=1,J2
        DO J=1,J2
          IF(SIGT2(J).EQ.0.0) THEN
            PIJW(IPIJ+(J-1)*J2+I)=WPR(I,J)
          ELSE
            PIJW(IPIJ+(J-1)*J2+I)=WPR(I,J)/SIGT2(J)
          ENDIF
        ENDDO
      ENDDO
      DO I=1,J2
        DO JC=1,J3
          PISW(IPIS+(JC-1)*J2+I)=WPR(I,-JC)
          IF(SIGT2(I).EQ.0.0) THEN
            PSJW(IPIS+(I-1)*J3+JC)=WPR(-JC,I)
          ELSE
            PSJW(IPIS+(I-1)*J3+JC)=WPR(-JC,I)/SIGT2(I)
          ENDIF
        ENDDO
      ENDDO
      DO IC=1,J3
        DO JC=1,J3
          PSSW(IPSS+(JC-1)*J3+IC)=WPR(-IC,-JC)
        ENDDO
      ENDDO
    ELSE
      ! COMPUTE THE SCATTERING-REDUCED COLLISION AND ESCAPE MATRICES.
      DO I=1,J2
        DO J=1,J2
          IF(SIGT2(J).EQ.0.0) THEN
            PIJW(IPIJ+(J-1)*J2+I)=0.0
          ELSE
            PIJW(IPIJ+(J-1)*J2+I)=-WPR(I,J)*SIGW2(J)/SIGT2(J)
          ENDIF
        ENDDO
        PIJW(IPIJ+(I-1)*J2+I)=1.0+PIJW(IPIJ+(I-1)*J2+I)
      ENDDO
      CALL ALINV(J2,PIJW(IPIJ+1),J2,IER)
      IF(IER.NE.0) CALL XABORT('MUSP: SINGULAR MATRIX.')
      ALLOCATE(WORK(J2))
      DO I=1,J2
        DO K=1,J2
          WORK(K)=PIJW(IPIJ+(K-1)*J2+I)
        ENDDO
        DO J=1,J2
          WGAR=0.0
          DO K=1,J2
            IF(SIGT2(J).EQ.0.0) THEN
              WGAR=WGAR+WORK(K)*WPR(K,J)
            ELSE
              WGAR=WGAR+WORK(K)*WPR(K,J)/SIGT2(J)
            ENDIF
          ENDDO
          PIJW(IPIJ+(J-1)*J2+I)=WGAR
        ENDDO
        DO JC=1,J3
          WGAR=0.0
          DO K=1,J2
            WGAR=WGAR+WORK(K)*WPR(K,-JC)
          ENDDO
          PISW(IPIS+(JC-1)*J2+I)=WGAR
        ENDDO
      ENDDO
      DEALLOCATE(WORK)
      !
      ! COMPUTE THE SCATTERING-REDUCED COLLISION PROBABILITY MATRIX
      ! FOR INCOMING NEUTRONS.
      DO IC=1,J3
        DO J=1,J2
          IF(SIGT2(J).EQ.0.0) THEN
            WGAR=WPR(-IC,J)
          ELSE
            WGAR=WPR(-IC,J)/SIGT2(J)
          ENDIF
          DO K=1,J2
            IF(SIGT2(K).NE.0.0) THEN
              WGAR=WGAR+WPR(-IC,K)*PIJW(IPIJ+(J-1)*J2+K)*SIGW2(K)/SIGT2(K)
            ENDIF
          ENDDO
          PSJW(IPIS+(J-1)*J3+IC)=WGAR
        ENDDO
      ENDDO
      !
      ! COMPUTE THE SCATTERING-REDUCED TRANSMISSION PROBABILITY MATRIX.
      DO IC=1,J3
        DO JC=1,J3
          WGAR=WPR(-IC,-JC)
          DO K=1,J2
            IF(SIGT2(K).NE.0.0) THEN
              WGAR=WGAR+WPR(-IC,K)*PISW(IPIS+(JC-1)*J2+K)*SIGW2(K)/SIGT2(K)
            ENDIF
          ENDDO
          PSSW(IPSS+(JC-1)*J3+IC)=WGAR
        ENDDO
      ENDDO
    ENDIF
    DEALLOCATE(DPROBX,DPROB,VOLSUR,MATALB)
    IF(IMPX.GE.8) THEN
      IF(LSKIP) THEN
        IN=1
      ELSE
        IN=2
      ENDIF
      CALL SYBPRX(IN,J3,J2,IMACRO,SIGT2,SIGW2,PIJW(IPIJ+1),PISW(IPIS+1), &
      & PSJW(IPIS+1),PSSW(IPSS+1))
    ENDIF
    DEALLOCATE(SIGW2,SIGT2,MAT2)
    J1=J1+J2
    IPIJ=IPIJ+J2*J2
    IPIS=IPIS+J2*J3
    IPSS=IPSS+J3*J3
  ENDDO
  ! end of SYB004 equivalent
  !----
  ! COMPUTE THE GLOBAL SCATTERING-REDUCED COLLISION PROBABILITY MATRIX
  !----
  ALLOCATE(IMAC(NREGIO),PP(NREGIO,NREGIO))
  CALL LCMGET(IPTRK,'MERGE_MACRO',IMAC)
  PP(:NREGIO,:NREGIO)=0.0
  IPIJ=0
  DO JKG=1,NGEN
    J2=NMC_NODE(JKG+1)-NMC_NODE(JKG)
    I1=0
    DO IKK=1,NMERGE
      IKG=IGEN(IKK)
      I2=NMC_NODE(IKG+1)-NMC_NODE(IKG)
      IF(IKG.EQ.JKG) THEN
        DO J=1,J2
          DO I=1,J2
            PP(IMAC(I1+I),IMAC(I1+J))=PIJW(IPIJ+(J-1)*J2+I)
          ENDDO
        ENDDO
      ENDIF
      I1=I1+I2
    ENDDO
    IPIJ=IPIJ+J2*J2
  ENDDO
  !----
  ! COMPUTE PSSB=A*(I-PSS*A)**-1
  !----
  ALLOCATE(IFR(NIFR),ALB(NIFR),MIX(NMIX),DVX(NMIX))
  CALL LCMGET(IPTRK,'IFR',IFR)
  CALL LCMGET(IPTRK,'ALB',ALB)
  CALL LCMGET(IPTRK,'MIX',MIX)
  CALL LCMGET(IPTRK,'DVX',DVX)
  IJAT=MAXVAL(MIX)
  ALLOCATE(PSSB(IJAT,2*IJAT))
  PSSB(:IJAT,:2*IJAT)=0.0
  DO I=1,IJAT
    PSSB(I,I)=1.0
  ENDDO
  DO ICEL=1,NMCEL
    IKK=INUM(ICEL)
    IKG=IGEN(IKK)
    J3=NMC_SURF(IKG+1)-NMC_SURF(IKG)
    IT=0
    DO IK=1,IKK-1
      IT=IT+(NMC_SURF(IGEN(IK)+1)-NMC_SURF(IGEN(IK)))
    ENDDO
    IS=0
    DO IK=1,ICEL-1
      IS=IS+(NMC_SURF(IGEN(INUM(IK))+1)-NMC_SURF(IGEN(INUM(IK))))
    ENDDO
    IPSS=0
    DO IK=1,IKG-1
      IPSS=IPSS+(NMC_SURF(IK+1)-NMC_SURF(IK))**2
    ENDDO
    DO JC=1,J3
      J1=IFR(IS+JC)
      J2=MIX(IT+JC)
      ALBEDO=ALB(IS+JC)
      PSSB(J1,IJAT+J2)=PSSB(J1,IJAT+J2)+ALBEDO*DVX(IT+JC)
      DO IC=1,J3
        J2=MIX(IT+IC)
        PSSB(J1,J2)=PSSB(J1,J2)-PSSW(IPSS+(JC-1)*J3+IC)*ALBEDO*DVX(IT+IC)
      ENDDO
    ENDDO
  ENDDO
  CALL ALSB(IJAT,IJAT,PSSB,IER,IJAT)
  IF(IER.NE.0) CALL XABORT('MUSP: SINGULAR MATRIX.')
  !----
  !  COMPUTATION OF PISW*PSSB*PSJW
  !----
  I1=0
  DO IKK=1,NMERGE
    IKG=IGEN(IKK)
    I2=NMC_NODE(IKG+1)-NMC_NODE(IKG)
    I3=NMC_SURF(IKG+1)-NMC_SURF(IKG)
    IT=0
    DO IK=1,IKK-1
      IT=IT+(NMC_SURF(IGEN(IK)+1)-NMC_SURF(IGEN(IK)))
    ENDDO
    IPIS=0
    DO IK=1,IKG-1
      IPIS=IPIS+(NMC_NODE(IK+1)-NMC_NODE(IK))*(NMC_SURF(IK+1)-NMC_SURF(IK))
    ENDDO
    DO I=1,I2
      DO IC=1,I3
        ICC=MIX(IT+IC)
        ZZZ=PISW(IPIS+(IC-1)*I2+I)*SIGN(1.0,DVX(IT+IC))
        J1=0
        DO JKK=1,NMERGE
          JKG=IGEN(JKK)
          J2=NMC_NODE(JKG+1)-NMC_NODE(JKG)
          J3=NMC_SURF(JKG+1)-NMC_SURF(JKG)
          JT=0
          DO IK=1,JKK-1
            JT=JT+(NMC_SURF(IGEN(IK)+1)-NMC_SURF(IGEN(IK)))
          ENDDO
          IPSJ=0
          DO IK=1,JKG-1
            IPSJ=IPSJ+(NMC_NODE(IK+1)-NMC_NODE(IK))*(NMC_SURF(IK+1)-NMC_SURF(IK))
          ENDDO
          DO J=1,J2
            DO JC=1,J3
              JCC=MIX(JT+JC)
              PBJ=PSJW(IPSJ+(J-1)*J3+JC)
              PP(IMAC(I1+I),IMAC(J1+J))=PP(IMAC(I1+I),IMAC(J1+J))+ZZZ*DVX(JT+JC)* &
              & PSSB(JCC,IJAT+ICC)*PBJ
            ENDDO
          ENDDO
          J1=J1+J2
        ENDDO
      ENDDO
    ENDDO
    I1=I1+I2
  ENDDO
  ! end of SYBRX3 equivalent
  DEALLOCATE(PSSB,DVX,MIX,ALB,IFR)
  DEALLOCATE(PSSW,PSJW,PISW,PIJW)
  DEALLOCATE(NMC_SURF,NMC_NODE)
  DEALLOCATE(INUM,IGEN)
  !
  IF(IMPX.GE.7) THEN
    WRITE (6,170) (J,J=1,NREGIO)
    DO I=1,NREGIO
      WRITE (6,180) I,(PP(I,J),J=1,NREGIO)
    ENDDO
    WRITE (6,'(//)')
  ENDIF
  IF((IMPX.GE.10).OR.(IMPX.LT.0)) THEN
    ! CHECK THE RECIPROCITY CONDITIONS.
    VOLTOT=0.0
    DO I=1,NREGIO
      VOLTOT=VOLTOT+VOL(I)
    ENDDO
    VOLTOT=VOLTOT/REAL(NREGIO)
    WRK=0.0
    DO I=1,NREGIO
      DO J=1,NREGIO
        AAA=PP(I,J)*VOL(I)
        BBB=PP(J,I)*VOL(J)
        WRK=MAX(WRK,ABS(AAA-BBB)/VOLTOT)
      ENDDO
    ENDDO
    IF(WRK.GE.EPS1) WRITE (6,150) WRK
    IF(WRK.GE.EPS1) CALL XABORT('MUSP: non symmetric matrices.')
    ! CHECK THE CONSERVATION CONDITIONS.
    IF(.NOT.ILK) THEN
      WRK=0.0
      DO I=1,NREGIO
        F1=1.0
        DO J=1,NREGIO
          AAA=PP(I,J)
          F1=F1-AAA*(SIGT0(MAT(J))-SIGW0(MAT(J)))
        ENDDO
        WRK=AMAX1(WRK,ABS(F1))
      ENDDO
      IF(WRK.GE.EPS1) WRITE (6,160) WRK
      IF(WRK.GE.EPS1) CALL XABORT('MUSP: non conservative matrices.')
    ENDIF
  ENDIF
  !
  IC=0
  DO IKK=1,NREGIO
    IOF=(IKK-1)*NREGIO
    DO JKK=1,IKK
      IC=IC+1
      PIJ(IC)=PP(JKK,IKK)*VOL(JKK)
    ENDDO
  ENDDO
  DEALLOCATE(IMAC,PP)
  RETURN
  !
  150 FORMAT (/56H MUSP: THE SCATTERING-REDUCED PIJ DO NOT MEET THE RECIPR, &
  & 25HOCITY CONDITIONS. RECIP =,1P,E10.3/)
  160 FORMAT (/56H MUSP: THE SCATTERING-REDUCED PIJ DO NOT MEET THE CONSER, &
  & 25HVATION CONDITIONS. LEAK =,1P,E10.3/)
  170 FORMAT (//47H MUSP: SCATTERING-REDUCED COLLISION PROBABILITY, &
  & 9H MATRIX ://(11X,2HJ=,I4,:,5X,2HJ=,I4,:,5X,2HJ=,I4,:,5X,2HJ=, &
  & I4,:,5X,2HJ=,I4,:,5X,2HJ=,I4,:,5X,2HJ=,I4,:,5X,2HJ=,I4,:,5X,2HJ=, &
  & I4,:,5X,2HJ=,I4,:,5X,2HJ=,I4))
  180 FORMAT (3H I=,I4,2H: ,1P,11E11.3/(9X,11E11.3))
END SUBROUTINE MUSP