Initial commit from Polytechnique Montreal

1 files changed, 706 insertions, 0 deletions

diff --git a/Dragon/src/PSOUR.f b/Dragon/src/PSOUR.f
new file mode 100644
index 0000000..688c5dd
--- /dev/null
+++ b/Dragon/src/PSOUR.f
@@ -0,0 +1,706 @@
+*DECK PSOUR
+      SUBROUTINE PSOUR(NENTRY,HENTRY,IENTRY,JENTRY,KENTRY)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Set the transition source from companion particles in a coupled-
+* particle problem or the isotropic/monodirectional boundary sources.
+*
+*Copyright:
+* Copyright (C) 2020 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 and C. Bienvenue
+*
+*Parameters: input/output
+* NENTRY  number of LCM objects or files used by the operator.
+* HENTRY  name of each LCM object or file:
+*         HENTRY(1): creation/modification type(L_SOURCE);
+*         HENTRY(2): read-only type(L_MACROLIB) for the source;
+*         HENTRY(3): read-only type(L_TRACKING);
+*         HENTRY(4): read-only type(L_FLUX) for companion particle 1;
+*         HENTRY(5): read-only type(L_FLUX) for companion particle 1;
+*         HENTRY(6): read-only type(L_FLUX) for companion particle 2;
+*         ...
+* IENTRY  type of each LCM object or file:
+*         =1 LCM memory object; =2 XSM file; =3 sequential binary file;
+*         =4 sequential ascii file.
+* JENTRY  access of each LCM object or file:
+*         =0 the LCM object or file is created;
+*         =1 the LCM object or file is open for modifications;
+*         =2 the LCM object or file is open in read-only mode.
+* KENTRY  LCM object address or file unit number.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER      NENTRY,IENTRY(NENTRY),JENTRY(NENTRY)
+      TYPE(C_PTR)  KENTRY(NENTRY)
+      CHARACTER    HENTRY(NENTRY)*12
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER (IOUT=6,NSTATE=40,MAXPAR=10)
+      TYPE(C_PTR) IPSOUR,JPSOUR,KPSOUR,IPTRK,IPGEOM,IPMAC,JPMAC,KPMAC,
+     1 IPFLX(MAXPAR),IPSOUR2(MAXPAR),JPFLX,JPTRK
+      TYPE(C_PTR) JPSOURA,JPSOURB,KPSOURA,KPSOURB
+      CHARACTER HSIGN*12,TEXT12*12,CMODUL*12,HPART*6,HSMG*131
+      CHARACTER(LEN=12) HPFLX(MAXPAR)
+      CHARACTER(LEN=1) HPARTC(MAXPAR)
+      INTEGER ISTATE(NSTATE),MESH_LEN(3),GN,STYPE,CPT
+      DOUBLE PRECISION DFLOTT
+      REAL DIR(3),ZNORM,BSNORM
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      REAL, ALLOCATABLE, DIMENSION(:) :: NORM,SOURTEMP
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: MATCOD
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: FUNKNO,SUNKNO,BS
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: ISISOUR,NBS1,NBS2,NBS
+      REAL, ALLOCATABLE, DIMENSION(:) :: XXX,YYY,ZZZ,ISOUR,XMIN,XMAX,
+     1 YMIN,YMAX,ZMIN,ZMAX
+      INTEGER, ALLOCATABLE, DIMENSION(:,:) :: BSINFO
+      INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: BSINFO1,BSINFO2
+      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: BS1,BS2
+*----
+*  PARAMETER VALIDATION.
+*----
+      IF((IENTRY(1).NE.1).AND.(IENTRY(1).NE.2)) CALL XABORT('PSOUR: LI'
+     1 //'NKED LIST OR XSM FILE EXPECTED AT LHS.')
+      IF(JENTRY(1).NE.0.AND.JENTRY(1).NE.1) CALL XABORT('PSOUR: ENTRY '
+     1 //'IN CREATE OR MODIFICATION MODE EXPECTED.')
+      IPSOUR=KENTRY(1)
+      IPMAC=C_NULL_PTR
+      IPTRK=C_NULL_PTR
+      JPTRK=C_NULL_PTR
+      IPGEOM=C_NULL_PTR
+      IPFLX(:)=C_NULL_PTR
+      NPART0=0
+      NSOUR2=0
+      NTRK=0
+      DO I=2,NENTRY
+        IF((IENTRY(I).NE.1).AND.(IENTRY(I).NE.2)) CALL XABORT('PSOUR: '
+     1  //'LINKED LIST OR XSM FILE EXPECTED AT RHS.')
+        IF(JENTRY(I).NE.2) CALL XABORT('PSOUR: ENTRY IN READ-ONLY MODE'
+     1  //' EXPECTED.')
+        CALL LCMGTC(KENTRY(I),'SIGNATURE',12,HSIGN)
+        IF(HSIGN.EQ.'L_MACROLIB') THEN
+          IPMAC=KENTRY(I)
+        ELSE IF(HSIGN.EQ.'L_LIBRARY') THEN
+          IPMAC=LCMGID(KENTRY(I),'MACROLIB')
+        ELSE IF(HSIGN.EQ.'L_TRACK') THEN
+          IF(NTRK.EQ.0) THEN
+            IPTRK=KENTRY(I)
+          ELSE IF(NTRK.EQ.1) THEN
+            JPTRK=KENTRY(I)
+          ELSE
+            CALL XABORT('PSOUR: TOO MUCH TRACKING OBJECTS DEFINED.')
+          ENDIF
+          NTRK=NTRK+1
+        ELSE IF(HSIGN.EQ.'L_GEOM') THEN
+          IPGEOM=KENTRY(I)
+        ELSE IF(HSIGN.EQ.'L_FLUX') THEN
+          NPART0=NPART0+1
+          IF(NPART0.GT.MAXPAR) CALL XABORT('PSOUR: MAXPAR OVERFLOW.')
+          IPFLX(NPART0)=KENTRY(I)
+          HPFLX(NPART0)=HENTRY(I)
+        ELSE IF(HSIGN.EQ.'L_SOURCE') THEN
+          NSOUR2=NSOUR2+1
+          IF(NSOUR2.GT.MAXPAR) CALL XABORT('PSOUR: MAXPAR OVERFLOW.')
+          IPSOUR2(NSOUR2)=KENTRY(I)
+        ELSE
+          TEXT12=HENTRY(I)
+          CALL XABORT('PSOUR: SIGNATURE OF '//TEXT12//' IS '//HSIGN//
+     1    '. NON EXPECTED TYPE.')
+        ENDIF
+      ENDDO
+*----
+*  RECOVER MACROLIB INFORMATION FOR THE SOURCE
+*----
+      IF(.NOT.C_ASSOCIATED(IPMAC)) CALL XABORT('PSOUR: NO MACROLIB.')
+      CALL LCMGET(IPMAC,'STATE-VECTOR',ISTATE)
+      NG2=ISTATE(1)
+      NMAT=ISTATE(2)
+      NANIS=ISTATE(3)-1
+      IADJ=ISTATE(13)
+      CALL LCMGTC(IPMAC,'PARTICLE',6,HPART)
+*----
+*  RECOVER TRACKING INFORMATION
+*----
+      IF(.NOT.C_ASSOCIATED(IPTRK)) CALL XABORT('PSOUR: NO TRACKING.')
+      CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATE)
+      NREG=ISTATE(1)
+      NUNS=ISTATE(2)
+      NDIM=ISTATE(9)
+      LX=ISTATE(12)
+      LY=ISTATE(13)
+      LZ=ISTATE(14)
+      NLF=ISTATE(15)
+      ALLOCATE(MATCOD(NREG))
+      CALL LCMGET(IPTRK,'MATCOD',MATCOD)
+      CALL LCMGTC(IPTRK,'TRACK-TYPE',12,CMODUL)
+*----
+*  RECOVER GEOMETRY INFORMATION
+*----
+      IF(C_ASSOCIATED(IPGEOM)) THEN
+        IF(NDIM.EQ.1) THEN
+          CALL LCMLEN(IPGEOM,'MESHX',MESH_LEN(1),ITYLCM)
+          ALLOCATE(XXX(MESH_LEN(1)),YYY(0),ZZZ(0))
+          CALL LCMGET(IPGEOM,'MESHX',XXX)
+        ELSE IF(NDIM.EQ.2) THEN
+          CALL LCMLEN(IPGEOM,'MESHX',MESH_LEN(1),ITYLCM)
+          CALL LCMLEN(IPGEOM,'MESHY',MESH_LEN(2),ITYLCM)
+          ALLOCATE(XXX(MESH_LEN(1)),YYY(MESH_LEN(2)),ZZZ(0))
+          CALL LCMGET(IPGEOM,'MESHX',XXX)
+          CALL LCMGET(IPGEOM,'MESHY',YYY)
+        ELSE IF(NDIM.EQ.3) THEN
+          CALL LCMLEN(IPGEOM,'MESHX',MESH_LEN(1),ITYLCM)
+          CALL LCMLEN(IPGEOM,'MESHY',MESH_LEN(2),ITYLCM)
+          CALL LCMLEN(IPGEOM,'MESHZ',MESH_LEN(3),ITYLCM)
+          ALLOCATE(XXX(MESH_LEN(1)),YYY(MESH_LEN(2)),ZZZ(MESH_LEN(3)))
+          CALL LCMGET(IPGEOM,'MESHX',XXX)
+          CALL LCMGET(IPGEOM,'MESHY',YYY)
+          CALL LCMGET(IPGEOM,'MESHZ',ZZZ)
+        ENDIF
+      ENDIF
+*----
+*  RECOVER ADDITIONAL SOURCE INFORMATION
+*----
+      ALLOCATE(SUNKNO(NUNS,NG2),SOURTEMP(NUNS))
+      SUNKNO(:NUNS,:NG2)=0.0
+      DO NSOURIN=1,NSOUR2
+        ! VOLUMIC SOURCES
+        IF(IADJ.EQ.0) THEN
+          JPSOUR=LCMGID(IPSOUR2(NSOURIN),'DSOUR') 
+        ELSE IF(IADJ.EQ.1) THEN
+          JPSOUR=LCMGID(IPSOUR2(NSOURIN),'ASOUR') 
+        ENDIF
+        KPSOUR=LCMGIL(JPSOUR,1)
+        DO IG=1,NG2
+          CALL LCMGDL(KPSOUR,IG,SOURTEMP(:NUNS))
+          SUNKNO(:NUNS,IG)=SUNKNO(:NUNS,IG)+SOURTEMP(:NUNS)
+        ENDDO
+      ENDDO
+*----
+* INITIALIZATION FLUX AND RECOVER INPUT SOURCE INFORMATION IF AVAILABLE
+*----
+      ALLOCATE(NBS1(NG2),BS1(1,1,1),BSINFO1(1,1,1))
+      NBS1(:)=0
+      BSNORM=-1.0
+      ISBS=0
+      BS1(:,:,:)=0.0
+      BSINFO1(:,:,:)=0
+
+      IF(NDIM.EQ.1) THEN
+        MAXL=1
+      ELSE IF(NDIM.EQ.2) THEN
+        MAXL=MAX(LX,LY)
+      ELSE
+        MAXL=MAX(LX*LY,LX*LZ,LY*LZ)
+      ENDIF
+
+      IF(JENTRY(1).EQ.1) THEN
+        CALL LCMGTC(IPSOUR,'SIGNATURE',12,HSIGN)
+        IF(HSIGN.NE.'L_SOURCE') THEN
+          TEXT12=HENTRY(1)
+          CALL XABORT('SOUR: SIGNATURE OF'//TEXT12//' IS '//HSIGN//
+     1    '. L_SOURCE EXPECTED.')
+        ENDIF
+        ! VOLUMIC SOURCES
+        IF(IADJ.EQ.0) THEN
+          JPSOUR=LCMGID(IPSOUR,'DSOUR') 
+        ELSE IF(IADJ.EQ.1) THEN
+          JPSOUR=LCMGID(IPSOUR,'ASOUR') 
+        ENDIF
+        KPSOUR=LCMGIL(JPSOUR,1)
+        DO IG=1,NG2
+          CALL LCMGDL(KPSOUR,IG,SOURTEMP(:NUNS))
+          SUNKNO(:NUNS,IG)=SUNKNO(:NUNS,IG)+SOURTEMP(:NUNS)  
+        ENDDO
+        ! BOUNDARY SOURCES
+        CALL LCMLEN(IPSOUR,'NBS',ILEN,ITYLCM)
+        IF(ILEN.GT.0) THEN
+          ISBS=1
+          CALL LCMGET(IPSOUR,'NBS',NBS1)
+          DEALLOCATE(BS1,BSINFO1)
+          ALLOCATE(BS1(MAXL,NG2,SUM(NBS1)),BSINFO1(2,IG,SUM(NBS1)))
+          JPSOURA=LCMGID(IPSOUR,'BS')
+          JPSOURB=LCMGID(IPSOUR,'BSINFO')
+          DO IG=1,NG2
+            IF(NBS1(IG).GT.0) THEN
+              KPSOURA=LCMGIL(JPSOURA,IG)
+              KPSOURB=LCMGIL(JPSOURB,IG)
+              DO N=1,NBS1(IG)
+                CALL LCMGDL(KPSOURA,N,BS1(1,IG,N))
+                CALL LCMGDL(KPSOURB,N,BSINFO1(1,IG,N))
+              ENDDO
+            ENDIF
+          ENDDO
+        ENDIF
+        CALL LCMLEN(IPSOUR,'NORM-FS',ILEN,ITYLCM)
+        IF(ILEN.NE.0) CALL LCMGET(IPSOUR,'NORM-FS',BSNORM)
+      ENDIF
+      DEALLOCATE(SOURTEMP)
+*----
+*  READ THE INPUT DATA
+*----
+      IMPX=1
+      NPART=0
+      STYPE=-1
+      NSOUR=-1
+      ALLOCATE(ISOUR(NG2),ISISOUR(NG2))
+      ISISOUR=0    
+      ISOUR=0.0
+      ISIDEF=0
+      ISXDEF=0
+      ISYDEF=0
+      ISZDEF=0
+      MONOP=0
+      ISDIRDEF=0
+   10 CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+   11 IF(INDIC.NE.3) CALL XABORT('PSOUR: CHARACTER DATA EXPECTED.')
+      IF(TEXT12.EQ.'EDIT') THEN
+*       READ THE PRINT INDEX.
+        CALL REDGET(INDIC,IMPX,FLOTT,TEXT12,DFLOTT)
+        IF(INDIC.NE.1) CALL XABORT('PSOUR: INTEGER DATA EXPECTED.')
+      ELSE IF(TEXT12.EQ.'PARTICLE') THEN
+*       READ THE PARTICLE TYPE ('N', 'G', 'B', 'C', 'P')
+        NPART=NPART+1
+        IF(NPART.GT.NPART0) CALL XABORT('PSOUR: NPART0 OVERFLOW.')
+        CALL REDGET(INDIC,NITMA,FLOTT,HPARTC(NPART),DFLOTT)
+        IF(INDIC.NE.3) CALL XABORT('PSOUR: CHARACTER DATA EXPECTED.')
+        TEXT12='GROUP-'//HPARTC(NPART)
+        CALL LCMLEN(IPMAC,TEXT12,ILENG,ITYLCM)
+        IF(ILENG.EQ.0) THEN
+          WRITE(HSMG,'(14HPSOUR: RECORD ,A12,22H IS NOT AVAILABLE IN T,
+     1    12HHE MACROLIB.)') TEXT12
+          CALL LCMLIB(IPMAC)
+          CALL XABORT(HSMG)
+        ENDIF
+      ELSE IF(TEXT12.EQ.'ISO') THEN
+        ! ISOTROPIC SOURCE KEYWORD
+        IF(STYPE.NE.-1) CALL XABORT('PSOUR: ONLY ONE SOURCE TYPE ALLOW'
+     1  //'ED.')
+        STYPE=0
+        CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+        IF(INDIC.NE.1) CALL XABORT('PSOUR: INTEGER DATA EXPECTED'
+     1  //' (NUMBER OF SOURCES).')
+        NSOUR=NITMA
+        IF(NSOUR.LE.0) CALL XABORT('PSOUR: INVALID NUMBER OF SOURCES.')
+        ALLOCATE(XMIN(NSOUR),XMAX(NSOUR),YMIN(NSOUR),YMAX(NSOUR),
+     1  ZMIN(NSOUR),ZMAX(NSOUR))
+      ELSE IF(TEXT12.EQ.'MONO') THEN
+        ! MONODIRECTIONNAL BOUNDARY SOURCE KEYWORD
+        IF(STYPE.NE.-1) CALL XABORT('PSOUR: ONLY ONE SOURCE TYPE ALLOW'
+     1  //'ED.')
+        STYPE=1
+        CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+        IF(INDIC.NE.1) CALL XABORT('PSOUR: INTEGER DATA EXPECTED'
+     1  //' (NUMBER OF SOURCES).')
+        NSOUR=NITMA
+        IF(NSOUR.LE.0) CALL XABORT('PSOUR: INVALID NUMBER OF SOURCES.')
+        ALLOCATE(XMIN(NSOUR),XMAX(NSOUR),YMIN(NSOUR),YMAX(NSOUR),
+     1  ZMIN(NSOUR),ZMAX(NSOUR))
+! SOURCE INTENSITY PER ENERGY GROUP KEYWORD
+      ELSE IF(TEXT12.EQ.'INTG') THEN
+        ISIDEF=1
+        CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+        IF(INDIC.EQ.2) THEN
+          ISOUR(1)=FLOTT
+          IF(ISOUR(1).LT.0.0) CALL XABORT('PSOUR: INVALID INTENSITY.')
+          DO I=2,NG2
+            CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+            IF(INDIC.NE.2) CALL XABORT('PSOUR: REAL DATA EXPECTED.')
+            ISOUR(I)=FLOTT
+            IF(ISOUR(I).LT.0.0) CALL XABORT('PSOUR: INVALID INTENSITY.')
+          ENDDO
+          ISISOUR(1:NG2)=1
+        ELSE IF(INDIC.EQ.1) THEN
+          GN=NITMA
+          IF(GN.GT.NG2.OR.GN.LT.1) CALL XABORT('PSOUR: INVALID GROUP'
+     1    //' NUMBER.')
+          ISISOUR(GN)=1
+          CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+          IF(INDIC.NE.2) CALL XABORT('PSOUR: REAL DATA EXPECTED.')
+          ISOUR(GN)=FLOTT      
+          IF(ISOUR(I).LT.0.0) CALL XABORT('PSOUR: INVALID INTENSITY.')
+  12      CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+          IF(INDIC.EQ.1) THEN
+            GN=NITMA
+            IF(GN.GT.NG2.OR.GN.LT.1) CALL XABORT('PSOUR: INVALID GROUP'
+     1      //' NUMBER.')
+            CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+            IF(INDIC.NE.2) CALL XABORT('PSOUR: REAL DATA EXPECTED.')
+            ISOUR(GN)=FLOTT      
+            IF(ISOUR(I).LT.0.0) CALL XABORT('PSOUR: INVALID INTENSITY.')
+            GO TO 12
+          ELSE
+            GO TO 11
+          ENDIF
+        ELSE
+          CALL XABORT('PSOUR: REAL OR INTEGER DATA EXPECTED.')
+        ENDIF
+! (X,Y,Z) LIMITS KEYWORDS
+      ELSE IF(TEXT12.EQ.'XLIM') THEN
+        DO I=1,NSOUR
+        CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+        IF(INDIC.NE.2) CALL XABORT('PSOUR: REAL DATA EXPECTED.')
+        XMIN(I)=FLOTT
+        CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+        IF(INDIC.NE.2) CALL XABORT('PSOUR: REAL DATA EXPECTED.')
+        XMAX(I)=FLOTT
+        ISXDEF=1
+        ENDDO
+      ELSE IF(TEXT12.EQ.'YLIM') THEN       
+        IF(NDIM.LT.2) CALL XABORT('PSOUR: INVALID USE OF YLIM, DIM < 2')
+        DO I=1,NSOUR
+        CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+        IF(INDIC.NE.2) CALL XABORT('PSOUR: REAL DATA EXPECTED.')
+        YMIN(I)=FLOTT
+        CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+        IF(INDIC.NE.2) CALL XABORT('PSOUR: REAL DATA EXPECTED.')
+        YMAX(I)=FLOTT
+        ISYDEF=1
+        ENDDO
+      ELSE IF(TEXT12.EQ.'ZLIM') THEN
+        IF(NDIM.LT.3) CALL XABORT('PSOUR: INVALID USE OF YLIM, DIM < 3')
+        DO I=1,NSOUR
+        CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+        IF(INDIC.NE.2) CALL XABORT('PSOUR: REAL DATA EXPECTED.')
+        ZMIN(I)=FLOTT
+        CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+        IF(INDIC.NE.2) CALL XABORT('PSOUR: REAL DATA EXPECTED.')
+        ZMAX(I)=FLOTT
+        ISZDEF=1
+        ENDDO
+! BOUNDARY SOURCE LOCATION KEYWORD
+      ELSE IF(TEXT12.EQ.'X-') THEN
+        IF(MONOP.NE.0) CALL XABORT('PSOUR: BOUNDARY SOURCE ALREADY'
+     1  //'DEFINED')
+        MONOP=-1
+      ELSE IF(TEXT12.EQ.'X+') THEN
+        IF(MONOP.NE.0) CALL XABORT('PSOUR: BOUNDARY SOURCE ALREADY'
+     1  //'DEFINED')
+        MONOP=1
+      ELSE IF(TEXT12.EQ.'Y-') THEN
+        IF(NDIM.LT.2) CALL XABORT('PSOUR: INVALID USE OF YLIM, DIM < 2')
+        IF(MONOP.NE.0) CALL XABORT('PSOUR: BOUNDARY SOURCE ALREADY'
+     1  //'DEFINED')
+        MONOP=-2
+      ELSE IF(TEXT12.EQ.'Y+') THEN
+        IF(NDIM.LT.2) CALL XABORT('PSOUR: INVALID USE OF YLIM, DIM < 2')
+        IF(MONOP.NE.0) CALL XABORT('PSOUR: BOUNDARY SOURCE ALREADY'
+     1  //'DEFINED')
+        MONOP=2
+      ELSE IF(TEXT12.EQ.'Z-') THEN
+        IF(NDIM.LT.3) CALL XABORT('PSOUR: INVALID USE OF YLIM, DIM < 3')
+        IF(MONOP.NE.0) CALL XABORT('PSOUR: BOUNDARY SOURCE ALREADY'
+     1  //'DEFINED')
+        MONOP=-3
+      ELSE IF(TEXT12.EQ.'Z+') THEN
+        IF(NDIM.LT.3) CALL XABORT('PSOUR: INVALID USE OF YLIM, DIM < 3')
+        IF(MONOP.NE.0) CALL XABORT('PSOUR: BOUNDARY SOURCE ALREADY'
+     1  //'DEFINED')
+        MONOP=3
+      ELSE IF(TEXT12.EQ.'DIR') THEN
+        ! MONODIRECTIONAL SOURCE DIRCTION (2 ANGLES)
+        ISDIRDEF=1 
+        CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+        IF(INDIC.NE.2) CALL XABORT('PSOUR: REAL DATA EXPECTED.')
+        TEMP1=FLOTT 
+        CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+        IF(INDIC.NE.2) CALL XABORT('PSOUR: REAL DATA EXPECTED.')
+        TEMP2=FLOTT 
+        CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+        IF(INDIC.EQ.2) THEN
+          DIR(1)=TEMP1 !MU
+          DIR(2)=TEMP2 !ETA
+          DIR(3)=FLOTT !XI
+          IF(DIR(1).GT.1.0.OR.DIR(1).LT.-1.0.OR.DIR(2).GT.1.0.OR.
+     1    DIR(2).LT.-1.0.OR.DIR(3).GT.1.0.OR.DIR(3).LT.-1.0.OR.
+     2    DIR(1)**2+DIR(2)**2+DIR(3)**2.GT.1.0) CALL XABORT('PSOUR:'
+     3    //' INVALID DIRECTION COSINES VALUES.')
+        ELSE
+          IF(TEMP1.GT.180.0.OR.TEMP1.LT.0.0.OR.TEMP2.GT.90.0.OR.
+     1    TEMP2.LT.0.0) CALL XABORT('PSOUR: INVALID POLAR OR AZIMUTAL'
+     2    //' ANGLE.')
+          DIR(1)=COS(TEMP1*3.1416/180)                   !MU
+          DIR(2)=SQRT(1-DIR(1)**2)*COS(TEMP2*3.1416/180) !ETA
+          DIR(3)=SQRT(1-DIR(1)**2)*SIN(TEMP2*3.1416/180) !XI
+          GO TO 11
+        ENDIF
+      ELSE IF(TEXT12.EQ.';') THEN
+        GO TO 20
+      ELSE
+        CALL XABORT('PSOUR: '//TEXT12//' IS AN INVALID KEYWORD.')
+      ENDIF
+      GO TO 10
+*----
+*  COMPUTE THE TRANSITION SOURCE
+*----
+   20 IF(NPART.GT.0) THEN
+        IF(NPART.NE.NPART0) THEN
+          WRITE(HSMG,'(30HPSOUR: INVALID NUMBER OF RHS (,I2,8H) NUMBER,
+     1    25H OF COMPANION PARTICLES =,I2,1H.)') NENTRY,NPART
+          CALL XABORT(HSMG)
+        ELSE IF(STYPE.GE.0) THEN
+          CALL XABORT('PSOUR: TRANSITION AND BOUNDARY SOURCES PRESENT.')
+        ENDIF
+        IF(IMPX.GT.0) WRITE(IOUT,100)
+        IF(.NOT.C_ASSOCIATED(JPTRK)) CALL XABORT('PSOUR: NO TRACKING'
+     1  //' ASSOCIATED WITH THE COMPANION PARTICLE.')
+        DO IPART=1,NPART
+          CALL LCMGET(IPFLX(IPART),'STATE-VECTOR',ISTATE)
+          NG1=ISTATE(1)
+          NUNF=ISTATE(2)
+          ALLOCATE(FUNKNO(NUNF,NG1))
+          IF(IMPX.GT.0) WRITE(IOUT,110) IPART,HPARTC(NPART),HPFLX(NPART)
+          TEXT12='GROUP-'//HPARTC(IPART)
+          JPMAC=LCMGID(IPMAC,TEXT12)
+          IF(IADJ.EQ.0) THEN
+            JPFLX=LCMGID(IPFLX(IPART),'FLUX')
+          ELSE IF(IADJ.EQ.1) THEN
+            JPFLX=LCMGID(IPFLX(IPART),'AFLUX')
+          ENDIF
+          DO IG=1,NG1
+            IF(IADJ.EQ.0) THEN
+              CALL LCMGDL(JPFLX,IG,FUNKNO(1,IG))
+            ELSE IF(IADJ.EQ.1) THEN
+              CALL LCMGDL(JPFLX,NG1-IG+1,FUNKNO(1,IG))
+            ENDIF
+          ENDDO
+          DO IG=1,NG2
+            KPMAC=LCMGIL(JPMAC,IG)
+            CALL PSOUSN(NUNF,NUNS,IG,IPTRK,JPTRK,KPMAC,NANIS,NREG,
+     1      NMAT,NG1,NG2,MATCOD,FUNKNO,SUNKNO)
+          ENDDO
+          DEALLOCATE(FUNKNO)
+        ENDDO
+        ALLOCATE(NORM(NPART))
+        DO IPART=1,NPART
+          CALL LCMLEN(IPFLX(IPART),'NORM-FS',ILEN,ITYLCM)
+          IF(ILEN.GT.0) CALL LCMGET(IPFLX(IPART),'NORM-FS',NORM(IPART))
+        ENDDO
+        CALL LCMPUT(IPSOUR,'NORM-FS',1,2,SUM(NORM))
+        DEALLOCATE(MATCOD,NORM)
+*----
+*  BOUNDARY SOURCES
+*----
+      ELSE IF(STYPE.GE.0) THEN
+        IF(ISIDEF.EQ.0) CALL XABORT('PSOUR: SOURCE INTENSITY NOT DEFIN'
+     1  //'ED.')
+        IF(MONOP.EQ.0.AND.STYPE.EQ.1) CALL XABORT('PSOUR: BOUNDARY SOU'
+     1  //'RCE POSITION NOT DEFINED.')
+        IF(ISDIRDEF.EQ.0.AND.STYPE.EQ.1) CALL XABORT('PSOUR:' 
+     1  //' MONODIRECTIONAL SOURCE DIRECTION NOT DEFINED.')
+        IF(IMPX.GT.0) THEN
+          IF(STYPE.EQ.0) THEN
+            WRITE(IOUT,'(/34H PSOUR: ISOTROPIC BOUNDARY SOURCE.)')
+          ELSE IF(STYPE.EQ.1) THEN
+            WRITE(IOUT,'(/40H PSOUR: MONODIRECTIONAL BOUNDARY SOURCE.)')
+          ENDIF
+        ENDIF
+
+        ! X-BOUNDARY
+        IF(ISXDEF.EQ.1) THEN
+          IF(ABS(MONOP).EQ.1) CALL XABORT('PSOUR: BOUNDARY SOURCE'
+     1    //' X- OR X+ : XLIM SHOULD NOT BE DEFINED.')
+          DO I=1,NSOUR
+            IF(XMIN(I).LT.XXX(1).OR.XMAX(I).GT.XXX(MESH_LEN(1)))
+     1      CALL XABORT('PSOUR: SOURCE X-LIMITS OUTSIDE THE GEOMETRY.')
+            IF(XMIN(I).GE.XMAX(I)) CALL XABORT('PSOUR: SOURCE X-LIMITS'
+     1      //' INVALID.')
+          ENDDO
+        ELSE
+          IF(STYPE.EQ.0) CALL XABORT('PSOUR: SOURCE X-LIMITS NOT DEFIN'
+     1    //'ED.')
+        ENDIF
+
+        ! Y-BOUNDARY
+        IF(ISYDEF.EQ.1) THEN
+          IF(ABS(MONOP).EQ.2) CALL XABORT('PSOUR: BOUNDARY SOURCE'
+     1    //' Y- OR Y+ : YLIM SHOULD NOT BE DEFINED.')
+          DO I=1,NSOUR
+            IF(YMIN(I).LT.YYY(1).OR.YMAX(I).GT.YYY(MESH_LEN(2)))
+     1      CALL XABORT('PSOUR: SOURCE Y-LIMITS OUTSIDE THE GEOMETRY.')
+            IF(YMIN(I).GE.YMAX(I)) CALL XABORT('PSOUR: SOURCE Y-LIMITS'
+     1      //' INVALID.')
+          ENDDO
+        ELSE
+          IF(NDIM.GE.2.AND.STYPE.EQ.0) CALL XABORT('PSOUR: SOURCE Y-LI'
+     1    //'MITS NOT DEFINED.')
+          IF(NDIM.GE.2.AND.STYPE.EQ.1.AND.(ABS(MONOP).NE.2))
+     1    CALL XABORT('PSOUR: SOURCE Y-LIMITS NOT DEFINED.')
+        ENDIF
+
+        ! Z-BOUNDARY
+        IF(ISZDEF.EQ.1) THEN
+          IF(ABS(MONOP).EQ.3) CALL XABORT('PSOUR: BOUNDARY SOURCE'
+     1    //' Z- OR Z+ : ZLIM SHOULD NOT BE DEFINED.')
+          DO I=1,NSOUR
+            IF(ZMIN(I).LT.ZZZ(1).OR.ZMAX(I).GT.ZZZ(MESH_LEN(3)))
+     1      CALL XABORT('PSOUR: SOURCE Z-LIMITS OUTSIDE THE GEOMETRY.')
+            IF(ZMIN(I).GE.ZMAX(I)) CALL XABORT('PSOUR: SOURCE Z-LIMITS'
+     1      //' INVALID.')
+          ENDDO
+        ELSE
+          IF(NDIM.GE.3.AND.STYPE.EQ.0) CALL XABORT('PSOUR: SOURCE Z-LI'
+     1    //'MITS NOT DEFINED.')
+          IF(NDIM.GE.3.AND.STYPE.EQ.1.AND.(ABS(MONOP).NE.3))
+     1    CALL XABORT('PSOUR: SOURCE Z-LIMITS NOT DEFINED.')
+        ENDIF
+*----
+*  COMPUTE THE FIXED EXTERNAL SOURCES
+*----
+        IF(NDIM.EQ.1) THEN
+          NBST=SUM(ISISOUR)
+        ELSEIF(NDIM.EQ.2) THEN
+          NBST=SUM(ISISOUR)*2
+        ELSE
+          NBST=SUM(ISISOUR)*4
+        ENDIF
+
+        ALLOCATE(BSINFO2(2,NG2,NBST),BS2(MAXL,NG2,NBST),NBS2(NG2))
+        NBS2=0
+
+        IF(STYPE.EQ.0) THEN
+          CALL PSOISO(IPTRK,IPGEOM,NREG,LX,LY,LZ,NG2,NUNS,NDIM,
+     1    NSOUR,ISOUR,XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,XXX,YYY,ZZZ,
+     2    MESH_LEN,SUNKNO,ZNORM)
+        ELSE IF(STYPE.EQ.1) THEN
+          CALL PSOMON(IPTRK,IPGEOM,NREG,LX,LY,LZ,NG2,NDIM,NSOUR,
+     1    ISOUR,ISISOUR,XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,XXX,YYY,
+     2    ZZZ,MESH_LEN,BSINFO2,BS2,MAXL,ZNORM,DIR,MONOP,NBST,NBS2)
+        ELSE
+          CALL XABORT('PSOUR: SOURCE TYPE EXPECTED.')
+        ENDIF
+
+        DEALLOCATE(XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,XXX,YYY,ZZZ,ISOUR,
+     1  ISISOUR)
+
+        ! SAVE BOUNDARY SOURCE
+        IF(STYPE.EQ.0.AND.ISBS.EQ.1) THEN
+          JPSOURA=LCMLID(IPSOUR,'BS',NG2)
+          JPSOURB=LCMLID(IPSOUR,'BSINFO',NG2)
+          DO IG=1,NG2
+            IF(NBS1(IG).GT.0) THEN
+              KPSOURA=LCMLIL(JPSOURA,IG,NBS1(IG))
+              KPSOURB=LCMLIL(JPSOURB,IG,NBS1(IG))
+              DO N=1,NBS1(IG)
+                CALL LCMPDL(KPSOURA,N,MAXL,2,BS1(1,IG,N))
+                CALL LCMPDL(KPSOURB,N,2,1,BSINFO1(1,IG,N))
+              ENDDO
+            ENDIF
+          ENDDO
+          CALL LCMPUT(IPSOUR,'NBS',NG2,1,NBS1)
+        ELSE IF(STYPE.EQ.1.AND.ISBS.EQ.0) THEN
+          JPSOURA=LCMLID(IPSOUR,'BS',NG2)
+          JPSOURB=LCMLID(IPSOUR,'BSINFO',NG2)
+          DO IG=1,NG2
+            IF(NBS2(IG).GT.0) THEN
+              KPSOURA=LCMLIL(JPSOURA,IG,NBS2(IG))
+              KPSOURB=LCMLIL(JPSOURB,IG,NBS2(IG))
+              DO N=1,NBS2(IG)
+                CALL LCMPDL(KPSOURA,N,MAXL,2,BS2(1,IG,N))
+                CALL LCMPDL(KPSOURB,N,2,1,BSINFO2(1,IG,N))
+              ENDDO
+            ENDIF
+          ENDDO
+          CALL LCMPUT(IPSOUR,'NBS',NG2,1,NBS2)
+          DEALLOCATE(BS2,BSINFO2)      
+        ELSE IF(STYPE.EQ.1.AND.ISBS.EQ.1) THEN
+          JPSOURA=LCMLID(IPSOUR,'BS',NG2)
+          JPSOURB=LCMLID(IPSOUR,'BSINFO',NG2)
+          ALLOCATE(NBS(NG2))
+          NBS(:NG2)=NBS1(:NG2)+NBS2(:NG2)
+          DO IG=1,NG2
+            ALLOCATE(BS(MAXL,NBS(IG)),BSINFO(2,NBS(IG)))
+            BSINFO(:2,:NBS1(IG))=BSINFO1(:2,IG,:NBS1(IG))
+            BS(:MAXL,:NBS1(IG))=BS1(:MAXL,IG,:NBS1(IG))
+            CPT=1
+            DO N2=1,NBS2(IG)
+              ISD=0
+              DO N1=1,NBS1(IG)
+                IF(BSINFO2(1,IG,N2).EQ.BSINFO(1,N1).AND.
+     1             BSINFO2(2,IG,N2).EQ.BSINFO(2,N1)) THEN
+                  BS(:MAXL,N1) = BS(:MAXL,N1) + BS2(:MAXL,IG,N2)
+                  NBS(IG)=NBS(IG)-1
+                  ISD=1
+                ENDIF
+              ENDDO
+              IF(ISD.EQ.0) THEN
+                BSINFO(:2,NBS1(IG)+CPT)=BSINFO2(:2,IG,N2)
+                BS(:MAXL,NBS1(IG)+CPT)=BS2(:MAXL,IG,N2)
+                CPT=CPT+1
+              ENDIF
+            ENDDO 
+            IF(NBS(IG).GT.0) THEN
+              KPSOURA=LCMLIL(JPSOURA,IG,NBS(IG))
+              KPSOURB=LCMLIL(JPSOURB,IG,NBS(IG))
+              DO N=1,NBS(IG)
+                CALL LCMPDL(KPSOURA,N,MAXL,2,BS(1,N))
+                CALL LCMPDL(KPSOURB,N,2,1,BSINFO(1,N))
+              ENDDO
+            ENDIF
+            DEALLOCATE(BS,BSINFO)
+          ENDDO
+          CALL LCMPUT(IPSOUR,'NBS',NG2,1,NBS)
+          DEALLOCATE(BS2,BSINFO2,NBS,NBS1,NBS2)
+        ENDIF
+
+        ! Save the source normalization factor
+        IF(BSNORM.LT.0) THEN
+          CALL LCMPUT(IPSOUR,'NORM-FS',1,2,ZNORM)
+        ELSE
+          CALL LCMPUT(IPSOUR,'NORM-FS',1,2,ZNORM+BSNORM)
+        ENDIF
+      ENDIF
+      DEALLOCATE(BS1,BSINFO1)
+*----
+*  SAVE THE TRANSITION SOURCE ON LCM
+*----
+      IOF=1
+      NDIR=0
+      NCST=0
+      IF(IADJ.EQ.0) THEN
+        NDIR=1
+        JPSOUR=LCMLID(IPSOUR,'DSOUR',NDIR)
+      ELSE IF(IADJ.EQ.1) THEN
+        NCST=1
+        JPSOUR=LCMLID(IPSOUR,'ASOUR',NCST)
+      ENDIF
+      KPSOUR=LCMLIL(JPSOUR,IOF,NG2)
+      DO IG=1,NG2
+        CALL LCMPDL(KPSOUR,IG,NUNS,2,SUNKNO(1,IG))
+      ENDDO
+      DEALLOCATE(SUNKNO)
+*----
+*  SAVE THE SIGNATURE AND STATE VECTOR
+*----
+      HSIGN='L_SOURCE'
+      CALL LCMPTC(IPSOUR,'SIGNATURE',12,HSIGN)
+      ISTATE(:NSTATE)=0
+      ISTATE(1)=NG2
+      ISTATE(2)=NUNS
+      ISTATE(3)=NDIR
+      ISTATE(4)=NCST
+      IF(IMPX.GT.0) WRITE(6,120) (ISTATE(I),I=1,4)
+      CALL LCMPUT(IPSOUR,'STATE-VECTOR',NSTATE,1,ISTATE)
+      RETURN
+*
+  100 FORMAT(/44H PSOUR: TRANSITIONS FROM COMPANION PARTICLES/4X,
+     1 14HPARTICLE......, 3X,13H PRIMARY_FLUX)
+  110 FORMAT(1X,I4,A13,4X,A12)
+  120 FORMAT(/8H OPTIONS/8H -------/
+     1 7H NG    ,I8,28H   (NUMBER OF ENERGY GROUPS)/
+     2 7H NUN   ,I8,40H   (NUMBER OF UNKNOWNS PER ENERGY GROUP)/
+     3 7H NDIR  ,I8,35H   (NUMBER OF DIRECT FIXED SOURCES)/
+     4 7H NCST  ,I8,36H   (NUMBER OF ADJOINT FIXED SOURCES))
+      END