*----
*  Name          : rep900_mphy.x2m
*  Type          : DONJON CLE2000 procedure
*
*  Local parameters:
*    densB1   boron density for the first calculation
*    densB2   boron density for the second calculation
*
*  Calculation options:
*    Splitx   Radial mesh-splitting of assemblies
*    Splitz   Axial mesh-splitting of assemblies
*    TypeCal  Boron concentration computation option: "critical" or
*             "fixed"
*
*----
*  STRUCTURES, MODULES and PROCEDURES definitions
*----
LINKED_LIST GeoCoeur Track Flux Cpo L0 Reseau Matex
            Lib Macro2 Fmap ;
MODULE      TRIVAT: GREP: FIND0: UTL: DELETE: ABORT: END: ;

PROCEDURE   GetDim GeoCo SetFuelMap SetParam ThermoParam ;
PROCEDURE   assertS ;
SEQ_ASCII   _ACompo ;
REAL        DELTA REFVALUE ;
*----
*  Local Parameters
*----
REAL    densB1 := 2000.0 ;
REAL    densB2 := 1000.0 ;
*----
*  Calculation options
*----
INTEGER Splitx Splitz := 1 1        ;
STRING  TypeCal       := "critical" ;
REAL    powi := 2775.0 ;

REAL Rcomb1 Rcomb2 Rcomb3 Rcomb4 
     R_int_TG R_ext_TG R_int_TI R_ext_TI R_int_G R_ext_G 
     Cote CoteL
     RmodeN1 RmodeN2 
     R_ext_G2 R_ext_TGm R_ext_TIm
     N_U5_UOX N_U5_MOX
     N_Pu8 N_Pu9 N_Pu0 N_Pu1 N_Pu2 N_Am1
     densU_UOX densU_MOX densPu_MOX
     teneur1_Pu teneur2_Pu teneur3_Pu
     temp_comb temp_mode dens_mode
     temp_comb_0 temp_mode_0 dens_mode_0 ;

*----
*  Recover the Multicompo
*----
Cpo := _ACompo ;
UTL: Cpo :: DIR ;

REAL ORF IRC ORC IRCGT ORCGT ORAG ORAN 
     IRCA ORCA ORDET ORDETC Lame LPitch Height ;

REAL pbore_0 := 500. ; ! natural boron concentration in ppm

  GetDim ::
    >>Rcomb1<< >>Rcomb2<< >>Rcomb3<< >>Rcomb4<<
    >>R_int_TG<< >>R_ext_TG<< >>R_int_TI<< >>R_ext_TI<<
    >>R_int_G<< >>R_ext_G<< >>Cote<< >>Lame<< >>CoteL<<
    >>RmodeN1<< >>RmodeN2<< >>R_ext_G2<< >>R_ext_TGm<< >>R_ext_TIm<<
    >>N_U5_UOX<< >>N_U5_MOX<< >>N_Pu8<< >>N_Pu9<< >>N_Pu0<<
    >>N_Pu1<< >>N_Pu2<< >>N_Am1<< >>densU_UOX<< >>densU_MOX<<
    >>densPu_MOX<< >>teneur1_Pu<< >>teneur2_Pu<< >>teneur3_Pu<<
    >>temp_comb_0<< >>temp_mode_0<< >>dens_mode_0<< ;

EVALUATE temp_comb temp_mode dens_mode := 800.0 600.0 0.659 ;

EVALUATE
     ORF IRC ORC IRCGT ORCGT
     ORAG ORAN IRCA ORCA
     ORDET ORDETC
     Lame LPitch Height
     :=
     Rcomb4 R_int_TG R_ext_TG R_int_TI R_ext_TI
     R_int_TG R_ext_TG R_int_TI R_ext_TI
     R_int_TG R_ext_TG
     Lame Cote 10.0
     ;

REAL rPavX := 17.0 ; ! assemblage 17 x 17
INTEGER iPavX := 17 ; ! assemblage 17 x 17

REAL    dx   := rPavX LPitch * 2.0 Lame * + ;
INTEGER MaxR := iPavX iPavX * 33 * Splitx * Splitx * Splitz * ;
INTEGER Iter := 1 ;
REAL keff1 keff2 keffn ;
REAL Y1 Y2 ROOT YNEW ;
REAL densB ;
LOGICAL CONV ;

