*****************************************************************
*                                                               *
* Procedure :  Steady.c2m                                       *
* Purpose   :  reactor steady state computation including       *
*              neutron physics, thermalhydraulics, local xenon  *
*              effects and critical boron search feature        *
* Author    :  V. Salino                                        *
*                                                               *
* CALL      :                                                   *
*  Power Fmap Matex Thermo MicroF Burn := Steady                *
*        Fmap Matex Thermo MicroF Burn ReflXS Track ::          *
*        <<CritCB>> <<CBinit>>                                  *
*        >>CB<< >>keff<< ;                                      *
*                                                               *
* Parameters: input                                             *
*  CritCB    set to True to search critical boron concentration *
*            computation, False otherwise (logical variable).   *
*  CBinit    initialization for boron concentration (ppm) when  *
*            CritCB.EQ.True, otherwise boron concentration to   *
*            be used.                                           *
*  Prel      relative power=requested power/nominal power.      *
*  InletTmod inlet coolant temperature (Celsius).               *
*                                                               *
* Parameters: output                                            *
*  CB        critical boron concentration if CritCB.EQ.True,    *
*            can be negative in which case CB is extrapolated   *
*            and should not considered physical but fictive.    *
*  keff      k-effective, differs from 1 only if CB is          *
*            negative or if CritCB.EQ.False.                    *
*                                                               *
*****************************************************************
PARAMETER Power Flux Thermo MicroF Burn Fmap Matex SapUOX SapMOX 
          MacroRefl Track ::
  ::: LINKED_LIST Power Flux Thermo MicroF Burn Fmap Matex SapUOX 
                  SapMOX MacroRefl Track ; ;
MODULE SCR: NCR: MACINI: RESINI: TRIVAA: FLUD: FLPOW: THM: EVO: SIM:
       GREP: DELETE: ABORT: END: ;
*--
* Local objects and procedure
*--
LINKED_LIST FluxAv MacroTot MacroF System PowerIter 
            OldPower
             ;
PROCEDURE ComparePdist ;
*--
* Procedure arguments
*--
LOGICAL CritCB evo pred ;
INTEGER BUindex ;
REAL CB ;
  :: >>CritCB<< >>CB<< >>evo<< >>pred<< ;
*--
* Local variables
*--
REAL CBmax := 2000.0 ;
REAL CBmin :=    0.0 ;
REAL ErrorRho ErrorCB ErrorTFuel ErrorDCool ErrorTCool ErrorPdis ;
REAL EpsRho EpsCB EpsTFuel EpsDCool EpsPdis :=
     0.5    0.05  0.5      5.0E-5   5.0E-4   ;
INTEGER Iter := 1 ;
REAL CBinterp CBp1 CBm1 keff DeltaRho Rho Rhom1 ;
LOGICAL CONV := $False_L ;
REAL DiffBorWorth ;
*--
* Thermalhydraulics parameters
*--
REAL dx       := 21.5 ;
REAL Tot_tub  := 6.6E-03 2.0 ** $Pi_R * 25.0  * ;
REAL Tot_pin  := 4.7E-03 2.0 ** $Pi_R * 264.0 * ;
REAL asssect  := dx dx * 1.E-04 * Tot_tub - Tot_pin - ;
REAL coresect := 157.0 asssect * ;
REAL sass     := dx dx * 1.E-04 * ;
*--
* Initialized with a THM computation using a flat power distribution
*--
REAL Ptot := 2750.0 ;
REAL BundPow := Ptot 157.0 / 29.0 / 1.0E+3 * (* MW to kW *) ;
Fmap := RESINI: Fmap ::
  BUNDLE-POW SAME <<BundPow>> ;
*--
* Iteration loop with feedback, flux, xenon and simplified
* thermalhydraulics
*--
IF pred THEN
FluxAv := Flux ;
ENDIF ;
*--
REPEAT
*--
* Determine CB to use for interpolation
*--
  IF CB CBmin < THEN
    EVALUATE CBinterp := CBmin ;
  ELSEIF CB CBmax > THEN
    EVALUATE CBinterp := CBmax ;
  ELSE
    EVALUATE CBinterp := CB ;
  ENDIF ;
*--
* Thermalhydraulics computation, xenon saturation
*--
  IF Iter 1 = THEN
    Fmap := SIM: Fmap ::
    EDIT 1
    SET-PARAM 'C-BORE' <<CBinterp>> ;
  ELSE
  System := DELETE: System ;
    Thermo Fmap := THM: Thermo Fmap :: EDIT 0
      SET-PARAM 'C-BORE' <<CBinterp>> ;
    ! Evolution jusqu'a l'equilibre seulement si demande
    ! Sinon on va rechercher la CB crit sans etre a l'equilibre Xe
    IF evo THEN
      ECHO "bugg1" ;
      Burn MicroF := EVO: Burn MicroF PowerIter ::
        EDIT 0 FLUX_POW DEPL 100. DAY KEEP ;
      ECHO "bugg2" ;
      Burn MicroF := EVO: Burn MicroF PowerIter ::
        EDIT 0 FLUX_POW SAVE 100. DAY KEEP ;
    ENDIF ;
  ENDIF ;
