*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