summaryrefslogtreecommitdiff
path: root/Donjon/src/THM.f
diff options
context:
space:
mode:
authorstainer_t <thomas.stainer@oecd-nea.org>2025-09-08 13:48:49 +0200
committerstainer_t <thomas.stainer@oecd-nea.org>2025-09-08 13:48:49 +0200
commit7dfcc480ba1e19bd3232349fc733caef94034292 (patch)
tree03ee104eb8846d5cc1a981d267687a729185d3f3 /Donjon/src/THM.f
Initial commit from Polytechnique Montreal
Diffstat (limited to 'Donjon/src/THM.f')
-rw-r--r--Donjon/src/THM.f1489
1 files changed, 1489 insertions, 0 deletions
diff --git a/Donjon/src/THM.f b/Donjon/src/THM.f
new file mode 100644
index 0000000..4917d83
--- /dev/null
+++ b/Donjon/src/THM.f
@@ -0,0 +1,1489 @@
+*DECK THM
+ SUBROUTINE THM(NENTRY,HENTRY,IENTRY,JENTRY,KENTRY)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Simplified thermal-hydraulics module.
+*
+*Copyright:
+* Copyright (C) 2013 Ecole Polytechnique de Montreal.
+*
+*Author(s):
+* A. Hebert, P. Gallet and V. Salino
+* 02/2025: C. HUET - Modifications to include pressure drop calculation
+* 08/2025: M.Bellier & R. Guasch modifified to include :
+* - drift-flux model
+* - variable axial properties
+*
+*
+*Parameters: input
+* NENTRY number of data structures transfered to this module.
+* HENTRY name of the data structures.
+* IENTRY data structure type where:
+* IENTRY=1 for LCM memory object;
+* IENTRY=2 for XSM file;
+* IENTRY=3 for sequential binary file;
+* IENTRY=4 for sequential ASCII file.
+* JENTRY access permission for the data structure where:
+* JENTRY=0 for a data structure in creation mode;
+* JENTRY=1 for a data structure in modifications mode;
+* JENTRY=2 for a data structure in read-only mode.
+* KENTRY data structure pointer.
+*
+*Comments:
+* The THM: module specification is:
+* THERMO MAPFL := THM: [ THERMO ] MAPFL :: (descthm) ;
+* where
+* THERMO : name of the \emph{thermo) object that will be created or updated
+* by the THM: module. Object \emph{thermo} contains thermal-hydraulics
+* information set or computed by THM: in transient or in permanent
+* conditions such as the distribution of the enthalpy, the pressure, the
+* velocity, the density and the temperatures of the coolant for all the
+* channels in the geometry. It also contains all the values of the fuel
+* temperatures in transient or in permanent conditions according to the
+* discretisation chosen for the fuel rods.
+* MAPFL : name of the \emph{map} object containing fuel regions description
+* and local parameter informations.
+* (descthm) : structure describing the input data to the THM: module.
+*
+*-----------------------------------------------------------------------
+*
+ USE GANLIB
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ INTEGER NENTRY,IENTRY(NENTRY),JENTRY(NENTRY)
+ TYPE(C_PTR) KENTRY(NENTRY)
+ CHARACTER HENTRY(NENTRY)*12
+*----
+* LOCAL VARIABLES
+*----
+ PARAMETER(NSTATE=40,IOUT=6,PI=3.141592654,ZKILO=1.0E3)
+ PARAMETER(IMAXO=1000,JMAXO=100,KMAXO=200,NMAXO=40,MAXRAD=10)
+ PARAMETER(DTEMPR=5.0,DTEMPT=40.0,DPRESS=4.0)
+ CHARACTER TEXT*40,TEXT12*12,HSIGN*12,PNAME*12,TXTDIR*12,HSMG*131,
+ > UCONDF*12,UCONDC*12,SNAME*32,SCOMP*32,FNAME*32,FCOMP*32
+ INTEGER ISTATE(NSTATE),TIMEIT,ITIME
+ REAL STATE(NSTATE),DTIME,KHGAP,KHCONV,WTEFF
+ REAL POULET,HX(IMAXO),HY(JMAXO),HZ(KMAXO)
+ REAL RPRAD(MAXRAD),FPRAD(MAXRAD),TERP(MAXRAD)
+ DOUBLE PRECISION DFLOT,DSUM
+ LOGICAL LPRAD
+ TYPE(C_PTR) IPTHM,IPMAP,JPMAP,KPMAP,JPTHM,KPTHM,LPTHM,MPTHM,
+ > KPTHMI,LPTHMI,MPTHMI
+*----
+* ALLOCATABLE ARRAYS
+*----
+ INTEGER, ALLOCATABLE, DIMENSION(:) :: NUM,IREFSC
+ REAL, ALLOCATABLE, DIMENSION(:) :: XX,YY,ZZ,BURN,BURN2,PW,FRO,
+ 1 FNFUCST,FNTGCST,FRACPU
+ REAL, ALLOCATABLE, DIMENSION(:) :: FPOWER,KCONDF,KCONDC,TIMESR,
+ 1 TPOWER,PFORM,DTERP
+ REAL, ALLOCATABLE, DIMENSION(:,:,:) :: XBURN,POW,TCOMB,DCOOL,
+ 1 TCOOL,TSURF,PCOOL,HCOOL
+ REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: RAD
+ DOUBLE PRECISION ARF,ARCI,ARCE,DARF,DARC
+ REAL, ALLOCATABLE, DIMENSION(:,:) :: VAL,RVAL
+ REAL, ALLOCATABLE, DIMENSION(:) :: ACOOL,PCH,HD,FFUEL,FCOOL
+ REAL, ALLOCATABLE, DIMENSION(:) :: RAPCOOL,RAPFUEL,PM
+ REAL, ALLOCATABLE, DIMENSION(:,:) :: FNFU,FNTG
+*----
+* PARAMETER VALIDATION
+*----
+ IF(NENTRY.NE.2)CALL XABORT('@THM: 2 PARAMETERS EXPECTED.')
+ IPTHM=KENTRY(1)
+ IPMAP=KENTRY(2)
+ IF((IENTRY(1).NE.1).AND.(IENTRY(1).NE.2))CALL XABORT('@THM:'
+ 1 //' LCM OBJECT EXPECTED AT FIRST LHS.')
+ IF((IENTRY(2).NE.1).AND.(IENTRY(2).NE.2))CALL XABORT('@THM:'
+ 1 //' LCM OBJECT EXPECTED AT SECOND LHS.')
+ CALL LCMGTC(IPMAP,'SIGNATURE',12,HSIGN)
+ IF(HSIGN.NE.'L_MAP')THEN
+ TEXT=HENTRY(2)
+ CALL XABORT('@THM: SIGNATURE OF '//TEXT//' IS '//HSIGN//
+ 1 '. L_MAP EXPECTED.')
+ ENDIF
+*----
+* RECOVER L_MAP STATE-VECTOR
+*----
+ CALL LCMGET(IPMAP,'STATE-VECTOR',ISTATE)
+ NB=ISTATE(1)
+ NCH=ISTATE(2)
+ NPARM=ISTATE(8)
+ NSIMS=ISTATE(13)
+ ALLOCATE(FNFUCST(NCH),FNTGCST(NCH),FRACPU(NCH))
+*----
+* READ DATA
+*----
+ IMPX=1
+ ITIME=0
+ DTIME=0.0
+ FPUISS=0.974
+ CFLUX=2.0E+6
+ SPEED=0.0
+ TINLET=0.0
+ POULET=0.0
+ POROS=0.05
+ ICONDF=0
+ ICONDC=0
+ IHGAP=0
+ IHCONV=0
+ IFRCDI=0
+ ISUBM=1
+ RC=0.0
+ RIG=0.0
+ RGG=0.0
+ RTG=0.0
+ PITCH=0.0
+ MAXIT1=50
+ MAXIT2=50
+ MAXIT3=50
+ IFLUID=0
+ IFUEL=0
+ IGAP=0
+ IPRES=0
+ IDFM=0
+ ERMAXT=1.0
+ ERMAXC=1.0E-3
+ NFD=5
+ NDTOT=8
+ NPRAD=0
+ TIMEIT=0
+ NPOWER=0
+ RELAX=1.0
+ RTIME=0.0
+ WTEFF=5.0/9.0 ! Rowlands weighting factor
+ EPSR=0.0
+ THETA=0.0
+ UCONDF='CELSIUS'
+ UCONDC='CELSIUS'
+ TPOW=0.0
+ FNFUCST(:NCH)=1.0
+ FNTGCST(:NCH)=0.0
+ FRACPU(:NCH)=0.0
+ IF(JENTRY(1).EQ.1) THEN
+ CALL LCMGET(IPTHM,'STATE-VECTOR',ISTATE)
+ IF(ISTATE(1).NE.NCH) CALL XABORT('THM: INVALID STATE VECTOR FO'
+ > //'R IPTHM OPJECT.')
+ MAXIT1=ISTATE(3)
+ MAXIT2=ISTATE(4)
+ MAXIT3=ISTATE(5)
+ NFD=ISTATE(6)
+ NDTOT=ISTATE(7)
+ ITIME=ISTATE(8)
+ TIMEIT=ISTATE(9)
+ IHGAP=ISTATE(10)
+ IHCONV=ISTATE(11)
+ ICONDF=ISTATE(12)
+ ICONDC=ISTATE(13)
+ IFRCDI=ISTATE(14)
+ ISUBM=ISTATE(15)
+ IF(ICONDF.EQ.1) NCONDF=ISTATE(16)
+ IF(ICONDC.EQ.1) NCONDC=ISTATE(17)
+ NPRAD=ISTATE(18)
+ IFLUID=ISTATE(20)
+ IGAP=ISTATE(21)
+ IPRES=ISTATE(22)
+ CALL LCMGET(IPTHM,'REAL-PARAM',STATE)
+ DTIME=STATE(1)
+ FPUISS=STATE(2)
+ CFLUX=STATE(3)
+ SPEED=STATE(4)
+ POULET=STATE(5)
+ TINLET=STATE(6)
+ POROS=STATE(7)
+ RC=STATE(8)
+ RIG=STATE(9)
+ RGG=STATE(10)
+ RTG=STATE(11)
+ PITCH=STATE(12)
+ ERMAXT=STATE(13)
+ ERMAXC=STATE(14)
+ RELAX=STATE(15)
+ RTIME=STATE(16)
+ IF(IHGAP.EQ.1) KHGAP=STATE(17)
+ IF(IHCONV.EQ.1) KHCONV=STATE(18)
+ WTEFF=STATE(19)
+ TPOW=STATE(20)
+ EPSR=STATE(22)
+ THETA=STATE(23)
+*----
+* RECOVER CELL-DEPENDENT DATA
+*----
+ CALL LCMGET(IPTHM,'NB-FUEL',FNFUCST)
+ CALL LCMGET(IPTHM,'NB-TUBE',FNTGCST)
+ CALL LCMGET(IPTHM,'FRACT-PU',FRACPU)
+*----
+* RECOVER CONDUCTIVITY INFORMATION ON LCM OBJECT THM
+*----
+ IF(ICONDF.EQ.1) THEN
+ ALLOCATE(KCONDF(NCONDF+3))
+ CALL LCMGET(IPTHM,'KCONDF',KCONDF)
+ CALL LCMGTC(IPTHM,'UCONDF',12,UCONDF)
+ ENDIF
+ IF(ICONDC.EQ.1) THEN
+ ALLOCATE(KCONDC(NCONDC+1))
+ CALL LCMGET(IPTHM,'KCONDC',KCONDC)
+ CALL LCMGTC(IPTHM,'UCONDC',12,UCONDC)
+ ENDIF
+ ENDIF
+*----
+* READ INPUT DATA
+*----
+ IPICK=0
+ LPRAD=.FALSE.
+ 10 CALL REDGET(ITYP,NITMA,FLOT,TEXT,DFLOT)
+ IF(ITYP.EQ.10)GO TO 60
+ 20 IF(ITYP.NE.3)CALL XABORT('@THM: CHARACTER DATA EXPECTED.')
+ IF(TEXT.EQ.'EDIT') THEN
+* Read printing index
+ CALL REDGET(ITYP,IMPX,FLOT,TEXT,DFLOT)
+ IF(ITYP.NE.1)CALL XABORT('@THM: INTEGER FOR EDIT EXPECTED.')
+ ELSE IF(TEXT.EQ.'TIME') THEN
+* Time at beginning of time-step (s).
+ CALL REDGET(ITYP,NITMA,RTIME,TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR RTIME EXPECTED.')
+ CALL REDGET(ITYP,NITMA,FLOT,TEXT,DFLOT)
+ IF(ITYP.EQ.2) THEN
+ DTIME=FLOT
+ ELSE IF(ITYP.EQ.3) THEN
+ GO TO 20
+ ELSE
+ CALL XABORT('@THM: REAL FOR DTIME EXPECTED.')
+ ENDIF
+ ITIME=1
+ ELSE IF(TEXT.EQ.'FLUID') THEN
+* Read fluid type
+ CALL REDGET(ITYP,NITMA,FLOT,TEXT,DFLOT)
+ IF(ITYP.NE.3)CALL XABORT('@THM: CHARACTER FOR FLUID EXPECTED.')
+ IF(TEXT.EQ.'H2O') THEN
+ IFLUID=0
+ ELSE IF(TEXT.EQ.'D2O') THEN
+ IFLUID=1
+ ELSE IF(TEXT.EQ.'SALT') THEN
+ IFLUID=2
+ CALL REDGET(ITYP,NITMA,FLOT,SNAME,DFLOT)
+ IF(ITYP.NE.3) THEN
+ CALL XABORT('@THM: CHARACTER FOR FLUID SALT NAME EXPECTED'
+ > //'.')
+ ENDIF
+ CALL REDGET(ITYP,NITMA,FLOT,SCOMP,DFLOT)
+ IF(ITYP.NE.3) THEN
+ CALL XABORT('@THM: CHARACTER FOR FLUID SALT COMPOSITION'
+ > //'EXPECTED.')
+ ENDIF
+ ELSE
+ CALL XABORT('@THM: INVALID FLUID TYPE.')
+ ENDIF
+ ELSE IF(TEXT.EQ.'FUEL') THEN
+* Read fuel type
+ CALL REDGET(ITYP,NITMA,FLOT,TEXT,DFLOT)
+ IF(ITYP.NE.3)CALL XABORT('@THM: CHARACTER FOR FLUID EXPECTED.')
+ IF(TEXT.EQ.'UO2') THEN
+ IFUEL=0
+ ELSE IF(TEXT.EQ.'SALT') THEN
+ IFUEL=1
+ CALL REDGET(ITYP,NITMA,FLOT,FNAME,DFLOT)
+ IF(ITYP.NE.3) THEN
+ CALL XABORT('@THM: CHARACTER FOR FLUID EXPECTED.')
+ ENDIF
+ CALL REDGET(ITYP,NITMA,FLOT,FCOMP,DFLOT)
+ IF(ITYP.NE.3) THEN
+ CALL XABORT('@THM: CHARACTER FOR FLUID EXPECTED.')
+ ENDIF
+ ELSE
+ CALL XABORT('@THM: INVALID FUEL TYPE.')
+ ENDIF
+ ELSE IF(TEXT.EQ.'FPUISS') THEN
+* Coolant power factor
+ CALL REDGET(ITYP,NITMA,FLOT,TEXT,DFLOT)
+ IF(ITYP.EQ.2) THEN
+ FPUISS=FLOT
+ ELSE
+ CALL XABORT('@THM: REAL FOR FPUISS EXPECTED.')
+ ENDIF
+ ELSE IF(TEXT.EQ.'CRITFL') THEN
+* Critical heat flux (W/m^2)
+ CALL REDGET(ITYP,NITMA,CFLUX,TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR CFLUX EXPECTED.')
+ ELSE IF(TEXT.EQ.'CWSECT') THEN
+* Core coolant section (m^2)
+ CALL REDGET(ITYP,NITMA,CWSECT,TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR CWSECT EXPECTED.')
+* Coolant flow (m^3/h)
+ CALL REDGET(ITYP,NITMA,FLOW,TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR FLOW EXPECTED.')
+ SPEED=FLOW/(3600.0*CWSECT)
+ ELSE IF(TEXT.EQ.'INLET-Q') THEN
+* Core coolant section (m^2)
+ IF((POULET.EQ.0.0).OR.(TINLET.EQ.0.0)) CALL XABORT('@THM: INLE'
+ > //'T INFORMATION NOT SET BEFORE USING INLET-Q.')
+ CALL REDGET(ITYP,NITMA,CWSECT,TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR CWSECT EXPECTED.')
+* Inlet mass flow rate (kg/s)
+ CALL REDGET(ITYP,NITMA,QFLUID,TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR QFLUID EXPECTED.')
+ IF(IFLUID.EQ.0) THEN
+ CALL THMPT(POULET,TINLET,RHOL,R2,R3,R4,R5)
+ ELSE IF(IFLUID.EQ.1) THEN
+ CALL THMHPT(POULET,TINLET,RHOL,R2,R3,R4,R5)
+ ELSE IF(IFLUID.EQ.2) THEN
+ CALL THMSPT(SNAME,SCOMP,TINLET,RHOL,R2,R3,R4,R5,IMPX)
+ ENDIF
+ SPEED=QFLUID/(CWSECT*RHOL)
+ ELSE IF(TEXT.EQ.'SPEED') THEN
+* Coolant velocity (m/s)
+ CALL REDGET(ITYP,NITMA,SPEED,TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR SPEED EXPECTED.')
+ ELSE IF(TEXT.EQ.'INLET') THEN
+* The POULET and TINLET informations are used to compute initial
+* enthalpy and water density.
+* Outlet pressure (Pa)
+ CALL REDGET(ITYP,NITMA,POULET,TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR POULET EXPECTED.')
+* Inlet temperature (K)
+ CALL REDGET(ITYP,NITMA,TINLET,TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR TINLET EXPECTED.')
+ ELSE IF(TEXT.EQ.'PUFR') THEN
+ ICONDF=0
+* Plutonium mass enrichment
+ CALL THMINP('PUFR',NCH,FRACPU)
+ ELSE IF(TEXT.EQ.'POROS') THEN
+ ICONDF=0
+* Oxyde porosity
+ CALL REDGET(ITYP,NITMA,POROS,TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR POROS EXPECTED.')
+ ELSE IF(TEXT.EQ.'CONDF') THEN
+ IF(ICONDF.EQ.1)DEALLOCATE(KCONDF)
+ ICONDF=1
+* Fuel conductivity expressed as a function of fuel temperature
+* (function = polynomial + inverse term)
+ CALL REDGET(ITYP,NCONDF,FLOT,TEXT,DFLOT)
+ IF(ITYP.NE.1)CALL XABORT('@THM: INTEGER FOR CONDF EXPECTED.')
+ IF(NCONDF.LT.0)CALL XABORT('@THM: NCONDF MUST BE LARGER OR '
+ > //'EQUAL TO 0.')
+ ALLOCATE(KCONDF(NCONDF+3))
+ DO I=1,NCONDF+1
+ CALL REDGET(ITYP,NITMA,KCONDF(I),TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR KCONDF EXPECTED.')
+ ENDDO
+ CALL REDGET(ITYP,NITMA,FLOT,TEXT12,DFLOT)
+ IF(ITYP.NE.3)CALL XABORT('@THM: CHARACTER DATA EXPECTED (INV, '
+ > //'CELSIUS OR KELVIN) IN CONDF STATEMENT.')
+ IF(TEXT12.EQ.'INV') THEN
+ CALL REDGET(ITYP,NITMA,KCONDF(NCONDF+2),TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR INV EXPECTED.')
+ CALL REDGET(ITYP,NITMA,KCONDF(NCONDF+3),TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR REF EXPECTED.')
+ CALL REDGET(ITYP,NITMA,FLOT,TEXT12,DFLOT)
+ ELSE
+ KCONDF(NCONDF+2)=0.0 ! Coefficient for the inverse term
+ KCONDF(NCONDF+3)=-273.15 ! Reference for the inverse term
+ ENDIF
+ IF((TEXT12.NE.'CELSIUS').AND.(TEXT12.NE.'KELVIN')) THEN
+ CALL XABORT('@THM: UNIT KEYWORD EXPECTED (CELSIUS OR '
+ > //'KELVIN) IN CONDF STATEMENT.')
+ ENDIF
+ UCONDF=TEXT12
+ ELSE IF(TEXT.EQ.'CONDC') THEN
+ IF(ICONDC.EQ.1)DEALLOCATE(KCONDC)
+ ICONDC=1
+* Clad conductivity expressed as a polynomial of clad temperature
+ CALL REDGET(ITYP,NCONDC,FLOT,TEXT,DFLOT)
+ IF(ITYP.NE.1)CALL XABORT('@THM: INTEGER FOR CONDC EXPECTED.')
+ IF(NCONDC.LT.0)CALL XABORT('@THM: NCONDC MUST BE LARGER OR '
+ > //'EQUAL TO 0.')
+ ALLOCATE(KCONDC(NCONDC+1))
+ DO I=1,NCONDC+1
+ CALL REDGET(ITYP,NITMA,KCONDC(I),TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR KCONDC EXPECTED.')
+ ENDDO
+ CALL REDGET(ITYP,NITMA,FLOT,UCONDC,DFLOT)
+ IF((ITYP.NE.3).OR.((UCONDC.NE.'CELSIUS').AND.
+ > (UCONDC.NE.'KELVIN'))) THEN
+ CALL XABORT('@THM: UNIT KEYWORD EXPECTED (CELSIUS OR '
+ > //'KELVIN) IN CONDC STATEMENT.')
+ ENDIF
+ ELSE IF(TEXT.EQ.'HGAP') THEN
+ IHGAP=1
+* Fixed, user-chosen value of the HGAP (heat exchange coefficient
+* of the gap) (W/m^2/K)
+ CALL REDGET(ITYP,NITMA,KHGAP,TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR HGAP EXPECTED.')
+ ELSE IF(TEXT.EQ.'HCONV') THEN
+ IHCONV=1
+* Fixed, user-chosen value of the HCONV (heat transfer coefficient
+* between clad and fluid) (W/m^2/K)
+ CALL REDGET(ITYP,NITMA,KHCONV,TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR HCONV EXPECTED.')
+ ELSE IF(TEXT.EQ.'TEFF') THEN
+* Surface temperature's weighting factor in effective fuel
+* temperature
+ CALL REDGET(ITYP,NITMA,WTEFF,TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR TEFF EXPECTED.')
+ ELSE IF(TEXT.EQ.'FORCEAVE') THEN
+* Force the use of the average value approximation for fuel
+* conductivity
+ IFRCDI=1
+ ELSE IF(TEXT.EQ.'MONO') THEN
+* one-phase flow model
+ ISUBM=0
+ ELSE IF(TEXT.EQ.'BOWR') THEN
+* Bowring's correlation
+ ISUBM=1
+ ELSE IF(TEXT.EQ.'SAHA') THEN
+* Saha-Zuber correlation
+ ISUBM=2
+ ELSE IF(TEXT.EQ.'RADIUS') THEN
+* Fuel pellet radius (m)
+ CALL REDGET(ITYP,NITMA,RC,TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR RC EXPECTED.')
+* Internal clad rod radius (m)
+ CALL REDGET(ITYP,NITMA,RIG,TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR RIG EXPECTED.')
+* External clad rod radius (m)
+ CALL REDGET(ITYP,NITMA,RGG,TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR RGG EXPECTED.')
+* Guide tube radius (m)
+ CALL REDGET(ITYP,NITMA,RTG,TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR RTG EXPECTED.')
+ ELSE IF(TEXT.EQ.'ASSMB') THEN
+* Number of active fuel rods
+ CALL THMINP('NB-FUEL',NCH,FNFUCST)
+* Number of guide tubes
+ CALL THMINP('NB-TUBE',NCH,FNTGCST)
+ ELSE IF(TEXT.EQ.'CLUSTER') THEN
+* Hexagonal pitch (m)
+ CALL REDGET(ITYP,NITMA,PITCH,TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR PITCH EXPECTED.')
+* Number of active fuel pins in cluster
+ CALL THMINP('NB-FUEL',NCH,FNFUCST)
+ ELSE IF(TEXT.EQ.'CONV') THEN
+* Number of conduction iterations
+ CALL REDGET(ITYP,MAXIT1,FLOT,TEXT,DFLOT)
+ IF(ITYP.NE.1)CALL XABORT('@THM: INTEGER FOR MAXIT1 EXPECTED.')
+* Number of center-pellet iterations
+ CALL REDGET(ITYP,MAXIT2,FLOT,TEXT,DFLOT)
+ IF(ITYP.NE.1)CALL XABORT('@THM: INTEGER FOR MAXIT2 EXPECTED.')
+* Number of flow iterations
+ CALL REDGET(ITYP,MAXIT3,FLOT,TEXT,DFLOT)
+ IF(ITYP.NE.1)CALL XABORT('@THM: INTEGER FOR MAXIT3 EXPECTED.')
+* Temperature maximum error (K)
+ CALL REDGET(ITYP,NITMA,ERMAXT,TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR ERMAXT EXPECTED.')
+* maximum relative error for the calculation of the properties
+* in the coolant (pressure, enthalpy, density, velocity,...)
+ CALL REDGET(ITYP,NITMA,ERMAXC,TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR ERMAXC EXPECTED.')
+ ELSE IF(TEXT.EQ.'RODMESH') THEN
+* Number of discretisation points in fuel
+ CALL REDGET(ITYP,NFD,FLOT,TEXT,DFLOT)
+ IF(ITYP.NE.1)CALL XABORT('@THM: INTEGER FOR NFD EXPECTED.')
+* Number of discretisation points in fuel rod (fuel+cladding)
+ CALL REDGET(ITYP,NDTOT,FLOT,TEXT,DFLOT)
+ IF(ITYP.NE.1)CALL XABORT('@THM: INTEGER FOR NDTOT EXPECTED.')
+ ELSE IF(TEXT.EQ.'RELAX') THEN
+* Relaxation parameter
+ CALL REDGET(ITYP,NITMA,RELAX,TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR RELAX EXPECTED.')
+ ELSE IF(TEXT.EQ.'RAD-PROF') THEN
+* Set radial power profile
+ NPRAD=0
+ LPRAD=.TRUE.
+ 30 CALL REDGET(ITYP,NITMA,FLOT,TEXT,DFLOT)
+ IF(ITYP.EQ.3)GO TO 20
+ NPRAD=NPRAD+1
+ RPRAD(NPRAD)=FLOT
+ IF(NPRAD.GT.MAXRAD) CALL XABORT('@THM: MAXRAD OVERFLOW.')
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR RAD-PROF-X EXPECTED.')
+ IF(RPRAD(NPRAD).LT.0.0)CALL XABORT('@THM: R TOO SMALL.')
+ IF(RPRAD(NPRAD).GT.RC)CALL XABORT('@THM: R TOO LARGE.')
+ CALL REDGET(ITYP,NITMA,FPRAD(NPRAD),TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR RAD-PROF-F EXPECTED.')
+ GO TO 30
+ ELSE IF(TEXT.EQ.'POWER-LAW') THEN
+* The total power in W generated in the fuel is defined as
+* T-POWER*TIME-LAW(t).
+ CALL REDGET(ITYP,NITMA,TPOW,TEXT12,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR T-POWER EXPECTED.')
+ CALL REDGET(ITYP,NPOWER,FLOT,TEXT12,DFLOT)
+ IF(ITYP.NE.1)CALL XABORT('@THM: INTEGER VALUE EXPECTED.')
+ ALLOCATE(TIMESR(NPOWER),TPOWER(NPOWER))
+ DO I=1,NPOWER
+ CALL REDGET(ITYP,NITMA,TIMESR(I),TEXT12,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR TIME EXPECTED.')
+ CALL REDGET(ITYP,NITMA,TPOWER(I),TEXT12,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR POWER EXPECTED.')
+ ENDDO
+ CALL LCMPUT(IPTHM,'TIME-SR1',NPOWER,2,TIMESR)
+ CALL LCMPUT(IPTHM,'POWER-SR1',NPOWER,2,TPOWER)
+ DEALLOCATE(TPOWER,TIMESR)
+ ELSE IF(TEXT.EQ.'F-RUG') THEN
+* Rugosity of the fuel rod
+ CALL REDGET(ITYP,NITMA,EPSR,TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR F-RUG EXPECTED.')
+ ELSE IF(TEXT.EQ.'THETA') THEN
+* Angle of the fuel channel
+ CALL REDGET(ITYP,NITMA,THETA,TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR THETA EXPECTED.')
+ ELSE IF(TEXT.EQ.'PDROP') THEN
+* Pressure drop identification
+ CALL REDGET(ITYP,IPRES,FLOT,TEXT,DFLOT)
+ IF(ITYP.NE.1)CALL XABORT('@THM: INTEGER FOR IPRES EXPECTED.')
+ ELSE IF(TEXT.EQ.'DFM') THEN
+* Drift Flux Model identification
+ CALL REDGET(ITYP,IDFM,FLOT,TEXT,DFLOT)
+ IF(ITYP.NE.1)CALL XABORT('@THM: INTEGER FOR IDFM EXPECTED.')
+ ELSE IF(TEXT.EQ.'SET-PARAM') THEN
+* Reset a global parameter
+ CALL REDGET(ITYP,NITMA,FLOT,TEXT,DFLOT)
+ IF(ITYP.NE.3)CALL XABORT('@THM: CHARACTER NAME EXPECTED.')
+ CALL REDGET(ITYP,NITMA,VALUE,TEXT,DFLOT)
+ IF(ITYP.NE.2)CALL XABORT('@THM: REAL FOR VALUE EXPECTED.')
+ JPMAP=LCMGID(IPMAP,'PARAM')
+ DO 40 IPAR=1,NPARM
+ KPMAP=LCMGIL(JPMAP,IPAR)
+ CALL LCMGTC(KPMAP,'P-NAME',12,PNAME)
+ CALL LCMGET(KPMAP,'P-TYPE',ITYPE)
+ IF(ITYPE.EQ.1) THEN
+ IF(PNAME.EQ.TEXT) THEN
+ CALL LCMPUT(KPMAP,'P-VALUE',1,2,VALUE)
+ IF(IMPX.GT.0) WRITE(6,500) PNAME,VALUE
+ GO TO 10
+ ELSE
+ GO TO 40
+ ENDIF
+ ELSE IF(ITYPE.EQ.2) THEN
+ CALL XABORT('@THM: CANNOT RESET LOCAL PARAMETER: '//TEXT)
+ ENDIF
+ 40 CONTINUE
+ CALL XABORT('@THM: GLOBAL PARAMETER NAME NOT FOUND: '//TEXT)
+ ELSE IF(TEXT.EQ.';') THEN
+ GO TO 60
+ ELSE IF(TEXT.EQ.'PICK') THEN
+ IPICK=1
+ GO TO 60
+ ELSE
+ CALL XABORT('@THM: INVALID KEYWORD: '//TEXT//'.')
+ ENDIF
+ GO TO 10
+*----
+* TEST DATA INPUT
+*----
+ 60 IF(TINLET.LE.273.15) CALL XABORT('@THM: INLET TEMPERATURE MUST BE'
+ > //' HIGHER THAN 273.15K.')
+ IF(SPEED.EQ.0.0) CALL XABORT('@THM: ZERO COOLANT SPEED.')
+ IF(POULET.EQ.0.0) CALL XABORT('@THM: ZERO OUTLET PRESSURE.')
+ IF(RC.EQ.0.0) CALL XABORT('@THM: ZERO FUEL PELLET RADIUS.')
+ IF(RIG.EQ.0.0) CALL XABORT('@THM: ZERO INTERNAL CLAD ROD RADIUS.')
+ IF(RGG.EQ.0.0) CALL XABORT('@THM: ZERO EXTERNAL CLAD ROD RADIUS.')
+ IF(NDTOT.GT.NMAXO) CALL XABORT('@THM: NFD OVERFLOW, TOO MANY FUE'
+ > //'L DOMAINS')
+ IF(NDTOT.LT.8) CALL XABORT('@THM: NDTOT MUST AT LEAST BE EQUAL T'
+ > //'O 8')
+ IF(NFD.LT.4) CALL XABORT('@THM: NFD MUST AT LEAST BE EQUAL TO 4')
+ IF(NFD.GE.NDTOT) CALL XABORT('@THM: NFD MUST BE LOWER THAN NDTO'
+ > //'T.')
+ IF((RELAX.LE.0.0).OR.(RELAX.GT.1.0)) CALL XABORT('@THM: RELAX '
+ > //'PARAMETER EXPECTED BETWEEN 0<RELAX<=1.')
+ IF((WTEFF.LT.0.0).OR.(WTEFF.GT.1.0)) CALL XABORT('@THM: WTEFF '
+ > //'PARAMETER EXPECTED BETWEEN 0<=WTEFF<=1.')
+ IF(ITIME.EQ.1) RELAX=1.0
+ IF((RC.NE.RIG).AND.(IFUEL.EQ.1)) CALL XABORT('@THM: WITH MOLTEN'
+ > //' SALT FUEL INNER CLAD RADIUS MUST BE EQUAL TO FUEL RADIUS')
+*----
+* PRINT CHANNEL-DEPENDENT DATA
+*----
+ IF(IMPX.GT.1) THEN
+ WRITE(6,'(/28H THM: CHANNEL-DEPENDENT DATA)')
+ I1=1
+ DO I=1,(NCH-1)/8+1
+ I2=I1+7
+ IF(I2.GT.NCH) I2=NCH
+ WRITE(6,'(//8H CHANNEL,8(I8,6X,1H|))') (J,J=I1,I2)
+ WRITE(6,'(8H NB-FUEL,8(F10.2,4X,1H|))') (FNFUCST(J),J=I1,I2)
+ WRITE(6,'(8H NB-TUBE,8(F10.2,4X,1H|))') (FNTGCST(J),J=I1,I2)
+ WRITE(6,'(8H PUFR ,8(1P,E13.4,2H |))') (FRACPU(J),J=I1,I2)
+ I1=I1+8
+ ENDDO
+ ENDIF
+*----
+* SET POWER DISTRIBUTION
+*----
+ ALLOCATE(FRO(NFD-1))
+ IF(NPRAD.EQ.0) THEN
+ FRO(:NFD-1)=1.0
+ ELSE
+ IF(.NOT.LPRAD) THEN
+ CALL LCMGET(IPTHM,'RAD-PROF_R',RPRAD)
+ CALL LCMGET(IPTHM,'RAD-PROF_F',FPRAD)
+ ELSE
+ CALL LCMPUT(IPTHM,'RAD-PROF_R',NPRAD,2,RPRAD)
+ CALL LCMPUT(IPTHM,'RAD-PROF_F',NPRAD,2,FPRAD)
+ ENDIF
+ DAR1=0.0
+ DELT=0.5*RC**2/REAL(NFD-1)
+ DO IM=1,NFD-1
+ DAR2=DAR1+DELT
+ RADM=SQRT(DAR1+DAR2)
+ CALL ALTERP(.FALSE.,NPRAD,RPRAD(1),RADM,.FALSE.,TERP(1))
+ DSUM=0.0D0
+ DO J=1,NPRAD
+ DSUM=DSUM+TERP(J)*FPRAD(J)
+ ENDDO
+ FRO(IM)=REAL(DSUM)
+ DAR1=DAR2
+ ENDDO
+ ENDIF
+ IF(IMPX.GT.1) WRITE(6,480) (FRO(IM),IM=1,NFD-1)
+*----
+* RECOVER GEOMAP STATE-VECTOR
+* ISTATE(1): 7 = XYZ, 9 = HEXZ
+* In 3d hexagonal, NY=0, but THM: expects a 3D geometry, so we set
+* NY=1 and continue.
+*----
+ JPMAP=LCMGID(IPMAP,'GEOMAP')
+ CALL LCMGET(JPMAP,'STATE-VECTOR',ISTATE)
+ NX=ISTATE(3)
+ NY=ISTATE(4)
+ NZ=ISTATE(5)
+ NEL=ISTATE(6)
+ IF((ISTATE(1).EQ.9).AND.(NY.EQ.0)) NY=1
+ IF(NX.GT.IMAXO) CALL XABORT('@THM: NX OVERFLOW.')
+ IF(NY.GT.JMAXO) CALL XABORT('@THM: NY OVERFLOW.')
+ IF(NZ.GT.KMAXO) CALL XABORT('@THM: NZ OVERFLOW.')
+ ALLOCATE(PCH(NZ),ACOOL(NZ),HD(NZ))
+ ALLOCATE(FNFU(NCH,NZ),FNTG(NCH,NZ),RAPCOOL(NZ),
+ > RAPFUEL(NZ),PM(NZ),FFUEL(NZ),FCOOL(NZ))
+*----
+* RECOVER REACTOR MESH IN METER
+* The arrays HX, HY, and HZ contain the mesh size in X-, Y-, and
+* Z-direction and are used to determine the volume of a mesh, i.e.
+* V(I,J,K)=HX(I)*HY(J)*HZ(K)
+* For 3D hexagonal, set HX and HY to the square root of the SA surface
+* SASS
+*----
+ ALLOCATE(XX(NX+1),YY(NY+1),ZZ(NZ+1))
+ IF(ISTATE(1).EQ.7) THEN
+ CALL LCMGET(JPMAP,'MESHX',XX)
+ CALL LCMGET(JPMAP,'MESHY',YY)
+ ENDIF
+ CALL LCMGET(JPMAP,'MESHZ',ZZ)
+ IF(ISTATE(1).EQ.9) THEN
+ CALL LCMGET(JPMAP,'SIDE',SIDE)
+ SASS=1.5*SQRT(3.0)*SIDE*SIDE/1.0E4
+ DO 70 I=1,NX
+ HX(I) = SQRT(SASS)
+ 70 CONTINUE
+ DO 80 I=1,NY
+ HY(I) = SQRT(SASS)
+ 80 CONTINUE
+ ELSE
+ DO 90 I=1,NX
+ HX(I)=(XX(I+1)-XX(I))/100.0
+ 90 CONTINUE
+ DO 100 I=1,NY
+ HY(I)=(YY(I+1)-YY(I))/100.0
+ 100 CONTINUE
+ ENDIF
+ DO 110 I=1,NZ
+ HZ(I)=(ZZ(I+1)-ZZ(I))/100.0
+ 110 CONTINUE
+ DO 120 I=1,NZ+1
+ ZZ(I)=ZZ(I)/100.0
+ 120 CONTINUE
+ CALL LCMPUT(IPTHM,'MESHZ',NZ+1,2,ZZ)
+ DEALLOCATE(ZZ,YY,XX)
+*----
+* RECOVER LOCAL PARAMETER INFORMATION FROM L_MAP OBJECT
+*----
+ ALLOCATE(NUM(NEL),BURN(NCH*NB),PW(NCH*NB))
+ CALL LCMGET(IPMAP,'BMIX',NUM)
+ CALL LCMLEN(IPMAP,'BURN-INST',ILONG,ITYLCM)
+ IF(ILONG.EQ.NCH*NB) THEN
+ CALL LCMGET(IPMAP,'BURN-INST',BURN)
+ ELSE
+ CALL LCMLEN(IPMAP,'BURN-BEG',ILONG,ITYLCM)
+ IF(ILONG.NE.NCH*NB) CALL XABORT('@THM: MISSING BURNUP INFO ON '
+ > //'FUELMAP.')
+ ALLOCATE(BURN2(NCH*NB))
+ CALL LCMGET(IPMAP,'BURN-BEG',BURN)
+ CALL LCMGET(IPMAP,'BURN-END',BURN2)
+ DO I=1,NCH*NB
+ BURN(I)=(BURN(I)+BURN2(I))/2.0
+ ENDDO
+ DEALLOCATE(BURN2)
+ ENDIF
+ CALL LCMLEN(IPTHM,'POWER-SR1',NPOWER,ITYLCM)
+ IF(NPOWER.NE.0) THEN
+* USE POWER TIME LAW
+ IF(IMPX.GT.0) WRITE(6,*) 'THM: T-POWER = ',TPOW,' W'
+ IF(TPOW.EQ.0.0) CALL XABORT('@THM: T-POWER NOT DEFINED.')
+ IF(NCH.NE.1) CALL XABORT('@THM: NCH=1 EXPECTED.')
+ ALLOCATE(TIMESR(NPOWER),TPOWER(NPOWER),DTERP(NPOWER))
+ CALL LCMGET(IPTHM,'TIME-SR1',TIMESR)
+ CALL LCMGET(IPTHM,'POWER-SR1',TPOWER)
+ IF(ITIME.EQ.0) THEN
+ CALL ALTERP(.FALSE.,NPOWER,TIMESR(1),RTIME,.FALSE.,DTERP(1))
+ ELSE
+ IF(DTIME.EQ.0.0) CALL XABORT('@THM: DTIME NOT DEFINED.')
+ CALL ALTERI(.FALSE.,NPOWER,TIMESR(1),RTIME,RTIME+DTIME,
+ > DTERP(1))
+ DO J=1,NPOWER
+ DTERP(J)=DTERP(J)/DTIME
+ ENDDO
+ ENDIF
+ DPOW=0.0D0
+ DO J=1,NPOWER
+ DPOW=DPOW+DTERP(J)*TPOWER(J)
+ ENDDO
+ DPOW=DPOW*TPOW
+ DEALLOCATE(DTERP,TPOWER,TIMESR)
+ CALL LCMLEN(IPMAP,'AXIAL-FPW',ILONG,ITYLCM)
+ IF(ILONG.NE.NB) CALL XABORT('THM: NO AXIAL-FPW ON THE FUELMAP')
+ ALLOCATE(PFORM(NB))
+ CALL LCMGET(IPMAP,'AXIAL-FPW',PFORM)
+ DO I=1,NB
+ PW(I)=DPOW*PFORM(I)*1.0E-3
+ ENDDO
+ DEALLOCATE(PFORM)
+ ELSE
+* RECOVER POWER FROM FUELMAP
+ CALL LCMGET(IPMAP,'BUND-PW',PW)
+ ENDIF
+ IF(IMPX.GT.2) THEN
+ PTOT=0.0
+ DO I=1,NCH*NB
+ PTOT=PTOT+PW(I)
+ ENDDO
+ ENDIF
+*----
+* REBUILD LOCAL PARAMETER INFORMATION FOR THM
+*----
+ ALLOCATE(IREFSC(NCH))
+ CALL LCMLEN(IPMAP,'REF-SCHEME',ILONG,ITYLCM)
+ IF(ILONG.EQ.NCH) THEN
+ CALL LCMGET(IPMAP,'REF-SCHEME',IREFSC)
+ ELSE
+ IREFSC(:NCH)=1
+ ENDIF
+ ALLOCATE(XBURN(NZ,NX,NY),POW(NZ,NX,NY))
+ XBURN(:NZ,:NX,:NY)=0.0
+ POW(:NZ,:NX,:NY)=0.0
+ ICH=0
+ DO 165 IY=1,NY
+ DO 160 IX=1,NX
+ IEL=(IY-1)*NX+IX
+ DO 130 IZ=1,NZ
+ IF(NUM((IZ-1)*NX*NY+IEL).NE.0) GO TO 140
+ 130 CONTINUE
+ GO TO 160
+ 140 ICH=ICH+1
+ IB=0
+ DO 150 IZ=1,NZ
+ IF(NUM((IZ-1)*NX*NY+IEL).EQ.0) GO TO 150
+ IB=IB+1
+ IMA=(IB-1)*NCH+ICH
+ IF(IREFSC(ICH).GT.0) THEN
+ XBURN(IZ,IX,IY)=BURN(IMA)
+ POW(IZ,IX,IY)=PW(IMA)*1.0E3
+ ELSE
+ XBURN(NZ-IZ+1,IX,IY)=BURN(IMA)
+ POW(NZ-IZ+1,IX,IY)=PW(IMA)*1.0E3
+ ENDIF
+ 150 CONTINUE
+ IF(IB.NE.NB) CALL XABORT('@THM: INVALID NUMBER OF BUNDLES.')
+ 160 CONTINUE
+ 165 CONTINUE
+ IF(ICH.NE.NCH) CALL XABORT('@THM: INVALID NUMBER OF CHANNELS.')
+ DEALLOCATE(PW,BURN)
+*----
+* RECOVER AVERAGE FUEL TEMPERATURE FIELD FROM THM OBJECT
+*----
+ ALLOCATE(TCOMB(NZ,NX,NY))
+ TCOMB(:NZ,:NX,:NY)=0.0
+ JPMAP=LCMGID(IPMAP,'PARAM')
+ DO 220 IPAR=1,NPARM
+ KPMAP=LCMGIL(JPMAP,IPAR)
+ CALL LCMGTC(KPMAP,'P-NAME',12,PNAME)
+ IF(PNAME.EQ.'T-FUEL') THEN
+ CALL LCMGET(KPMAP,'P-TYPE',ITYPE)
+ ALLOCATE(VAL(NCH,NB))
+ IF(ITYPE.EQ.1) THEN
+ IF(IMPX.GT.0) WRITE(6,510) 'GLOBAL',PNAME
+ CALL LCMGET(KPMAP,'P-VALUE',FLOT)
+ DO 175 ICH=1,NCH
+ DO 170 IB=1,NB
+ VAL(ICH,IB)=FLOT
+ 170 CONTINUE
+ 175 CONTINUE
+ ELSE IF(ITYPE.EQ.2) THEN
+ IF(IMPX.GT.0) WRITE(6,510) 'LOCAL',PNAME
+ CALL LCMGET(KPMAP,'P-VALUE',VAL)
+ ENDIF
+ ICH=0
+ DO 215 IY=1,NY
+ DO 210 IX=1,NX
+ IEL=(IY-1)*NX+IX
+ DO 180 IZ=1,NZ
+ IF(NUM((IZ-1)*NX*NY+IEL).NE.0) GO TO 190
+ 180 CONTINUE
+ GO TO 210
+ 190 ICH=ICH+1
+ IB=0
+ DO 200 IZ=1,NZ
+ IF(NUM((IZ-1)*NX*NY+IEL).EQ.0) GO TO 200
+ IB=IB+1
+ IF(IREFSC(ICH).GT.0) THEN
+ TCOMB(IZ,IX,IY)=VAL(ICH,IB)
+ ELSE
+ TCOMB(NZ-IZ+1,IX,IY)=VAL(ICH,IB)
+ ENDIF
+ 200 CONTINUE
+ 210 CONTINUE
+ 215 CONTINUE
+ DEALLOCATE(VAL)
+ ENDIF
+ 220 CONTINUE
+ DEALLOCATE(IREFSC)
+*----
+* TEST TO COMPUTE STEADY-STATE OR TRANSIENT CALCULATION
+*----
+ IF(ITIME.EQ.0) THEN
+ GO TO 230
+ ELSE IF(ITIME.EQ.1) THEN
+ GO TO 310
+ ELSE
+ CALL XABORT('@THM: UNEXPECTED VALUE FOR ITIME.')
+ ENDIF
+*----
+* CALL DRIVER FOR STEADY-STATE CALCULATION
+*----
+* memory allocation for the steady-state calculation
+ 230 ALLOCATE(DCOOL(NZ,NX,NY),TCOOL(NZ,NX,NY),TSURF(NZ,NX,NY),
+ > PCOOL(NZ,NX,NY),HCOOL(NZ,NX,NY),RAD((NDTOT-1),NZ,NX,NY))
+ DCOOL(:NZ,:NX,:NY)=0.0
+ TCOOL(:NZ,:NX,:NY)=0.0
+ TSURF(:NZ,:NX,:NY)=0.0
+ PCOOL(:NZ,:NX,:NY)=0.0
+ HCOOL(:NZ,:NX,:NY)=0.0
+ RAD(:NDTOT-1,:NZ,:NX,:NY)=0.0
+*----
+* COMPUTE FUEL RADII
+*----
+ ALLOCATE(RVAL((NDTOT-1),NZ))
+ IF(JENTRY(1).EQ.0) THEN
+ WRITE(6,*)'RC,RIG=',RC,RIG
+*CGT THERE IS GAP
+ IF(RC.NE.RIG) THEN
+ IGAP=0
+ ARF=0.5*RC**2 ! at fuel radius
+ ARCI=0.5*RIG**2 ! at internal clad radius
+ ARCE=0.5*RGG**2 ! at external clad radius
+ DARF=ARF/REAL(NFD-1)
+ DARC=(ARCE-ARCI)/REAL(NDTOT-NFD-2)
+ DO IEL=1,NZ
+ RVAL(1,IEL)=0.0
+ DO I=1,NFD-1
+ RVAL(I+1,IEL)=REAL(SQRT(2.0D0*REAL(I)*DARF))
+ ENDDO
+ DO I=NFD+1,NDTOT-1
+ RVAL(I,IEL)=REAL(SQRT(2.0D0*(ARCI+REAL(I-NFD-1)*DARC)))
+ ENDDO
+ ENDDO
+ ELSE
+*CGT NO GAP
+ IGAP=1
+ ARF=0.5*RC**2 ! at fuel radius
+ ARCE=0.5*RGG**2 ! at external clad radius
+ DARF=ARF/REAL(NFD-1)
+ DARC=(ARCE-ARF)/REAL(NDTOT-NFD-1)
+ DO IEL=1,NZ
+ RVAL(1,IEL)=0.0
+ DO I=1,NFD
+ RVAL(I+1,IEL)=REAL(SQRT(2.0D0*REAL(I)*DARF))
+ ENDDO
+ DO I=NFD+1,NDTOT-1
+ RVAL(I,IEL)=REAL(SQRT(2.0D0*(ARF+REAL(I-NFD)*DARC)))
+ ENDDO
+ ENDDO
+ ENDIF
+ CALL LCMPUT(IPTHM,'REF-RAD',(NDTOT-1)*NZ,2,RVAL)
+ JPTHM=LCMDID(IPTHM,'HISTORY-DATA')
+ KPTHM=LCMDID(JPTHM,'TIMESTEP0000')
+ LPTHM=LCMLID(KPTHM,'CHANNEL',NCH)
+ ELSE
+ JPTHM=LCMGID(IPTHM,'HISTORY-DATA')
+ KPTHM=LCMGID(JPTHM,'TIMESTEP0000')
+ LPTHM=LCMGID(KPTHM,'CHANNEL')
+ ENDIF
+*----
+* LOOP OVER REACTOR CHANNELS
+*----
+ ICH=0
+ SUMSEC=0.0
+ DO 265 IY=1,NY
+ DO 260 IX=1,NX
+ IEL=IX+(IY-1)*NX
+ DO 240 IZ=1,NZ
+ IF(NUM((IZ-1)*NX*NY+IEL).NE.0) GO TO 250
+ 240 CONTINUE
+ GO TO 260
+ 250 ICH=ICH+1
+*----
+* COMPUTE HYDRAULICS CONSTANTS
+* SASS: assembly cross section in m^2
+* RC: fuel pellet radius in m
+* RTG: guide tube radius in m
+* ACOOL: coolant cross section per assembly in m^2
+* RAPCOOL: assembly over coolant volumic ratio
+* RAPFUEL: assembly over fuel volumic ratio
+* FCOOL: power density fraction in coolant.
+* FFUEL: power density fraction in fuel.
+* PCH: heating perimeter in m
+* PM: perimeter in contact with flow in m
+* HD: hydraulic diameter of one assembly in m
+* SPEED: inlet flow velocity in m/s
+*----
+ IF(FNFUCST(ICH).EQ.0.0) GO TO 260
+ SASS=HX(IX)*HY(IY)
+ DO K=1, NZ
+ FNFU(ICH,K)=FNFUCST(ICH)
+ FNTG(ICH,K)=FNTGCST(ICH)
+ IF(PITCH.EQ.0.0) THEN
+* PWR ASSEMBLY
+ ACOOL(K)=SASS-FNFU(ICH,K)*PI*RGG*RGG-FNTG(ICH,K)*PI*RTG*RTG
+ RAPCOOL(K)=SASS/ACOOL(K)
+ PCH(K)=FNFU(ICH,K)*2.0*PI*RGG
+ PM=PCH(K)+FNTG(ICH,K)*2.0*PI*RTG
+ SUMSEC=SUMSEC+ACOOL(K)
+ ELSE
+* CANDU CLUSTER
+ ATOTHEX=3.0*PITCH**2.0*(3.0)**0.5/2.0
+ ATIGEHEX=3.0*PI*RGG*RGG
+ ACOOL(K)=ATOTHEX-ATIGEHEX
+ PM(K)=6.0*PI*RGG
+ PCH(K)=PM(K)
+ RAPCOOL(K)=3.0*SASS/(FNFU(ICH,K)*ACOOL(K))
+ SUMSEC=SUMSEC+FNFU(ICH,K)*ACOOL(K)/3.0
+ ENDIF
+ RAPFUEL(K)=SASS/(FNFU(ICH,K)*PI*RC*RC)
+ FCOOL(K)=(1.0-FPUISS)*RAPCOOL(K)
+ FFUEL(K)=FPUISS*RAPFUEL(K)
+ HD(K)=4.0*ACOOL(K)/PM(K)
+ IF(HD(K).LE.0.) CALL XABORT('THM: NEGATIVE HYDRAULIC
+ >DIAMETER(1).')
+ ENDDO
+*----
+* RECOVER STEADY-STATE RADII
+*----
+ IF(JENTRY(1).EQ.0) THEN
+ RAD(:,:,IX,IY)=RVAL(:,:)
+ ELSE IF(JENTRY(1).EQ.1) THEN
+ MPTHM=LCMGIL(LPTHM,ICH)
+ CALL LCMGET(MPTHM,'RADII',RAD(1,1,IX,IY))
+ ENDIF
+*----
+* EXECUTION OF THE STEADY-STATE DRIVER PROGRAM
+*----
+ MPTHM=LCMDIL(LPTHM,ICH)
+ CALL THMDRV(MPTHM,IMPX,IX,IY,NZ,XBURN(1,IX,IY),SASS,HZ,CFLUX,
+ > POROS,FNFU(ICH,:),NFD,NDTOT,IFLUID,SNAME,SCOMP,IGAP,IFUEL,FNAME,
+ > FCOMP,FCOOL,FFUEL,ACOOL,
+ > HD,PCH,RAD(1,1,IX,IY),MAXIT1,MAXIT2,ERMAXT,SPEED,TINLET,POULET,
+ > FRACPU(ICH),ICONDF,NCONDF,KCONDF,UCONDF,ICONDC,NCONDC,KCONDC,
+ > UCONDC,IHGAP,KHGAP,IHCONV,KHCONV,WTEFF,IFRCDI,ISUBM,FRO,
+ > POW(1,IX,IY),IPRES,IDFM,TCOMB(1,IX,IY),DCOOL(1,IX,IY),
+ > TCOOL(1,IX,IY),TSURF(1,IX,IY),HCOOL(1,IX,IY),PCOOL(1,IX,IY))
+ 260 CONTINUE
+ 265 CONTINUE
+ IF(IMPX.GT.1) WRITE(6,610) SUMSEC,CWSECT
+ DEALLOCATE(RVAL)
+ CALL LCMPUT(KPTHM,'TIME',1,2,RTIME)
+ IF(IMPX.GT.1) WRITE(6,470) 'TIMESTEP0000',RTIME
+*----
+* PRINT AVERAGED THERMALHYDRAULICS PROPERTIES OVER THE CORE MAP
+*----
+ IF(IMPX.GT.1) THEN
+ CALL THMAVG(IPMAP,IMPX,NX,NY,NZ,NCH,TCOMB,TSURF,DCOOL,TCOOL,
+ > PCOOL,HCOOL,POW,NSIMS)
+ ENDIF
+ DEALLOCATE(RAD,HCOOL)
+ GO TO 400
+*----
+* CALL DRIVER FOR TRANSIENT CALCULATION
+*----
+* memory allocation for the transient calculation
+ 310 ALLOCATE(TSURF(NZ,NX,NY),TCOOL(NZ,NX,NY),DCOOL(NZ,NX,NY),
+ > PCOOL(NZ,NX,NY))
+ TSURF(:NZ,:NX,:NY)=0.0
+ TCOOL(:NZ,:NX,:NY)=0.0
+ DCOOL(:NZ,:NX,:NY)=0.0
+ PCOOL(:NZ,:NX,:NY)=0.0
+*----
+* RECOVER TIME INDEX AT INITIAL CONDITIONS
+*----
+ JPTHM=LCMDID(IPTHM,'HISTORY-DATA')
+ KPTHMI=LCMGID(JPTHM,'TIMESTEP0000')
+ CALL LCMGET(KPTHMI,'TIME',TIMEPR)
+ IF(ABS(RTIME-TIMEPR).LE.1.0E-3*DTIME) THEN
+ TIMEIT=1
+ ELSE
+ DO I=1,TIMEIT
+ WRITE(TXTDIR,'(8HTIMESTEP,I4.4)') I
+ KPTHMI=LCMGID(JPTHM,TXTDIR)
+ CALL LCMGET(KPTHMI,'TIME',TIMEPR)
+ IF(ABS(RTIME-TIMEPR).LE.1.0E-3*DTIME) THEN
+ TIMEIT=I+1
+ GO TO 315
+ ENDIF
+ ENDDO
+ WRITE(HSMG,'(45H@THM: UNABLE TO FIND INITIAL CONDITIONS AT T=,
+ > 1P,E14.4,3H S.)') RTIME
+ CALL XABORT(HSMG)
+ ENDIF
+ 315 LPTHMI=LCMGID(KPTHMI,'CHANNEL')
+ WRITE(TXTDIR,'(8HTIMESTEP,I4.4)') TIMEIT
+ KPTHM=LCMDID(JPTHM,TXTDIR)
+ LPTHM=LCMLID(KPTHM,'CHANNEL',NCH)
+ IF(IMPX.GT.1) WRITE(6,530) TIMEIT,RTIME,RTIME+DTIME
+*----
+* LOOP OVER REACTOR CHANNELS
+*----
+ ICH=0
+ DO 355 IY=1,NY
+ DO 350 IX=1,NX
+ IEL=IX+(IY-1)*NX
+ DO 320 IZ=1,NZ
+ IF(NUM((IZ-1)*NX*NY+IEL).NE.0) GO TO 330
+ 320 CONTINUE
+ GO TO 350
+ 330 ICH=ICH+1
+*----
+* COMPUTE HYDRAULICS CONSTANTS
+*----
+ IF(FNFUCST(ICH).EQ.0.0) GO TO 350
+ SASS=HX(IX)*HY(IY)
+ DO K=1, NZ
+ FNFU(ICH,K)=FNFUCST(ICH)
+ FNTG(ICH,K)=FNTGCST(ICH)
+ IF(PITCH.EQ.0.0) THEN
+* PWR ASSEMBLY
+ ACOOL(K)=SASS-FNFU(ICH,K)*PI*RGG*RGG-FNTG(ICH,K)*PI*RTG*RTG
+ RAPCOOL(K)=SASS/ACOOL(K)
+ PCH(K)=FNFU(ICH,K)*2.0*PI*RGG
+ PM=PCH(K)+FNTG(ICH,K)*2.0*PI*RTG
+ ELSE
+* CANDU CLUSTER
+ ATOTHEX=3.0*PITCH**2.0*(3.0)**0.5/2.0
+ ATIGEHEX=3.0*PI*RGG*RGG
+ ACOOL(K)=ATOTHEX-ATIGEHEX
+ PM(K)=6.0*PI*RGG
+ PCH(K)=PM(K)
+ RAPCOOL(K)=3.0*SASS/(FNFU(ICH,K)*ACOOL(K))
+ ENDIF
+ RAPFUEL(K)=SASS/(FNFU(ICH,K)*PI*RC*RC)
+ FCOOL(K)=(1.0-FPUISS)*RAPCOOL(K)
+ FFUEL(K)=FPUISS*RAPFUEL(K)
+ HD(K)=4.0*ACOOL(K)/PM(K)
+ IF(HD(K).LE.0.) CALL XABORT('THM: NEGATIVE HYDRAULIC
+ >DIAMETER(2).')
+ ENDDO
+*----
+* EXECUTION OF THE TRANSIENT DRIVER PROGRAM
+*----
+ MPTHMI=LCMGIL(LPTHMI,ICH)
+ MPTHM=LCMDIL(LPTHM,ICH)
+ CALL THMTRS(MPTHMI,MPTHM,IMPX,IX,IY,NZ,XBURN(1,IX,IY),SASS,HZ,
+ > DTIME,CFLUX,POROS,FNFU(ICH,:),NFD,NDTOT,IFLUID,SNAME,SCOMP,
+ > IGAP,IFUEL,FNAME,FCOMP,
+ > FCOOL,FFUEL,ACOOL,HD,PCH,MAXIT3,MAXIT1,MAXIT2,ERMAXT,ERMAXC,
+ > SPEED,TINLET,POULET,FRACPU(ICH),ICONDF,NCONDF,KCONDF,UCONDF,
+ > ICONDC,NCONDC,KCONDC,UCONDC,IHGAP,KHGAP,IHCONV,KHCONV,WTEFF,
+ > IFRCDI,ISUBM,FRO,POW(1,IX,IY),TCOMB(1,IX,IY),DCOOL(1,IX,IY),
+ > TCOOL(1,IX,IY),TSURF(1,IX,IY))
+ 350 CONTINUE
+ 355 CONTINUE
+ CALL LCMPUT(KPTHM,'TIME',1,2,RTIME+DTIME)
+ IF(IMPX.GT.1) WRITE(6,470) TXTDIR,RTIME+DTIME
+*----
+* RECOVER LOCAL PARAMETER INFORMATION COMPUTED BY THMDRV OR THMTRS
+*----
+ 400 ERRA1=0.0
+ ERRA2=0.0
+ ERRA3=0.0
+ ERRBB=0.0
+ ZMINA1=1.0E10
+ ZMINA2=1.0E10
+ ZMINA3=1.0E10
+ ZMINBB=1.0E10
+ ZMAXA1=0.0
+ ZMAXA2=0.0
+ ZMAXA3=0.0
+ ZMAXBB=0.0
+ RATIOX=0.0
+ ALLOCATE(IREFSC(NCH))
+ CALL LCMLEN(IPMAP,'REF-SCHEME',ILONG,ITYLCM)
+ IF(ILONG.EQ.NCH) THEN
+ CALL LCMGET(IPMAP,'REF-SCHEME',IREFSC)
+ ELSE
+ IREFSC(:NCH)=1
+ ENDIF
+ JPMAP=LCMGID(IPMAP,'PARAM')
+ DO 460 IPAR=1,NPARM
+ KPMAP=LCMGIL(JPMAP,IPAR)
+ CALL LCMGTC(KPMAP,'P-NAME',12,PNAME)
+ IF((PNAME.EQ.'T-FUEL').OR.(PNAME.EQ.'D-COOL').OR.
+ 1 (PNAME.EQ.'T-COOL').OR.(PNAME.EQ.'T-SURF').OR.
+ 2 (PNAME.EQ.'P-COOL')) THEN
+ CALL LCMGET(KPMAP,'P-TYPE',ITYPE)
+ ALLOCATE(VAL(NCH,NB))
+ RELAX0=1.0
+ IF(ITYPE.EQ.1) THEN
+ IF(IMPX.GT.0) WRITE(6,510) 'GLOBAL',PNAME
+ CALL LCMGET(KPMAP,'P-VALUE',FLOT)
+ DO 415 ICH=1,NCH
+ DO 410 IB=1,NB
+ VAL(ICH,IB)=FLOT
+ 410 CONTINUE
+ 415 CONTINUE
+ ELSE IF(ITYPE.EQ.2) THEN
+ RELAX0=RELAX
+ IF(IMPX.GT.0) WRITE(6,510) 'LOCAL',PNAME
+ CALL LCMGET(KPMAP,'P-VALUE',VAL)
+ ENDIF
+ ICH=0
+ DO 455 IY=1,NY
+ DO 450 IX=1,NX
+ IEL=(IY-1)*NX+IX
+ DO 420 IZ=1,NZ
+ IF(NUM((IZ-1)*NX*NY+IEL).NE.0) GO TO 430
+ 420 CONTINUE
+ GO TO 450
+ 430 ICH=ICH+1
+ IB=0
+ DO 440 IZ=1,NZ
+ IF(NUM((IZ-1)*NX*NY+IEL).EQ.0) GO TO 440
+ IB=IB+1
+ FLOT=0.0
+ IF(PNAME.EQ.'T-FUEL') THEN
+ IF(IREFSC(ICH).GT.0) THEN
+ FLOT=RELAX0*TCOMB(IZ,IX,IY)+(1.0-RELAX0)*VAL(ICH,IB)
+ ELSE
+ FLOT=RELAX0*TCOMB(NZ-IZ+1,IX,IY)+(1.0-RELAX0)*VAL(ICH,IB)
+ ENDIF
+ IF(ITIME.EQ.0) ERRA1=MAX(ERRA1,ABS(VAL(ICH,IB)-FLOT))
+ ZMINA1=MIN(ZMINA1,FLOT)
+ ZMAXA1=MAX(ZMAXA1,FLOT)
+ ELSE IF(PNAME.EQ.'D-COOL') THEN
+ IF(IREFSC(ICH).GT.0) THEN
+ FLOT=RELAX0*DCOOL(IZ,IX,IY)/ZKILO+(1.0-RELAX0)*VAL(ICH,IB)
+ ELSE
+ FLOT=RELAX0*DCOOL(NZ-IZ+1,IX,IY)/ZKILO+(1.0-RELAX0)
+ > *VAL(ICH,IB)
+ ENDIF
+ IF(ITIME.EQ.0) ERRA2=MAX(ERRA2,ABS(VAL(ICH,IB)-FLOT))
+ ZMINA2=MIN(ZMINA2,FLOT)
+ ZMAXA2=MAX(ZMAXA2,FLOT)
+ ELSE IF(PNAME.EQ.'T-COOL') THEN
+ IF(IREFSC(ICH).GT.0) THEN
+ FLOT=RELAX0*TCOOL(IZ,IX,IY)+(1.0-RELAX0)*VAL(ICH,IB)
+ ELSE
+ FLOT=RELAX0*TCOOL(NZ-IZ+1,IX,IY)+(1.0-RELAX0)*VAL(ICH,IB)
+ ENDIF
+ IF(ITIME.EQ.0) ERRA3=MAX(ERRA3,ABS(VAL(ICH,IB)-FLOT))
+ ZMINA3=MIN(ZMINA3,FLOT)
+ ZMAXA3=MAX(ZMAXA3,FLOT)
+ ELSE IF(PNAME.EQ.'T-SURF') THEN
+ IF(IREFSC(ICH).GT.0) THEN
+ FLOT=RELAX0*TSURF(IZ,IX,IY)+(1.0-RELAX0)*VAL(ICH,IB)
+ ELSE
+ FLOT=RELAX0*TSURF(NZ-IZ+1,IX,IY)+(1.0-RELAX0)*VAL(ICH,IB)
+ ENDIF
+ IF(ITIME.EQ.0) ERRBB=MAX(ERRBB,ABS(VAL(ICH,IB)-FLOT))
+ ZMINBB=MIN(ZMINBB,FLOT)
+ ZMAXBB=MAX(ZMAXBB,FLOT)
+ ELSE IF(PNAME.EQ.'P-COOL') THEN
+ IF(IREFSC(ICH).GT.0) THEN
+ FLOT=RELAX0*PCOOL(IZ,IX,IY)+(1.0-RELAX0)*VAL(ICH,IB)
+ ELSE
+ FLOT=RELAX0*PCOOL(NZ-IZ+1,IX,IY)+(1.0-RELAX0)*VAL(ICH,IB)
+ ENDIF
+ IF(ITIME.EQ.0) ERRA4=MAX(ERRBB,ABS(VAL(ICH,IB)-FLOT))
+ ZMINA4=MIN(ZMINBB,FLOT)
+ ZMAXA4=MAX(ZMAXBB,FLOT)
+ ELSE
+ CALL XABORT('@THM: INVALID PARAMETER TYPE: '// PNAME//'.')
+ ENDIF
+ VAL(ICH,IB)=FLOT
+ 440 CONTINUE
+ 450 CONTINUE
+ 455 CONTINUE
+ ITYPE=2
+ CALL LCMPUT(KPMAP,'P-TYPE',1,1,ITYPE)
+ CALL LCMPUT(KPMAP,'P-VALUE',NCH*NB,2,VAL)
+ CALL LCMLEN(IPMAP,'AXIAL-FPW',JLONG,ITYLCM)
+ DD1=0.0
+ DD2=0.0
+ IF(JLONG.NE.0) THEN
+ ALLOCATE(FPOWER(NB))
+ IF(JLONG.NE.NB) CALL XABORT('THM: UNABLE TO FIND RECORD AXIA'
+ 1 //'L-FPW IN THE FUELMAP.')
+ CALL LCMGET(IPMAP,'AXIAL-FPW',FPOWER)
+ DO ICH=1,NCH
+ DO IB=1,NB
+ DD1=DD1+VAL(ICH,IB)*FPOWER(IB)**2
+ DD2=DD2+FPOWER(IB)**2
+ ENDDO
+ ENDDO
+ DEALLOCATE(FPOWER)
+ ELSE
+ ALLOCATE(PW(NCH*NB))
+ CALL LCMGET(IPMAP,'BUND-PW',PW)
+ ITOT=0
+ DO IB=1,NB
+ DO ICH=1,NCH
+ ITOT=ITOT+1
+ DD1=DD1+VAL(ICH,IB)*PW(ITOT)**2
+ DD2=DD2+PW(ITOT)**2
+ ENDDO
+ ENDDO
+ DEALLOCATE(PW)
+ ENDIF
+ TMOY0=DD1/DD2
+ TEXT12='AVG-'//PNAME(:8)
+ CALL LCMLEN(IPTHM,TEXT12,KLONG,ITYLCM)
+ IF(((PNAME.EQ.'T-FUEL').OR.(PNAME.EQ.'T-COOL').OR.
+ 1 (PNAME.EQ.'P-COOL')).AND.(KLONG.GT.0)) THEN
+ CALL LCMGET(IPTHM,TEXT12,TMOY0I)
+ IF(PNAME.EQ.'T-FUEL') THEN
+ RATIO=ABS(TMOY0/DTEMPR-TMOY0I/DTEMPR)
+ IF(IMPX.GT.0) WRITE(6,490) TEXT12,TMOY0I,TMOY0,RATIO
+ RATIOX=MAX(RATIOX,RATIO)
+ ELSE IF(PNAME.EQ.'T-COOL') THEN
+ RATIO=ABS(TMOY0/DTEMPT-TMOY0I/DTEMPT)
+ IF(IMPX.GT.0) WRITE(6,490) TEXT12,TMOY0I,TMOY0,RATIO
+ RATIOX=MAX(RATIOX,RATIO)
+ ELSE IF(PNAME.EQ.'P-COOL') THEN
+ RATIO=ABS(TMOY0/DPRESS-TMOY0I/DPRESS)
+ IF(IMPX.GT.0) WRITE(6,490) TEXT12,TMOY0I,TMOY0,RATIO
+ RATIOX=MAX(RATIOX,RATIO)
+ ENDIF
+ ENDIF
+ CALL LCMPUT(IPTHM,TEXT12,1,2,TMOY0)
+ DEALLOCATE(VAL)
+ ENDIF
+ IF(PNAME.EQ.'T-FUEL') THEN
+ IF(ITIME.EQ.0) CALL LCMPUT(IPTHM,'ERROR-T-FUEL',1,2,ERRA1)
+ CALL LCMPUT(IPTHM,'MIN-T-FUEL',1,2,ZMINA1)
+ CALL LCMPUT(IPTHM,'MAX-T-FUEL',1,2,ZMAXA1)
+ IF(IMPX.GT.0) WRITE(6,520) 'FUEL TEMPERATURE',ERRA1,'K',
+ 1 ZMINA1,'K',ZMAXA1,'K'
+ ELSE IF(PNAME.EQ.'D-COOL') THEN
+ IF(ITIME.EQ.0) CALL LCMPUT(IPTHM,'ERROR-D-COOL',1,2,ERRA2)
+ CALL LCMPUT(IPTHM,'MIN-D-COOL',1,2,ZMINA2)
+ CALL LCMPUT(IPTHM,'MAX-D-COOL',1,2,ZMAXA2)
+ IF(IMPX.GT.0) WRITE(6,520) 'COOLANT DENSITY',ERRA2,'g/cc',
+ 1 ZMINA2,'g/cc',ZMAXA2,'g/cc'
+ ELSE IF(PNAME.EQ.'T-COOL') THEN
+ IF(ITIME.EQ.0) CALL LCMPUT(IPTHM,'ERROR-T-COOL',1,2,ERRA3)
+ CALL LCMPUT(IPTHM,'MIN-T-COOL',1,2,ZMINA3)
+ CALL LCMPUT(IPTHM,'MAX-T-COOL',1,2,ZMAXA3)
+ IF(IMPX.GT.0) WRITE(6,520) 'COOLANT TEMPERATURE',ERRA3,'K',
+ 1 ZMINA3,'K',ZMAXA3,'K'
+ ELSE IF(PNAME.EQ.'T-SURF') THEN
+ IF(ITIME.EQ.0) CALL LCMPUT(IPTHM,'ERROR-T-SURF',1,2,ERRBB)
+ IF(IMPX.GT.0) WRITE(6,520) 'FUEL SURFACE TEMPERATURE',ERRBB,
+ 1 'K',ZMINBB,'K',ZMAXBB,'K'
+ ELSE IF(PNAME.EQ.'P-COOL') THEN
+ IF(ITIME.EQ.0) CALL LCMPUT(IPTHM,'ERROR-P-COOL',1,2,ERRA4)
+ CALL LCMPUT(IPTHM,'MIN-P-COOL',1,2,ZMINA4)
+ CALL LCMPUT(IPTHM,'MAX-P-COOL',1,2,ZMAXA4)
+ IF(IMPX.GT.0) WRITE(6,520) 'COOLANT PRESSURE',ERRA4,'Pa',
+ 1 ZMINA4,'Pa',ZMAXA4,'Pa'
+ ENDIF
+ 460 CONTINUE
+ DEALLOCATE(IREFSC)
+*----
+* SAVE CONDUCTIVITY INFORMATION ON LCM OBJECT THM
+*----
+ IF(ICONDF.EQ.1) THEN
+ CALL LCMPUT(IPTHM,'KCONDF',NCONDF+3,2,KCONDF)
+ CALL LCMPTC(IPTHM,'UCONDF',12,UCONDF)
+ ENDIF
+ IF(ICONDC.EQ.1) THEN
+ CALL LCMPUT(IPTHM,'KCONDC',NCONDC+1,2,KCONDC)
+ CALL LCMPTC(IPTHM,'UCONDC',12,UCONDC)
+ ENDIF
+*----
+* RELEASE MEMORY
+*----
+ DEALLOCATE(ACOOL,PCH,HD)
+ DEALLOCATE(RAPCOOL,RAPFUEL,PM,FNFU,FNTG,FFUEL,FCOOL)
+ DEALLOCATE(PCOOL,DCOOL,TCOOL,TSURF,TCOMB)
+ DEALLOCATE(NUM)
+ DEALLOCATE(POW,XBURN)
+ DEALLOCATE(FRO)
+ IF(ICONDF.EQ.1)DEALLOCATE(KCONDF)
+ IF(ICONDC.EQ.1)DEALLOCATE(KCONDC)
+*----
+* STATE-VECTOR FOR THM
+*----
+ HSIGN='L_THM'
+ CALL LCMPTC(IPTHM,'SIGNATURE',12,HSIGN)
+ ISTATE(:NSTATE)=0
+ ISTATE(1)=NCH
+ ISTATE(2)=NZ
+ ISTATE(3)=MAXIT1
+ ISTATE(4)=MAXIT2
+ ISTATE(5)=MAXIT3
+ ISTATE(6)=NFD
+ ISTATE(7)=NDTOT
+ ISTATE(8)=ITIME
+ ISTATE(9)=TIMEIT
+ ISTATE(10)=IHGAP
+ ISTATE(11)=IHCONV
+ ISTATE(12)=ICONDF
+ ISTATE(13)=ICONDC
+ ISTATE(14)=IFRCDI
+ ISTATE(15)=ISUBM
+ IF(ICONDF.EQ.1) ISTATE(16)=NCONDF
+ IF(ICONDC.EQ.1) ISTATE(17)=NCONDC
+ ISTATE(18)=NPRAD
+ ISTATE(19)=NPOWER
+ ISTATE(20)=IFLUID
+ ISTATE(21)=IGAP
+ ISTATE(22)=IPRES
+ ISTATE(23)=IDFM
+ CALL LCMPUT(IPTHM,'STATE-VECTOR',NSTATE,1,ISTATE)
+ STATE(:NSTATE)=0.0
+ STATE(1)=DTIME
+ STATE(2)=FPUISS
+ STATE(3)=CFLUX
+ STATE(4)=SPEED
+ STATE(5)=POULET
+ STATE(6)=TINLET
+ STATE(7)=POROS
+ STATE(8)=RC
+ STATE(9)=RIG
+ STATE(10)=RGG
+ STATE(11)=RTG
+ STATE(12)=PITCH
+ STATE(13)=ERMAXT
+ STATE(14)=ERMAXC
+ STATE(15)=RELAX
+ STATE(16)=RTIME
+ IF(IHGAP.EQ.1) STATE(17)=KHGAP
+ IF(IHCONV.EQ.1) STATE(18)=KHCONV
+ STATE(19)=WTEFF
+ STATE(20)=TPOW
+ STATE(21)=RATIOX
+ STATE(22)=EPSR
+ STATE(23)=THETA
+ CALL LCMPUT(IPTHM,'REAL-PARAM',NSTATE,2,STATE)
+ IF(IMPX.GT.0) THEN
+ WRITE(6,540) ISTATE(:15),ISTATE(18:23)
+ IF(ISTATE(10).EQ.1) WRITE(6,550) (ISTATE(16))
+ IF(ISTATE(11).EQ.1) WRITE(6,560) (ISTATE(17))
+ WRITE(6,570) STATE(:16),STATE(19),STATE(21:23)
+ IF(ISTATE(10).EQ.1) WRITE(6,580) (STATE(17))
+ IF(ISTATE(11).EQ.1) WRITE(6,590) (STATE(18))
+ IF(ISTATE(19).GT.0) WRITE(6,600) (STATE(20))
+ ENDIF
+ IF(IMPX.GT.4) CALL LCMLIB(IPTHM)
+*----
+* SAVE CELL-DEPENDENT DATA
+*----
+ CALL LCMPUT(IPTHM,'NB-FUEL',NCH,2,FNFUCST)
+ CALL LCMPUT(IPTHM,'NB-TUBE',NCH,2,FNTGCST)
+ CALL LCMPUT(IPTHM,'FRACT-PU',NCH,2,FRACPU)
+ DEALLOCATE(FRACPU,FNTGCST,FNFUCST)
+*----
+* RECOVER THE VARIATION RATIO AND SAVE IT IN A CLE-2000 VARIABLE
+*----
+ IF(IPICK.EQ.1) THEN
+ CALL REDGET(ITYP,NITMA,FLOT,TEXT12,DFLOT)
+ IF(ITYP.NE.-2) CALL XABORT('THM: OUTPUT REAL EXPECTED.')
+ ITYP=2
+ CALL REDPUT(ITYP,NITMA,RATIOX,TEXT12,DFLOT)
+ CALL REDGET(ITYP,NITMA,FLOT,TEXT12,DFLOT)
+ IF((ITYP.NE.3).OR.(TEXT12.NE.';')) THEN
+ CALL XABORT('THM: ; CHARACTER EXPECTED.')
+ ENDIF
+ ENDIF
+ RETURN
+*
+ 470 FORMAT(/11H THM: SAVE ,A,9H AT TIME=,1P,E12.4,3H S.)
+ 480 FORMAT(/31H THM: RADIAL POWER FORM FACTORS/(1P,10E12.4))
+ 490 FORMAT(/18H THM: PARAMETER = ,A,1P,E12.4,3H ->,E12.4,7H RATIO=,
+ 1 E12.4)
+ 500 FORMAT(/27H THM: SET GLOBAL PARAMETER ,A,2H =,1P E12.4)
+ 510 FORMAT(/14H THM: RECOVER ,A,13H PARAMETER = ,A,1H.)
+ 520 FORMAT(/15H THM: ERROR ON ,A,2H =,F12.3,1X,A,13H MIN VALUE =,
+ 1 F12.3,1X,A,13H MAX VALUE =,F12.3,1X,A)
+ 530 FORMAT(/28H THM: PERFORM TRANSIENT STEP,I5,9H BETWEEN ,1P,E14.4,
+ 1 4H AND,E14.4,3H S.)
+ 540 FORMAT(/
+ 1 14H STATE VECTOR:/
+ 2 7H NZ ,I9,27H (NUMBER OF AXIAL MESHES)/
+ 3 7H NCH ,I9,43H (NUMBER OF CHANNELS IN THE RADIAL PLANE)/
+ 4 7H MAXIT1,I9,36H (NUMBER OF CONDUCTION ITERATIONS)/
+ 5 7H MAXIT2,I9,39H (NUMBER OF CENTER-PELLET ITERATIONS)/
+ 6 7H MAXIT3,I9,30H (NUMBER OF FLOW ITERATIONS)/
+ 7 7H NFD ,I9,32H (NUMBER OF FUEL RADIAL ZONES)/
+ 8 7H NDTOT ,I9,36H (NUMBER OF DISCRETISATION POINTS)/
+ 9 7H ITIME ,I9,21H (CALCULATION TYPE)/
+ 1 7H TIMEIT,I9,30H (TRANSIENT ITERATION INDEX)/
+ 2 7H IHGAP ,I9,34H (HGAP FLAG (0=DEFAULT/1=FIXED))/
+ 3 7H IHCONV,I9,42H (HCONV FLAG (0=DITTUS-BOELTER/1=FIXED))/
+ 4 7H ICONDF,I9,46H (FUEL CONDUCTIVITY FLAG (0=STORA-CHENEBAULT,
+ 5 54H (UOX), COMETHE (MOX)/1=USER-PROVIDED FUNCTION OF FUEL,
+ 6 14H TEMPERATURE))/
+ 7 7H ICONDC,I9,39H (CLAD CONDUCTIVITY FLAG (0=DEFAULT/1,
+ 8 47H=USER-PROVIDED POLYNOMIAL OF CLAD TEMPERATURE))/
+ 9 7H IFRCDI,I9,40H (FUEL CONDUCTIVITY APPROXIMATION FLAG,
+ 1 44H (0=DEFAULT/1=AVERAGE APPROXIMATION FORCED))/
+ 2 7H ISUBM ,I9,47H (BOILING MODEL FLAG (0=ONE-PHASE/1=BOWRING C,
+ 3 37HORRELATION/2=SAHA-ZUBER CORRELATION))/
+ 4 7H NPRAD ,I9,47H (RADIAL POWER FORM FACTOR (0=FLAT/NUMBER OF ,
+ 5 8HPOINTS))/
+ 6 7H NPOWER,I9,36H (NUMBER OF POINTS IN POWER-TABLE)/
+ 7 7H IFLUID,I9,32H (TYPE OF FLUID (0=H2O/1=D2O))/
+ 8 7H IGAP ,I9,39H (GAP IS CONSIDERED (0=GAP/1=NO GAP))/
+ 9 7H IPRES ,I9,46H (PRESSURE DROP (0=CONSTANT/1=NON CONSTANT))/
+ 1 7H IDFM ,I9,47H (DRIFT FLUX MODEL (0=HEM1/1=EPRI/2=MODEBSTIO,
+ 2 21HN/3=GERAMP/4=CHEXAL)))
+ 550 FORMAT(
+ 1 7H NCONDF,I9,43H (DEGREE OF FUEL CONDUCTIVITY POLYNOMIAL))
+ 560 FORMAT(
+ 1 7H NCONDC,I9,43H (DEGREE OF CLAD CONDUCTIVITY POLYNOMIAL))
+ 570 FORMAT(/
+ 1 12H REAL PARAM:,1P/
+ 2 7H DTIME ,E12.4,19H (TIME STEP IN S)/
+ 3 7H FPUISS,E12.4,25H (COOLANT POWER FACTOR)/
+ 4 7H CFLUX ,E12.4,32H (CRITICAL HEAT FLUX IN W/M^2)/
+ 5 7H SPEED ,E12.4,28H (COOLANT VELOCITY IN M/S)/
+ 6 7H POULET,E12.4,34H (OUTLET COOLANT PRESSURE IN PA)/
+ 7 7H TINLET,E12.4,35H (INLET COOLANT TEMPERATURE IN K)/
+ 8 7H POROS ,E12.4,19H (OXYDE POROSITY)/
+ 9 7H RC ,E12.4,28H (FUEL PELLET RADIUS IN M)/
+ 1 7H RIG ,E12.4,34H (INTERNAL CLAD ROD RADIUS IN M)/
+ 2 7H RGG ,E12.4,34H (EXTERNAL CLAD ROD RADIUS IN M)/
+ 3 7H RTG ,E12.4,27H (GUIDE TUBE RADIUS IN M)/
+ 4 7H PITCH ,E12.4,24H (HEXAGONAL SIDE IN M)/
+ 5 7H ERMAXT,E12.4,35H (TEMPERATURE MAXIMUM ERROR IN K)/
+ 6 7H ERMAXC,E12.4,32H (FLOW MAXIMUM RELATIVE ERROR)/
+ 7 7H RELAX ,E12.4,25H (RELAXATION PARAMETER)/
+ 8 7H RTIME ,E12.4,20H (TIME VALUE IN S)/
+ 9 7H WTEFF ,E12.4,44H (SURFACE TEMPERATURE WEIGHTING FACTOR IN ,
+ 1 27HEFFECTIVE FUEL TEMPERATURE)/
+ 2 7H RATIOX,E12.4,35H (MAXIMUM OF VARIABLE VARIATIONS)/
+ 3 7H EPSR ,E12.4,34H (RUGOSITY IN M OF THE FUEL ROD)/
+ 4 7H THETA ,E12.4,41H (ANGLE IN RADIANS OF THE FUEL CHANNEL))
+ 580 FORMAT(7H HGAP ,1P,E12.4,20H (HGAP IN W/m^2/K))
+ 590 FORMAT(7H HCONV ,1P,E12.4,21H (HCONV IN W/m^2/K))
+ 600 FORMAT(7H HCONV ,1P,E12.4,22H (POWER FACTOR IN W))
+ 610 FORMAT(/37H THM: CORE COOLANT SECTION. COMPUTED=,1P,E9.2,
+ 1 7H GIVEN=,E9.2,4H m2.)
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 18:52:37 +0000