*DECK PKINI SUBROUTINE PKINI(NENTRY,HENTRY,IENTRY,JENTRY,KENTRY) * *----------------------------------------------------------------------- * *Purpose: * Collect input information for the point kinetic module. * *Copyright: * Copyright (C) 2017 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 * NENTRY number of data structures transfered to this module. * HENTRY name of the data structures. * IENTRY data structure type where: * IENTRY=1 for LCM memory object; * IENTRY=2 for XSM file; * IENTRY=3 for sequential binary file; * IENTRY=4 for sequential ASCII file. * JENTRY access permission for the data structure where: * JENTRY=0 for a data structure in creation mode; * JENTRY=1 for a data structure in modifications mode; * JENTRY=2 for a data structure in read-only mode. * KENTRY data structure pointer. * *Comments: * The PKINI: module specification is: * MAPFL := PKINI: MAPFL :: (descpkini) ; * where * MAPFL : name of the \emph{map} object containing fuel regions description * and global parameter informations. * (descpkini) : structure describing the input data to the PKINI: module. * *----------------------------------------------------------------------- * USE GANLIB *---- * SUBROUTINE ARGUMENTS *---- INTEGER NENTRY,IENTRY(NENTRY),JENTRY(NENTRY) TYPE(C_PTR) KENTRY(NENTRY) CHARACTER HENTRY(NENTRY)*12 *---- * LOCAL VARIABLES *---- PARAMETER (NSTATE=40,MAXALP=10,MAXTIM=2) INTEGER ISTATE(NSTATE) TYPE(C_PTR) IPMAP,JPMAP,KPMAP,JPPAR,KPPAR REAL LAMBDA,TIMES(MAXTIM) DOUBLE PRECISION DFLOT LOGICAL LCUBIC CHARACTER TEXT12*12,HSIGN*12,HPNAME*12,HSMG*131,HPARAM(MAXALP)*12 *---- * ALLOCATABLE ARRAYS *---- REAL, ALLOCATABLE, DIMENSION(:) :: BETAI,LAMBDAI,X,Y,PARAMI,FPOWER REAL, ALLOCATABLE, DIMENSION(:,:) :: VAL DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: YINIT *---- * RECOVER THE FUELMAP *---- IF(NENTRY.NE.1) CALL XABORT('PKINI: ONE PARAMETER EXPECTED.') IF((IENTRY(1).NE.1).AND.(IENTRY(1).NE.2))CALL XABORT('PKINI:' 1 //' LCM OBJECT EXPECTED.') IF(JENTRY(1).NE.1) CALL XABORT('PKINI: SECOND ENTRY IN MODIFICATI' 1 //'ON MODE EXPECTED.') IPMAP=KENTRY(1) CALL LCMGTC(IPMAP,'SIGNATURE',12,HSIGN) IF(HSIGN.NE.'L_MAP') THEN TEXT12=HENTRY(1) CALL XABORT('PKINI: SIGNATURE OF '//TEXT12//' IS '//HSIGN// 1 '. L_MAP EXPECTED.') ENDIF CALL LCMGET(IPMAP,'STATE-VECTOR',ISTATE) NB=ISTATE(1) NCH=ISTATE(2) NPARM=ISTATE(8) IF(NPARM.GT.0) JPPAR=LCMGID(IPMAP,'PARAM') IF(NCH.NE.1) CALL XABORT('PKINI: ONE CHANNEL EXPECTED.') CALL LCMSIX(IPMAP,'P-KINETIC',1) *---- * READ INPUT DATA *---- IMPX=1 NGROUP=0 NALPHA=0 NPTIME=0 EPSILON=1.0E-2 POW0=0.0 LAMBDA=0.0 10 CALL REDGET(ITYP,NITMA,FLOT,TEXT12,DFLOT) IF(ITYP.NE.3)CALL XABORT('PKINI: CHARACTER DATA EXPECTED.') IF(TEXT12.EQ.'EDIT') THEN * READ PRINTING INDEX CALL REDGET(ITYP,IMPX,FLOT,TEXT12,DFLOT) IF(ITYP.NE.1)CALL XABORT('PKINI: INTEGER FOR EDIT EXPECTED.') ELSE IF(TEXT12.EQ.'POWER') THEN * Initial power (MW) CALL REDGET(ITYP,NITMA,POW0,TEXT12,DFLOT) IF(ITYP.NE.2)CALL XABORT('PKINI: REAL FOR P0 EXPECTED.') ELSE IF(TEXT12.EQ.'LAMBDA') THEN * Prompt neutron generation time (s) CALL REDGET(ITYP,NITMA,LAMBDA,TEXT12,DFLOT) IF(ITYP.NE.2)CALL XABORT('PKINI: REAL FOR LAMBDA EXPECTED.') CALL LCMPUT(IPMAP,'LAMBDA',1,2,LAMBDA) ELSE IF(TEXT12.EQ.'EPSILON') THEN * Rugge-Kutta EPSILON CALL REDGET(ITYP,NITMA,FLOT,TEXT12,DFLOT) IF(ITYP.NE.2)CALL XABORT('PKINI: REAL FOR EPSILON EXPECTED.') CALL LCMPUT(IPMAP,'EPSILON',1,2,FLOT) ELSE IF(TEXT12.EQ.'TIME') THEN * Set initial time and stage length (s) CALL REDGET(ITYP,NITMA,T,TEXT12,DFLOT) IF(ITYP.NE.2)CALL XABORT('PKINI: REAL FOR T EXPECTED.') CALL REDGET(ITYP,NITMA,DT,TEXT12,DFLOT) IF(ITYP.NE.2)CALL XABORT('PKINI: REAL FOR DT EXPECTED.') CALL LCMPUT(IPMAP,'T-VALUE_INIT',1,2,T) IF(IMPX.GT.0) WRITE(6,100) T,DT ELSE IF(TEXT12.EQ.'NGROUP') THEN * Read printing index CALL REDGET(ITYP,NGROUP,FLOT,TEXT12,DFLOT) IF(ITYP.NE.1)CALL XABORT('PKINI: INTEGER FOR NGROUP EXPECTED.') ELSE IF(TEXT12.EQ.'BETAI') THEN * Delayed neutron fraction IF(NGROUP.EQ.0)CALL XABORT('PKINI: NGROUP NOT DEFINED.') ALLOCATE(BETAI(NGROUP)) DO IG=1,NGROUP CALL REDGET(ITYP,NITMA,BETAI(IG),TEXT12,DFLOT) IF(ITYP.NE.2)CALL XABORT('PKINI: REAL FOR BETAI EXPECTED.') ENDDO CALL LCMPUT(IPMAP,'BETAI',NGROUP,2,BETAI) ELSE IF(TEXT12.EQ.'LAMBDAI') THEN * Delayed neutron time constant IF(NGROUP.EQ.0)CALL XABORT('PKINI: NGROUP NOT DEFINED.') ALLOCATE(LAMBDAI(NGROUP)) DO IG=1,NGROUP CALL REDGET(ITYP,NITMA,LAMBDAI(IG),TEXT12,DFLOT) IF(ITYP.NE.2)CALL XABORT('PKINI: REAL FOR LAMBDAI EXPECTED.') ENDDO CALL LCMPUT(IPMAP,'LAMBDAI',NGROUP,2,LAMBDAI) ELSE IF(TEXT12.EQ.'ALPHA') THEN IF(NPARM.EQ.0)CALL XABORT('PKINI: NPARM NOT DEFINED.') JPMAP=LCMLID(IPMAP,'ALPHA',MAXALP) 20 CALL REDGET(ITYP,NITMA,FLOT,TEXT12,DFLOT) IF(ITYP.NE.3)CALL XABORT('PKINI: CHARACTER DATA EXPECTED.') IF(TEXT12.EQ.'ENDA') GO TO 10 NALPHA=NALPHA+1 IF(NALPHA.GT.MAXALP) CALL XABORT('PKINI: MAXALP OVERFLOW.') DO IPAR=1,NPARM KPPAR=LCMGIL(JPPAR,IPAR) CALL LCMGTC(KPPAR,'P-NAME',12,HPNAME) IF(HPNAME.EQ.TEXT12) GO TO 25 ENDDO WRITE(HSMG,'(24HPKINI: GLOBAL PARAMETER ,A,16H IS NOT DEFINED , 1 15HIN THE FUELMAP.)') TEXT12 CALL XABORT(HSMG) 25 KPMAP=LCMDIL(JPMAP,NALPHA) CALL LCMPTC(KPMAP,'P-NAME',12,TEXT12) HPARAM(NALPHA)=TEXT12 CALL REDGET(ITYP,NITMA,FLOT,TEXT12,DFLOT) IF(ITYP.NE.3)CALL XABORT('PKINI: CHARACTER DATA EXPECTED.') IF(TEXT12.EQ.'DIRECT') THEN ITYPE=1 ELSE IF(TEXT12.EQ.'DERIV') THEN ITYPE=2 ELSE IF(TEXT12.EQ.'SQDERIV') THEN ITYPE=3 ELSE CALL XABORT('PKINI: DIRECT OR DERIV EXPECTED.') ENDIF CALL REDGET(ITYP,NXY,FLOT,TEXT12,DFLOT) LCUBIC=.FALSE. IF(ITYP.EQ.3) THEN IF(TEXT12.EQ.'LINEAR') THEN LCUBIC=.FALSE. ELSE IF(TEXT12.EQ.'CUBIC') THEN LCUBIC=.TRUE. ELSE CALL XABORT('PKINI: LINEAR OR CUBIC EXPECTED.') ENDIF CALL REDGET(ITYP,NXY,FLOT,TEXT12,DFLOT) ENDIF IF(ITYP.NE.1)CALL XABORT('PKINI: INTEGER FOR NXY EXPECTED(1).') ALLOCATE(X(NXY),Y(NXY)) DO I=1,NXY CALL REDGET(ITYP,NITMA,X(I),TEXT12,DFLOT) IF(ITYP.NE.2)CALL XABORT('PKINI: REAL FOR X EXPECTED(1).') CALL REDGET(ITYP,NITMA,Y(I),TEXT12,DFLOT) IF(ITYP.NE.2)CALL XABORT('PKINI: REAL FOR Y EXPECTED(1).') ENDDO CALL LCMPUT(KPMAP,'ALPHA-LAW-P',NXY,2,X) CALL LCMPUT(KPMAP,'ALPHA-LAW-R',NXY,2,Y) CALL LCMPUT(KPMAP,'ALPHA-LAW-T',1,1,ITYPE) CALL LCMPUT(KPMAP,'ALPHA-LAW-I',1,5,LCUBIC) DEALLOCATE(Y,X) GO TO 20 ELSE IF(TEXT12.EQ.'PTIME') THEN IF(NALPHA.EQ.0)CALL XABORT('PKINI: NO FEEDBACK PARAMETERS.') JPMAP=LCMGID(IPMAP,'ALPHA') CALL REDGET(ITYP,NITMA,FLOT,TEXT12,DFLOT) IF(ITYP.NE.3)CALL XABORT('PKINI: CHARACTER DATA EXPECTED.') 30 IF(TEXT12.EQ.'ENDP') GO TO 10 NPTIME=NPTIME+1 IALP=0 DO IAL=1,NALPHA IALP=IAL IF(TEXT12.EQ.HPARAM(IAL)) GO TO 40 ENDDO WRITE(HSMG,'(24HPKINI: GLOBAL PARAMETER ,A,16H IS NOT A FEEDBA, 1 13HCK PARAMETER.)') TEXT12 CALL XABORT(HSMG) 40 KPMAP=LCMDIL(JPMAP,IALP) LCUBIC=.FALSE. 50 CALL REDGET(ITYP,NITMA,FLOT,TEXT12,DFLOT) IF(ITYP.EQ.1) THEN NXY=NITMA ALLOCATE(X(NXY),Y(NXY)) DO I=1,NXY CALL REDGET(ITYP,NITMA,X(I),TEXT12,DFLOT) IF(ITYP.NE.2)CALL XABORT('PKINI: REAL FOR X EXPECTED(2).') CALL REDGET(ITYP,NITMA,Y(I),TEXT12,DFLOT) IF(ITYP.NE.2)CALL XABORT('PKINI: REAL FOR Y EXPECTED(2).') ENDDO CALL REDGET(ITYP,NITMA,FLOT,TEXT12,DFLOT) IF(ITYP.NE.3)CALL XABORT('PKINI: CHARACTER DATA EXPECTED.') ELSE IF(ITYP.EQ.3) THEN IF(TEXT12.EQ.'LINEAR') THEN LCUBIC=.FALSE. GO TO 50 ELSE IF(TEXT12.EQ.'CUBIC') THEN LCUBIC=.TRUE. GO TO 50 ELSE IF(TEXT12.EQ.'T-DELT') THEN GO TO 60 ELSE CALL XABORT('PKINI: LINEAR, CUBIC OR T-DELT EXPECTED.') ENDIF 60 CALL REDGET(ITYP,NITMA,FLOT,TEXT12,DFLOT) IF(ITYP.EQ.1) THEN NXY=NITMA CALL REDGET(ITYP,NITMA,T1,TEXT12,DFLOT) IF(ITYP.NE.2) CALL XABORT('PKINI: REAL DATA EXPECTED(1).') ELSE IF (ITYP.EQ.2) THEN NXY=1001 T1=FLOT ELSE CALL XABORT('PKINI: INTEGER OR REAL DATA EXPECTED.') ENDIF ALLOCATE(X(NXY),Y(NXY)) CALL REDGET(ITYP,NITMA,T2,TEXT12,DFLOT) IF(ITYP.NE.2) CALL XABORT('PKINI: REAL DATA EXPECTED(2).') IF(T2.LE.T1) CALL XABORT('PKINI: T2 > T1 EXPECTED.') DELT=(T2-T1)/REAL(NXY-1) TT=T1 CALL REDGET(ITYP,NITMA,FLOT,TEXT12,DFLOT) IF(ITYP.NE.3)CALL XABORT('PKINI: CHARACTER DATA EXPECTED.') IF(TEXT12.NE.'P-VALV') CALL XABORT('PKINI: P-VALV EXPECTED.') CALL REDGET(ITYP,NITMA,GAMMA,TEXT12,DFLOT) IF(ITYP.NE.2) CALL XABORT('PKINI: REAL DATA EXPECTED(3).') CALL REDGET(ITYP,NITMA,PINIT,TEXT12,DFLOT) IF(ITYP.NE.2) CALL XABORT('PKINI: REAL DATA EXPECTED(4).') CALL REDGET(ITYP,NITMA,PFINAL,TEXT12,DFLOT) IF(ITYP.NE.2) CALL XABORT('PKINI: REAL DATA EXPECTED(5).') CALL REDGET(ITYP,NITMA,TB1,TEXT12,DFLOT) IF(ITYP.NE.2) CALL XABORT('PKINI: REAL DATA EXPECTED(6).') IF(TB1.LT.T1) CALL XABORT('PKINI: INVALID VALUE OF TB1.') CALL REDGET(ITYP,NITMA,B1,TEXT12,DFLOT) IF(ITYP.NE.2) CALL XABORT('PKINI: REAL DATA EXPECTED(7).') CALL REDGET(ITYP,NITMA,FLOT,TEXT12,DFLOT) IF(ITYP.NE.3)CALL XABORT('PKINI: CHARACTER DATA EXPECTED.') TB2=T2 PFINA2=PFINAL IF(TEXT12.EQ.'RESET') THEN CALL REDGET(ITYP,NITMA,PFINA2,TEXT12,DFLOT) IF(ITYP.NE.2) CALL XABORT('PKINI: REAL DATA EXPECTED(8).') IF(PFINA2.GT.PFINAL) CALL XABORT('PKINI: INVALID VALUE OF' > //' PFINA2.') CALL REDGET(ITYP,NITMA,TB2,TEXT12,DFLOT) IF(ITYP.NE.2) CALL XABORT('PKINI: REAL DATA EXPECTED(9).') IF(TB2.LE.TB1) CALL XABORT('PKINI: INVALID VALUE OF TB2.') CALL REDGET(ITYP,NITMA,B2,TEXT12,DFLOT) IF(ITYP.NE.2) CALL XABORT('PKINI: REAL DATA EXPECTED(10).') CALL REDGET(ITYP,NITMA,FLOT,TEXT12,DFLOT) IF(ITYP.NE.3)CALL XABORT('PKINI: CHARACTER DATA EXPECTED.') ENDIF * ALPHA=2.0*GAMMA/(GAMMA-1.0) YGAR=PINIT DO I=1,NXY X(I)=TT IF(TT.LE.TB1) THEN Y(I)=PINIT ELSE IF(TT.LE.TB2) THEN Y(I)=MAX(PINIT/((B1*(TT-TB1)+1.0))**ALPHA,PFINAL) YGAR=Y(I) ELSE Y(I)=MAX(YGAR/((B2*(TT-TB2)+1.0))**ALPHA,PFINA2) ENDIF TT=TT+DELT ENDDO ELSE CALL XABORT('PKINI: INTEGER OR CHARACTER DATA EXPECTED.') ENDIF CALL LCMPUT(KPMAP,'TIME-LAW-T',NXY,2,X) CALL LCMPUT(KPMAP,'TIME-LAW-P',NXY,2,Y) CALL LCMPUT(KPMAP,'TIME-LAW-I',1,5,LCUBIC) DEALLOCATE(Y,X) GO TO 30 ELSE IF(TEXT12.EQ.';') THEN GO TO 70 ELSE CALL XABORT('PKINI: INVALID KEYWORD: '//TEXT12//'.') ENDIF GO TO 10 *---- * RECOVER THE INITIAL GLOBAL PARAMETER VALUES *---- 70 ALLOCATE(PARAMI(NALPHA)) IF(IMPX.GT.0) WRITE(6,110) DO IAL=1,NALPHA DO IPAR=1,NPARM KPPAR=LCMGIL(JPPAR,IPAR) CALL LCMGTC(KPPAR,'P-NAME',12,HPNAME) IF(HPNAME.EQ.HPARAM(IAL)) GO TO 80 ENDDO CALL XABORT('PKINI: GLOBAL PARAMETER NOT FOUND.') 80 CALL LCMLEN(KPPAR,'P-VALUE',ILONG,ITYLCM) IF(ILONG.EQ.0) THEN WRITE(HSMG,'(33HPKINI: VALUE OF GLOBAL PARAMETER ,A,6H IS NO, 1 20H SET IN THE FUELMAP.)') HPNAME CALL XABORT(HSMG) ENDIF CALL LCMGET(KPPAR,'P-VALUE',PARAMI(IAL)) IF(IMPX.GT.0) WRITE(6,120) IAL,HPNAME,PARAMI(IAL) ENDDO CALL LCMPUT(IPMAP,'P-VALUE-INIT',NALPHA,2,PARAMI) CALL LCMPTC(IPMAP,'P-NAME',12,NALPHA,HPARAM) *---- * SET INITIAL SOLUTION OF POINT KINETIC EQUATIONS *---- IF(POW0.EQ.0.0) CALL XABORT('PKINI: INITIAL POWER NOT DEFINED.') IF(NGROUP.EQ.0) CALL XABORT('PKINI: NGROUP NOT DEFINED.') ALLOCATE(YINIT(NGROUP+1)) YINIT(1)=POW0 DO I=2,NGROUP+1 YINIT(I)=POW0*BETAI(I-1)/(LAMBDAI(I-1)*LAMBDA) ENDDO *---- * SAVE INITIAL CONDITIONS *---- ISTAGE=1 TIMES(:MAXTIM)=-1.0E30 TIMES(1)=T TEXT12='PKIN-DAT0001' IF(IMPX.GT.0) WRITE(6,130) T,TEXT12 CALL LCMSIX(IPMAP,TEXT12,1) CALL LCMPUT(IPMAP,'P-VALUE',NALPHA,2,PARAMI) CALL LCMPUT(IPMAP,'Y-VALUE',NGROUP+1,4,YINIT) CALL LCMPUT(IPMAP,'T-VALUE',1,2,T) CALL LCMPUT(IPMAP,'DT-VALUE',1,2,DT) CALL LCMPUT(IPMAP,'I-VALUE',1,1,ISTAGE) CALL LCMSIX(IPMAP,' ',2) DEALLOCATE(PARAMI,YINIT,LAMBDAI,BETAI) CALL LCMPUT(IPMAP,'PKIN-TIMES',MAXTIM,2,TIMES) *---- * CREATE STATE VECTOR AND SAVE POWER INFORMATION *---- ISTATE(:NSTATE)=0 ISTATE(1)=ISTAGE ISTATE(2)=NGROUP ISTATE(3)=NALPHA ISTATE(4)=NPTIME CALL LCMPUT(IPMAP,'STATE-VECTOR',NSTATE,1,ISTATE) IF(IMPX.GT.0) WRITE(6,140) ISTATE(2:4) IF(IMPX.GT.1) CALL LCMLIB(IPMAP) CALL LCMSIX(IPMAP,' ',2) IF(IMPX.GT.0) WRITE(6,150) POW0 CALL LCMPUT(IPMAP,'REACTOR-PW',1,2,POW0) CALL LCMLEN(IPMAP,'BUND-PW',ILONG,ITYLCM) CALL LCMLEN(IPMAP,'AXIAL-FPW',JLONG,ITYLCM) IF((ILONG.EQ.NCH*NB).AND.(JLONG.EQ.NB)) THEN ALLOCATE(VAL(NCH,NB),FPOWER(ILONG)) CALL LCMGET(IPMAP,'AXIAL-FPW',FPOWER) DSUM=0.0 DO IB=1,NB DSUM=DSUM+FPOWER(IB) ENDDO DO ICH=1,NCH DO IB=1,NB VAL(ICH,IB)=FPOWER(IB)*POW0*1.0E3/(DSUM*REAL(NCH)) ENDDO ENDDO CALL LCMPUT(IPMAP,'BUND-PW',NCH*NB,2,VAL) DEALLOCATE(FPOWER,VAL) ENDIF RETURN * 100 FORMAT(/34H PKINI: DEFINE INITIAL STAGE TIME=,1P,E12.4,2H S, 1 10H DURATION=,E12.4,2H S) 110 FORMAT(43H PKINI: GLOBAL PARAMETER -- INITIAL VALUES:) 120 FORMAT(1X,I8,A14,3H = ,1P,E12.4) 130 FORMAT(/40H PKINI: SAVE INFORMATION RELATED TO TIME,1P,E12.4, 1 27H S IN LCM DIRECTORY NAMED ',A12,2H'.) 140 FORMAT(/ 1 14H STATE VECTOR:/ 2 7H NGROUP,I9,39H (NUMBER OF DELAYED PRECURSOR GROUPS)/ 3 7H NALPHA,I9,34H (NUMBER OF FEEDBACK PARAMETERS)/ 4 7H NPTIME,I9,40H (NUMBER OF TIME-DEPENDENT PARAMETERS)) 150 FORMAT(/22H PKINI: INITIAL POWER=,1P,E12.4,3H MW) END