Initial commit from Polytechnique Montreal

1 files changed, 182 insertions, 0 deletions

diff --git a/Donjon/src/RODMOD.f b/Donjon/src/RODMOD.f
new file mode 100644
index 0000000..25cd8e9
--- /dev/null
+++ b/Donjon/src/RODMOD.f
@@ -0,0 +1,182 @@
+*DECK RODMOD
+      SUBROUTINE RODMOD(IPMAP,NGRP,MAXMIX,NCH,NB,RNAME,INS,INSS,
+     > HMIX,RMIX,NREB,RODSIZE,INTROD,INI,NUMMIX,NCALL)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Modification of the rod data stored in the PARAM folder of a fuel map
+*
+*Copyright:
+* Copyright (C) 2017 Ecole Polytechnique de Montreal.
+*
+*Author(s): 
+* G. Tixier
+*
+*Parameters: input
+* IPMAP   pointer to the fuel map
+* NGRP    number of rod groups
+* MAXMIX  maximum number of rod mix
+* NCH     number of fuel channels
+* NB      number of fuel bundles per channel
+* RNAME   name of rod group
+* INS     rod insertion for each rod group
+* INSS    rod insertion step (in cm)
+* HMIX    height of rod mix (if more than 2 rod mix are defined)
+* RMIX    number associated to rod mix
+* NREB    number of bottom-reflective meshes
+* RODSIZE total number of meshes for the fuel zone (=NCH*NB)
+* INTROD  final value of data stored in the PARAM folder after 
+*         calculation
+* INI     initial value of data stored in the PARAM folder (no rod
+*         inserted)
+* NUMMIX  number of rod mix for each rod group
+* NCALL   number to distinguish first or other calls to the ROD: 
+*         module
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER NGRP,MAXMIX,NREB,NCH,NB,RODSIZE,NUMMIX(NGRP)
+      INTEGER INSM,NCALL
+      REAL HMIX(NGRP*MAXMIX),RMIX(NGRP*MAXMIX),INTROD(RODSIZE),INS(NGRP)
+      REAL INI,INSS
+      CHARACTER(LEN=3) RNAME(NGRP),RASS(NCH)
+      TYPE(C_PTR) IPMAP
+*----
+*  LOCAL VARIABLES
+*----
+      CHARACTER TEXT*3
+      INTEGER I,J,K,L,M
+      REAL FLOT,F1,F2,ICM,HMB,HMT,ZMI1,ZMI2
+      DOUBLE PRECISION DFLOT
+      TYPE(C_PTR) JPMAP,MPMAP
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      REAL, ALLOCATABLE, DIMENSION(:) :: GEOZZ
+
+* Recover axial meshing
+      JPMAP=LCMGID(IPMAP,'GEOMAP')
+      CALL LCMLEN(JPMAP,'MESHZ',ILONG,ITYLCM)
+      ALLOCATE(GEOZZ(ILONG))
+      CALL LCMGET(JPMAP,'MESHZ',GEOZZ)
+
+* Recover rod parameters
+      MPMAP=LCMGID(IPMAP,'ROD-INFO')
+      CALL LCMGET(MPMAP,'ROD-INIT',INI)
+      CALL LCMGET(MPMAP,'INS-MAX',INSM)
+      CALL LCMGET(MPMAP,'STEP-CM',INSS)
+      CALL LCMGET(MPMAP,'REFL-BOTTOM',NREB)      
+      CALL LCMGTC(MPMAP,'ROD-NAME',3,NGRP,RNAME)
+      CALL LCMGET(MPMAP,'ROD-INSERT',INS)
+      CALL LCMGET(MPMAP,'ROD-RIN',RMIX)
+      CALL LCMGET(MPMAP,'ROD-NBZONE',NUMMIX)
+      CALL LCMGET(MPMAP,'ROD-HEIGHT',HMIX)
+      IF(NB+NREB+1.GT.ILONG) CALL XABORT('RODMOD: GEOZZ OVERFLOW.')
+
+      IF(NCALL.EQ.1) THEN
+* If it is the first call to the ROD: module, the rod map is stored in
+* the fuel map
+      N=1
+      CALL REDGET(ITYP,NITMA,FLOT,TEXT,DFLOT)
+      IF(ITYP.NE.3)CALL XABORT('@RODMOD: CHARACTER DATA EXPECTED.')
+      DO WHILE(N.LE.NCH)
+        RASS(N)=TEXT
+        CALL REDGET(ITYP,NITMA,FLOT,TEXT,DFLOT)
+        N=N+1
+      END DO
+      CALL LCMPTC(MPMAP,'ROD-MAP',3,NCH,RASS)
+      ENDIF
+* Recover rod map
+      CALL LCMGTC(MPMAP,'ROD-MAP',3,NCH,RASS)
+      I=1
+      DO WHILE(I.LE.NCH)
+        K=NB
+        DO WHILE(K.GE.1)
+        IF((RASS(I).EQ.'|').OR.(RASS(I).EQ.'-').OR.(RASS(I).EQ.'-|-'
+     >  )) THEN
+* If no control rod is defined
+          INTROD((K-1)*NCH+I)=INI
+          GO TO 10
+        ELSE
+* Recover control rod information
+          J=1
+          DO WHILE(J.LE.NGRP)
+            IF(RASS(I).EQ.RNAME(J)) THEN
+              EXIT
+            ELSE
+              J=J+1
+            ENDIF 
+          END DO
+          M=NUMMIX(J)
+          DO WHILE(M.GE.1)
+            IF(M.EQ.1) THEN
+              ICM=INS(J)*INSS
+            ELSE
+              ICM=INS(J)*INSS-HMIX(J+(M-2)*NGRP)
+            ENDIF
+* Mesh size calculations
+            HMB=GEOZZ(NB+NREB+1)-GEOZZ(K+NREB)
+            HMT=GEOZZ(NB+NREB+1)-GEOZZ(K+NREB+1)
+            IF(ICM.LT.0.0) THEN
+              M=M-1
+            ELSE
+              ZMI1=HMB-ICM
+              IF(ZMI1.LE.0.0) THEN
+                INTROD((K-1)*NCH+I)=RMIX(J+(M-1)*NGRP)
+                EXIT
+              ELSE
+                ZMI2=HMT-ICM
+                IF(ZMI2.GE.0.0) THEN
+                  IF(M.EQ.1) THEN
+                    INTROD((K-1)*NCH+I)=INI
+                    GO TO 10
+                  ELSE
+                    M=M-1
+                  ENDIF
+                ELSE
+* Calculation of the proportion of control rod for the mesh considered
+                    F1=HMB-ICM
+                    F2=GEOZZ(K+NREB+1)-GEOZZ(K+NREB)-F1
+                  IF(M.EQ.1) THEN
+                    INTROD((K-1)*NCH+I)=(F2*RMIX(J+(M-1)*NGRP)
+     >              +F1*INI)/(F1+F2)
+                    GO TO 10
+                  ELSE
+                    IF(M.GE.3) THEN
+* It is not possible to have two interfaces in one mesh
+                    IF(((INS(J)*INSS-HMIX(J+(M-3)*NGRP)).GE.HMT).AND.
+     >              ((INS(J)*INSS-HMIX(J+(M-3)*NGRP)).LE.HMB))
+     >              CALL XABORT('@RODMOD: IT IS NOT POSSIBLE TO HAVE'
+     1   //' A MESH WITH MORE THAN TWO INTERFACES. HMIX HAS TO BE'
+     1   //' BIGGER')
+                    ENDIF
+* If two mixtures fill one mesh, we have to compute a fraction of
+* insertion for each rod mixture
+                    INTROD((K-1)*NCH+I)=(F2*RMIX(J+(M-1)*NGRP)+
+     >              F1*RMIX(J+(M-2)*NGRP))/(F1+F2)
+                    EXIT
+                  ENDIF
+                ENDIF
+              ENDIF
+            ENDIF
+          END DO
+        ENDIF
+        K=K-1
+        END DO
+   10   L=K-1
+* If a control rod does not fill one mesh entirely, the meshes located
+* at the bottom of this mesh do not contain control rod
+        DO WHILE(L.GE.1)
+          INTROD((L-1)*NCH+I)=INI
+          L=L-1
+        END DO
+      I=I+1
+      END DO
+      DEALLOCATE(GEOZZ)
+      RETURN
+      END