path: root/Donjon/src/THM.f

*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