summaryrefslogtreecommitdiff
path: root/Trivac/src/FLDSMB.f
diff options
context:
space:
mode:
Diffstat (limited to 'Trivac/src/FLDSMB.f')
-rwxr-xr-xTrivac/src/FLDSMB.f475
1 files changed, 475 insertions, 0 deletions
diff --git a/Trivac/src/FLDSMB.f b/Trivac/src/FLDSMB.f
new file mode 100755
index 0000000..f479a9a
--- /dev/null
+++ b/Trivac/src/FLDSMB.f
@@ -0,0 +1,475 @@
+*DECK FLDSMB
+ SUBROUTINE FLDSMB (IPTRK,IPSYS,IPFLUX,LL4,ITY,NUN,NGRP,ICL1,ICL2,
+ 1 IMPX,IMPH,TITR,EPS2,MAXOUT,MAXINR,EPSINR,EVECT,FKEFF)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Solution of a multigroup eigenvalue system for the calculation of the
+* direct neutron flux in BIVAC. Use the preconditionned power method
+* with a two-parameter SVAT acceleration technique.
+*
+*Copyright:
+* Copyright (C) 2002 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
+* IPTRK L_TRACK pointer to the BIVAC tracking information.
+* IPSYS L_SYSTEM pointer to system matrices.
+* IPFLUX L_FLUX pointer to the solution.
+* LL4 order of the system matrices.
+* ITY type of algorithm: 1: Diffusion theory; 11: Simplified PN
+* approximation.
+* NUN number of unknowns in each energy group.
+* NGRP number of energy groups.
+* ICL1 number of free iterations in one cycle of the inverse power
+* method.
+* ICL2 number of accelerated iterations in one cycle.
+* IMPX print parameter: =0: no print; =1: minimum printing;
+* =2: iteration history is printed.
+* IMPH type of histogram processing:
+* =0: no action is taken;
+* =1: the flux is compared to a reference flux stored on LCM;
+* =2: the convergence histogram is printed;
+* =3: the convergence histogram is printed with axis and
+* titles. The plotting file is completed;
+* =4: the convergence histogram is printed with axis, acce-
+* leration factors and titles. The plotting file is
+* completed.
+* TITR title.
+* EPS2 convergence criteria for the flux.
+* MAXOUT maximum number of outer iterations.
+* MAXINR maximum number of thermal iterations.
+* EPSINR thermal iteration epsilon.
+* EVECT initial estimate of the unknown vector.
+*
+*Parameters: output
+* EVECT converged unknown vector.
+* FKEFF effective multiplication factor.
+*
+*Reference:
+* A. H\'ebert, 'Preconditioning the power method for reactor
+* calculations', Nucl. Sci. Eng., 94, 1 (1986).
+*
+*-----------------------------------------------------------------------
+*
+ USE GANLIB
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ TYPE(C_PTR) IPTRK,IPSYS,IPFLUX
+ CHARACTER*72 TITR
+ INTEGER LL4,ITY,NUN,NGRP,ICL1,ICL2,IMPX,IMPH,MAXOUT,MAXINR
+ REAL FKEFF,EPS2,EPSINR,EVECT(NUN,NGRP)
+*----
+* LOCAL VARIABLES
+*----
+ PARAMETER (MMAXX=250,EPS1=1.0E-5)
+ CHARACTER*12 TEXT12
+ LOGICAL LOGTES
+ DOUBLE PRECISION AEAE,AEAG,AEAH,AGAG,AGAH,AHAH,BEBE,BEBG,BEBH,
+ 1 BGBG,BGBH,BHBH,AEBE,AEBG,AEBH,AGBE,AGBG,AGBH,AHBE,AHBG,AHBH,
+ 2 X,DXDA,DXDB,Y,DYDA,DYDB,Z,DZDA,DZDB,F,D2F(2,3),EVAL,ALP,BET,
+ 3 FMIN
+ INTEGER ITITR(18)
+ REAL ERR(MMAXX),ALPH(MMAXX),BETA(MMAXX)
+ DOUBLE PRECISION, PARAMETER :: ALP_TAB(24) = (/ 0.2, 0.4, 0.6,
+ 1 0.8, 1.0, 1.2, 1.5, 2.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0,
+ 2 40.0, 45.0, 50.0, 55.0, 60.0, 65.0, 70.0, 75.0, 80.0, 85.0 /)
+ DOUBLE PRECISION, PARAMETER :: BET_TAB(11) = (/ -1.0, -0.8, -0.6,
+ 1 -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0 /)
+ REAL, DIMENSION(:), ALLOCATABLE :: WORK
+ REAL, DIMENSION(:,:), ALLOCATABLE :: GRAD1,GRAD2,GAR1,GAR2,GAR3,
+ 1 GAF1,GAF2,GAF3
+*----
+* SCRATCH STORAGE ALLOCATION
+*----
+ ALLOCATE(GRAD1(NUN,NGRP),GRAD2(NUN,NGRP),GAR1(NUN,NGRP),
+ 1 GAR2(NUN,NGRP),GAR3(NUN,NGRP),GAF1(NUN,NGRP),GAF2(NUN,NGRP),
+ 2 GAF3(NUN,NGRP),WORK(NUN))
+*
+* TKT : CPU TIME FOR THE SOLUTION OF LINEAR SYSTEMS.
+* TKB : CPU TIME FOR BILINEAR PRODUCT EVALUATIONS.
+ TKT=0.0
+ TKB=0.0
+ CALL KDRCPU(TK1)
+*----
+* PRECONDITIONED POWER METHOD
+*----
+ EVAL=1.0
+ VVV=0.0
+ ISTART=1
+ TEST=0.0
+ IF(IMPX.GE.1) WRITE (6,600)
+ IF(IMPX.GE.2) WRITE (6,610)
+ DO 25 IGR=1,NGRP
+ WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR
+ CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,EVECT(1,IGR),GAR1(1,IGR))
+ DO 20 JGR=1,NGRP
+ IF(JGR.EQ.IGR) GO TO 20
+ WRITE(TEXT12,'(1HA,2I3.3)') IGR,JGR
+ CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+ IF(ILONG.EQ.0) GO TO 20
+ CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,EVECT(1,JGR),WORK(1))
+ DO 10 I=1,LL4
+ GAR1(I,IGR)=GAR1(I,IGR)-WORK(I)
+ 10 CONTINUE
+ 20 CONTINUE
+ 25 CONTINUE
+ CALL KDRCPU(TK2)
+ TKB=TKB+(TK2-TK1)
+*
+ M=0
+ 30 M=M+1
+*----
+* EIGENVALUE EVALUATION
+*----
+ CALL KDRCPU(TK1)
+ AEBE=0.0D0
+ BEBE=0.0D0
+ DO 75 IGR=1,NGRP
+ DO 40 I=1,LL4
+ GAF1(I,IGR)=0.0
+ 40 CONTINUE
+ DO 60 JGR=1,NGRP
+ WRITE(TEXT12,'(1HB,2I3.3)') IGR,JGR
+ CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+ IF(ILONG.EQ.0) GO TO 60
+ CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,EVECT(1,JGR),WORK(1))
+ DO 50 I=1,LL4
+ GAF1(I,IGR)=GAF1(I,IGR)+WORK(I)
+ 50 CONTINUE
+ 60 CONTINUE
+ DO 70 I=1,LL4
+ AEBE=AEBE+GAR1(I,IGR)*GAF1(I,IGR)
+ BEBE=BEBE+GAF1(I,IGR)**2
+ 70 CONTINUE
+ 75 CONTINUE
+ EVAL=AEBE/BEBE
+ CALL KDRCPU(TK2)
+ TKB=TKB+(TK2-TK1)
+*----
+* DIRECTION EVALUATION
+*----
+ DO 110 IGR=1,NGRP
+ CALL KDRCPU(TK1)
+ DO 80 I=1,LL4
+ GRAD1(I,IGR)=REAL(EVAL)*GAF1(I,IGR)-GAR1(I,IGR)
+ 80 CONTINUE
+ DO 100 JGR=1,IGR-1
+ WRITE(TEXT12,'(1HA,2I3.3)') IGR,JGR
+ CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+ IF(ILONG.EQ.0) GO TO 100
+ CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD1(1,JGR),WORK(1))
+ DO 90 I=1,LL4
+ GRAD1(I,IGR)=GRAD1(I,IGR)+WORK(I)
+ 90 CONTINUE
+ 100 CONTINUE
+ CALL KDRCPU(TK2)
+ TKB=TKB+(TK2-TK1)
+*
+ CALL KDRCPU(TK1)
+ WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR
+ CALL MTLDLS(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD1(1,IGR))
+ CALL KDRCPU(TK2)
+ TKT=TKT+(TK2-TK1)
+ 110 CONTINUE
+*----
+* PERFORM THERMAL (UP-SCATTERING) ITERATIONS
+*----
+ KTER=0
+ NADI=5 ! used with SPN approximations
+ IF(MAXINR.GT.1) THEN
+ CALL FLDBHR(IPTRK,IPSYS,.FALSE.,LL4,ITY,NUN,NGRP,ICL1,ICL2,
+ 1 IMPX,NADI,MAXINR,EPSINR,KTER,TKT,TKB,GRAD1)
+ ENDIF
+*----
+* DISPLACEMENT EVALUATION
+*----
+ F=0.0D0
+ DELS=ABS(REAL((EVAL-VVV)/EVAL))
+ VVV=REAL(EVAL)
+ CALL KDRCPU(TK1)
+*----
+* EVALUATION OF THE TWO ACCELERATION PARAMETERS ALP AND BET
+*----
+ ALP=1.0D0
+ BET=0.0D0
+ N=0
+ AEAE=0.0D0
+ AEAG=0.0D0
+ AEAH=0.0D0
+ AGAG=0.0D0
+ AGAH=0.0D0
+ AHAH=0.0D0
+ BEBG=0.0D0
+ BEBH=0.0D0
+ BGBG=0.0D0
+ BGBH=0.0D0
+ BHBH=0.0D0
+ AEBG=0.0D0
+ AEBH=0.0D0
+ AGBE=0.0D0
+ AGBG=0.0D0
+ AGBH=0.0D0
+ AHBE=0.0D0
+ AHBG=0.0D0
+ AHBH=0.0D0
+ DO 165 IGR=1,NGRP
+ WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR
+ CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD1(1,IGR),GAR2(1,IGR))
+ DO 120 I=1,LL4
+ GAF2(I,IGR)=0.0
+ 120 CONTINUE
+ DO 160 JGR=1,NGRP
+ IF(JGR.EQ.IGR) GO TO 140
+ WRITE(TEXT12,'(1HA,2I3.3)') IGR,JGR
+ CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+ IF(ILONG.EQ.0) GO TO 140
+ CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD1(1,JGR),WORK(1))
+ DO 130 I=1,LL4
+ GAR2(I,IGR)=GAR2(I,IGR)-WORK(I)
+ 130 CONTINUE
+ 140 WRITE(TEXT12,'(1HB,2I3.3)') IGR,JGR
+ CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+ IF(ILONG.EQ.0) GO TO 160
+ CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD1(1,JGR),WORK(1))
+ DO 150 I=1,LL4
+ GAF2(I,IGR)=GAF2(I,IGR)+WORK(I)
+ 150 CONTINUE
+ 160 CONTINUE
+ 165 CONTINUE
+ IF(1+MOD(M-ISTART,ICL1+ICL2).GT.ICL1) THEN
+ DO 175 IGR=1,NGRP
+ DO 170 I=1,LL4
+* COMPUTE (A ,A )
+ AEAE=AEAE+GAR1(I,IGR)**2
+ AEAG=AEAG+GAR1(I,IGR)*GAR2(I,IGR)
+ AEAH=AEAH+GAR1(I,IGR)*GAR3(I,IGR)
+ AGAG=AGAG+GAR2(I,IGR)**2
+ AGAH=AGAH+GAR2(I,IGR)*GAR3(I,IGR)
+ AHAH=AHAH+GAR3(I,IGR)**2
+* COMPUTE (B ,B )
+ BEBG=BEBG+GAF1(I,IGR)*GAF2(I,IGR)
+ BEBH=BEBH+GAF1(I,IGR)*GAF3(I,IGR)
+ BGBG=BGBG+GAF2(I,IGR)**2
+ BGBH=BGBH+GAF2(I,IGR)*GAF3(I,IGR)
+ BHBH=BHBH+GAF3(I,IGR)**2
+* COMPUTE (A ,B )
+ AEBG=AEBG+GAR1(I,IGR)*GAF2(I,IGR)
+ AEBH=AEBH+GAR1(I,IGR)*GAF3(I,IGR)
+ AGBE=AGBE+GAR2(I,IGR)*GAF1(I,IGR)
+ AGBG=AGBG+GAR2(I,IGR)*GAF2(I,IGR)
+ AGBH=AGBH+GAR2(I,IGR)*GAF3(I,IGR)
+ AHBE=AHBE+GAR3(I,IGR)*GAF1(I,IGR)
+ AHBG=AHBG+GAR3(I,IGR)*GAF2(I,IGR)
+ AHBH=AHBH+GAR3(I,IGR)*GAF3(I,IGR)
+ 170 CONTINUE
+ 175 CONTINUE
+*
+ 180 N=N+1
+ IF(N.GT.10) GO TO 185
+* COMPUTE X(M+1)
+ X=BEBE+ALP*ALP*BGBG+BET*BET*BHBH+2.0D0*(ALP*BEBG+BET*BEBH
+ 1 +ALP*BET*BGBH)
+ DXDA=2.0D0*(BEBG+ALP*BGBG+BET*BGBH)
+ DXDB=2.0D0*(BEBH+ALP*BGBH+BET*BHBH)
+* COMPUTE Y(M+1)
+ Y=AEAE+ALP*ALP*AGAG+BET*BET*AHAH+2.0D0*(ALP*AEAG+BET*AEAH
+ 1 +ALP*BET*AGAH)
+ DYDA=2.0D0*(AEAG+ALP*AGAG+BET*AGAH)
+ DYDB=2.0D0*(AEAH+ALP*AGAH+BET*AHAH)
+* COMPUTE Z(M+1)
+ Z=AEBE+ALP*ALP*AGBG+BET*BET*AHBH+ALP*(AEBG+AGBE)
+ 1 +BET*(AEBH+AHBE)+ALP*BET*(AGBH+AHBG)
+ DZDA=AEBG+AGBE+2.0D0*ALP*AGBG+BET*(AGBH+AHBG)
+ DZDB=AEBH+AHBE+ALP*(AGBH+AHBG)+2.0D0*BET*AHBH
+* COMPUTE F(M+1)
+ F=X*Y-Z*Z
+ D2F(1,1)=2.0D0*(BGBG*Y+DXDA*DYDA+X*AGAG-DZDA**2-2.0D0*Z*AGBG)
+ D2F(1,2)=2.0D0*BGBH*Y+DXDA*DYDB+DXDB*DYDA+2.0D0*X*AGAH
+ 1 -2.0D0*DZDA*DZDB-2.0D0*Z*(AGBH+AHBG)
+ D2F(2,2)=2.0D0*(BHBH*Y+DXDB*DYDB+X*AHAH-DZDB**2-2.0D0*Z*AHBH)
+ D2F(2,1)=D2F(1,2)
+ D2F(1,3)=DXDA*Y+X*DYDA-2.0D0*Z*DZDA
+ D2F(2,3)=DXDB*Y+X*DYDB-2.0D0*Z*DZDB
+* SOLUTION OF A LINEAR SYSTEM.
+ CALL ALSBD(2,1,D2F,IER,2)
+ IF(IER.NE.0) GO TO 185
+ ALP=ALP-D2F(1,3)
+ BET=BET-D2F(2,3)
+ IF(ALP.GT.100.0) GO TO 185
+ IF((ABS(D2F(1,3)).LE.1.0D-4).AND.(ABS(D2F(2,3)).LE.1.0D-4))
+ 1 GO TO 190
+ GO TO 180
+*
+* alternative algorithm in case of Newton-Raphton failure
+ 185 IF(IMPX.GT.0) WRITE(6,'(/30H FLDSMB: FAILURE OF THE NEWTON,
+ 1 55H-RAPHTON ALGORIHTHM FOR COMPUTING THE OVERRELAXATION PA,
+ 2 9HRAMETERS.)')
+ IAMIN=999
+ IBMIN=999
+ FMIN=HUGE(FMIN)
+ DO IA=1,SIZE(ALP_TAB)
+ ALP=ALP_TAB(IA)
+ DO IB=1,SIZE(BET_TAB)
+ BET=BET_TAB(IB)
+* COMPUTE X
+ X=BEBE+ALP*ALP*BGBG+BET*BET*BHBH+2.0D0*(ALP*BEBG+BET*BEBH
+ 1 +ALP*BET*BGBH)
+* COMPUTE Y
+ Y=AEAE+ALP*ALP*AGAG+BET*BET*AHAH+2.0D0*(ALP*AEAG+BET*AEAH
+ 1 +ALP*BET*AGAH)
+* COMPUTE Z
+ Z=AEBE+ALP*ALP*AGBG+BET*BET*AHBH+ALP*(AEBG+AGBE)
+ 1 +BET*(AEBH+AHBE)+ALP*BET*(AGBH+AHBG)
+* COMPUTE F
+ F=X*Y-Z*Z
+ IF(F.LT.FMIN) THEN
+ IAMIN=IA
+ IBMIN=IB
+ FMIN=F
+ ENDIF
+ ENDDO
+ ENDDO
+ ALP=ALP_TAB(IAMIN)
+ BET=BET_TAB(IBMIN)
+ 190 BET=BET/ALP
+ IF((ALP.LT.1.0D0).AND.(ALP.GT.0.0D0)) THEN
+ ALP=1.0D0
+ BET=0.0D0
+ ELSE IF(ALP.LE.0.0D0) THEN
+ ISTART=M+1
+ ALP=1.0D0
+ BET=0.0D0
+ ENDIF
+ DO 205 IGR=1,NGRP
+ DO 200 I=1,LL4
+ GRAD1(I,IGR)=REAL(ALP)*(GRAD1(I,IGR)+REAL(BET)*GRAD2(I,IGR))
+ GAR2(I,IGR)=REAL(ALP)*(GAR2(I,IGR)+REAL(BET)*GAR3(I,IGR))
+ GAF2(I,IGR)=REAL(ALP)*(GAF2(I,IGR)+REAL(BET)*GAF3(I,IGR))
+ 200 CONTINUE
+ 205 CONTINUE
+ ENDIF
+ CALL KDRCPU(TK2)
+ TKB=TKB+(TK2-TK1)
+*
+ LOGTES=(M.LT.ICL1).OR.(MOD(M-ISTART,ICL1+ICL2).EQ.ICL1-1)
+ IF(LOGTES.AND.(DELS.LE.EPS1)) THEN
+ DELT=0.0
+ DO 220 IGR=1,NGRP
+ DELN=0.0
+ DELD=0.0
+ DO 210 I=1,LL4
+ EVECT(I,IGR)=EVECT(I,IGR)+GRAD1(I,IGR)
+ GAR1(I,IGR)=GAR1(I,IGR)+GAR2(I,IGR)
+ GAF1(I,IGR)=GAF1(I,IGR)+GAF2(I,IGR)
+ GRAD2(I,IGR)=GRAD1(I,IGR)
+ GAR3(I,IGR)=GAR2(I,IGR)
+ GAF3(I,IGR)=GAF2(I,IGR)
+ DELN=MAX(DELN,ABS(GAF2(I,IGR)))
+ DELD=MAX(DELD,ABS(GAF1(I,IGR)))
+ 210 CONTINUE
+ IF(DELD.NE.0.0) DELT=MAX(DELT,DELN/DELD)
+ 220 CONTINUE
+ IF(IMPX.GE.2) WRITE (6,615) M,AEAE,AEAG,AEAH,AGAG,AGAH,AHAH,
+ 1 BEBE,ALP,BET,EVAL,F,DELS,DELT,N,BEBG,BEBH,BGBG,BGBH,BHBH,
+ 2 AEBE,AEBG,AEBH,AGBE,AGBG,AGBH,AHBE,AHBG,AHBH
+* COMPUTE THE CONVERGENCE HISTOGRAM.
+ IF((IMPH.GE.1).AND.(M.LE.MMAXX)) THEN
+ CALL FLDXCO(IPFLUX,LL4,NUN,EVECT(1,NGRP),.TRUE.,ERR(M))
+ ALPH(M)=REAL(ALP)
+ BETA(M)=REAL(BET)
+ ENDIF
+ IF(DELT.LE.EPS2) GO TO 240
+ ELSE
+ DO 235 IGR=1,NGRP
+ DO 230 I=1,LL4
+ EVECT(I,IGR)=EVECT(I,IGR)+GRAD1(I,IGR)
+ GAR1(I,IGR)=GAR1(I,IGR)+GAR2(I,IGR)
+ GAF1(I,IGR)=GAF1(I,IGR)+GAF2(I,IGR)
+ GRAD2(I,IGR)=GRAD1(I,IGR)
+ GAR3(I,IGR)=GAR2(I,IGR)
+ GAF3(I,IGR)=GAF2(I,IGR)
+ 230 CONTINUE
+ 235 CONTINUE
+ IF(IMPX.GE.2) WRITE (6,620) M,AEAE,AEAG,AEAH,AGAG,AGAH,AHAH,
+ 1 BEBE,ALP,BET,EVAL,F,DELS,N,BEBG,BEBH,BGBG,BGBH,BHBH,AEBE,
+ 2 AEBG,AEBH,AGBE,AGBG,AGBH,AHBE,AHBG,AHBH
+* COMPUTE THE CONVERGENCE HISTOGRAM.
+ IF((IMPH.GE.1).AND.(M.LE.MMAXX)) THEN
+ CALL FLDXCO(IPFLUX,LL4,NUN,EVECT(1,NGRP),.TRUE.,ERR(M))
+ ALPH(M)=REAL(ALP)
+ BETA(M)=REAL(BET)
+ ENDIF
+ ENDIF
+ IF(M.EQ.1) TEST=DELS
+ IF((M.GT.5).AND.(DELS.GT.TEST)) CALL XABORT('FLDSMB: CONVERGENCE'
+ 1 //' FAILURE.')
+ IF(M.GE.MAXOUT) CALL XABORT('FLDSMB: MAXIMUM NUMBER OF ITERATION'
+ 1 //'S REACHED.')
+ GO TO 30
+*----
+* SOLUTION EDITION
+*----
+ 240 FKEFF=REAL(1.0D0/EVAL)
+ IF(IMPX.EQ.1) WRITE (6,640) M
+ IF(IMPX.GE.1) THEN
+ WRITE (6,650) TKT,TKB,TKT+TKB
+ WRITE (6,670) FKEFF
+ ENDIF
+ IF(IMPX.EQ.3) THEN
+ DO 250 IGR=1,NGRP
+ WRITE (6,680) IGR,(EVECT(I,IGR),I=1,LL4)
+ 250 CONTINUE
+ ENDIF
+ IF(IMPH.EQ.1) THEN
+ CALL LCMLEN(IPFLUX,'REF',ILONG,ITYLCM)
+ IF(ILONG.EQ.0) THEN
+ WRITE(6,'(40H FLDSMB: STORE A REFERENCE THERMAL FLUX.)')
+ CALL LCMPUT(IPFLUX,'REF',NUN,2,EVECT(1,NGRP))
+ ENDIF
+ ELSE IF(IMPH.GE.2) THEN
+ IGRAPH=0
+ 260 IGRAPH=IGRAPH+1
+ WRITE (TEXT12,'(5HHISTO,I3)') IGRAPH
+ CALL LCMLEN (IPFLUX,TEXT12,ILENG,ITYLCM)
+ IF(ILENG.EQ.0) THEN
+ MM=MIN(M,MMAXX)
+ READ (TITR,'(18A4)') ITITR
+ CALL LCMSIX (IPFLUX,TEXT12,1)
+ CALL LCMPUT (IPFLUX,'HTITLE',18,3,ITITR)
+ CALL LCMPUT (IPFLUX,'ALPHA',MM,2,ALPH)
+ CALL LCMPUT (IPFLUX,'BETA',MM,2,BETA)
+ CALL LCMPUT (IPFLUX,'ERROR',MM,2,ERR)
+ CALL LCMPUT (IPFLUX,'IMPH',1,1,IMPH)
+ CALL LCMSIX (IPFLUX,' ',2)
+ ELSE
+ GO TO 260
+ ENDIF
+ ENDIF
+*----
+* SCRATCH STORAGE DEALLOCATION
+*----
+ DEALLOCATE(GRAD1,GRAD2,GAR1,GAR2,GAR3,GAF1,GAF2,GAF3,WORK)
+ RETURN
+*
+ 600 FORMAT(1H1/50H FLDSMB: ITERATIVE PROCEDURE BASED ON PRECONDITION,
+ 1 16HED POWER METHOD./9X,16HDIRECT EQUATION.)
+ 610 FORMAT(//5X,17HBILINEAR PRODUCTS,48X,5HALPHA,3X,4HBETA,3X,
+ 1 12HEIGENVALUE..,12X,8HACCURACY,11(1H.),2X,1HN)
+ 615 FORMAT(1X,I3,1P,7E9.1,0P,2F8.3,E14.6,3E10.2,I4/(4X,1P,7E9.1))
+ 620 FORMAT(1X,I3,1P,7E9.1,0P,2F8.3,E14.6,2E10.2,10X,I4/(4X,1P,7E9.1))
+ 640 FORMAT(/23H FLDSMB: CONVERGENCE IN,I4,12H ITERATIONS.)
+ 650 FORMAT(/53H FLDSMB: CPU TIME USED TO SOLVE THE TRIANGULAR LINEAR,
+ 1 10H SYSTEMS =,F10.3/23X,34HTO COMPUTE THE BILINEAR PRODUCTS =,
+ 2 F10.3,20X,16HTOTAL CPU TIME =,F10.3)
+ 670 FORMAT(//42H FLDSMB: EFFECTIVE MULTIPLICATION FACTOR =,1P,E17.10/)
+ 680 FORMAT(//47H FLDSMB: EIGENVECTOR CORRESPONDING TO THE GROUP,I4
+ 1 //(5X,1P,8E14.5))
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 11:26:31 +0000