Initial commit from Polytechnique Montreal

1 files changed, 628 insertions, 0 deletions

diff --git a/Donjon/src/THMDRV.f b/Donjon/src/THMDRV.f
new file mode 100644
index 0000000..ff642e3
--- /dev/null
+++ b/Donjon/src/THMDRV.f
@@ -0,0 +1,628 @@
+*DECK THMDRV
+      SUBROUTINE THMDRV(MPTHM,IMPX,IX,IY,NZ,XBURN,VOLXY,HZ,CFLUX,POROS,
+     > FNFU,NFD,NDTOT,IFLUID,SNAME,SCOMP,IGAP,IFUEL,FNAME,FCOMP,FCOOL,
+     > FFUEL,ACOOL,HD,PCH,RAD,
+     > MAXIT1,MAXITL,ERMAXT,SPEED,TINLET,POUTLET,
+     > FRACPU,ICONDF,NCONDF,KCONDF,UCONDF,ICONDC,NCONDC,KCONDC,UCONDC,
+     > IHGAP,KHGAP,IHCONV,KHCONV,WTEFF,IFRCDI,ISUBM,FRO,POW,IPRES,IDFM,
+     > TCOMB, DCOOL,TCOOL,TSURF,HCOOL,PCOOL)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Driver of the steady-state thermal-hydraulics calculation.
+*
+*Copyright:
+* Copyright (C) 2012 Ecole Polytechnique de Montreal.
+*
+*Author(s): 
+* A. Hebert
+* C. Garrido 
+* 08/2023: Modifications to include Molten Salt heat transfer in coolant
+* C. Garrido 
+* 07/2024: Modifications to include Molten Salt heat transfer in static
+*          fuel
+* C. Huet
+* 02/2025: Modifications to include pressure drop calculation
+* R. Guasch & M. Bellier
+* 08/2025: Modifications to include mass+momentum+energy conservation equation
+*  solution using a Drift-Flux Model.
+*
+*Parameters: input
+* MPTHM   directory of the THM object containing steady-state
+*         thermohydraulics data.
+* IMPX    printing index (=0 for no print).
+* IX      position of mesh along X direction.
+* IY      position of mesh along Y direction.
+* NZ      number of meshes along Z direction (channel direction).
+* XBURN   burnup distribution in MWday/tonne.
+* VOLXY   mesh area in the radial plane.
+* HZ      Z-directed mesh widths.
+* CFLUX   critical heat flux in W/m^2.
+* POROS   oxyde porosity.
+* FNFU    number of active fuel rods in the fuel bundle.
+* NFD     number of discretization points in fuel region.
+* NDTOT   number of total discretization points in the the fuel
+*         pellet and the cladding.
+* IFLUID  type of fluid (0=H2O; 1=D2O; 2=SALT).
+* SNAME   Name of the molten salt (e.g. "LiF-BeF2")
+* SCOMP   Composition of the molten salt (e.g. "0.66-0.34")
+* FCOOL   power density fraction in coolant.
+* FFUEL   power density fraction in fuel.
+* ACOOL   coolant cross section area in m^2.
+* HD      hydraulic diameter of one assembly in m.
+* PCH     heating perimeter in m.
+* RAD     fuel and clad radii in m.
+* MAXIT1  maximum number of conduction iterations.
+* MAXITL  maximum number of center-pellet iterations.
+* ERMAXT  convergence criterion.
+* SPEED   inlet flow velocity in m/s.
+* TINLET  inlet temperature in K.
+* POUTLET outlet pressure in Pa.
+* FRACPU  plutonium fraction in fuel.
+* ICONDF  fuel conductivity flag (0=Stora-Chenebault or COMETHE/
+*         1=user-provided polynomial + inverse term).
+* NCONDF  degree of user-provided fuel conductivity polynomial.
+* KCONDF  polynomial coefficients for fuel conductivity in W/m/K^(k+1)
+*         (except for the two last coefficients which belongs to the
+*         inverse term).
+* UCONDF  required unit of temperature in polynomial for fuel
+*         conductivity (KELVIN or CELSIUS).
+* ICONDC  clad conductivity flag (0=default/1=user-provided
+*         polynomial).
+* NCONDC  degree of user-provided clad conductivity polynomial.
+* KCONDC  polynomial coefficients for clad conductivity in W/m/K^(k+1).
+* UCONDC  required unit of temperature in polynomial for clad
+*         conductivity (KELVIN or CELSIUS).
+* IHGAP   flag indicating HGAP chosen (0=default/1=user-provided).
+* KHGAP   fixed user-provided HGAP value in W/m^2/K.
+* IHCONV  flag indicating HCONV chosen (0=default/1=user-provided).
+* KHCONV  fixed user-provided HCONV value in W/m^2/K.
+* WTEFF   surface temperature's weighting factor in effective fuel
+*         temperature.
+* IFRCDI  flag indicating if average approximation is forced during
+*         fuel conductivity evaluation (0=default/1=average
+*         approximation forced).
+* ISUBM   subcooling model (0: one-phase; 1: Jens-Lottes model;
+*         2: Saha- Zuber model).
+* FRO     radial power form factors.
+* POW     power distribution in W.
+* IGAP    Flag indicating if the gap is considered (0=gap/1=no gap)
+* IFUEL   type of fuel (0=UO2/MOX; 1=SALT).
+* FNAME   Name of the molten salt (e.g. "LiF-BeF2")
+* FCOMP   Composition of the molten salt (e.g. "0.66-0.34")
+* IPRES   flag indicating if pressure is to be computed (0=nonstant/
+*         1=variable).
+* IDFM    flag indicating if the drift flux model is to be used 
+*         (0=Without modifications(Chexal correlation for epsilon, no drift flux model in the Navier-Stokes equations)
+*           /1=EPRI/2=MODEBSTION/3=GERAMP/4=HEM1(VGJ=0)) 
+*
+*Parameters: output
+* TCOMB   averaged fuel temperature distribution in K.
+* DCOOL   coolant density distribution in g/cc.
+* TCOOL   coolant temperature distribution in K.
+* TSURF   surface fuel temperature distribution in K.
+* HCOOL   coolant enthalpty distribution in J/kg.
+* PCOOL   coolant pressure distribution in Pa.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+      USE t_saltdata
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) MPTHM
+      INTEGER IMPX,IX,IY,NZ,NFD,NDTOT,IFLUID,MAXIT1,MAXITL,IHGAP,IGAP,
+     > IFUEL,IPRES, IDFM
+      REAL XBURN(NZ),VOLXY,CFLUX,POROS,FRACPU,ERMAXT,
+     > SPEED,TINLET,POUTLET,
+     > FFUEL(NZ),ACOOL(NZ),RAD(NDTOT-1,NZ),FNFU(NZ),FCOOL(NZ),HZ(NZ),
+     > KCONDF(NCONDF+3),KCONDC(NCONDC+1),KHGAP,KHCONV,WTEFF,FRO(NFD-1),
+     > POW(NZ),TCOMB(NZ),DCOOL(NZ),TCOOL(NZ),TSURF(NZ),HCOOL(NZ),
+     > PCOOL(NZ),MUT(NZ), HD(NZ), PCH(NZ)
+      CHARACTER UCONDF*12,UCONDC*12
+*----
+*  LOCAL VARIABLES
+*----
+      TYPE(tpdata) STP,FTP
+      PARAMETER (KMAXO=100,MAXNPO=40)
+      REAL TRE11(MAXNPO),RADD(MAXNPO),ENT(4),MFLOW,TLC(NZ)
+      CHARACTER HSMG*131,SNAME*32,SCOMP*32,FNAME*32,FCOMP*32
+      REAL XS(4),TC1,PC(NZ),TP(NZ),RHOL,XFL(NZ),EPS(NZ),HINLET,
+     > TCLAD(NZ),ENTH(NZ),SLIP(NZ),AGM(NZ),QFUEL(NZ),QCOOL(NZ),K11,
+     > VLIQ(NZ),VVAP(NZ)
+      INTEGER KWA(NZ)
+      REAL XX2(MAXNPO),XX3(MAXNPO),ZF(2)
+      DATA XS/-0.861136,-0.339981,0.339981,0.861136/
+
+      REAL TBUL(NZ),VGJprime(NZ),HLV(NZ),DGCOOL(NZ),DLCOOL(NZ)
+
+      INTEGER I
+      REAL PINLET, ERRV, ERRP, ERRD, NORMV, NORMP, NORMD
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      REAL, ALLOCATABLE, DIMENSION(:) :: VCOOL,TCENT
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: TEMPT
+
+      REAL, ALLOCATABLE, DIMENSION(:) :: PTEMP, VTEMP, DTEMP
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(VCOOL(NZ),TEMPT(NDTOT,NZ),TCENT(NZ))
+      ALLOCATE(PTEMP(NZ), VTEMP(NZ), DTEMP(NZ))
+*----
+*  COMPUTE THE INLET FLOW ENTHALPY AND VELOCITY
+*   INITIALIZE PINLET TO POUTLET, WILL BE UPDATED IF IPRES=1
+*   ELSE PINLET = POUTLET
+*----
+      PINLET = POUTLET
+      IF(NDTOT.GT.MAXNPO) CALL XABORT('THMDRV: MAXNPO OVERFLOW.')
+      IF(IFLUID.EQ.0) THEN
+         CALL THMSAT(PINLET,TSAT)
+      ELSE IF(IFLUID.EQ.1) THEN
+         CALL THMHST(PINLET,TSAT)
+*CGT TODO: GET ALSO FREEZING??
+      ELSE IF(IFLUID.EQ.2) THEN
+         CALL THMSGT(SNAME,SCOMP,STP,IMPX)
+         CALL THMSST(STP,TSAT,IMPX)
+*CGT
+      ENDIF
+      IF (IFUEL.EQ.1) THEN
+         CALL THMSGT(FNAME,FCOMP,FTP,IMPX)
+      ENDIF
+
+      IF(TINLET.GT.TSAT) THEN
+         WRITE(HSMG,'(27HTHMDRV: INLET TEMPERATURE (,1P,E12.4,
+     >   40H K) GREATER THAN SATURATION TEMPERATURE.)') TINLET
+         CALL XABORT(HSMG)
+      ENDIF
+      IF(IFLUID.EQ.0) THEN
+        CALL THMPT(PINLET,TINLET,RHOIN,HINLET,R3,R4,R5)
+      ELSE IF(IFLUID.EQ.1) THEN
+        CALL THMHPT(PINLET,TINLET,RHOIN,HINLET,R3,R4,R5)
+      ELSE IF(IFLUID.EQ.2) THEN
+        CALL THMSPT(STP,TINLET,RHOIN,HINLET,R3,R4,R5,IMPX)
+      ENDIF
+      MFLOW=SPEED*RHOIN
+      HMSUP=HINLET
+*----
+*  INITIALIZE VALUES OF STEAM QUALITIES, VOID FRACTION AND DENSITY
+*  PRESSURE, VELOCITY AND TEMPERATURE OF THE COOLANT ALONG THE CHANNEL.
+*---
+      DO K=1,NZ
+         EPS(K)=0.0
+         XFL(K)=0.0
+         SLIP(K)=1.0
+         KWA(K)=0
+         MUT(K)=0.0
+         QFUEL(K)=0.0
+         VGJprime(K)=0.0
+         HLV(K)=0.0
+
+         PCOOL(K)=PINLET
+         VCOOL(K)=MFLOW/RHOIN
+         DCOOL(K)=RHOIN
+         DLCOOL(K)=RHOIN
+         DGCOOL(K)=0.0 
+         TCOOL(K)=TINLET
+         HCOOL(K)=HINLET
+*----
+*  COMPUTE THE SATURATION TEMPERATURE AND THE THERMODYNAMIC PROPERTIES
+*  IF THE PRESSURE DROP IS COMPUTED
+*---
+
+        IF (IPRES.EQ.1) THEN
+          IF(POW(K).EQ.0.0) CYCLE
+          IF(IFLUID.EQ.0) THEN
+            CALL THMSAT(PCOOL(K),TSAT)
+          ELSE IF(IFLUID.EQ.1) THEN
+            CALL THMHST(PCOOL(K),TSAT)
+          ENDIF
+
+          TB=TSAT-0.1
+          IF(TCOOL(K).LT.TB) THEN
+            IF(IFLUID.EQ.0) THEN
+              CALL THMPT(PCOOL(K),TCOOL(K),RHOIN,H11,K11,MUT(K),C11)
+            ELSE IF(IFLUID.EQ.1) THEN
+              CALL THMHPT(PCOOL(K),TCOOL(K),RHOIN,H11,K11,MUT(K),C11)
+            ELSE IF(IFLUID.EQ.2) THEN
+              CALL THMSPT(STP,TCOOL(K),R11,H11,K11,MUT(K),C11,IMPX)
+            ENDIF
+          ELSE
+            IF(IFLUID.EQ.0) THEN
+              CALL THMPT(PCOOL(K),TB,R11,H11,K11,MUT(K),C11)
+            ELSE IF(IFLUID.EQ.1) THEN
+              CALL THMHPT(PCOOL(K),TB,R11,H11,K11,MUT(K),C11)
+            ELSE IF(IFLUID.EQ.2) THEN
+              CALL THMSPT(STP,TB,R11,H11,K11,MUT(K),C11,IMPX)
+            ENDIF
+          ENDIF
+        ENDIF
+      ENDDO
+*----
+*  MAIN LOOP ALONG THE 1D CHANNEL.
+*----
+      ERRV = 1.0
+      ERRP = 1.0
+      ERRD = 1.0
+      NORMP = PINLET
+      NORMV = SPEED
+      NORMD = RHOIN
+      I=0
+      IF (IPRES .EQ. 0) GOTO 30
+   10 CONTINUE
+*----
+*  UPDATE HINLET FUNCTION OF INLET PRESSURE AND TEMPERATURE
+*----
+       HMSUP=HINLET
+       SPEED=MFLOW/DCOOL(1)
+*----
+*  WHILE LOOP FOR PRESSURE AND VELOCITY CONVERGENCE
+*  CHECK FOR CONVERGENCE
+*----
+        IF (I .GT. 1000) GOTO 20
+        IF ((ERRP.LT.5E-4).AND.(ERRV.LT.5E-4).AND.(IDFM.EQ.0)) GOTO 20
+
+        IF ((IDFM.GT.0).AND.(I.GT.3)) THEN
+        IF ((ERRP.LT.5E-4).AND.(ERRV.LT.5E-4).AND.(ERRD.LT.5E-4)) THEN
+        GOTO 20
+        ENDIF
+        ENDIF
+
+          I = I + 1
+
+          PTEMP = PCOOL
+          VTEMP = VCOOL
+          DTEMP = DCOOL
+
+          SPEED = MFLOW/DCOOL(1)
+          CALL THMPV(SPEED, PCOOL(NZ), VCOOL, DCOOL, 
+     >              PCOOL, TCOOL, MUT, XFL, HD, NZ,
+     >              HZ, EPS, DLCOOL,DGCOOL, VGJprime, IDFM, ACOOL)
+* Extrapolate from first two values of PCOOL to get PINLET at first face.
+* This ensures that computed HINLET is not HCOOL(1)
+          PINLET = (3.0/2.0)*PCOOL(1) - (1.0/2.0)*PCOOL(2)
+          IF (IFLUID.EQ.0) THEN
+            CALL THMPT(PINLET, TINLET, RHOIN, HINLET, R3, R4, R5)
+          ELSE IF(IFLUID.EQ.1) THEN
+            CALL THMHPT(PINLET,TINLET,RHOIN,HINLET,R3,R4,R5)
+          ELSE IF(IFLUID.EQ.2) THEN
+            CALL THMSPT(STP,TINLET,RHOIN,HINLET,R3,R4,R5,IMPX)
+          ENDIF
+*   Update inlet enthalpy based on computed inlet pressure.
+          HMSUP = HINLET
+   30 CONTINUE
+*----
+*  MAIN LOOP ALONG THE 1D CHANNEL.
+*----
+      K0=0 ! onset of nucleate boiling point
+      DO K=1,NZ
+        IF(POW(K).EQ.0.0) CYCLE
+        IF(IFLUID.EQ.0) THEN
+          CALL THMSAT(PCOOL(K),TSAT)
+        ELSE IF(IFLUID.EQ.1) THEN
+          CALL THMHST(PCOOL(K),TSAT)
+        ENDIF
+        TBUL(K)=TSAT
+*----
+*  COMPUTE THE LINEAR POWER, THE VOLUMIC POWER AND THE THERMAL EXCHANGE
+*  COEFFICIENT OF THE GAP
+*----
+        DV=VOLXY*HZ(K)
+*       linear power in W/m
+        POWLIN=(POW(K)/DV)*VOLXY/FNFU(K)
+*       volumic power in W/m^3
+        QFUEL(K)=POW(K)*FFUEL(K)/DV
+        QCOOL(K)=POW(K)*FCOOL(K)/DV
+*----
+*  INITIALIZATION OF PINCELL TEMPERATURES
+*----
+        IF(POW(K).EQ.0.0) CYCLE
+        IF(IMPX.GT.4) WRITE(6,190) K
+        DO L=1,NDTOT
+           TRE11(L)=TCOMB(K)
+        ENDDO
+        DO L=1,NDTOT-1
+           RADD(L)=RAD(L,K)
+        ENDDO
+*----
+*  COMPUTE THE POWER DENSITY AND HEAT FLOW ALONG THE CHANNEL
+*----
+*       out-of-clad heat flow in W/m2
+        IF(IMPX.GT.5) WRITE(6,'(15H THMDRV: QFUEL(,I5,2H)=,1P,E12.4,
+     >  6H W/m2.)') K,QFUEL(K)
+        PHI2=0.5*QFUEL(K)*RAD(NFD,K)**2/RAD(NDTOT-1,K)
+        IF(PHI2.GT.CFLUX) THEN
+          WRITE(HSMG,'(23HTHMDRV: THE HEAT FLUX (,1P,E12.4,5H) IS ,
+     >    37HGREATER THAN THE CRITICAL HEAT FLUX (,E12.4,2H).)')
+     >    PHI2,CFLUX
+          CALL XABORT(HSMG)
+        ENDIF
+*----
+*  COMPUTE FOUR VALUES OF ENTHALPY IN J/KG TO BE USED IN GAUSSIAN
+*  INTEGRATION. DELTH1 IS THE ENTHALPY INCREASE IN EACH AXIAL MESH.
+*---- 
+        IF (IDFM.EQ.0) THEN
+          DELTH1=(PCH(K)/ACOOL(K)*PHI2+QCOOL(K))*HZ(K)/MFLOW
+        ELSE
+          DELTH1= (PCH(K)/ACOOL(K)*PHI2+QCOOL(K))*HZ(K)*ACOOL(K)
+        ENDIF
+        IF ((K.GT.1).AND.(IDFM.GT.0)) THEN
+          DELTH1= (PCH(K)/ACOOL(K)*PHI2+QCOOL(K))*HZ(K)*ACOOL(K)
+          DELTH1 = DELTH1 + ((VCOOL(K-1) + EPS(K-1)*(DLCOOL(K-1)-
+     >      DGCOOL(K-1))/DCOOL(K-1)*VGJprime(K-1))
+     >      + (VCOOL(K) + EPS(K)*(DLCOOL(K)-DGCOOL(K))/
+     >      DCOOL(K)*VGJprime(K)))/2*(PCOOL(K-1)*ACOOL(K-1)-PCOOL(K)
+     >      *ACOOL(K))
+          DELTH1 = DELTH1 +(EPS(K-1)*DGCOOL(K-1)*(DLCOOL(K-1)/
+     >      DCOOL(K-1))*HLV(K-1)*VGJprime(K-1)*ACOOL(K-1))-(EPS(K)*
+     >      DGCOOL(K)*(DLCOOL(K)/DCOOL(K))*HLV(K)*VGJprime(K)*ACOOL(K))   
+          DELTH1 = DELTH1/MFLOW/ACOOL(K)
+        ENDIF
+        DO I1=1,4
+          POINT=(1.0+XS(I1))/2.0
+          ENT(I1)=HMSUP+POINT*DELTH1
+        ENDDO
+        HMSUP=HMSUP+DELTH1
+*----
+*  COMPUTE THE VALUE OF THE DENSITY AND THE CLAD-COOLANT HEAT TRANSFER
+*  COEFFICIENT
+*----
+        IF(K.GT.1) THEN
+          XFL(K)=XFL(K-1)
+          EPS(K)=EPS(K-1)
+          SLIP(K)=SLIP(K-1)
+        ENDIF
+*CGT
+        IF ((IFLUID.EQ.0).OR.(IFLUID.EQ.1)) THEN
+          CALL THMH2O(0,IX,IY,K,K0,PCOOL(K),MFLOW,HMSUP,ENT,HD(K),
+     >    IFLUID,IHCONV,KHCONV,ISUBM,RAD(NDTOT-1,K),ZF,VCOOL(K),
+     >    IDFM,PHI2,XFL(K),EPS(K),SLIP(K),ACOOL(K),PCH(K),HZ(K),TCALO,
+     >    RHO,RHOL,RHOG,TRE11(NDTOT),KWA(K),VGJprime(K),HLV(K))
+        ELSEIF (IFLUID.EQ.2) THEN
+          CALL THMSAL(IMPX,0,IX,IY,K,K0,MFLOW,HMSUP,ENT,HD(K),STP,
+     >    IHCONV,KHCONV,ISUBM,RAD(NDTOT-1,K),ZF,PHI2,XFL(K),
+     >    EPS(K),SLIP(K),HZ(K),TCALO,RHO,RHOL,TRE11(NDTOT),
+     >    KWA(K))
+        ENDIF
+*CGT
+*----
+*  STEADY-STATE SOLUTION OF THE CONDUCTION EQUATIONS IN A FUEL PIN.
+*----
+        DTINV=0.0
+        IF(IGAP.EQ.0) THEN
+          CALL THMROD(IMPX,NFD,NDTOT-1,MAXIT1,MAXITL,ERMAXT,DTINV,RADD,
+     >    TRE11,TRE11,QFUEL(K),FRO,TRE11(NDTOT),POWLIN,XBURN(K),
+     >    POROS,FRACPU,ICONDF,NCONDF,KCONDF,UCONDF,ICONDC,NCONDC,
+     >    KCONDC,UCONDC,IHGAP,KHGAP,IFRCDI,TC1,XX2,XX3,ZF)
+        ELSE
+          CALL THMRNG(IMPX,NFD,NDTOT-1,MAXIT1,MAXITL,ERMAXT,DTINV,RADD,
+     >    TRE11,TRE11,QFUEL(K),FRO,TRE11(NDTOT),XBURN(K),
+     >    POROS,FRACPU,ICONDF,NCONDF,KCONDF,UCONDF,ICONDC,NCONDC,
+     >    KCONDC,UCONDC,IFRCDI,IFUEL,FTP,TC1,XX2,XX3,ZF)
+        ENDIF
+*
+        DO K1=1,NDTOT-1
+          TRE11(K1)=XX2(K1)+TRE11(NDTOT)*XX3(K1)
+        ENDDO
+*----
+*  RECOVER MESHWISE TEMPERATURES AND FLUID DENSITY. BY DEFAULT, USE THE
+*  ROWLANDS FORMULA TO COMPUTE THE EFFECTIVE FUEL TEMPERATURE, OTHERWISE
+*  USE USER-SPECIFIED WEIGHTING FACTOR.
+*----  
+        TCOMB(K)=(1.0-WTEFF)*TC1+WTEFF*TRE11(NFD)
+        TCENT(K)=TC1
+        TSURF(K)=TRE11(NFD)
+        TCLAD(K)=TRE11(NDTOT)
+        TCOOL(K)=TCALO
+        DCOOL(K)=RHO
+        DLCOOL(K)=RHOL
+        HCOOL(K)=HMSUP
+        PC(K)=PINLET
+        TP(K)=TCLAD(K)
+        TLC(K)=TCOOL(K)
+        ENTH(K)=HCOOL(K)
+        AGM(K)=MFLOW ! constant flow rate
+        DO K2=1,NDTOT
+          TEMPT(K2,K)=TRE11(K2)
+        ENDDO
+        IF (IPRES .EQ. 0) THEN
+          PCOOL(K)=PINLET
+          VCOOL(K)=MFLOW/DCOOL(K)
+        ENDIF
+*----
+*  COMPUTE THE SATURATION TEMPERATURE AND THE THERMODYNAMIC PROPERTIES
+*  IF THE PRESSURE DROP IS COMPUTED
+*---
+        IF (IPRES.EQ.1) THEN
+          IF(POW(K).EQ.0.0) CYCLE
+          IF(IFLUID.EQ.0) THEN
+            CALL THMSAT(PCOOL(K),TSAT)
+          ELSE IF(IFLUID.EQ.1) THEN
+            CALL THMHST(PCOOL(K),TSAT)
+          ENDIF
+ 
+          TB=TSAT-0.1
+          IF(TCOOL(K).LT.TB) THEN
+            IF(IFLUID.EQ.0) THEN
+             CALL THMPT(PCOOL(K),TCOOL(K),RHOIN,H11,K11,MUT(K),C11)
+            ELSE IF(IFLUID.EQ.1) THEN
+             CALL THMHPT(PCOOL(K),TCOOL(K),RHOIN,H11,K11,MUT(K),C11)
+            ELSE IF(IFLUID.EQ.2) THEN
+             CALL THMSPT(STP,TCOOL(K),R11,H11,K11,MUT(K),C11,IMPX)
+            ENDIF
+          ELSE
+            IF(IFLUID.EQ.0) THEN
+             CALL THMPT(PCOOL(K),TB,R11,H11,K11,MUT(K),C11)
+            ELSE IF(IFLUID.EQ.1) THEN
+             CALL THMHPT(PCOOL(K),TB,R11,H11,K11,MUT(K),C11)
+            ELSE IF(IFLUID.EQ.2) THEN
+             CALL THMSPT(STP,TB,R11,H11,K11,MUT(K),C11,IMPX)
+            ENDIF
+          ENDIF
+         ENDIF
+      ENDDO
+*----
+* IF THE PRESSURE DROP IS COMPUTED, COMPUTE THE 
+* THE PRESSURE AND VELOCITY RESIDUALS
+* IF DFM IS ACTIVATED, COMPUTE DCOOL RESIDUALS
+*----
+      IF (IPRES .EQ. 0) GOTO 20
+      ERRV = 0.0
+      ERRP = 0.0
+      ERRD = 0.0
+      NORMV = 0.0
+      NORMP = 0.0
+      NORMD = 0.0
+
+      DO K=1,NZ
+*         Under relaxation of coolant pressure and velocity.
+          VCOOL(K) = 0.40*VCOOL(K) + (1.00-0.40)*VTEMP(K)
+          PCOOL(K) = 0.40*PCOOL(K) + (1.00-0.40)*PTEMP(K)
+          ERRV = ERRV + (VCOOL(K)-VTEMP(K))**2
+          NORMV = NORMV + VCOOL(K)**2
+          ERRP = ERRP + (PCOOL(K)-PTEMP(K))**2
+          NORMP = NORMP + PCOOL(K)**2
+        IF (IDFM.GT.0) THEN
+*         Under relaxation of coolant density.
+          DCOOL(K) = 0.40*DCOOL(K) + (1.00-0.40)*DTEMP(K)
+          ERRD = ERRD + (DCOOL(K) - DTEMP(K))**2
+          NORMD = NORMD + DCOOL(K)**2
+        ENDIF
+      ENDDO
+      NORMV = SQRT(NORMV)
+      NORMP = SQRT(NORMP)
+      ERRV = SQRT(ERRV) / NORMV
+      ERRP = SQRT(ERRP) / NORMP
+      IF (IDFM.GT.0) THEN
+        NORMD = SQRT(NORMD)
+        ERRD = SQRT(ERRD) / NORMD
+      ENDIF
+      GO TO 10
+
+   20 CONTINUE
+
+      IF (I.GE.1000) THEN
+        PRINT *, 'ERRV =', ERRV
+        PRINT *, 'ERRP =', ERRP
+        PRINT *, 'ERRD =', ERRD
+        CALL XABORT('THMDRV: MAXIMUM NB OF ITERATIONS REACHED.')
+      ELSE IF(IMPX.GT.0) THEN
+        WRITE(6,'(37H THMDRV: CONVERGENCE REACHED AT ITER=,I5,1H.)') I
+      ENDIF
+
+*----
+* RECONSTRUCT THE PHASE VELOCITIES FROM VCOOL, EPS and VGJ
+*----
+      DO K=1,NZ
+        IF (IDFM .GT. 0) THEN
+          VLIQ(K) = VCOOL(K) - (1.0/(1.0- EPS(K)) - DLCOOL(K)/DCOOL(K))
+     >     * VGJprime(K)
+          VVAP(K) = VCOOL(K) + DLCOOL(K)/DCOOL(K) *VGJprime(K)
+        ELSE
+          VLIQ(K) = VCOOL(K)
+          VVAP(K) = VCOOL(K) 
+        ENDIF
+      ENDDO     
+*----
+* PRINT THE THERMOHYDRAULICAL PARAMETERS
+*----
+      IF(IMPX.GT.4) THEN
+        WRITE(6,250) 'POW', POW(:NZ)
+        WRITE(6,250) 'PCOOL', PCOOL(:NZ)
+        WRITE(6,250) 'VCOOL', VCOOL(:NZ)
+        WRITE(6,250) 'DCOOL', DCOOL(:NZ)
+        WRITE(6,250) 'TCOOL', TCOOL(:NZ)
+        WRITE(6,250) 'EPS', EPS(:NZ)
+        WRITE(6,250) 'XFL', XFL(:NZ)
+        WRITE(6,250) 'TSAT', TBUL(:NZ)
+        WRITE(6,250) 'MUT', MUT(:NZ)
+      ENDIF
+*----
+*  PRINT THE OUTLET THERMOHYDRAULICAL PARAMETERS
+*----
+      IF(IMPX.GT.3) THEN
+        WRITE(6,'(/16H THMDRV: CHANNEL,2I6/1X,27(1H-))') IX,IY
+        WRITE(6,210) ' ____________________________________________',
+     >          '_____________________________________________________',
+     >          '_____________________________________________________',
+     >          '______________'
+        WRITE(6,210) '|     |   TCOMB    |   TSURF    |    DCOOL ',
+     >          '   |    TCOOL    |    PCOOL    |    HCOOL    |    ',
+     >          'QFUEL    |    QCOOL    |    VOID   |     QUAL    |',
+     >          '     SLIP    |  FLOW  |',
+     >          '|     |     K      |     K      |    Kg/m3    |   ',
+     >          '   K      |     Pa      |    J/Kg     |    W/m3   ',
+     >          '  |    W/m3     |           |             |       ',
+     >          '      | REGIME |'
+        WRITE(6,210) '|_____|____________|____________|____________',
+     >          '_|_____________|_____________|_____________|_________',
+     >          '____|_____________|___________|_____________|________',
+     >          '_____|________|'
+        DO L=NZ,1,-1
+          IF(L.EQ.1) THEN
+            WRITE(6,220) '| BOT |',TCOMB(L),' |',TSURF(L),
+     >            ' |',DCOOL(L),' |',TCOOL(L),' |',PCOOL(L),
+     >            ' |',HCOOL(L),' |',QFUEL(L),' |',QCOOL(L),' |',
+     >            EPS(L),' |',XFL(L),' |',SLIP(L),' |',KWA(L),' |'
+          ELSEIF(L.EQ.NZ) THEN
+            WRITE(6,220) '| TOP |',TCOMB(L),' |',TSURF(L),
+     >            ' |',DCOOL(L),' |',TCOOL(L),' |',PCOOL(L),
+     >            ' |',HCOOL(L),' |',QFUEL(L),' |',QCOOL(L),' |',
+     >            EPS(L),' |',XFL(L),' |',SLIP(L),' |',KWA(L),' |'
+          ELSE
+            WRITE(6,230) '| ',L,' |',TCOMB(L),' |',TSURF(L),
+     >            ' |',DCOOL(L),' |',TCOOL(L),' |',PCOOL(L),
+     >            ' |',HCOOL(L),' |',QFUEL(L),' |',QCOOL(L),' |',
+     >            EPS(L),' |',XFL(L),' |',SLIP(L),' |',KWA(L),' |'
+          ENDIF
+        ENDDO
+        WRITE(6,210) '|_____|____________|____________|____________',
+     >          '_|_____________|_____________|_____________|_________',
+     >          '____|_____________|___________|_____________|________',
+     >          '_____|________|'
+        WRITE(6,240) MFLOW
+      ENDIF
+*----
+*  MODIFICATION OF THE VECTORS TO FIT THE GEOMETRY OF THE CHANNELS AND
+*  THE BUNDLES AND WRITE THE DATA IN LCM OBJECT THM 
+*----
+      CALL LCMPUT(MPTHM,'PRESSURE',NZ,2,PCOOL)
+      CALL LCMPUT(MPTHM,'DENSITY',NZ,2,DCOOL)
+      CALL LCMPUT(MPTHM,'LIQUID-DENS',NZ,2,DLCOOL)
+      CALL LCMPUT(MPTHM,'ENTHALPY',NZ,2,HCOOL)
+      CALL LCMPUT(MPTHM,'VELOCITIES',NZ,2,VCOOL)
+      CALL LCMPUT(MPTHM,'V-LIQ',NZ,2,VLIQ)
+      CALL LCMPUT(MPTHM,'V-VAP',NZ,2,VVAP)
+      CALL LCMPUT(MPTHM,'EPSILON',NZ,2,EPS)
+      CALL LCMPUT(MPTHM,'EPSOUT',1,2,EPS(NZ))
+      CALL LCMPUT(MPTHM,'XFL',NZ,2,XFL)
+      CALL LCMPUT(MPTHM,'CENTER-TEMPS',NZ,2,TCENT)
+      CALL LCMPUT(MPTHM,'COOLANT-TEMP',NZ,2,TCOOL)
+      CALL LCMPUT(MPTHM,'POWER',NZ,2,POW)
+      CALL LCMPUT(MPTHM,'TEMPERATURES',NDTOT*NZ,2,TEMPT)
+      CALL LCMPUT(MPTHM,'PINLET',1,2,PINLET)
+      CALL LCMPUT(MPTHM,'TINLET',1,2,TINLET)
+      CALL LCMPUT(MPTHM,'VINLET',1,2,SPEED)
+      CALL LCMPUT(MPTHM,'POULET',1,2,POUTLET)
+      CALL LCMPUT(MPTHM,'RADII',(NDTOT-1)*NZ,2,RAD)
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(TCENT,TEMPT,VCOOL)
+      DEALLOCATE(PTEMP, VTEMP, DTEMP)
+      RETURN
+*
+  190 FORMAT(/21H THMDRV: AXIAL SLICE=,I5)
+  210 FORMAT(1X,A,A,A,A)
+  220 FORMAT(1X,A,F11.2,A,F11.2,A,F12.4,A,F12.2,A,3P,E12.4,
+     >       A,1P,E12.4,A,1P,E12.4,A,1P,E12.4,A,0P,F10.4,A,E12.4,A,
+     >       E12.4,A,I5,2X,A)
+  230 FORMAT(1X,A,I3,A,F11.2,A,F11.2,A,F12.4,A,F12.2,A,3P,E12.4,
+     >       A,1P,E12.4,A,1P,E12.4,A,1P,E12.4,A,0P,F10.4,A,E12.4,A,
+     >       E12.4,A,I5,2X,A)
+  240 FORMAT(7H MFLXT=,1P,E12.4,8H Kg/m2/s)
+  250 FORMAT(9H THMDRV: ,A6,1H:,1P,11E12.4/(4X,12E12.4))
+      END