1 files changed, 109 insertions, 0 deletions
diff --git a/Dragon/src/MCGPJJ.f b/Dragon/src/MCGPJJ.f
new file mode 100644
index 0000000..c837c56
--- /dev/null
+++ b/Dragon/src/MCGPJJ.f
@@ -0,0 +1,109 @@
+*DECK MCGPJJ
+      SUBROUTINE MCGPJJ(IPTRK,IPRINT,NDIM,NANI,MAXNU,NPJJM,KEYANI)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute the effective number of pjj intermodes to be stored in 2D or
+* 3D and the corresponding index when an expansion up to order L 
+* of the scattering cross-section is considered in order to
+* construct the source term of the scalar flux moments for
+* a method of characteristics iteration.
+*
+*Copyright:
+* Copyright (C) 2002 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): R. Le Tellier
+*
+*Parameters: input
+* IPTRK  pointer to the tracking (L_TRACK signature).
+* IPRINT print flag (> 1 for print).
+* NDIM   number of dimensions for the geometry.
+* NANI   scattering anisotropy (=1 for isotropic scattering).
+* MAXNU  number of angular modes nu=(l,m).
+* KEYANI 'mode to l' index: l=KEYANI(nu):
+*        Cartesian 2D KEYANI(NU)=INT(0.5*(SQRT(REAL(1+8*(NU-1)))-1.0));
+*        Cartesian 3D KEYANI(NU)=INT(SQRT(REAL(NU-1))).
+*
+*Parameters: output
+* NPJJM  number of non-vanishing pjj modes.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+      IMPLICIT NONE
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPTRK
+      INTEGER IPRINT,NDIM,NANI,MAXNU,NPJJM,KEYANI(MAXNU)
+*----
+*  LOCAL VARIABLES
+*----
+      INTEGER MAXPJJM,IL,INU,ILP,INUP,L,LT,DT,NMODE,NMODO,NPJJM0
+      LOGICAL EVEN
+      CHARACTER CDIM*2
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: IPJJM
+*--- 
+*  Compute the number of effective modes
+*---
+      MAXPJJM=MAXNU*(MAXNU+1)/2
+      L=NANI-1
+      LT=L/2
+      DT=L-2*LT
+      IF (NDIM.EQ.2) THEN
+         CDIM='2D'
+         NMODE=(LT+1)**2        ! number of 'even l' modes
+         NMODO=(LT+1)*(LT+2*DT) ! number of 'odd l' modes
+      ELSE ! NDIM.EQ.3
+         CDIM='3D'
+         NMODE=(L+1-DT)*(LT+1)
+         NMODO=(L+1+DT)*(LT+DT)
+      ENDIF
+*     'even l' and 'odd l' modes are unconnected
+*            (i.e. pjj('even l' <- 'odd l') = 0)
+*     and pjj(nu <- nu') = pjj(nu' <- nu)
+*     so the effective number of pjj is
+      NPJJM=NMODE*(NMODE+1)/2+NMODO*(NMODO+1)/2
+*
+      IF (IPRINT.GT.1) THEN
+         WRITE(*,*) '--------------------'
+         WRITE(*,*) 'ANISOTROPY INDEX FOR L=',NANI-1,' IN ',CDIM
+         CALL PRINIM('NU->L ',KEYANI(1),MAXNU)
+         WRITE(*,*) '--------------------'
+         WRITE(*,*) NPJJM,
+     1      ' PJJ(NU <- NU'') MODES OUT OF',MAXPJJM,' TO BE STORED'
+      ENDIF
+*---
+*  Compute and store the indexes for these modes
+*---
+      ALLOCATE(IPJJM(2*NPJJM))
+      NPJJM0=-1
+      DO INU=1,MAXNU
+         IL=KEYANI(INU)
+         DO INUP=1,INU
+            ILP=KEYANI(INUP)
+            EVEN=(MOD((IL+ILP),2).EQ.0)
+            IF (EVEN) THEN
+               NPJJM0=NPJJM0+1
+               IPJJM(NPJJM0+1)=INU
+               IPJJM(NPJJM+NPJJM0+1)=INUP
+            ENDIF
+         ENDDO
+      ENDDO
+      NPJJM0=NPJJM0+1
+      IF(NPJJM0.NE.NPJJM) CALL XABORT('MCGPIJ: bug.')
+      CALL LCMPUT(IPTRK,'PJJIND$MCCG',2*NPJJM,1,IPJJM)
+      IF (IPRINT.GT.1) THEN
+         WRITE(*,*) 'INDEXES FOR THE CORRESPONDING',NPJJM,' MODES:'
+         CALL PRINIM('-> NU ',IPJJM(1),NPJJM)
+         CALL PRINIM('NU''-> ',IPJJM(NPJJM+1),NPJJM)
+      ENDIF
+      DEALLOCATE(IPJJM)
+*
+      RETURN
+      END