diff options
Diffstat (limited to 'Trivac/src/KINSRC.f')
| -rwxr-xr-x | Trivac/src/KINSRC.f | 257 |
1 files changed, 257 insertions, 0 deletions
diff --git a/Trivac/src/KINSRC.f b/Trivac/src/KINSRC.f new file mode 100755 index 0000000..1af7f00 --- /dev/null +++ b/Trivac/src/KINSRC.f @@ -0,0 +1,257 @@ +*DECK KINSRC + SUBROUTINE KINSRC(IPTRK,IPSYS,CMOD,IMPX,IFL,IPR,IEXP,NGR,NBM, + 1 NBFIS,NDG,ITY,LL4,NUN,NUP,PDC,TTF,TTP,DT,ADJ,OVR,CHI,CHD,SGF, + 2 SGD,OMEGA,EVECT,PC,SRC) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Compute the space-time kinetics source for a known neutron flux. +* +*Copyright: +* Copyright (C) 2010 Ecole Polytechnique de Montreal. +* +*Author(s): A. Hebert +* +*Parameters: input +* IPTRK pointer to L_TRACK object. +* IPSYS pointer to L_SYSTEM object. +* CMOD name of the assembly door (BIVAC or TRIVAC). +* IMPX print parameter (equal to zero for no print). +* IFL integration scheme for fluxes: =1 implicit; +* =2 Crank-Nicholson; =3 theta. +* IPR integration scheme for precursors: =1 implicit; +* =2 Crank-Nicholson; =3 theta; =4 exponential. +* IEXP exponential transformation flag (=1 to activate). +* NGR number of energy groups. +* NBM number of material mixtures. +* NBFIS number of fissile isotopes. +* NDG number of delayed-neutron groups. +* ITY type of solution: =1: classical Bivac/diffusion; +* =2: classical Trivac/diffusion; =3 Raviart-Thomas in +* Trivac/diffusion; =11: Bivac/SPN; =13 Trivac/SPN. +* LL4 order of the system matrices. +* NUN total number of unknowns per energy group. +* NUP total number of precursor unknowns per precursor group. +* PDC precursor decay constants. +* TTF value of theta-parameter for fluxes. +* TTP value of theta-parameter for precursors. +* DT current time increment. +* ADJ flag for adjoint space-time kinetics calculation +* OVR reciprocal neutron velocities/DT. +* CHI steady-state fission spectrum. +* CHD delayed fission spectrum +* SGF nu*fission macroscopic x-sections/keff. +* SGD delayed nu*fission macroscopic x-sections/keff. +* OMEGA exponential transformation parameter. +* EVECT neutron flux. +* PC precursor concentrations. +* +*Parameters: output +* SRC space-time kinetics source. +* +*----------------------------------------------------------------------- +* + USE GANLIB +*---- +* SUBROUTINE ARGUMENTS +*---- + TYPE(C_PTR) IPTRK,IPSYS + INTEGER IMPX,IFL,IPR,IEXP,NGR,NBM,NBFIS,NDG,ITY,LL4,NUN,NUP + REAL PDC(NDG),TTF,TTP,DT,OVR(NBM,NGR),CHI(NBM,NBFIS,NGR), + 1 CHD(NBM,NBFIS,NGR,NDG),SGF(NBM,NBFIS,NGR),SGD(NBM,NBFIS,NGR,NDG), + 2 OMEGA(NBM,NGR),EVECT(NUN,NGR),PC(NUP,NDG,NBFIS) + DOUBLE PRECISION SRC(NUN,NGR) + CHARACTER CMOD*12 + LOGICAL ADJ +*---- +* LOCAL VARIABLES +*---- + PARAMETER(IOS=6) + DOUBLE PRECISION DTF,DTP,DARG,DK,DSM + LOGICAL LFIS + CHARACTER TEXT12*12 + REAL, DIMENSION(:), ALLOCATABLE :: WORK1,WORK2,CHEXP + REAL, DIMENSION(:), POINTER :: AGAR + TYPE(C_PTR) AGAR_PTR +*---- +* SCRATCH STORAGE ALLOCATION +*---- + ALLOCATE(WORK1(LL4),WORK2(NBM),CHEXP(NBM)) +* + DTF=9999.0D0 + DTP=9999.0D0 + IF(IFL.EQ.1)THEN + DTF=1.0D0 + ELSEIF(IFL.EQ.2)THEN + DTF=0.5D0 + ELSEIF(IFL.EQ.3)THEN + DTF=DBLE(TTF) + ENDIF + IF(IPR.EQ.2)THEN + DTP=0.5D0 + ELSEIF(IPR.EQ.3)THEN + DTP=DBLE(TTP) + ENDIF +* + IF(IMPX.GT.0) WRITE(IOS,1001) CMOD + SRC(:NUN,:NGR)=0.0D0 + DO 200 IGR=1,NGR + WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR + CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,EVECT(1,IGR),WORK1) + DO 10 IND=1,LL4 + SRC(IND,IGR)=-(1.0D0-DTF)*WORK1(IND) + 10 CONTINUE + DO 40 JGR=1,NGR + IF(JGR.EQ.IGR) GO TO 40 + IF(.NOT.ADJ) THEN + WRITE(TEXT12,'(1HA,2I3.3)') IGR,JGR + ELSE + WRITE(TEXT12,'(1HA,2I3.3)') JGR,IGR + ENDIF + CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM) + IF(ILONG.EQ.0) GO TO 40 + IF((CMOD.EQ.'BIVAC').OR.(ITY.EQ.13))THEN + CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,EVECT(1,JGR),WORK1) + DO 20 IND=1,LL4 + SRC(IND,IGR)=SRC(IND,IGR)+(1.0D0-DTF)*WORK1(IND) + 20 CONTINUE + ELSEIF(CMOD.EQ.'TRIVAC')THEN + CALL LCMGPD(IPSYS,TEXT12,AGAR_PTR) + CALL C_F_POINTER(AGAR_PTR,AGAR,(/ ILONG /)) + DO 30 IND=1,ILONG + SRC(IND,IGR)=SRC(IND,IGR)+(1.0D0-DTF)*AGAR(IND)*EVECT(IND,JGR) + 30 CONTINUE + ENDIF + 40 CONTINUE +*---- +* PRECURSOR CONTRIBUTION +*---- + DO 180 IFIS=1,NBFIS + DO 90 IDG=1,NDG + DARG=PDC(IDG)*DT + IF(IPR.EQ.1)THEN + DK=1.0D0/(1.0D0+DARG) + ELSEIF(IPR.EQ.4)THEN + DK=DEXP(-DARG) + ELSE + DK=(1.0D0-(1.0D0-DTP)*DARG)/(1.0D0+DTP*DARG) + ENDIF + DSM=1.0D0-DTF+DTF*DK + LFIS=.FALSE. + DO 50 IBM=1,NBM + LFIS=LFIS.OR.(CHD(IBM,IFIS,IGR,IDG).NE.0.0) + 50 CONTINUE + IF(LFIS) THEN + IF(.NOT.ADJ) THEN + DO 60 IBM=1,NBM + IF(IEXP.EQ.0) THEN + CHEXP(IBM)=CHD(IBM,IFIS,IGR,IDG) + ELSE +* exponential transformation + CHEXP(IBM)=CHD(IBM,IFIS,IGR,IDG)*EXP(-OMEGA(IBM,IGR)*DT) + ENDIF + 60 CONTINUE + ELSE + DO 70 IBM=1,NBM + IF(IEXP.EQ.0) THEN + CHEXP(IBM)=SGD(IBM,IFIS,IGR,IDG) + ELSE +* exponential transformation + CHEXP(IBM)=SGD(IBM,IFIS,IGR,IDG)*EXP(-OMEGA(IBM,IGR)*DT) + ENDIF + 70 CONTINUE + ENDIF + IF(CMOD.EQ.'BIVAC')THEN + CALL KINBLM(IPTRK,NBM,LL4,CHEXP(1),PC(1,IDG,IFIS),WORK1) + ELSEIF(CMOD.EQ.'TRIVAC')THEN + CALL KINTLM(IPTRK,NBM,LL4,CHEXP(1),PC(1,IDG,IFIS),WORK1) + ENDIF + DO 80 IND=1,LL4 + SRC(IND,IGR)=SRC(IND,IGR)+PDC(IDG)*DSM*WORK1(IND) + 80 CONTINUE + ENDIF + 90 CONTINUE +*---- +* FISSION CONTRIBUTION +*---- + DO 170 JGR=1,NGR + IF(.NOT.ADJ) THEN + DO 100 IBM=1,NBM + WORK2(IBM)=CHI(IBM,IFIS,IGR)*SGF(IBM,IFIS,JGR) + 100 CONTINUE + ELSE + DO 110 IBM=1,NBM + WORK2(IBM)=CHI(IBM,IFIS,JGR)*SGF(IBM,IFIS,IGR) + 110 CONTINUE + ENDIF + IF(CMOD.EQ.'BIVAC')THEN + CALL KINBLM(IPTRK,NBM,LL4,WORK2(1),EVECT(1,JGR),WORK1) + ELSEIF(CMOD.EQ.'TRIVAC')THEN + CALL KINTLM(IPTRK,NBM,LL4,WORK2(1),EVECT(1,JGR),WORK1) + ENDIF + DO 120 IND=1,LL4 + SRC(IND,IGR)=SRC(IND,IGR)+(1.0D0-DTF)*WORK1(IND) + 120 CONTINUE + DO 160 IDG=1,NDG + DARG=PDC(IDG)*DT + IF(IPR.EQ.1)THEN + DK=0.0D0 + ELSEIF(IPR.EQ.4)THEN + DK=(1.0D0-DEXP(-DARG))/DARG-DEXP(-DARG) + ELSE + DK=(1.0D0-DTP)*DARG/(1.0D0+DTP*DARG) + ENDIF + DSM=1.0D0-DTF-DTF*DK + IF(.NOT.ADJ) THEN + DO 130 IBM=1,NBM + WORK2(IBM)=CHD(IBM,IFIS,IGR,IDG)*SGD(IBM,IFIS,JGR,IDG) + 130 CONTINUE + ELSE + DO 140 IBM=1,NBM + WORK2(IBM)=CHD(IBM,IFIS,JGR,IDG)*SGD(IBM,IFIS,IGR,IDG) + 140 CONTINUE + ENDIF + IF(CMOD.EQ.'BIVAC')THEN + CALL KINBLM(IPTRK,NBM,LL4,WORK2(1),EVECT(1,JGR),WORK1) + ELSEIF(CMOD.EQ.'TRIVAC')THEN + CALL KINTLM(IPTRK,NBM,LL4,WORK2(1),EVECT(1,JGR),WORK1) + ENDIF + DO 150 IND=1,LL4 + SRC(IND,IGR)=SRC(IND,IGR)-DSM*WORK1(IND) + 150 CONTINUE + 160 CONTINUE + 170 CONTINUE + 180 CONTINUE +*---- +* 1/V CONTRIBUTION +*---- + IF(CMOD.EQ.'BIVAC')THEN + CALL KINBLM(IPTRK,NBM,LL4,OVR(1,IGR),EVECT(1,IGR),WORK1) + ELSEIF(CMOD.EQ.'TRIVAC')THEN + CALL KINTLM(IPTRK,NBM,LL4,OVR(1,IGR),EVECT(1,IGR),WORK1) + ENDIF + DO 190 IND=1,LL4 + SRC(IND,IGR)=SRC(IND,IGR)+WORK1(IND) + 190 CONTINUE + 200 CONTINUE +*---- +* EDITION +*---- + IF(IMPX.GT.5) THEN + WRITE(IOS,1002) + DO 210 IGR=1,NGR + WRITE(IOS,1003) IGR,(SRC(IND,IGR),IND=1,LL4) + 210 CONTINUE + ENDIF +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(CHEXP,WORK2,WORK1) + RETURN +* + 1001 FORMAT(/1X,'COMPUTING THE SPACE-TIME KINETICS SOURCE VECTOR', + 1 1X,'ACCORDING TO THE TRACKING TYPE: ',A6/) + 1002 FORMAT(/1X,'=> COMPUTED SPACE-TIME KINETICS SOURCE VECTOR') + 1003 FORMAT(/15H NEUTRON GROUP=,I5/(1P,8D14.5)) + END |