*--
* Saphyb Interpolation
* Historic effects are neglected for simplicity's sake
*--
  MicroF := SCR: MicroF SapUOX SapMOX Fmap ::
      EDIT 0 RES
      MICRO LINEAR
      TABLE SapUOX 'BURN'
      MIX 4
         INST-BURN
        SET TF  526.85 ! Celsius
        SET TCA 326.85 ! Celsius
        SET DCA 0.659 ! g/cm3
        SET ppmB 500.0 ! ppm
        ADD 'TF' 526.85 MAP
          REF 'BURN' SAMEASREF ENDREF
        ADD 'TCA' 326.85 MAP
          REF 'BURN' SAMEASREF ENDREF
        ADD 'DCA' 0.659 MAP
          REF 'BURN' SAMEASREF ENDREF
        ADD 'ppmB' 500.0 MAP
          REF 'BURN' SAMEASREF ENDREF
        MICRO ALL
* Recherche poids Xenon : commenter/decommenter
!             'Xe135'   *
!             'I135'    0.
              'B10'     *
              'B11'     *
* Temporairement, Sm/Pm sont remis a leurs valeurs de Saphyb
* (equivalent a ne pas les particulariser). Apres premier
* appel a Steady, on pourrait enlever cela pour les laisser
* evoluer librement.
!             'Nd147'   0.
!             'Pm147'   0.
!             'Pm148'   0.
!             'Pm148m'  0.
!             'Pm149'   0.
!             'Sm149'   0.
             '*MAC*RES' 1.
      ENDMIX
      TABLE SapMOX 'BURN'
      MIX 5
         INST-BURN
        SET TF  526.85 ! Celsius
        SET TCA 326.85 ! Celsius
        SET DCA 0.659 ! g/cm3
        SET ppmB 500.0 ! ppm
        ADD 'TF' 526.85 MAP
          REF 'BURN' SAMEASREF ENDREF
        ADD 'TCA' 326.85 MAP
          REF 'BURN' SAMEASREF ENDREF
        ADD 'DCA' 0.659 MAP
          REF 'BURN' SAMEASREF ENDREF
        ADD 'ppmB' 500.0 MAP
          REF 'BURN' SAMEASREF ENDREF
        MICRO ALL
* Recherche poids Xenon : commenter/decommenter
!             'Xe135'   *
!             'I135'    0.
              'B10'     *
              'B11'     *
* Temporairement, Sm/Pm sont remis a leurs valeurs de Saphyb
* (equivalent a ne pas les particulariser). Apres premier
* appel a Steady, on pourrait enlever cela pour les laisser
* evoluer librement.
!             'Nd147'   0.
!             'Pm147'   0.
!             'Pm148'   0.
!             'Pm148m'  0.
!             'Pm149'   0.
!             'Sm149'   0.
             '*MAC*RES' 1.
      ENDMIX
      CHAIN
        B10      STABLE
        B11      STABLE
        I135     NG    0.
        Xe135    NG    0. FROM DECAY 1.0E+00 I135
        Nd147   STABLE
        Pm147   STABLE
        Pm148   STABLE
        Pm148m  STABLE
        Pm149   STABLE
        Sm149   STABLE
*--------- With Samarium -----------------------------------
             Nd147    NG    0.
             Pm147    NG    0. FROM DECAY 1.0E+00 Nd147
             Pm148    NG    0. FROM NG    5.3E-01 Pm147
             Pm148m   NG    0. FROM NG    4.7E-01 Pm147
             Pm149    NG    0. FROM NG    1.0E+00 Pm148
                                    NG    1.0E+00 Pm148m
             Sm149    NG    0. FROM DECAY 1.0E+00 Pm149
        MACT     NFTOT 0.
      ENDCHAIN
  ;
  MacroF := MicroF :: STEP UP "MACROLIB" ;
*--
* Flux and power computation
*--
  MacroTot Matex := MACINI: Matex MacroRefl MacroF ;
  MacroF := DELETE: MacroF ;
  System := TRIVAA: MacroTot Track :: EDIT 0 ;
  MacroTot := DELETE: MacroTot ;
  IF Iter 1 = THEN
    Flux := FLUD: Flux System Track :: EDIT 0
      ADI 4 EXTE 1.0E-5 ACCE 5 3 ;
  ELSE
    Flux := FLUD: Flux System Track :: EDIT 0
      RELAX 0.15 ;
  ENDIF ;
  GREP: Flux :: GETVAL 'K-EFFECTIVE ' 1  >>keff<<  ;
  EVALUATE Rho := 1. 1. keff / - 1.0E5 * ;
  IF Iter 1 > THEN
    PowerIter := DELETE: PowerIter ;
  ENDIF ;
  PowerIter Fmap := FLPOW: Fmap Flux Track Matex ::
    EDIT 0 PTOT <<Ptot>> PRINT DISTR POWER ;
