Initial commit from Polytechnique Montreal

1 files changed, 304 insertions, 0 deletions

diff --git a/Donjon/src/RESGEO.f b/Donjon/src/RESGEO.f
new file mode 100644
index 0000000..0c168f7
--- /dev/null
+++ b/Donjon/src/RESGEO.f
@@ -0,0 +1,304 @@
+*DECK RESGEO
+      SUBROUTINE RESGEO(IPMAP,IPMTX,LX,LY,LZ,NFUEL,IMPX,IGEO,NX,NY,NZ,
+     1 NCH,NB,NTOT,LNAP,IPCPO)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Create and check the fuel-map geometry.
+*
+*Copyright:
+* Copyright (C) 2007 Ecole Polytechnique de Montreal.
+*
+*Author(s): 
+* E. Varin, D. Sekki and V. Descotes
+*
+*Update(s):
+* R. Chambon 2014
+*
+*Parameters: input
+* IPMAP  pointer to fuel-map information.
+* IPMTX  pointer to matex information.
+* LX     number of elements along x-axis in geometry.
+* LY     number of elements along y-axis in geometry.
+* LZ     number of elements along z-axis in geometry.
+* NFUEL  number of fuel types.
+* IMPX   printing index (=0 for no print).
+* IGEO   type of geometry (=7 or =9)
+*
+*Parameters: output
+* NX     number of elements along x-axis in fuel map.
+* NY     number of elements along y-axis in fuel map.
+* NZ     number of elements along z-axis in fuel map.
+* NCH    number of reactor channels.
+* NB     number of fuel bundles per channel.
+* NTOT   total number of fuel bundles.
+* LNAP   Flag to call NAP: module to unfold geometry at assembly level
+* IPCPO  pointer to multicompo information
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPMAP,IPMTX,IPCPO,IPGNW
+      INTEGER LX,LY,LZ,NFUEL,IGEO,NX,NY,NZ,NCH,NB,NTOT
+      LOGICAL LNAP
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER(NSTATE=40,IOUT=6,EPSI=1.0E-4)
+      INTEGER ISTATE(NSTATE),JENT(1),IENT(1),JENT2(3),IENT2(3),NCODE(6),
+     1 ICODE(6)
+      TYPE(C_PTR) KENT(1),KENT2(3)
+      REAL    GEOXX(LX+1),GEOYY(LY+1),GEOZZ(LZ+1),GEOSI,GMAPSI,ZCODE(6)
+      CHARACTER HENT(1)*12,HENT2(3)*12,TEXT*12
+      DOUBLE PRECISION DFLOT
+*----
+*  ALLOCATABLE STATEMENTS
+*----
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: ISPLX,ISPLY,ISPLZ,MAT
+      REAL, ALLOCATABLE, DIMENSION(:) :: GMAPX,GMAPY,GMAPZ
+*----
+*  FUEL-MAP GEOMETRY
+*----
+      IF(IMPX.GT.1)WRITE(IOUT,*)'** CREATING FUEL-MAP GEOMETRY **'
+      CALL LCMSIX(IPMAP,'GEOMAP',1)
+      NENT=1
+      JENT(1)=0
+      HENT(1)='GEOMAP'
+      IENT(1)=1
+      KENT(1)=IPMAP
+      CALL GEOD(NENT,HENT,IENT,JENT,KENT)
+      IF(IMPX.GT.3)CALL LCMLIB(IPMAP)
+*
+      IF(LNAP) THEN
+*----
+*  FUEL-MAP GEOMETRY UNFOLDING WITH NAP:
+*----
+        IF(.NOT.C_ASSOCIATED(IPCPO)) THEN
+          CALL XABORT('RESGEO: COMPO LCM OBJECT MISSING AT RHS.')
+        ENDIF
+        CALL REDGET(ITYP,NITMA,FLOT,TEXT,DFLOT)
+        IF(ITYP.NE.3)CALL XABORT('@RESGEO: CHARACTER DATA EXPECTED.')
+        IF(TEXT.NE.':::') CALL XABORT('@RESGEO: ::: keyword EXPECTED.')
+        CALL REDGET(ITYP,NITMA,FLOT,TEXT,DFLOT)
+        IF(ITYP.NE.3)CALL XABORT('@RESGEO: CHARACTER DATA EXPECTED.')
+        IF(TEXT.NE.'NAP:') CALL XABORT('@RESGEO: NAP: keyword '
+     1   //'EXPECTED.')
+        CALL LCMOP(IPGNW,'GEONEW',0,1,0)
+        CALL LCMSIX(IPMAP,' ',0)
+        CALL LCMSIX(IPMAP,'GEOMAP',1)
+        NENT2=3
+        JENT2(1)=0
+        JENT2(2)=2
+        JENT2(3)=2
+        HENT2(1)='GEONEW'
+        HENT2(1)='GEOOLD'
+        HENT2(1)='COMPO'
+        IENT2(1)=1
+        IENT2(2)=1
+        IENT2(3)=1
+        KENT2(1)=IPGNW
+        KENT2(2)=IPMAP
+        KENT2(3)=IPCPO
+        CALL NAP(NENT2,HENT2,IENT2,JENT2,KENT2)
+        CALL LCMSIX(IPMAP,' ',0)
+        IF(IMPX.GT.3)CALL LCMLIB(IPMAP)
+        CALL LCMDEL(IPMAP,'GEOMAP')
+        IF(IMPX.GT.3)CALL LCMLIB(IPMAP)
+        CALL LCMSIX(IPMAP,'GEOMAP',1)
+        IF(IMPX.GT.3)CALL LCMLIB(IPMAP)
+        CALL LCMEQU(IPGNW,IPMAP)
+        IF(IMPX.GT.3)CALL LCMLIB(IPMAP)
+        CALL LCMCL(IPGNW,1)
+      ENDIF
+****
+      CALL LCMSIX(IPMAP,' ',0)
+      IF(IMPX.GT.3)CALL LCMLIB(IPMAP)
+      CALL LCMSIX(IPMAP,'GEOMAP',1)
+      IF(IMPX.GT.3)CALL LCMLIB(IPMAP)
+****
+      ISTATE(:NSTATE)=0
+      CALL LCMGET(IPMAP,'STATE-VECTOR',ISTATE)
+      IF(ISTATE(1).NE.IGEO) CALL XABORT('@RESGEO: THE GEOMETRY '
+     1 // 'IN FUEL-MAP MUST HAVE THE SAME TYPE AS IN THE MATEX-OBJECT')
+      IGEO=ISTATE(1)
+      NX=ISTATE(3)
+      NY=ISTATE(4)
+      NZ=ISTATE(5)
+*----
+*  READ FUEL-MAP GEOMETRY AND PERFORM MESH-SPLITTING
+*----
+      IMPX0=MAX(0,IMPX-1)
+      NX2=NX
+      NY2=NY
+      IF(IGEO.GE.8) NY2=1
+      NZ2=NZ
+      ALLOCATE(ISPLX(NX2),ISPLY(NY2),ISPLZ(NZ2))
+      ISPLTL=0
+      ISPLTH=0
+      IHEX=0
+      CALL LCMLEN(IPMAP,'SPLITL',ILEN,ITYLCM)
+      IF(ILEN.GT.0) CALL LCMGET(IPMAP,'SPLITL',ISPLTL)
+      CALL LCMLEN(IPMAP,'SPLITH',ILEN,ITYLCM)
+      IF(ILEN.GT.0) CALL LCMGET(IPMAP,'SPLITH',ISPLTH)
+      IF(IGEO.LT.8) THEN
+        CALL LCMLEN(IPMAP,'SPLITX',ILONG,ITYLCM)
+        IF(ILONG.GT.0) THEN
+          CALL LCMGET(IPMAP,'SPLITX',ISPLX)
+          NX2=0
+          DO IOLD=1,NX
+            NX2=NX2+ISPLX(IOLD)
+          ENDDO
+        ENDIF
+        CALL LCMLEN(IPMAP,'SPLITY',ILONG,ITYLCM)
+        IF(ILONG.GT.0) THEN
+          CALL LCMGET(IPMAP,'SPLITY',ISPLY)
+          NY2=0
+          DO IOLD=1,NY
+            NY2=NY2+ISPLY(IOLD)
+          ENDDO
+        ENDIF
+      ELSEIF((ISPLTH.NE.0).AND.((IGEO.EQ.8).OR.(IGEO.EQ.9))) THEN
+        NX2=NX*6*(ISPLTH**2)
+        CALL LCMGET(IPMAP,'IHEX',IHEX)
+      ELSEIF((ISPLTL.NE.0).AND.((IGEO.EQ.8).OR.(IGEO.EQ.9))) THEN
+        NX2=NX*3*(ISPLTL**2)
+        CALL LCMGET(IPMAP,'IHEX',IHEX)
+      ENDIF
+      CALL LCMLEN(IPMAP,'SPLITZ',ILONG,ITYLCM)
+      IF(ILONG.GT.0) THEN
+        CALL LCMGET(IPMAP,'SPLITZ',ISPLZ)
+        NZ2=0
+        DO IOLD=1,NZ
+          NZ2=NZ2+ISPLZ(IOLD)
+        ENDDO
+      ENDIF
+      MAXPTS=NX2*NY2*NZ2
+      MAXX=NX2
+      IF(IHEX.EQ.1) THEN
+        MAXPTS=12*MAXPTS
+        MAXX=12*MAXX
+      ELSE IF((IHEX.EQ.2).OR.(IHEX.EQ.3)) THEN
+        MAXPTS=6*MAXPTS
+        MAXX=6*MAXX
+      ELSE IF(IHEX.EQ.4) THEN
+        MAXPTS=4*MAXPTS
+        MAXX=4*MAXX
+      ELSE IF(IHEX.EQ.5) THEN
+        MAXPTS=3*MAXPTS
+        MAXX=3*MAXX
+      ELSE IF((IHEX.GE.6).AND.(IHEX.LE.8)) THEN
+        MAXPTS=2*MAXPTS
+        MAXX=2*MAXX
+      ENDIF
+      ALLOCATE(MAT(MAXPTS),GMAPX(MAXX+1),GMAPY(NY2+1),GMAPZ(NZ2+1))
+      CALL READ3D(MAXX,NY2,NZ2,MAXPTS,IPMAP,IHEX,IR,ILK,SIDE,GMAPX,
+     1 GMAPY,GMAPZ,IMPX0,NX2,NY2,NZ2,MAT,NEL,NCODE,ICODE,ZCODE,ISPLX,
+     2 ISPLY,ISPLZ,ISPLH,ISPLL)
+      IF((NEL.NE.NX2*NY2*NZ2).AND.(IHEX.EQ.0))CALL XABORT('@RESGEO: WR'
+     1 // 'ONG GEOMETRY.')
+      IF((NEL.NE.NX2*NZ2).AND.(IHEX.NE.0))CALL XABORT('@RESGEO: WRONG ' 
+     1 // 'HEXAGONAL GEOMETRY, WRONG NUMBER OF ELEMENTS.')
+      DEALLOCATE(MAT,ISPLZ,ISPLY,ISPLX)
+      IF(IMPX.GT.2)WRITE(IOUT,*)'CHECKING FUEL-MAP GEOMETRY'
+      IF((IGEO.NE.7).AND.(IGEO.NE.9))CALL XABORT('@RESGEO: ONLY '
+     1 //'3D-CARTESIAN OR 3D-HEXAGONAL GEOMETRY ALLOWED.')
+      IF(IHEX.EQ.0) THEN
+        IF((LX.LT.NX).OR.(LY.LT.NY).OR.(LZ.LT.NZ)) THEN
+        WRITE(IOUT,*) 'Geometry LX=',LX,', LY=',LY,' and LZ=',LZ, 
+     1   ' must be greater or equal to map ',
+     2   'NX=',NX,' NY=',NY,' and NZ=',NZ
+          CALL XABORT('@RESGEO: WRONG GEOMETRY DEFINITION.')
+        ENDIF  
+      ELSE
+        IF((LX.LT.NX).OR.(LZ.LT.NZ)) THEN
+        WRITE(IOUT,*) 'Geometry LX=',LX,' and LZ=',LZ, 
+     1   ' must be greater or equal to map ',
+     2   'NX=',NX,' and NZ=',NZ
+          CALL XABORT('@RESGEO: WRONG GEOMETRY DEFINITION.')
+        ENDIF  
+      ENDIF 
+      IF(NZ.LT.NB)THEN
+        WRITE(IOUT,*)'@RESGEO: FOUND NZ=',NZ,' LESS THAN NB=',NB
+        CALL XABORT('@RESGEO: WRONG FUEL-MAP GEOMETRY DEFINITION.')
+      ENDIF
+*----
+*  CHECK MESHX OR SIDE
+*----
+      IF(IGEO.EQ.7) THEN
+        GEOXX(:LX+1)=0.0
+        CALL LCMGET(IPMTX,'MESHX',GEOXX)
+        DO 10 IMP=1,NX+1
+        DO IGM=1,LX+1
+          IF(ABS(GMAPX(IMP)-GEOXX(IGM)).LT.EPSI)THEN
+            GEOXX(IGM)=GMAPX(IMP)
+            GOTO 10
+          ENDIF
+        ENDDO
+        WRITE(IOUT,*)'@RESGEO: MESHX IN L_MAP ',GMAPX(IMP)
+        CALL XABORT('@RESGEO: UNABLE TO FIND THIS MESHX IN L_GEOM.')
+   10   CONTINUE
+        CALL LCMPUT(IPMTX,'MESHX',LX+1,2,GEOXX)
+      ELSE IF(IGEO.EQ.9) THEN
+        ISPLTL=0
+        NY=1
+        CALL LCMGET(IPMAP,'SIDE',GMAPSI)
+        CALL LCMLEN(IPMAP,'SPLITL',ILONG,ITYLCM)
+        IF(ILONG.GT.0) CALL LCMGET(IPMAP,'SPLITL',ISPLTL)
+        IF(ISPLTL.EQ.0) ISPLTL=1
+        GMAPSI=GMAPSI/REAL(ISPLTL)
+        CALL LCMGET(IPMTX,'SIDE',GEOSI)
+        IF(ABS(GMAPSI-GEOSI).LT.EPSI)THEN
+          GEOSI=GMAPSI
+          GOTO 20
+        ENDIF
+        WRITE(IOUT,*)'@RESGEO: SIDE IN L_MAP ',GMAPSI, GEOSI
+        CALL XABORT('@RESGEO: UNABLE TO FIND THIS SIDE IN L_GEOM.')
+   20   CONTINUE
+        CALL LCMPUT(IPMTX,'SIDE',1,2,GEOSI)
+      ENDIF
+*----
+*  CHECK MESHY (ONLY IF 3D-CARTESIAN GEOMETRY)
+*----
+      IF(IGEO.EQ.7) THEN
+        GEOYY(:LY+1)=0.0
+        CALL LCMGET(IPMTX,'MESHY',GEOYY)
+        DO 30 IMP=1,NY+1
+        DO IGM=1,LY+1
+          IF(ABS(GMAPY(IMP)-GEOYY(IGM)).LT.EPSI)THEN
+            GEOYY(IGM)=GMAPY(IMP)
+            GOTO 30
+          ENDIF
+        ENDDO
+        WRITE(IOUT,*)'@RESGEO: MESHY IN FUEL MAP ',GMAPY(IMP)
+        CALL XABORT('@RESGEO: UNABLE TO FIND THIS MESHY IN L_GEOM.')
+   30   CONTINUE
+        CALL LCMPUT(IPMTX,'MESHY',LY+1,2,GEOYY)
+      ENDIF
+*----
+*  CHECK MESHZ
+*----
+      GEOZZ(:LZ+1)=0.0
+      CALL LCMGET(IPMTX,'MESHZ',GEOZZ)
+      DO 50 IMP=1,NZ+1
+      DO IGM=1,LZ+1
+        IF(ABS(GMAPZ(IMP)-GEOZZ(IGM)).LT.EPSI)THEN
+          GEOZZ(IGM)=GMAPZ(IMP)
+          GOTO 50
+        ENDIF
+      ENDDO
+      WRITE(IOUT,*)'@RESGEO: MESHZ IN FUEL MAP ',GMAPZ(IMP)
+      CALL XABORT('@RESGEO: UNABLE TO FIND THIS MESHZ IN L_GEOM.')
+   50 CONTINUE
+      CALL LCMPUT(IPMTX,'MESHZ',LZ+1,2,GEOZZ)
+      DEALLOCATE(GMAPZ,GMAPY,GMAPX)
+*----
+*  CHECK FUEL MIXTURES
+*----
+      CALL RESPFM(IPMAP,IPMTX,NX,NY,NZ,LX,LY,LZ,NFUEL,IMPX,IGEO,NCH,NB,
+     1 NTOT)
+      RETURN
+      END