GeoCoeur Matex := GeoCo  ::
   <<Splitx>> <<Splitz>> <<MaxR>> <<dx>> ;

Track := TRIVAT: GeoCoeur ::
   EDIT 1 MAXR <<MaxR>> DUAL 2 3 ;

*--
* Fuel map definition
*--
Fmap Matex := SetFuelMap Matex :: <<dx>> ;

Fmap := SetParam Fmap :: <<densB1>> <<temp_comb_0>> <<temp_mode_0>>
        <<dens_mode_0>> ;

*--
* Calculation
*--
Fmap Matex Flux := ThermoParam Fmap Matex Cpo Track
 :: <<powi>> <<densB1>> <<pbore_0>> <<temp_comb_0>> <<temp_mode_0>>
    <<dens_mode_0>> ;

GREP: Flux   ::
  GETVAL 'K-EFFECTIVE ' 1  >>keff1<<  ;

ECHO "densB1=" densB1 "K-effective 1 = " keff1   ;
IF Splitx Splitz * 1 = THEN
  assertS Flux :: 'K-EFFECTIVE' 1 0.9939373 ;
ELSEIF Splitx Splitz * 4 = THEN
  assertS Flux :: 'K-EFFECTIVE' 1 0.9940094 ;
ENDIF ;

IF TypeCal "critical" = THEN

  Flux DELETE: Flux ;

  Fmap Matex Flux := ThermoParam Fmap Matex Cpo Track
   :: <<powi>> <<densB2>> <<pbore_0>> <<temp_comb_0>> <<temp_mode_0>>
      <<dens_mode_0>> ;

  GREP: Flux   ::
    GETVAL 'K-EFFECTIVE ' 1  >>keff2<<   ;
  ECHO "densB2=" densB2 "K-effective    2 = " keff2  ;
  IF Splitx Splitz * 1 = THEN
    assertS Flux :: 'K-EFFECTIVE' 1 1.024981 ;
  ELSEIF Splitx Splitz * 4 = THEN
    assertS Flux :: 'K-EFFECTIVE' 1 1.025050 ;
  ENDIF ;

  EVALUATE Y1 := 1.0 1.0 keff1 / - ;
  EVALUATE Y2 := 1.0 1.0 keff2 / - ;
  L0 := FIND0: :: ITMAX 25 TOL 1.0E-4
                  POINT X <<densB1>> Y <<Y1>>
                  POINT X <<densB2>> Y <<Y2>>
                  >>CONV<< >>ROOT<< ;
  REPEAT

    Flux DELETE: Flux ;

    EVALUATE densB := ROOT ;

    Fmap Matex Flux := ThermoParam Fmap Matex Cpo Track
     :: <<powi>> <<densB>> <<pbore_0>> <<temp_comb_0>> <<temp_mode_0>>
        <<dens_mode_0>> ;

    GREP: Flux   ::
      GETVAL 'K-EFFECTIVE ' 1  >>keffn<<   ;
    ECHO "K-effective iter " Iter "et bore " ROOT " = " keffn   ;

    EVALUATE YNEW  := 1.0 1.0 keffn / - ;
    L0 := FIND0: L0 :: Y <<YNEW>>
                       >>CONV<< >>ROOT<< ;

    EVALUATE Iter := Iter 1 + ;

  UNTIL CONV ;
  ECHO "Zero is =" ROOT ;
  IF Splitx Splitz * 1 = THEN
    EVALUATE REFVALUE := 1.804631E+03 ;
  ELSEIF Splitx Splitz * 4 = THEN
    EVALUATE REFVALUE := 1.807573E+03 ;
  ENDIF ;
  EVALUATE DELTA := ROOT REFVALUE - REFVALUE / ABS ;
  IF DELTA 5.0E-3 < THEN
    ECHO "TEST SUCCESSFUL; DELTA=" DELTA ;
  ELSE
    ECHO "------------" ;
    ECHO "TEST FAILURE" ;
    ECHO "------------" ;
    ECHO "REFERENCE=" REFVALUE " CALCULATED=" ROOT ;
    ABORT: ;
  ENDIF ;

ENDIF ;

ECHO "test rep900_mphy.x2m completed" ;
END: ;
QUIT "LIST" .