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/NUMER3.f | 745 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 745 insertions(+) create mode 100644 Dragon/src/NUMER3.f (limited to 'Dragon/src/NUMER3.f') diff --git a/Dragon/src/NUMER3.f b/Dragon/src/NUMER3.f new file mode 100644 index 0000000..0d73db0 --- /dev/null +++ b/Dragon/src/NUMER3.f @@ -0,0 +1,745 @@ +*DECK NUMER3 + SUBROUTINE NUMER3 (NCOUR,MULTC,NCODE,ZCODE,LX,LY,LZ,IORI,ISM, + 1 POURCE,IMPX,NMBLK,IFR,ALB,SUR,NMERGE,INUM,MIX,DVX,NGEN,IGEN, + 2 XX,YY,ZZ) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Surface renumbering algorithm for Cartesian geometry. +* The 3-D DP-1 approximation is not implemented. +* +*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/output +* NCOUR number of surfaces per block (input); number of out-currents +* per block (output). +* MULTC type of multicell approximation: +* =1 Roth; =2 Roth X ncour; =3 DP-0; =4 DP-1. +* NCODE type of boundary condition on each side of the domain: +* =0 not used; =1 VOID; =2 REFL; +* =3 DIAG; =4 TRAN; =5 SYME. +* ZCODE value of the albedo on each side of the domain. +* LX number of blocks along the X-axis. +* LY number of blocks along the Y-axis. +* LZ number of blocks along the Z-axis. +* IORI orientation of the blocks. +* ISM permutation index corresponding to each orientation +* (ISM(I,N)=I is the natural orientation). +* POURCE weight associated with each merged block. +* IMPX print flag (equal to 0 for no print). +* NMBLK total number of blocks in the domain. +* IFR index-number of in-currents. +* ALB transmission/albedo associated with each in-current. +* SUR surface associated with each in-current. +* NMERGE total number of merged cells for which specific values +* of the neutron flux and reactions rates are required. +* Many cells with different position in the domain can +* be merged before the neutron flux calculation if they +* own the same generating cell (NMERGE.le.NMBLK). +* INUM index-number of the merged cell associated to each cell. +* MIX index-number of out-currents. +* DVX weight associated with each out-current. +* Note: IFR, ALB, MIX and DVX contains information to rebuild +* the geometrical 'A' matrix. +* NGEN total number of generating blocks in the cartesian domain. +* IGEN index-number of the generating block associated with each +* merged block. +* XX X-thickness of the generating blocks. +* YY Y-thickness of the generating blocks. +* ZZ Z-thickness of the generating blocks. +* +*----------------------------------------------------------------------- +* +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER NCOUR,MULTC,NCODE(6),LX,LY,LZ,IORI(NMBLK),ISM(6,8), + 1 IMPX,NMBLK,IFR(12*NMBLK),NMERGE,INUM(NMBLK),MIX(12*NMERGE), + 2 NGEN,IGEN(NMERGE) + REAL ZCODE(6),POURCE(NMERGE),ALB(12*NMBLK),SUR(12*NMBLK), + 1 DVX(12*NMERGE),XX(NGEN),YY(NGEN),ZZ(NGEN) +*---- +* LOCAL VARIABLES +*---- + PARAMETER (EPS=1.0E-5) + LOGICAL LL1,LL2,LOG1,LOG2,LOG3 + CHARACTER DIRR(6)*2,DIRZ(12)*2,HSMG*131 + INTEGER IDDD(6),ISMZ(12) + REAL DDD(6) + INTEGER, ALLOCATABLE, DIMENSION(:) :: JF2 + REAL, ALLOCATABLE, DIMENSION(:) :: GG3 + SAVE DIRR + DATA DIRR/'X-','X+','Y-','Y+','Z-','Z+'/ +*---- +* SCRATCH STORAGE ALLOCATION +*---- + ALLOCATE(JF2(6*NMBLK),GG3(NMBLK)) +* + DO 100 I=1,NCOUR*NMERGE + MIX(I)=I + DVX(I)=1.0 +100 CONTINUE + IS1=0 + IS2=0 + LXY=LX*LY + LL1=((NCODE(2).EQ.3).AND.(NCODE(3).EQ.3)) + LL2=((NCODE(1).EQ.3).AND.(NCODE(4).EQ.3)) + IF (LL1) THEN + IS1=1 + LXY=LX*(LX+1)/2 + ELSE IF (LL2) THEN + IS2=1 + LXY=LX*(LX+1)/2 + ENDIF + IBLK=0 + DO 280 K0=1,LZ + DO 275 K1=1,LY + LXM=1 + LXP=LX + IF (LL1) LXP=K1 + IF (LL2) LXM=K1 + DO 270 K2=LXM,LXP + IBLK=IBLK+1 + IKK=INUM(IBLK) + FRX=1.0 + FRY=1.0 + FRZ=1.0 + IF (IKK.EQ.0) GO TO 265 + IS=NCOUR*(IBLK-1) + IT=NCOUR*(IKK-1) + II=IORI(IBLK) + DO 110 IC=1,6 + IDDD(IC)=-1 +110 CONTINUE + IF (K2.GT.1) IDDD(1)=IBLK-1 + IF (K2.LT.LX) IDDD(2)=IBLK+1 + IF (K1.GT.1) IDDD(3)=IBLK-(LXP-LXM+1)+IS1 + IF (K1.LT.LY) IDDD(4)=IBLK+(LXP-LXM+1)-IS2 + IF (K0.GT.1) IDDD(5)=IBLK-LXY + IF (K0.LT.LZ) IDDD(6)=IBLK+LXY +* + DO 120 IC=1,NCOUR + ALB(IS+IC)=1.0 + SUR(IS+IC)=0.0 + JBLK=IDDD(IC) + IF (JBLK.GT.0) THEN + JKK=INUM(JBLK) + JT=NCOUR*(JKK-1) + IF ((MOD(IC,2).EQ.1).AND.(JKK.GT.0)) THEN + IFR(IS+ISM(IC,II))=JT+ISM(IC+1,IORI(JBLK)) + ELSE IF ((MOD(IC,2).EQ.0).AND.(JKK.GT.0)) THEN + IFR(IS+ISM(IC,II))=JT+ISM(IC-1,IORI(JBLK)) + ELSE + IFR(IS+ISM(IC,II))=0 + ENDIF + IDDD(IC)=JKK + ELSE + IFR(IS+ISM(IC,II))=0 + ENDIF +120 CONTINUE +*---- +* VOID OR REFL BOUNDARY CONDITIONS +*---- + IKG=IGEN(IKK) + LOG1=(K2.EQ.1).OR.(IDDD(1).EQ.0) + LOG2=(NCODE(1).EQ.1).OR.(LL2.AND.(NCODE(3).EQ.1)) + LOG3=(NCODE(1).EQ.2).OR.(LL2.AND.(NCODE(3).EQ.2)) + IF (LOG1.AND.LOG2) THEN + ALB(IS+ISM(1,II))=-ZCODE(1) + IFR(IS+ISM(1,II))=IT+ISM(1,II) + ELSE IF (LOG1.AND.LOG3) THEN + ALB(IS+ISM(1,II))=-1.0 + IFR(IS+ISM(1,II))=IT+ISM(1,II) + ENDIF + IF(LOG1.AND.(NCODE(1).EQ.1)) SUR(IS+ISM(1,II))=YY(IKG)*ZZ(IKG) + IF(LOG1.AND.(NCODE(1).EQ.2)) SUR(IS+ISM(1,II))=YY(IKG)*ZZ(IKG) + LOG1=(K2.EQ.LX).OR.(IDDD(2).EQ.0) + LOG2=(NCODE(2).EQ.1).OR.(LL1.AND.(NCODE(4).EQ.1)) + LOG3=(NCODE(2).EQ.2).OR.(LL1.AND.(NCODE(4).EQ.2)) + IF (LOG1.AND.LOG2) THEN + ALB(IS+ISM(2,II))=-ZCODE(2) + IFR(IS+ISM(2,II))=IT+ISM(2,II) + ELSE IF (LOG1.AND.LOG3) THEN + ALB(IS+ISM(2,II))=-1.0 + IFR(IS+ISM(2,II))=IT+ISM(2,II) + ENDIF + IF(LOG1.AND.(NCODE(2).EQ.1)) SUR(IS+ISM(2,II))=YY(IKG)*ZZ(IKG) + IF(LOG1.AND.(NCODE(2).EQ.2)) SUR(IS+ISM(2,II))=YY(IKG)*ZZ(IKG) + LOG1=(K1.EQ.1).OR.(IDDD(3).EQ.0) + LOG2=(NCODE(3).EQ.1).OR.(LL1.AND.(NCODE(1).EQ.1)) + LOG3=(NCODE(3).EQ.2).OR.(LL1.AND.(NCODE(1).EQ.2)) + IF (LOG1.AND.LOG2) THEN + ALB(IS+ISM(3,II))=-ZCODE(3) + IFR(IS+ISM(3,II))=IT+ISM(3,II) + ELSE IF (LOG1.AND.LOG3) THEN + ALB(IS+ISM(3,II))=-1.0 + IFR(IS+ISM(3,II))=IT+ISM(3,II) + ENDIF + IF(LOG1.AND.(NCODE(3).EQ.1)) SUR(IS+ISM(3,II))=XX(IKG)*ZZ(IKG) + IF(LOG1.AND.(NCODE(3).EQ.2)) SUR(IS+ISM(3,II))=XX(IKG)*ZZ(IKG) + LOG1=(K1.EQ.LY).OR.(IDDD(4).EQ.0) + LOG2=(NCODE(4).EQ.1).OR.(LL2.AND.(NCODE(2).EQ.1)) + LOG3=(NCODE(4).EQ.2).OR.(LL2.AND.(NCODE(2).EQ.2)) + IF (LOG1.AND.LOG2) THEN + ALB(IS+ISM(4,II))=-ZCODE(4) + IFR(IS+ISM(4,II))=IT+ISM(4,II) + ELSE IF (LOG1.AND.LOG3) THEN + ALB(IS+ISM(4,II))=-1.0 + IFR(IS+ISM(4,II))=IT+ISM(4,II) + ENDIF + IF(LOG1.AND.(NCODE(4).EQ.1)) SUR(IS+ISM(4,II))=XX(IKG)*ZZ(IKG) + IF(LOG1.AND.(NCODE(4).EQ.2)) SUR(IS+ISM(4,II))=XX(IKG)*ZZ(IKG) + LOG1=(K0.EQ.1).OR.(IDDD(5).EQ.0) + IF (LOG1.AND.(NCODE(5).EQ.1)) THEN + ALB(IS+ISM(5,II))=-ZCODE(5) + IFR(IS+ISM(5,II))=IT+ISM(5,II) + ELSE IF (LOG1.AND.(NCODE(5).EQ.2)) THEN + ALB(IS+ISM(5,II))=-1.0 + IFR(IS+ISM(5,II))=IT+ISM(5,II) + ENDIF + IF(LOG1.AND.(NCODE(5).EQ.1)) SUR(IS+ISM(5,II))=XX(IKG)*YY(IKG) + IF(LOG1.AND.(NCODE(5).EQ.2)) SUR(IS+ISM(5,II))=XX(IKG)*YY(IKG) + LOG1=(K0.EQ.LZ).OR.(IDDD(6).EQ.0) + IF (LOG1.AND.(NCODE(6).EQ.1)) THEN + ALB(IS+ISM(6,II))=-ZCODE(6) + IFR(IS+ISM(6,II))=IT+ISM(6,II) + ELSE IF (LOG1.AND.(NCODE(6).EQ.2)) THEN + ALB(IS+ISM(6,II))=-1.0 + IFR(IS+ISM(6,II))=IT+ISM(6,II) + ENDIF + IF(LOG1.AND.(NCODE(6).EQ.1)) SUR(IS+ISM(6,II))=XX(IKG)*YY(IKG) + IF(LOG1.AND.(NCODE(6).EQ.2)) SUR(IS+ISM(6,II))=XX(IKG)*YY(IKG) +*---- +* CORRECT THE PARITY OF THE INTERFACE CURRENTS FOR DP-1 CASES WITH +* 'MIRROR' ORIENTATION +*---- + DO 125 IC=1,NCOUR + IF(II.GE.5) ALB(IS+IC)=-ALB(IS+IC) + JBLK=IBLK + IF((K2.GT.1).AND.(IC.EQ.1)) JBLK=IBLK-1 + IF((K2.LT.LX).AND.(IC.EQ.2)) JBLK=IBLK+1 + IF((K1.GT.1).AND.(IC.EQ.3)) JBLK=IBLK-(LXP-LXM+1)+IS1 + IF((K1.LT.LY).AND.(IC.EQ.4)) JBLK=IBLK+(LXP-LXM+1)-IS2 + IF((K0.GT.1).AND.(IC.EQ.5)) JBLK=IBLK-LXY + IF((K0.LT.LZ).AND.(IC.EQ.6)) JBLK=IBLK+LXY + IF(IORI(JBLK).GE.5) ALB(IS+ISM(IC,II))=-ALB(IS+ISM(IC,II)) +125 CONTINUE +*---- +* DIAG BOUNDARY CONDITION +*---- + IF (K1.EQ.K2) THEN + IF(LL1.OR.LL2) THEN + IKG=IGEN(IKK) + IF(XX(IKG).NE.YY(IKG)) CALL XABORT('NUMER3: A CELL ON THE ' + 1 //'DIAGONAL SYMMETRY AXIS IS NOT SQUARE.') + ENDIF + IF ((K1.EQ.1).AND.(NCODE(1).EQ.3).AND.(NCODE(3).EQ.5)) THEN + FRX=0.25 + ALB(IS+ISM(1,II))=-ALB(IS+ISM(2,II)) + IFR(IS+ISM(1,II))=IFR(IS+ISM(2,II)) + DVX(IT+ISM(1,II))=-DVX(IT+ISM(2,II)) + MIXNEW=MIX(IT+ISM(2,II)) + MIXOLD=MIX(IT+ISM(1,II)) + DO 130 KC=1,NCOUR + IF(MIX(IT+KC).EQ.MIXOLD) MIX(IT+KC)=MIXNEW +130 CONTINUE + ELSE IF (NCODE(1).EQ.3) THEN + FRX=0.5 + ALB(IS+ISM(1,II))=-ALB(IS+ISM(3,II)) + IFR(IS+ISM(1,II))=IFR(IS+ISM(3,II)) + DVX(IT+ISM(1,II))=-DVX(IT+ISM(3,II)) + MIXNEW=MIX(IT+ISM(3,II)) + MIXOLD=MIX(IT+ISM(1,II)) + DO 140 KC=1,NCOUR + IF(MIX(IT+KC).EQ.MIXOLD) MIX(IT+KC)=MIXNEW +140 CONTINUE + ENDIF + IF ((K1.EQ.LY).AND.(NCODE(2).EQ.3).AND.(NCODE(4).EQ.5)) THEN + FRX=0.25 + ALB(IS+ISM(2,II))=-ALB(IS+ISM(1,II)) + IFR(IS+ISM(2,II))=IFR(IS+ISM(1,II)) + DVX(IT+ISM(2,II))=-DVX(IT+ISM(1,II)) + MIXNEW=MIX(IT+ISM(1,II)) + MIXOLD=MIX(IT+ISM(2,II)) + DO 150 KC=1,NCOUR + IF(MIX(IT+KC).EQ.MIXOLD) MIX(IT+KC)=MIXNEW +150 CONTINUE + ELSE IF (NCODE(2).EQ.3) THEN + ALB(IS+ISM(2,II))=-ALB(IS+ISM(4,II)) + IFR(IS+ISM(2,II))=IFR(IS+ISM(4,II)) + DVX(IT+ISM(2,II))=-DVX(IT+ISM(4,II)) + MIXNEW=MIX(IT+ISM(4,II)) + MIXOLD=MIX(IT+ISM(2,II)) + DO 160 KC=1,NCOUR + IF(MIX(IT+KC).EQ.MIXOLD) MIX(IT+KC)=MIXNEW +160 CONTINUE + ENDIF + IF ((K1.EQ.1).AND.(NCODE(3).EQ.3).AND.(NCODE(1).EQ.5)) THEN + FRY=0.25 + ALB(IS+ISM(3,II))=-ALB(IS+ISM(4,II)) + IFR(IS+ISM(3,II))=IFR(IS+ISM(4,II)) + DVX(IT+ISM(3,II))=-DVX(IT+ISM(4,II)) + MIXNEW=MIX(IT+ISM(4,II)) + MIXOLD=MIX(IT+ISM(3,II)) + DO 170 KC=1,NCOUR + IF(MIX(IT+KC).EQ.MIXOLD) MIX(IT+KC)=MIXNEW +170 CONTINUE + ELSE IF (NCODE(3).EQ.3) THEN + FRY=0.5 + ALB(IS+ISM(3,II))=-ALB(IS+ISM(1,II)) + IFR(IS+ISM(3,II))=IFR(IS+ISM(1,II)) + DVX(IT+ISM(3,II))=-DVX(IT+ISM(1,II)) + MIXNEW=MIX(IT+ISM(1,II)) + MIXOLD=MIX(IT+ISM(3,II)) + DO 180 KC=1,NCOUR + IF(MIX(IT+KC).EQ.MIXOLD) MIX(IT+KC)=MIXNEW +180 CONTINUE + ENDIF + IF ((K1.EQ.LY).AND.(NCODE(4).EQ.3).AND.(NCODE(2).EQ.5)) THEN + FRY=0.25 + ALB(IS+ISM(4,II))=-ALB(IS+ISM(3,II)) + IFR(IS+ISM(4,II))=IFR(IS+ISM(3,II)) + DVX(IT+ISM(4,II))=-DVX(IT+ISM(3,II)) + MIXNEW=MIX(IT+ISM(3,II)) + MIXOLD=MIX(IT+ISM(4,II)) + DO 190 KC=1,NCOUR + IF(MIX(IT+KC).EQ.MIXOLD) MIX(IT+KC)=MIXNEW +190 CONTINUE + ELSE IF (NCODE(4).EQ.3) THEN + ALB(IS+ISM(4,II))=-ALB(IS+ISM(2,II)) + IFR(IS+ISM(4,II))=IFR(IS+ISM(2,II)) + DVX(IT+ISM(4,II))=-DVX(IT+ISM(2,II)) + MIXNEW=MIX(IT+ISM(2,II)) + MIXOLD=MIX(IT+ISM(4,II)) + DO 200 KC=1,NCOUR + IF(MIX(IT+KC).EQ.MIXOLD) MIX(IT+KC)=MIXNEW +200 CONTINUE + ENDIF + ENDIF +*---- +* TRAN BOUNDARY CONDITION +*---- + IF ((K2.EQ.1).AND.(NCODE(1).EQ.4)) THEN + JBLK=IBLK+LXP-LXM + IFR(IS+ISM(1,II))=NCOUR*(INUM(JBLK)-1)+ISM(2,IORI(JBLK)) + ENDIF + IF ((K2.EQ.LX).AND.(NCODE(2).EQ.4)) THEN + JBLK=IBLK+LXM-LXP + IFR(IS+ISM(2,II))=NCOUR*(INUM(JBLK)-1)+ISM(1,IORI(JBLK)) + ENDIF + IF ((K1.EQ.1).AND.(NCODE(3).EQ.4)) THEN + JBLK=IBLK+(LY-1)*LX + IFR(IS+ISM(3,II))=NCOUR*(INUM(JBLK)-1)+ISM(4,IORI(JBLK)) + ENDIF + IF ((K1.EQ.LY).AND.(NCODE(4).EQ.4)) THEN + JBLK=IBLK-(LY-1)*LX + IFR(IS+ISM(4,II))=NCOUR*(INUM(JBLK)-1)+ISM(3,IORI(JBLK)) + ENDIF + IF ((K0.EQ.1).AND.(NCODE(5).EQ.4)) THEN + JBLK=IBLK+(LZ-1)*LXY + IFR(IS+ISM(5,II))=NCOUR*(INUM(JBLK)-1)+ISM(6,IORI(JBLK)) + ENDIF + IF ((K0.EQ.LZ).AND.(NCODE(6).EQ.4)) THEN + JBLK=IBLK-(LZ-1)*LXY + IFR(IS+ISM(6,II))=NCOUR*(INUM(JBLK)-1)+ISM(5,IORI(JBLK)) + ENDIF +*---- +* SYME BOUNDARY CONDITION +*---- + IF ((K2.EQ.1).AND.(NCODE(1).EQ.5)) THEN + FRX=0.5 + ALB(IS+ISM(1,II))=-ALB(IS+ISM(2,II)) + IFR(IS+ISM(1,II))=IFR(IS+ISM(2,II)) + SUR(IS+ISM(3,II))=0.5*SUR(IS+ISM(3,II)) + SUR(IS+ISM(4,II))=0.5*SUR(IS+ISM(4,II)) + IF(ISM(5,II).NE.0) SUR(IS+ISM(5,II))=0.5*SUR(IS+ISM(5,II)) + IF(ISM(6,II).NE.0) SUR(IS+ISM(6,II))=0.5*SUR(IS+ISM(6,II)) + DVX(IT+ISM(1,II))=-DVX(IT+ISM(2,II)) + MIXNEW=MIX(IT+ISM(2,II)) + MIXOLD=MIX(IT+ISM(1,II)) + DO 210 KC=1,NCOUR + IF(MIX(IT+KC).EQ.MIXOLD) MIX(IT+KC)=MIXNEW +210 CONTINUE + ELSE IF ((K2.EQ.LX).AND.(NCODE(2).EQ.5)) THEN + FRX=0.5 + ALB(IS+ISM(2,II))=-ALB(IS+ISM(1,II)) + IFR(IS+ISM(2,II))=IFR(IS+ISM(1,II)) + SUR(IS+ISM(3,II))=0.5*SUR(IS+ISM(3,II)) + SUR(IS+ISM(4,II))=0.5*SUR(IS+ISM(4,II)) + IF(ISM(5,II).NE.0) SUR(IS+ISM(5,II))=0.5*SUR(IS+ISM(5,II)) + IF(ISM(6,II).NE.0) SUR(IS+ISM(6,II))=0.5*SUR(IS+ISM(6,II)) + DVX(IT+ISM(2,II))=-DVX(IT+ISM(1,II)) + MIXNEW=MIX(IT+ISM(1,II)) + MIXOLD=MIX(IT+ISM(2,II)) + DO 220 KC=1,NCOUR + IF(MIX(IT+KC).EQ.MIXOLD) MIX(IT+KC)=MIXNEW +220 CONTINUE + ENDIF + IF ((K1.EQ.1).AND.(NCODE(3).EQ.5)) THEN + FRY=0.5 + ALB(IS+ISM(3,II))=-ALB(IS+ISM(4,II)) + IFR(IS+ISM(3,II))=IFR(IS+ISM(4,II)) + SUR(IS+ISM(1,II))=0.5*SUR(IS+ISM(1,II)) + SUR(IS+ISM(2,II))=0.5*SUR(IS+ISM(2,II)) + IF(ISM(5,II).NE.0) SUR(IS+ISM(5,II))=0.5*SUR(IS+ISM(5,II)) + IF(ISM(6,II).NE.0) SUR(IS+ISM(6,II))=0.5*SUR(IS+ISM(6,II)) + DVX(IT+ISM(3,II))=-DVX(IT+ISM(4,II)) + MIXNEW=MIX(IT+ISM(4,II)) + MIXOLD=MIX(IT+ISM(3,II)) + DO 230 KC=1,NCOUR + IF(MIX(IT+KC).EQ.MIXOLD) MIX(IT+KC)=MIXNEW +230 CONTINUE + ELSE IF ((K1.EQ.LY).AND.(NCODE(4).EQ.5)) THEN + FRY=0.5 + ALB(IS+ISM(4,II))=-ALB(IS+ISM(3,II)) + IFR(IS+ISM(4,II))=IFR(IS+ISM(3,II)) + SUR(IS+ISM(1,II))=0.5*SUR(IS+ISM(1,II)) + SUR(IS+ISM(2,II))=0.5*SUR(IS+ISM(2,II)) + IF(ISM(5,II).NE.0) SUR(IS+ISM(5,II))=0.5*SUR(IS+ISM(5,II)) + IF(ISM(6,II).NE.0) SUR(IS+ISM(6,II))=0.5*SUR(IS+ISM(6,II)) + DVX(IT+ISM(4,II))=-DVX(IT+ISM(3,II)) + MIXNEW=MIX(IT+ISM(3,II)) + MIXOLD=MIX(IT+ISM(4,II)) + DO 240 KC=1,NCOUR + IF(MIX(IT+KC).EQ.MIXOLD) MIX(IT+KC)=MIXNEW +240 CONTINUE + ENDIF + IF ((K0.EQ.1).AND.(NCODE(5).EQ.5)) THEN + FRZ=0.5 + ALB(IS+ISM(5,II))=-ALB(IS+ISM(6,II)) + IFR(IS+ISM(5,II))=IFR(IS+ISM(6,II)) + SUR(IS+ISM(1,II))=0.5*SUR(IS+ISM(1,II)) + SUR(IS+ISM(2,II))=0.5*SUR(IS+ISM(2,II)) + SUR(IS+ISM(3,II))=0.5*SUR(IS+ISM(3,II)) + SUR(IS+ISM(4,II))=0.5*SUR(IS+ISM(4,II)) + DVX(IT+ISM(5,II))=-DVX(IT+ISM(6,II)) + MIXNEW=MIX(IT+ISM(6,II)) + MIXOLD=MIX(IT+ISM(5,II)) + DO 250 KC=1,NCOUR + IF(MIX(IT+KC).EQ.MIXOLD) MIX(IT+KC)=MIXNEW +250 CONTINUE + ELSE IF ((K0.EQ.LZ).AND.(NCODE(6).EQ.5)) THEN + FRZ=0.5 + ALB(IS+ISM(6,II))=-ALB(IS+ISM(5,II)) + IFR(IS+ISM(6,II))=IFR(IS+ISM(5,II)) + SUR(IS+ISM(1,II))=0.5*SUR(IS+ISM(1,II)) + SUR(IS+ISM(2,II))=0.5*SUR(IS+ISM(2,II)) + SUR(IS+ISM(3,II))=0.5*SUR(IS+ISM(3,II)) + SUR(IS+ISM(4,II))=0.5*SUR(IS+ISM(4,II)) + DVX(IT+ISM(6,II))=-DVX(IT+ISM(5,II)) + MIXNEW=MIX(IT+ISM(5,II)) + MIXOLD=MIX(IT+ISM(6,II)) + DO 260 KC=1,NCOUR + IF(MIX(IT+KC).EQ.MIXOLD) MIX(IT+KC)=MIXNEW +260 CONTINUE + ENDIF +* +265 GG3(IBLK)=FRX*FRY*FRZ +270 CONTINUE +275 CONTINUE +280 CONTINUE + DO 285 I=1,NCOUR*NMBLK + IFR(I)=MIX(IFR(I)) +285 CONTINUE +*---- +* ELIMINATION OF THE BLOCKS OUTSIDE THE DOMAIN +*---- + JBLK=0 + DO 300 IBLK=1,NMBLK + IKK=INUM(IBLK) + IF (IKK.GT.0) THEN + JBLK=JBLK+1 + INUM(JBLK)=IKK + IORI(JBLK)=IORI(IBLK) + GG3(JBLK)=GG3(IBLK) + IS=NCOUR*(IBLK-1) + JS=NCOUR*(JBLK-1) + DO 290 IC=1,NCOUR + IFR(JS+IC)=IFR(IS+IC) + ALB(JS+IC)=ALB(IS+IC) + SUR(JS+IC)=SUR(IS+IC) +290 CONTINUE + ENDIF +300 CONTINUE + NMBLK=JBLK +* + DO 310 IKK=1,NMERGE + POURCE(IKK)=0.0 +310 CONTINUE + DO 330 IBLK=1,NMBLK + IKK=INUM(IBLK) + POURCE(IKK)=POURCE(IKK)+GG3(IBLK) +330 CONTINUE +*---- +* VALIDATION OF VECTOR IFR +*---- + DO 345 IBLK=1,NMBLK + IS=NCOUR*(IBLK-1) + DO 340 IC=1,NCOUR + ISURF=IFR(IS+IC) + IF (ISURF.EQ.0) THEN + WRITE (HSMG,'(44HNUMER3: FAILURE OF THE SURFACE RENUMBERING A, + 1 12HLGORITHM(1).)') + GO TO 570 + ENDIF + JC=1+MOD(ISURF-1,NCOUR) + JT=NCOUR*((ISURF-1)/NCOUR) + IF (MIX(JT+JC).NE.ISURF) THEN + WRITE (HSMG,'(44HNUMER3: FAILURE OF THE SURFACE RENUMBERING A, + 1 12HLGORITHM(2).)') + GO TO 570 + ENDIF +340 CONTINUE +345 CONTINUE +*---- +* VALIDATION OF THE GEOMETRICAL RECIPROCITY AND MODIFICATION OF ALBEDOS +*---- + DO 355 IBLK=1,NMBLK + IKK=INUM(IBLK) + IKG=IGEN(IKK) + IS=NCOUR*(IBLK-1) + IT=NCOUR*(IKK-1) + DO 350 IC=1,NCOUR + ALB(IS+IC)=ALB(IS+IC)*GG3(IBLK)/POURCE(IKK) + IP=1+MOD(MIX(IT+IC)-1,NCOUR) + FR1=YY(IKG)*ZZ(IKG) + IF ((IP.EQ.3).OR.(IP.EQ.4)) FR1=XX(IKG)*ZZ(IKG) + IF ((IP.EQ.5).OR.(IP.EQ.6)) FR1=XX(IKG)*YY(IKG) + JP=1+MOD(IFR(IS+IC)-1,NCOUR) + JKG=IGEN(1+(IFR(IS+IC)-1)/NCOUR) + FR2=YY(JKG)*ZZ(JKG) + IF ((JP.EQ.3).OR.(JP.EQ.4)) FR2=XX(JKG)*ZZ(JKG) + IF ((JP.EQ.5).OR.(JP.EQ.6)) FR2=XX(JKG)*YY(JKG) + DELTA=ABS(FR1-FR2) + IF (ABS(FR1-FR2).GT.EPS) THEN + WRITE (HSMG,680) DIRR(IP),IKG,DIRR(JP),JKG + GO TO 570 + ENDIF +350 CONTINUE +355 CONTINUE +*---- +* COMPUTE VECTOR DVX +*---- + DO 395 IKK=1,NMERGE + IKG=IGEN(IKK) + IF (NCOUR.EQ.2) THEN + DDD(1)=0.5 + DDD(2)=0.5 + ELSE IF (NCOUR.EQ.4) THEN + SURFA=2.0*(XX(IKG)+YY(IKG)) + DO 360 IC=1,NCOUR + FR1=YY(IKG) + IF ((IC.EQ.3).OR.(IC.EQ.4)) FR1=XX(IKG) + DDD(IC)=FR1/SURFA +360 CONTINUE + ELSE IF (NCOUR.EQ.6) THEN + SURFA=2.0*(XX(IKG)*ZZ(IKG)+YY(IKG)*ZZ(IKG)+XX(IKG)*YY(IKG)) + DO 370 IC=1,NCOUR + FR1=YY(IKG)*ZZ(IKG) + IF ((IC.EQ.3).OR.(IC.EQ.4)) FR1=XX(IKG)*ZZ(IKG) + IF ((IC.EQ.5).OR.(IC.EQ.6)) FR1=XX(IKG)*YY(IKG) + DDD(IC)=FR1/SURFA +370 CONTINUE + ENDIF + IT=NCOUR*(IKK-1) + DO 390 IC=1,NCOUR + IF (MULTC.EQ.1) THEN +* ROTH APPROXIMATION. + DVX(IT+IC)=DDD(IC) + MIX(IT+IC)=IKK + ELSE + DELTA=0.0 + I1=MIX(IT+IC) + DO 380 JC=1,NCOUR + IF (MIX(IT+JC).EQ.I1) DELTA=DELTA+DDD(JC) +380 CONTINUE + ZSIGN=SIGN(1.0,DVX(IT+IC)) + DVX(IT+IC)=ZSIGN*DDD(IC)/DELTA + ENDIF +390 CONTINUE +395 CONTINUE + IJAS=NCOUR*NMBLK + IJAR=NCOUR*NMERGE +*---- +* RECOMPUTE VECTOR IFR FOR ROTH APPROXIMATION +*---- + IF (MULTC.EQ.1) THEN + DO 400 I=1,IJAS + IFR(I)=1+(IFR(I)-1)/NCOUR +400 CONTINUE + ENDIF +*---- +* REMOVE THE UNUSED SURFACE NUMBERS +*---- + DO 410 I=1,IJAS + JF2(I)=0 +410 CONTINUE + IJAT=0 + DO 420 I=1,IJAR + J=MIX(I) + IF (J.GT.IJAS) THEN + WRITE (HSMG,'(44HNUMER3: FAILURE OF THE SURFACE RENUMBERING A, + 1 12HLGORITHM(3).)') + GO TO 570 + ENDIF + IF (JF2(J).EQ.0) THEN + IJAT=IJAT+1 + JF2(J)=IJAT + ENDIF +420 CONTINUE + DO 430 I=1,IJAR + MIX(I)=JF2(MIX(I)) +430 CONTINUE + DO 440 I=1,IJAS + IFR(I)=JF2(IFR(I)) +440 CONTINUE +*---- +* INCLUDE THE DP-1 APPROXIMATION +*---- + IF ((MULTC.EQ.4).AND.(NCOUR.EQ.2)) THEN +* DP-1 APPROXIMATION IN 1-D. + DO 455 I1=IJAR,1,-1 + FR1=ABS(DVX(I1)) + JND=(MIX(I1)-1)*2 + DO 450 JCOUR=1,2 + JSURF=(I1-1)*2+JCOUR + DVX(JSURF)=FR1 + MIX(JSURF)=JND+JCOUR +450 CONTINUE +455 CONTINUE + DO 465 I1=IJAS,1,-1 + FR1=ABS(ALB(I1)) + FR2=SUR(I1) + JND=(IFR(I1)-1)*2 + DO 460 JCOUR=1,2 + JSURF=(I1-1)*2+JCOUR + ALB(JSURF)=FR1 + SUR(JSURF)=FR2 + IFR(JSURF)=JND+JCOUR +460 CONTINUE +465 CONTINUE + NCOUR=4 + ELSE IF ((MULTC.EQ.4).AND.(NCOUR.EQ.4)) THEN +* DP-1 APPROXIMATION IN 2-D. + DO 480 I1=IJAR,1,-1 + ZSIGN=SIGN(1.0,DVX(I1)) + FR1=ABS(DVX(I1)) + FR2=SUR(I1) + JND=(MIX(I1)-1)*3 + DO 470 JCOUR=1,3 + JSURF=(I1-1)*3+JCOUR + DVX(JSURF)=FR1 + MIX(JSURF)=JND+JCOUR +470 CONTINUE + DVX(JSURF)=ZSIGN*FR1 +480 CONTINUE + DO 500 I1=IJAS,1,-1 + ZSIGN=SIGN(1.0,ALB(I1)) + FR1=ABS(ALB(I1)) + FR2=SUR(I1) + JND=(IFR(I1)-1)*3 + DO 490 JCOUR=1,3 + JSURF=(I1-1)*3+JCOUR + ALB(JSURF)=FR1 + SUR(JSURF)=FR2 + IFR(JSURF)=JND+JCOUR +490 CONTINUE + ALB(JSURF)=ZSIGN*FR1 +500 CONTINUE + NCOUR=12 + ELSE IF ((MULTC.EQ.4).AND.(NCOUR.EQ.6)) THEN + CALL XABORT('NUMER3: INVALID OPTION.') + ELSE + DO 510 I=1,IJAS + ALB(I)=ABS(ALB(I)) +510 CONTINUE + DO 520 I=1,IJAR + DVX(I)=ABS(DVX(I)) +520 CONTINUE + ENDIF + IJAS=NCOUR*NMBLK + IJAR=NCOUR*NMERGE +*---- +* PRINT THE SURFACE NUMBERS AFTER MERGING +*---- + IF (IMPX.GT.2) THEN + WRITE (6,620) + MIN6=MIN(6,NCOUR) + WRITE (6,650) ('----------------',I=1,MIN6) + DO 560 IBLK=1,NMBLK + IKK=INUM(IBLK) + WRITE (6,630) IBLK,IKK,IGEN(IKK) + I1=IORI(IBLK) + IF ((MULTC.EQ.4).AND.(NCOUR.EQ.12)) THEN + DO 530 I=1,12 + ISMZ(I)=3*ISM(1+(I-1)/3,I1)+MOD(I-1,3)-2 + DIRZ(I)=DIRR(1+(I-1)/3) +530 CONTINUE + ELSE IF ((MULTC.EQ.4).AND.(NCOUR.EQ.4)) THEN + DO 540 I=1,4 + ISMZ(I)=2*ISM(1+(I-1)/2,I1)+MOD(I-1,2)-1 + DIRZ(I)=DIRR(1+(I-1)/2) +540 CONTINUE + ELSE + DO 550 I=1,NCOUR + ISMZ(I)=ISM(I,I1) + DIRZ(I)=DIRR(I) +550 CONTINUE + ENDIF + IT0=NCOUR*(IBLK-1) + IT1=NCOUR*(IKK-1) + WRITE (6,660) (DIRZ(I),I=1,MIN6) + WRITE (6,635) (MIX(IT1+ISMZ(I)),IFR(IT0+ISMZ(I)),I=1,MIN6) + WRITE (6,640) (ALB(IT0+ISMZ(I)),I=1,MIN6) + WRITE (6,645) (DVX(IT1+ISMZ(I)),I=1,MIN6) + IF (NCOUR.EQ.12) THEN + WRITE (6,660) (DIRZ(I),I=7,12) + WRITE (6,635) (MIX(IT1+ISMZ(I)),IFR(IT0+ISMZ(I)),I=7,12) + WRITE (6,640) (ALB(IT0+ISMZ(I)),I=7,12) + WRITE (6,645) (DVX(IT1+ISMZ(I)),I=7,12) + ENDIF + WRITE (6,650) ('----------------',I=1,MIN6) +560 CONTINUE + ENDIF +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(GG3,JF2) + RETURN +* +570 WRITE (6,620) + WRITE (6,650) ('----------------',I=1,NCOUR) + DO 580 IBLK=1,NMBLK + IKK=INUM(IBLK) + WRITE (6,630) IBLK,IKK,IGEN(IKK) + I1=IORI(IBLK) + IT0=NCOUR*(IBLK-1) + IT1=NCOUR*(IKK-1) + WRITE (6,660) (DIRR(I),I=1,NCOUR) + WRITE (6,635) (MIX(IT1+ISM(I,I1)),IFR(IT0+ISM(I,I1)),I=1,NCOUR) + WRITE (6,640) (ALB(IT0+ISM(I,I1)),I=1,NCOUR) + WRITE (6,645) (DVX(IT1+ISM(I,I1)),I=1,NCOUR) + WRITE (6,650) ('----------------',I=1,NCOUR) +580 CONTINUE + CALL XABORT(HSMG) +* +620 FORMAT (///31H SURFACE NUMBERS AFTER MERGING./) +630 FORMAT (7H BLOCK=,I5,5X,13HMERGED BLOCK=,I5,5X,12HGENERATING B, + 1 5HLOCK=,I5) +635 FORMAT (8H IN/OUT:,6(I6,2H /,I5,3H I)) +640 FORMAT (8H ALBEDO:,1P,6(E13.5,3H I)) +645 FORMAT (8H DVX:,1P,6(E13.5,3H I)) +650 FORMAT (8H -------,6(A16)) +660 FORMAT (/8H SIDE:,6(A9,6X,1HI)) +680 FORMAT (49HNUMER3: GEOMETRICAL RECIPROCITY CONDITION IS VIOL, + 1 10HATED (SIDE,A3,20H OF GENERATING BLOCK,I5,8H VS SIDE,A3, + 2 20H OF GENERATING BLOCK,I5,3H ).) + END -- cgit v1.2.3