*--
* Compute convergence errors (absolute values) on :
* - reactivity (pcm)
* - boron concentration (ppm)
* - 3D power distribution (relative discrepancy, maximum)
* - fuel temperature distribution (K, maximum)
* - moderator density distribution (g/cm3, maximum)
*--
  IF Iter 1 > THEN
    EVALUATE ErrorRho := Rho Rhom1 - ABS ;
    EVALUATE ErrorCB  := CB CBm1 - ABS ;
    Fmap OldPower PowerIter := ComparePdist :: >>ErrorPdis<< ;
    OldPower := DELETE: OldPower ;
    GREP: Thermo :: GETVAL 'ERROR-T-FUEL' 1  >>ErrorTFuel<<  ;
    GREP: Thermo :: GETVAL 'ERROR-D-COOL' 1  >>ErrorDCool<<  ;
  ENDIF ;
  OldPower := PowerIter ;
*--
* Reestimate critical CB. 3rd iteration is waited, in order to have a
* not-too-false k-effective, in particular regarding thermalhydraulics
* and neutron physics convergence. Otherwise predicted CB would be too
* far from the truth.
*--
  IF CritCB Iter 2 > * THEN
* To evaluate DiffBorWorth, we have 3 other possibilities :
* - Evaluate properly DiffBorWorth (internship ?), may be costly
* - Use the 2 previous computations to find (1 line with 2 points)
*   (may be instable during the first iterations but effective
*   near the end)
* - Use module FIND0 (may be long)
* Plot reactivity=f(CB), see which is going to be good
    EVALUATE DiffBorWorth := -7.0 ; ! pcm/ppm
    EVALUATE CBp1 := CBinterp Rho DiffBorWorth / - ;
  ENDIF ;
*--
* Prints
*--
  ECHO "#---" ;
  ECHO "Iter" Iter ;
  ECHO "keff=" keff "CB=" CB "Rho=" Rho ;
  ECHO "CBinterp=" CBinterp ;
  IF CritCB Iter 2 > * THEN
    ECHO "DiffBorWorth=" DiffBorWorth "CBp1" CBp1 ;
  ENDIF ;
  ECHO "---" ;
*--
* Check convergence
*--
  IF Iter 3 > THEN
    ECHO "Discrepancies between the previous iteration and current"
         "iteration" ;
    ECHO " Rho (pcm)     CB (ppm)      TFuel (K)     DCool (K)    "
         "Pdis (relat)" ;
    ECHO ErrorRho ErrorCB ErrorTFuel ErrorDCool ErrorPdis ;
    ECHO "" ;
    ECHO "Convergence criterion (epsilon)" ;
    ECHO " Rho (pcm)     CB (ppm)      TFuel (K)     DCool (K)    "
         "Pdis (relat)" ;
    ECHO EpsRho EpsCB EpsTFuel EpsDCool EpsPdis ;
    ECHO "" ;
    ECHO "Convergence reached ?" ;
    ECHO " Rho (pcm)     CB (ppm)      TFuel (K)     DCool (K)    "
         "Pdis (relat)" ;
    ECHO "" ErrorRho EpsRho < "      " ErrorCB EpsCB < "      "
         ErrorTFuel EpsTFuel < "       " ErrorDCool EpsDCool <
         "      " ErrorPdis  EpsPdis  < ;
    EVALUATE CONV := ErrorRho   EpsRho   <
                     ErrorCB    EpsCB    < *
                     ErrorTFuel EpsTFuel < *
                     ErrorDCool EpsDCool < *
                     ErrorPdis  EpsPdis  < * ;
  ENDIF ;
*--
* Prepare next iteration
*--
  IF CONV NOT THEN
    EVALUATE CBm1 Rhom1 := CB Rho ;
    IF CritCB Iter 2 > * THEN
      EVALUATE CB := CBp1 ;
    ENDIF ;
    EVALUATE Iter := Iter 1 + ;
  ENDIF ;
*--
  IF Iter 100 > THEN
    ECHO "Steady.c2m: maximum iteration reached (50)." ;
    ABORT: ;
  ENDIF ;
UNTIL Iter 4 = ;
*--
* Predicteur-correcteur
*--
IF pred THEN
  PowerIter := DELETE: PowerIter ;
  FluxAv := FLUD: FluxAv System Track :: EDIT 0
      RELAX 0.5 ;
  PowerIter Fmap := FLPOW: Fmap FluxAv Track Matex ::
    EDIT 0 PTOT <<Ptot>> PRINT DISTR POWER ;
  FluxAv := DELETE: FluxAv ;
ENDIF ;

*--
* Print Xe and I average fission yields. This average depends
* on burnup, thermalhydraulics and xenon conditions.
*--
REAL YieldI YieldXe ;
GREP: MicroF :: STEP UP 'DEPL-CHAIN'
                   GETVAL 'FISSIONYIELD' 1 >>YieldI<<
                   GETVAL 'FISSIONYIELD' 2 >>YieldXe<< ;
ECHO "Mean fission yield for I =" YieldI  ;
ECHO "Mean fission yield for Xe=" YieldXe ;

Power := PowerIter ;
PowerIter := DELETE: PowerIter ;
System := DELETE: System ;

  :: <<CB>> ;
  :: <<keff>> ; 

END: ;