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/SPHEQU.f | 757 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 757 insertions(+) create mode 100644 Dragon/src/SPHEQU.f (limited to 'Dragon/src/SPHEQU.f') diff --git a/Dragon/src/SPHEQU.f b/Dragon/src/SPHEQU.f new file mode 100644 index 0000000..b4b42dc --- /dev/null +++ b/Dragon/src/SPHEQU.f @@ -0,0 +1,757 @@ +*DECK SPHEQU + SUBROUTINE SPHEQU(NBMIX2,IPTRK2,IFTRAK,IPMACR,IPFLX,CDOOR,NSPH, + 1 KSPH,MAXIT,MAXNBI,EPSPH,IPRINT,IMC,NGCOND,NMERGE,NALBP,ISCAT, + 2 NREG2,NUN2,MAT2,VOL2,KEY2,MERG2,ILK,CTITRE,IGRMIN,IGRMAX,SPH) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Calculation of the SPH factors for the homogenization of any geometry +* using a transport-transport or transport-diffusion equivalence +* technique. The macro-calculation can be performed using a standard +* solution door of Dragon. +* +*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 +* NBMIX2 number of macro-mixtures. Equal to MAX(MAT2(IREG)) for +* (IREG .le. NREG2). +* IPTRK2 pointer to the tracking of the macro-geometry (L_TRACK +* signature). +* IFTRAK unit number of the sequential binary tracking file. +* IPMACR pointer to the reference macrolib (L_MACROLIB signature). +* IPFLX pointer towards an initialization flux (L_FLUX signature). +* CDOOR tracking operator used to track the macro-geometry. +* NSPH type of SPH algorithm: +* =2 homogeneous macro-calculation (non-iterative procedure +* or Hebert-Benoist SPH-5 procedure); +* =3 any type of pij macro-calculation; +* =4 any type of diffusion, SN, PN or SPN macro-calculation; +* KSPH type of SPH factor normalization: +* <0 asymptotic normatization with respect to mixture -KSPH; +* =1 average flux normalization; +* =2 Selengut normalization using ALBS information; +* =3 Selengut normalization using FD_B boundary fluxes; +* =4 generalized Selengut normalization (EDF type); +* =5 Selengut normalization with surface leakage; +* =6 Selengut with water gap normalization; +* =7 average flux normalization in fissile zones. +* The Hebert-Benoist procedure is used if NSPH=2 and KSPH=5. +* MAXIT maximum number of SPH iterations. +* MAXNBI acceptable number of SPH iterations with an increase in +* convergence error before aborting. +* EPSPH convergence criterion for stopping the SPH iterations. +* IPRINT print flag (equal to 0 for no print). +* IMC type of macro-calculation (=1 diffusion or SPN; +* =2 other options; +* =3 type PIJ with Bell acceleration). +* NGCOND number of condensed groups. +* NMERGE number of merged regions (equal to 1 or to the number of +* different flux components in the macro-calculation). +* NALBP number of physical albedos. +* ISCAT scattering anisotropy in the reference set of cross sections +* (=1 isotropic in LAB; =2 linearly-anisotropic in LAB). +* NREG2 number of macro-regions (in the macro-calculation). +* NUN2 number of unknowns in a one-group macro-calculation. +* MAT2 mixture index per macro-region. +* VOL2 volume of macro-regions. +* KEY2 pointer to flux values in unknown vector. +* MERG2 index of merged macro-regions per macro-mixture. +* ILK leakage switch. +* CTITRE title. +* IGRMIN first group to process. +* IGRMAX last group to process. +* +*Parameters: output +* SPH SPH homogenization factors. +* +*References(s): +* A. Hebert, 'A Consistent Technique for the Pin-by-Pin Homogenization +* of a Pressurized Water Reactor Assembly', Nucl. Sci. Eng., 113, 227 +* (1993). +* +* A. Hebert and P. Benoist, 'A Consistent Technique for the Global +* Homogenization of a Pressurized Water Reactor Assembly', Nucl. Sci. +* Eng., 109, 360 (1991). +* +* A. Hebert and G. Mathonniere, 'Development of a Third-Generation +* Superhomogeneisation Method for the Homogenization of a +* Pressurized Water Reactor Assembly', Nucl. Sci. Eng., 115, 129 +* (1993). +* +* T. Courau, M. Cometto, E. Girardi, D. Couyras and N. Schwartz, +* 'Elements of Validation of Pin-by-Pin Calculations with the Future +* EDF Calculation Scheme Based on APOLLO2 and COCAGNE Codes', +* Proceedings of ICAPP '08 Anaheim, CA USA, June 8-12, 2008. +* +*----------------------------------------------------------------------- +* + USE GANLIB + USE DOORS_MOD +*---- +* SUBROUTINE ARGUMENTS +*---- + TYPE(C_PTR) IPTRK2,IPMACR,IPFLX + INTEGER NBMIX2,IFTRAK,NSPH,KSPH,MAXIT,MAXNBI,IPRINT,IMC,NGCOND, + 1 NMERGE,NALBP,ISCAT,NREG2,NUN2,MAT2(NREG2),KEY2(NREG2), + 2 MERG2(NBMIX2),IGRMIN,IGRMAX + REAL VOL2(NREG2),SPH(NMERGE+NALBP,NGCOND) + LOGICAL ILK + CHARACTER CDOOR*12,CTITRE*72 +*---- +* LOCAL VARIABLES +*---- + PARAMETER(NSTATE=40) + INTEGER IPAR(NSTATE) + LOGICAL LNORM,LEXAC,LDIFF,REBFLG + DOUBLE PRECISION FLXTOT,VOLTOT + CHARACTER TEXT12*12,HSMG*131 + TYPE(C_PTR) IPADF,IPSYS2,JPSYS2,JPFLX,KPSYS2,IPSOU +*---- +* ALLOCATABLE ARRAYS +*---- + INTEGER, ALLOCATABLE, DIMENSION(:) :: NPSYS + REAL, ALLOCATABLE, DIMENSION(:) :: SIGMD,VOLMER,FACTOR,OUTG1, + 1 OUTG2,COURIN,COUROW,SNORM,SPHNEW,SIGG + REAL, ALLOCATABLE, DIMENSION(:,:) :: SIGMA,SIGMS,DIFF,FUNKNO, + 1 SUNKNO,FLXMER,ZLEAK,ALB1,ALB2,WORK + REAL, ALLOCATABLE, DIMENSION(:,:,:) :: SIGT,SIGW + REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: SUNMER + LOGICAL, ALLOCATABLE, DIMENSION(:) :: LFISS +*---- +* SCRATCH STORAGE ALLOCATION +*---- + ALLOCATE(NPSYS(NGCOND)) + ALLOCATE(SIGMA(0:NBMIX2,ISCAT+1),SIGMD(0:NBMIX2),COURIN(NGCOND), + 1 COUROW(NGCOND),VOLMER(NMERGE),SIGW(NMERGE,NGCOND,ISCAT+1), + 2 FLXMER(NMERGE,NGCOND),DIFF(NMERGE,NGCOND),ZLEAK(NMERGE,NGCOND), + 3 SUNMER(NMERGE,NGCOND,NGCOND,ISCAT),FUNKNO(NUN2,NGCOND), + 4 SUNKNO(NUN2,NGCOND),FACTOR(NMERGE),OUTG1(NGCOND),OUTG2(NGCOND), + 5 SNORM(NGCOND),ALB1(NALBP,NGCOND),ALB2(NALBP,NGCOND), + 6 LFISS(NMERGE)) +*---- +* CALCULATION OF THE REFERENCE MERGED/CONDENSED SET OF CROSS SECTIONS +*---- + IF(NALBP.GT.1) CALL XABORT('SPHEQU: NALBP<=1 EXPECTED.') + CALL LCMGET(IPMACR,'STATE-VECTOR',IPAR) + NL=IPAR(3) + NIFISS=IPAR(4) + ILEAKS=IPAR(9) + NW=MAX(1,IPAR(10)) + IF(IPAR(2).NE.NMERGE) THEN + CALL XABORT('SPHEQU: INVALID VALUE OF NMERGE.') + ELSE IF(IPAR(8).NE.NALBP) THEN + CALL XABORT('SPHEQU: INVALID VALUE OF NALBP.') + ENDIF + ALLOCATE(SIGT(NMERGE,NGCOND,NW+1),SIGMS(0:NBMIX2,NW+1)) + CALL SPHMAC(IPMACR,IPRINT,NMERGE,NALBP,NGCOND,ISCAT,NW,NIFISS, + 1 ILEAKS,VOLMER,FLXMER,SUNMER,SIGT,SIGW,DIFF,ZLEAK,OUTG2,ALB2, + 2 LFISS) +* + DO 30 INL=1,ISCAT + DO 20 IGR=1,NGCOND + DO 10 IKK=1,NMERGE + SUNMER(IKK,IGR,IGR,INL)=SUNMER(IKK,IGR,IGR,INL)-SIGW(IKK,IGR,INL) + 10 CONTINUE + 20 CONTINUE + 30 CONTINUE + LDIFF=.FALSE. + NLF=0 + NANI=0 + IF((NSPH.EQ.4).AND.((CDOOR.EQ.'BIVAC').OR.(CDOOR.EQ.'TRIVAC'))) + > THEN +* TRANSPORT-DIFFUSION EQUIVALENCE. + IF(.NOT.C_ASSOCIATED(IPTRK2)) THEN + CALL XABORT('SPHEQU: MACRO-TRACKING NOT DEFINED(1).') + ENDIF + CALL LCMGET(IPTRK2,'STATE-VECTOR',IPAR) + IF(CDOOR.EQ.'BIVAC') THEN + NLF=IPAR(14) + NANI=IPAR(16) + ELSE IF(CDOOR.EQ.'TRIVAC') THEN + NLF=IPAR(30) + NANI=IPAR(32) + ENDIF + LDIFF=(NLF.EQ.0).OR.(NANI.LT.0) + NANI=ABS(NANI) + IF(NANI.GT.ISCAT) CALL XABORT('SPHEQU: ISCAT OVERFLOW.') + IF(LDIFF) THEN + IF(ILEAKS.EQ.0) CALL XABORT('SPHEQU: UNABLE TO COMPUTE DIFF' + > //'USION COEFFICIENTS.') + ENDIF + ENDIF + IF(NLF.EQ.0) NANI=1 +*---- +* RECOVER THE AVERAGED GAP AND ROW FLUXES FOR EDF-TYPE NORMALIZATION. +*---- + IF(KSPH.GE.2) THEN + CALL LCMLEN(IPMACR,'ADF',ILONG,ITYLCM) + IF(ILONG.EQ.0) THEN + CALL LCMLIB(IPMACR) + CALL XABORT('SPHEQU: NO ADF DIRECTORY IN THE MACROLIB.') + ENDIF + IPADF=LCMDID(IPMACR,'ADF') + ENDIF + IF((KSPH.EQ.2).OR.(KSPH.EQ.5)) THEN + CALL LCMLEN(IPADF,'ALBS00',ILONG,ITYLCM) + IF(ILONG.NE.2*NGCOND) THEN + WRITE(HSMG,'(26HSPHEQU: BAD ALBS00 LENGTH=,I5,10H EXPECTED=, + > I5,1H.)') ILONG,2*NGCOND + CALL XABORT(HSMG) + ENDIF + ALLOCATE(WORK(NGCOND,2)) + CALL LCMGET(IPADF,'ALBS00',WORK) + COURIN(:NGCOND)=WORK(:NGCOND,1) + DEALLOCATE(WORK) + IF(IPRINT.GT.3) THEN + WRITE(6,'(/45H SPHEQU: THE VALUES OF ALBS PER MACRO-GROUPS , + > 3HARE)') + WRITE(6,'(1X,1P,10E13.5)') (COURIN(IGR),IGR=1,NGCOND) + ENDIF + ELSE IF((KSPH.EQ.3).OR.(KSPH.EQ.6)) THEN + CALL LCMLEN(IPADF,'FD_B',ILONG,ITYLCM) + IF(ILONG.NE.NMERGE*NGCOND) THEN + CALL LCMLIB(IPADF) + WRITE(HSMG,'(24HSPHEQU: BAD FD_B LENGTH=,I5,10H EXPECTED=, + > I5,5H (1).)') ILONG,NMERGE*NGCOND + CALL XABORT(HSMG) + ENDIF + ALLOCATE(WORK(NMERGE,NGCOND)) + CALL LCMGET(IPADF,'FD_B',WORK) + DO IGR=1,NGCOND + VOLTOT=0.0D0 + FLXTOT=0.0D0 + DO IKK=1,NMERGE + VOLTOT=VOLTOT+VOLMER(IKK) + FLXTOT=FLXTOT+WORK(IKK,IGR)*VOLMER(IKK) + ENDDO + COURIN(IGR)=REAL(FLXTOT/VOLTOT) + ENDDO + DEALLOCATE(WORK) + IF(IPRINT.GT.3) THEN + WRITE(6,'(/45H SPHEQU: THE VALUES OF FD_B PER MACRO-GROUPS , + > 3HARE)') + WRITE(6,'(1X,1P,10E13.5)') (COURIN(IGR),IGR=1,NGCOND) + ENDIF + ELSE IF(KSPH.EQ.4) THEN + CALL LCMLEN(IPADF,'FD_B',ILONG,ITYLCM) + IF(ILONG.NE.NMERGE*NGCOND) THEN + WRITE(HSMG,'(24HSPHEQU: BAD FD_B LENGTH=,I5,10H EXPECTED=, + > I5,5H (2).)') ILONG,NMERGE*NGCOND + CALL XABORT(HSMG) + ENDIF + ALLOCATE(WORK(NMERGE,NGCOND)) + CALL LCMGET(IPADF,'FD_B',WORK) + DO IGR=1,NGCOND + VOLTOT=0.0D0 + FLXTOT=0.0D0 + DO IKK=1,NMERGE + VOLTOT=VOLTOT+VOLMER(IKK) + FLXTOT=FLXTOT+WORK(IKK,IGR)*VOLMER(IKK) + ENDDO + COURIN(IGR)=REAL(FLXTOT/VOLTOT) + ENDDO + CALL LCMLEN(IPADF,'FD_H',ILONG,ITYLCM) + IF(ILONG.NE.NMERGE*NGCOND) THEN + WRITE(HSMG,'(24HSPHEQU: BAD FD_H LENGTH=,I5,10H EXPECTED=, + > I5,1H.)') ILONG,NMERGE*NGCOND + CALL XABORT(HSMG) + ENDIF + CALL LCMGET(IPADF,'FD_H',WORK) + DO IGR=1,NGCOND + VOLTOT=0.0D0 + FLXTOT=0.0D0 + DO IKK=1,NMERGE + VOLTOT=VOLTOT+VOLMER(IKK) + FLXTOT=FLXTOT+WORK(IKK,IGR)*VOLMER(IKK) + ENDDO + COUROW(IGR)=REAL(FLXTOT/VOLTOT) + ENDDO + DEALLOCATE(WORK) + IF(IPRINT.GT.3) THEN + WRITE(6,'(/45H SPHEQU: THE VALUES OF FD_B PER MACRO-GROUPS , + > 3HARE)') + WRITE(6,'(1X,1P,10E13.5)') (COURIN(IGR),IGR=1,NGCOND) + WRITE(6,'(/45H SPHEQU: THE VALUES OF FD_H PER MACRO-GROUPS , + > 3HARE)') + WRITE(6,'(1X,1P,10E13.5)') (COUROW(IGR),IGR=1,NGCOND) + ENDIF + ENDIF +*---- +* ITERATIVE STRATEGY USED TO COMPUTE THE SPH FACTORS +*---- + IPRIN2=MAX(0,IPRINT-5) + CALL LCMOP(IPSYS2,'SPH$SYS',0,1,0) + JPSYS2=LCMLID(IPSYS2,'GROUP',NGCOND) +*---- +* REMOVE DB2 LEAKAGE FROM SOURCES +*---- + IF(ILEAKS.NE.0) THEN + DO 50 IGR=1,NGCOND + DO 40 IKK=1,NMERGE + SUNMER(IKK,IGR,IGR,1)=SUNMER(IKK,IGR,IGR,1)-ZLEAK(IKK,IGR) + 40 CONTINUE + 50 CONTINUE + ENDIF +*---- +* SET SPH NORMALIZATION INFORMATION +*---- + IF(KSPH.LT.0) THEN +* ASYMPTOTIC NORMALIZATION WITH RESPECT TO MIXTURE -KSPH. + DO 70 IGR=1,NGCOND + VOLTOT=0.0D0 + FLXTOT=0.0D0 + DO 60 IKK=1,NMERGE + IF(IKK.EQ.-KSPH) THEN + VOLTOT=VOLTOT+VOLMER(IKK) + FLXTOT=FLXTOT+FLXMER(IKK,IGR)*VOLMER(IKK) + ENDIF + 60 CONTINUE + IF(VOLTOT.EQ.0.0) CALL XABORT('SPHEQU: ASYMPTOTIC NORMALIZATI' + > //'ON FAILURE.') + COURIN(IGR)=REAL(FLXTOT/VOLTOT) + SPH(:NMERGE,IGR)=1.0 + 70 CONTINUE + IF(IPRINT.GT.3) THEN + WRITE(6,910) (COURIN(IGR),IGR=1,NGCOND) + WRITE(6,'(/)') + ENDIF + ELSE IF(KSPH.EQ.1) THEN +* AVERAGE FLUX NORMALIZATION. + DO 90 IGR=1,NGCOND + VOLTOT=0.0D0 + FLXTOT=0.0D0 + DO 80 IKK=1,NMERGE + VOLTOT=VOLTOT+VOLMER(IKK) + FLXTOT=FLXTOT+FLXMER(IKK,IGR)*VOLMER(IKK) + 80 CONTINUE + COURIN(IGR)=REAL(FLXTOT/VOLTOT) + 90 CONTINUE + IF(IPRINT.GT.3) THEN + WRITE(6,910) (COURIN(IGR),IGR=1,NGCOND) + WRITE(6,'(/)') + ENDIF + ELSE IF((KSPH.EQ.2).OR.(KSPH.EQ.3).OR.(KSPH.EQ.5)) THEN +* SELENGUT NORMALIZATION. + DO 120 IGR=1,NGCOND + VOLTOT=0.0D0 + FLXTOT=0.0D0 + DO 100 IKK=1,NMERGE + VOLTOT=VOLTOT+VOLMER(IKK) + FLXTOT=FLXTOT+FLXMER(IKK,IGR)*VOLMER(IKK) + 100 CONTINUE + FLXTOT=FLXTOT/VOLTOT + DO 110 IKK=1,NMERGE + SPH(IKK,IGR)=REAL(FLXTOT)/COURIN(IGR) + 110 CONTINUE + 120 CONTINUE + ELSE IF(KSPH.EQ.4) THEN +* GENERALIZED SELENGUT NORMALIZATION (EDF-TYPE). + DO 150 IGR=1,NGCOND + VOLTOT=0.0D0 + FLXTOT=0.0D0 + DO 130 IKK=1,NMERGE + VOLTOT=VOLTOT+VOLMER(IKK) + FLXTOT=FLXTOT+FLXMER(IKK,IGR)*VOLMER(IKK) + 130 CONTINUE + FLXTOT=FLXTOT/VOLTOT + DO 140 IKK=1,NMERGE + SPH(IKK,IGR)=COUROW(IGR)/COURIN(IGR) + 140 CONTINUE + COURIN(IGR)=COURIN(IGR)*REAL(FLXTOT)/COUROW(IGR) + 150 CONTINUE + ELSE IF(KSPH.EQ.6) THEN +* MACRO-CALCULATION WATER GAP NORMALIZATION. +* COUROW = flux in gap in diffusion, initialized to flux in gap in transport +* FUNKNO = FLXMER : flux in diffusion, initialized to flux in transport + DO 160 IGR=1,NGCOND + COUROW(IGR)=COURIN(IGR) + SPH(:NMERGE,IGR)=1.0 + 160 CONTINUE + ELSE IF(KSPH.EQ.7) THEN +* AVERAGE FLUX NORMALIZATION IN FISSILE ZONES. + DO 180 IGR=1,NGCOND + VOLTOT=0.0D0 + FLXTOT=0.0D0 + DO 170 IKK=1,NMERGE + IF(LFISS(IKK)) THEN + VOLTOT=VOLTOT+VOLMER(IKK) + FLXTOT=FLXTOT+FLXMER(IKK,IGR)*VOLMER(IKK) + ENDIF + 170 CONTINUE + COURIN(IGR)=REAL(FLXTOT/VOLTOT) + 180 CONTINUE + IF(IPRINT.GT.3) THEN + WRITE(6,910) (COURIN(IGR),IGR=1,NGCOND) + WRITE(6,'(/)') + ENDIF + ENDIF + IF((NSPH.EQ.2).AND.(KSPH.NE.5)) GO TO 470 +*---- +* SPH ITERATIONS. +*---- + ITER=0 + IF(C_ASSOCIATED(IPFLX)) THEN + CALL LCMGET(IPFLX,'STATE-VECTOR',IPAR) + IF(IPAR(1).NE.NGCOND) CALL XABORT('SPHEQU: INVALID NB OF GROUP' + 1 //'S IN THE INITIALIZATION FLUX.') + IF(IPAR(2).NE.NUN2) CALL XABORT('SPHEQU: INVALID NB OF UNKNOWN' + 1 //'S IN THE INITIALIZATION FLUX.') + JPFLX=LCMGID(IPFLX,'FLUX') + DO 190 IGR=1,NGCOND + CALL LCMGDL(JPFLX,IGR,FUNKNO(1,IGR)) + 190 CONTINUE + ITER=1 + ELSE IF((CDOOR.EQ.'BIVAC').OR.(CDOOR.EQ.'TRIVAC')) THEN + DO 200 IGR=1,NGCOND + FUNKNO(:NUN2,IGR)=1.0 + 200 CONTINUE + ELSE + DO 220 IGR=1,NGCOND + FUNKNO(:NUN2,IGR)=0.0 + DO 210 IREG=1,NREG2 + IMAT=MAT2(IREG) + IF(IMAT.GT.0) FUNKNO(KEY2(IREG),IGR)=FLXMER(MERG2(IMAT),IGR) + 210 CONTINUE + 220 CONTINUE + ENDIF + OLDERR=1.0 + NBIERR=0 + 230 ITER=ITER+1 + IF(ITER.GE.MAXIT) THEN + WRITE(6,'(/46H SPHEQU: MAX. NUMBER OF ITERATIONS IS REACHED.)') + GO TO 440 + ENDIF + ERROR=0.0 + ERR2=0.0 + DO 240 IREG=1,NREG2 + IF(MAT2(IREG).GT.NBMIX2) THEN + CALL XABORT('SPHEQU: INVALID MACRO-MIXTURE INDEX.') + ENDIF + 240 CONTINUE +*---- +* SET MACROSCOPIC CROSS SECTIONS IN THE IPSYS2 LCM OBJECT. +*---- + NPSYS(:NGCOND)=0 + DO 310 IGR=IGRMIN,IGRMAX + SIGMS(0:NBMIX2,:NW+1)=0.0 + SIGMA(0:NBMIX2,:ISCAT+1)=0.0 + SIGMD(0:NBMIX2)=0.0 + DO 280 IREG=1,NREG2 + IMAT=MAT2(IREG) + IF(IMAT.EQ.0) GO TO 280 + IMERG=MERG2(IMAT) + IF(LDIFF) SIGMD(IMAT)=DIFF(IMERG,IGR)*SPH(IMERG,IGR) + IF(IMC.EQ.1) THEN + SIGMA(IMAT,1)=SIGW(IMERG,IGR,1)*SPH(IMERG,IGR) + DO 250 IW=1,NW+1 + IF(MOD(IW-1,2).EQ.0) THEN + SIGMS(IMAT,IW)=SIGT(IMERG,IGR,IW)*SPH(IMERG,IGR) + ELSE IF(MOD(IW-1,2).EQ.1) THEN + SIGMS(IMAT,IW)=SIGT(IMERG,IGR,IW)/SPH(IMERG,IGR) + ENDIF + 250 CONTINUE + ELSE IF(IMC.EQ.2) THEN +* TRANSPORT-PIJ EQUIVALENCE WITHOUT BELL FACTOR ACCELERATION. + SIGMA(IMAT,1)=SIGW(IMERG,IGR,1)*SPH(IMERG,IGR)+SIGT(IMERG,IGR,1) + > *(1.0-SPH(IMERG,IGR)) + DO 260 IW=1,NW+1 + SIGMS(IMAT,IW)=SIGT(IMERG,IGR,IW) + 260 CONTINUE + ELSE IF(IMC.EQ.3) THEN +* TRANSPORT-PIJ EQUIVALENCE WITH BELL FACTOR ACCELERATION. + SIGMA(IMAT,1)=0.0 + DO 270 IW=1,NW+1 + SIGMS(IMAT,IW)=SIGT(IMERG,IGR,IW) + 270 CONTINUE + ENDIF + 280 CONTINUE + NPSYS(IGR)=IGR + KPSYS2=LCMDIL(JPSYS2,IGR) + DO 290 IW=1,MIN(NW+1,10) + IF(IW.EQ.1) THEN + TEXT12='DRAGON-TXSC' + ELSE + WRITE(TEXT12,'(8HDRAGON-T,I1,3HXSC)') IW-1 + ENDIF + CALL LCMPUT(KPSYS2,TEXT12,NBMIX2+1,2,SIGMS(0,IW)) + 290 CONTINUE + CALL LCMPUT(KPSYS2,'DRAGON-S0XSC',NBMIX2+1,2,SIGMA(0,1)) + IF(LDIFF) THEN + SIGMD(0)=1.0E10 + CALL LCMPUT(KPSYS2,'DRAGON-DIFF',NBMIX2+1,2,SIGMD(0)) + ENDIF +*---- +* SPH CORRECTION OF PHYSICAL ALBEDOS +*---- + IF(NALBP.GT.0) THEN + DO 300 IAL=1,NALBP + FACT=0.5*(1.0-ALB2(IAL,IGR))/(1.0+ALB2(IAL,IGR))* + 1 SPH(NMERGE+IAL,IGR) + ALB1(IAL,IGR)=(1.0-2.0*FACT)/(1.0+2.0*FACT) + 300 CONTINUE + CALL LCMPUT(KPSYS2,'ALBEDO',NALBP,2,ALB1(1,IGR)) + ENDIF + 310 CONTINUE +*---- +* ASSEMBLY OF PIJ OR SYSTEM MATRICES AT ITERATION ITER. +*---- + IF(.NOT.C_ASSOCIATED(IPTRK2)) THEN + CALL XABORT('SPHEQU: MACRO-TRACKING NOT DEFINED(2).') + ENDIF + ISTRM=1 + KNORM=1 + IPHASE=2 + IF(NSPH.EQ.4) IPHASE=1 + IF(IPHASE.EQ.2) THEN + IPIJK=1 + ITPIJ=1 + LNORM=.FALSE. + CALL DOORPV(CDOOR,JPSYS2,NPSYS,IPTRK2,IFTRAK,IPRIN2,NGCOND, + > NREG2,NBMIX2,NANI,MAT2,VOL2,KNORM,IPIJK,ILK,ITPIJ,LNORM, + > CTITRE,NALBP) + ELSE + CALL DOORAV(CDOOR,JPSYS2,NPSYS,IPTRK2,IFTRAK,IPRIN2,NGCOND, + > NREG2,NBMIX2,NANI,NW,MAT2,VOL2,KNORM,ILK,CTITRE,NALBP,ISTRM) + ENDIF +*---- +* TRANSPORT-PIJ EQUIVALENCE WITH BELL FACTOR ACCELERATION. +*---- + IF(IMC.EQ.3) THEN + CALL SPHTRA(JPSYS2,ITER,NPSYS,KSPH,NREG2,NUN2,NMERGE,NALBP, + 1 NGCOND,SUNMER(1,1,1,1),FLXMER,NBMIX2,MAT2,VOL2,KEY2,MERG2,SPH, + 2 SIGW(1,1,1),SIGT(1,1,1),COURIN,FUNKNO) + SNORM(:NGCOND)=1.0 + GO TO 390 + ENDIF +*---- +* MACRO-FLUX CALCULATION AT ITERATION ITER. +*---- + ALLOCATE(SIGG(0:NBMIX2)) + DO IGR=1,NGCOND + SUNKNO(:NUN2,IGR)=0.0 + IF(NPSYS(IGR).EQ.0) CYCLE + IF(ITER.EQ.1) THEN + SIGG(0)=0.0 + DO IBM=1,NBMIX2 + PV=0.0 + DO JGR=1,NGCOND + PV=PV+SUNMER(MERG2(IBM),JGR,IGR,1)*FLXMER(MERG2(IBM),JGR) + ENDDO ! JGR + SIGG(IBM)=PV + ENDDO ! IBM + CALL DOORS(CDOOR,IPTRK2,NBMIX2,0,NUN2,SIGG,SUNKNO(1,IGR)) + ELSE + DO JGR=1,NGCOND + SIGG(0)=0.0 + DO IBM=1,NBMIX2 + SIGG(IBM)=SUNMER(MERG2(IBM),JGR,IGR,1)*SPH(MERG2(IBM),JGR) + ENDDO ! IBM + CALL DOORS(CDOOR,IPTRK2,NBMIX2,0,NUN2,SIGG,SUNKNO(1,IGR), + 1 FUNKNO(1,JGR)) + ENDDO + ENDIF + ENDDO ! IGR + DEALLOCATE(SIGG) +*---- +* COMPUTE THE MACRO-FLUX USING THE VECTORIAL DOOR. +*---- + IDIR=0 + IPSOU=C_NULL_PTR + LEXAC=.FALSE. + REBFLG=.FALSE. + CALL DOORFV(CDOOR,JPSYS2,NPSYS,IPTRK2,IFTRAK,IPRIN2,NGCOND, + 1 NBMIX2,IDIR,NREG2,NUN2,IPHASE,LEXAC,MAT2,VOL2,KEY2,CTITRE, + 2 SUNKNO,FUNKNO,IPMACR,IPSOU,REBFLG) +*---- +* COMPUTE MACRO-CALCULATION LEAKAGE RATES IF NALBP.GT.0 +*---- + IF(NALBP.GT.0) THEN + DO 340 IGR=1,NGCOND + OUTG1(IGR)=0.0 + DO 330 K=1,NREG2 + L=MAT2(K) + IF(L.EQ.0) GO TO 330 + IUN=KEY2(K) + IF(VOL2(K).EQ.0.0) GO TO 330 + IKK=MERG2(L) + OUTG1(IGR)=OUTG1(IGR)+(SIGW(IKK,IGR,1)-SIGT(IKK,IGR,1))* + > FUNKNO(IUN,IGR)*VOLMER(IKK)*SPH(IKK,IGR) + DO 320 JGR=1,NGCOND + OUTG1(IGR)=OUTG1(IGR)+SUNMER(IKK,JGR,IGR,1)*FUNKNO(IUN,JGR)* + > VOLMER(IKK)*SPH(IKK,JGR) + 320 CONTINUE + 330 CONTINUE + IF(IPRIN2.GT.0) WRITE(6,920) IGR,OUTG1(IGR) + 340 CONTINUE + ENDIF +*---- +* MACRO-FLUX NORMALIZATION. +*---- + IF(ILK.AND.(NALBP.EQ.0)) GO TO 390 + SNORM(:NGCOND)=1.0 + IF(KSPH.LT.0) THEN +* ASYMTTOTIC NORMALIZATION WITH RESPECT TO MIXTURE -KSPH. + IF(-KSPH.GT.NMERGE) CALL XABORT('SPHEQU: INVALID ASYMPTOTIC M' + > //'IXTURE SET.') + DO 360 IGR=1,NGCOND + IF(NPSYS(IGR).EQ.0) GO TO 360 + VOLTOT=0.0D0 + FLXTOT=0.0D0 + DO 350 IREG=1,NREG2 + IMAT=MAT2(IREG) + IF(IMAT.EQ.0) GO TO 350 + IF(MERG2(IMAT).EQ.-KSPH) THEN + VOLTOT=VOLTOT+VOL2(IREG) + FLXTOT=FLXTOT+FUNKNO(KEY2(IREG),IGR)*VOL2(IREG) + ENDIF + 350 CONTINUE + IF(VOLTOT.GT.0) THEN + FLXTOT=FLXTOT/VOLTOT + FFF=COURIN(IGR)/REAL(FLXTOT) + ELSE + FFF=1.0 + ENDIF + SNORM(IGR)=FFF + IF(IPRIN2.GT.0) WRITE(6,960) IGR,FFF + 360 CONTINUE + ELSE +* AVERAGE OR SELENGUT FLUX NORMALIZATION. + DO 380 IGR=1,NGCOND + IF(NPSYS(IGR).EQ.0) GO TO 380 + VOLTOT=0.0D0 + FLXTOT=0.0D0 + DO 370 IREG=1,NREG2 + IMAT=MAT2(IREG) + IF(IMAT.EQ.0) GO TO 370 + IMERG=MERG2(IMAT) + IF((KSPH.NE.7).OR.LFISS(IMERG)) THEN + VOLTOT=VOLTOT+VOL2(IREG) + FLXTOT=FLXTOT+FUNKNO(KEY2(IREG),IGR)*VOL2(IREG) + ENDIF + 370 CONTINUE + FLXTOT=FLXTOT/VOLTOT + IF(IPRINT.GE.100) THEN + WRITE(6,*)'FLXTOT: =',FLXTOT,'VOLTOT: =',VOLTOT + ENDIF + IF(KSPH.NE.6) THEN + SNORM(IGR)=COURIN(IGR)/REAL(FLXTOT) + ELSE +* Compute diffusion flux in water gap => COUROW(IGR) + IF(CDOOR.NE.'TRIVAC') CALL XABORT('SPHEQU: TRIVAC expected' + 1 //' as tracking module with SELE-GAP') + CALL SPHGAP(IPTRK2,IPRINT,NREG2,NUN2,MAT2,KEY2,FUNKNO(1,IGR), + 1 COUROW(IGR)) + IF(IPRINT.GE.100) THEN + WRITE(6,*)'COURIN: =',COURIN(IGR),'COUROW: =',COUROW(IGR) + ENDIF + SNORM(IGR)=COURIN(IGR)/COUROW(IGR) + ENDIF + IF(IPRIN2.GT.0) WRITE(6,960) IGR,SNORM(IGR) + 380 CONTINUE + ENDIF +*---- +* COMPUTE THE IMPROVED SPH FACTORS. +*---- + 390 ALLOCATE(SPHNEW(NMERGE+NALBP)) + DO 430 IGR=1,NGCOND + IF(NPSYS(IGR).EQ.0) GO TO 430 + VOLMER(:NMERGE)=0.0 + FACTOR(:NMERGE)=0.0 + SPHNEW(:NMERGE+NALBP)=1.0 + DO 400 IREG=1,NREG2 + IMAT=MAT2(IREG) + IF(IMAT.EQ.0) GO TO 400 + IKK=MERG2(IMAT) + VOLMER(IKK)=VOLMER(IKK)+VOL2(IREG) + FACTOR(IKK)=FACTOR(IKK)+FUNKNO(KEY2(IREG),IGR)*VOL2(IREG) + 400 CONTINUE + DO 410 IKK=1,NMERGE + IF(VOLMER(IKK).EQ.0.0) GO TO 410 + FACTOR(IKK)=FACTOR(IKK)/VOLMER(IKK) + SPHNEW(IKK)=FLXMER(IKK,IGR)/(SNORM(IGR)*FACTOR(IKK)) + IF(SPHNEW(IKK).LT.0.0) THEN + WRITE(6,980) IGR,IKK + SPHNEW(IKK)=1.0 + ENDIF + 410 CONTINUE + IF(NALBP.EQ.1) SPHNEW(NMERGE+1)=OUTG1(IGR)/(OUTG2(IGR)*SNORM(IGR)) + DO 420 IKK=1,NMERGE+NALBP + ERRT=ABS((SPHNEW(IKK)-SPH(IKK,IGR))/SPHNEW(IKK)) + ERR2=ERR2+ERRT*ERRT + ERROR=MAX(ERROR,ERRT) + SPH(IKK,IGR)=SPHNEW(IKK) + 420 CONTINUE + IF(IPRINT.GT.4) THEN + WRITE(6,930) 'NSPH',IGR,(SPH(IKK,IGR),IKK=1,NMERGE+NALBP) + ENDIF + IF(IPRINT.GT.5) THEN + WRITE(6,930) 'FUNKNO',IGR,(FUNKNO(IUNK,IGR),IUNK=1,NUN2) + ENDIF + 430 CONTINUE + DEALLOCATE(SPHNEW) + ERR2=SQRT(ERR2/(NMERGE*NGCOND)) + IF(IPRINT.GT.1) WRITE(6,935) ITER,ERROR,ERR2 + IF(IPRINT.GT.2) THEN + IF(ERROR.GE.EPSPH) WRITE(6,940) ((IKK,IGR,SPH(IKK,IGR), + > IKK=1,NMERGE+NALBP),IGR=1,NGCOND) + ENDIF + IF(ERR2.LT.EPSPH) GO TO 440 + IF((ITER.GT.1).AND.(ERR2.GT.OLDERR)) THEN + WRITE(6,970) ITER + NBIERR=NBIERR+1 + IF(NBIERR.GE.MAXNBI) THEN + WRITE(6,990) ITER + GO TO 440 + ENDIF + ENDIF + OLDERR=ERR2 + GO TO 230 + 440 WRITE(6,950) ITER +*---- +* RESET SOURCES TO NO DB2 LEAKAGE +*---- + IF(ILEAKS.NE.0) THEN + DO 460 IGR=1,NGCOND + DO 450 IKK=1,NMERGE + SUNMER(IKK,IGR,IGR,1)=SUNMER(IKK,IGR,IGR,1)+ZLEAK(IKK,IGR) + 450 CONTINUE + 460 CONTINUE + ENDIF + CALL LCMCL(IPSYS2,2) +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + 470 DEALLOCATE(LFISS,ALB2,ALB1,SNORM,OUTG2,OUTG1,FACTOR,SUNKNO, + 1 FUNKNO,SUNMER,ZLEAK,DIFF,FLXMER,SIGW,SIGT,VOLMER,COUROW,COURIN, + 2 SIGMD,SIGMS,SIGMA) + DEALLOCATE(NPSYS) + RETURN +* + 910 FORMAT(/44H SPHEQU: AVERAGE FLUXES PER MACRO-GROUPS ARE/ + > (1X,1P,10E13.5)) + 920 FORMAT(/8H SPHEQU:,5X,6HGROUP=,I4,15H MACRO LEAKAGE=,1P,E12.4) + 930 FORMAT(/26H SPHEQU: VALUES OF VECTOR ,A,9H IN GROUP,I5,4H ARE/ + > (1X,1P,10E13.5)) + 935 FORMAT(/14H SPHEQU: ITER=,I3,4X,6HERROR=,1P,E10.3,1X,6HERR 2=, + > E10.3) + 940 FORMAT(4X,4HSPH(,I3,1H,,I3,2H)=,F9.5,:,4X,4HSPH(,I3,1H,,I3,2H)=, + > F9.5,:,4X,4HSPH(,I3,1H,,I3,2H)=,F9.5,:,4X,4HSPH(,I3,1H,,I3,2H)=, + > F9.5,:,4X,4HSPH(,I3,1H,,I3,2H)=,F9.5) + 950 FORMAT(/40H SPHEQU: ENDING OF SPH CONVERGENCE AFTER,I5, + > 12H ITERATIONS.) + 960 FORMAT(/43H SPHEQU: FLUX NORMALIZATION FACTOR IN GROUP,I4,1H=,1P, + 1 E13.5) + 970 FORMAT(1X,'Warning: oscillations in SPH at iteration ',i10) + 980 FORMAT(1X,'Warning: negative SPH factor in group ',i5, + > ' and region ',i5,' set to 1.0') + 990 FORMAT(1X,'Warning: maximum of 3 error oscillations ', + >'in SPH convergence reached at iteration ',i10) + END -- cgit v1.2.3