diff options
Diffstat (limited to 'Donjon/src/PKINI.f')
| -rw-r--r-- | Donjon/src/PKINI.f | 417 |
1 files changed, 417 insertions, 0 deletions
diff --git a/Donjon/src/PKINI.f b/Donjon/src/PKINI.f new file mode 100644 index 0000000..48fac0a --- /dev/null +++ b/Donjon/src/PKINI.f @@ -0,0 +1,417 @@ +*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 |