From 7dfcc480ba1e19bd3232349fc733caef94034292 Mon Sep 17 00:00:00 2001 From: stainer_t Date: Mon, 8 Sep 2025 13:48:49 +0200 Subject: Initial commit from Polytechnique Montreal --- Dragon/src/PSOUR.f | 706 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 706 insertions(+) create mode 100644 Dragon/src/PSOUR.f (limited to 'Dragon/src/PSOUR.f') 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 -- cgit v1.2.3