Initial commit from Polytechnique Montreal

1 files changed, 825 insertions, 0 deletions

diff --git a/Dragon/src/g2s_generatingMC.f90 b/Dragon/src/g2s_generatingMC.f90
new file mode 100644
index 0000000..286fe2e
--- /dev/null
+++ b/Dragon/src/g2s_generatingMC.f90
@@ -0,0 +1,825 @@
+!
+!-----------------------------------------------------------------------
+!
+!Purpose:
+! Generate a dataset for a Monte Carlo code.
+!
+!Copyright:
+! Copyright (C) 2001 Ecole Polytechnique de Montreal.
+! This library is free software; you can redistribute it and/or
+! modify it under the terms of the GNU Lesser General Public
+! License as published by the Free Software Foundation; either
+! version 2.1 of the License, or (at your option) any later version.
+!
+!Author(s):
+! G. Civario (CS-SI)
+!
+!Comments:
+! fonctions:
+!  - initMC : initialisation des tableaux
+!  - destroyMC : liberation des pointeurs
+!  - createElemGene : constructeur d'un element generique
+!  - isSameEG : teste l'equivalence de deux elements generiques
+!  - isSameWayEG : teste si deux elements generiques ont meme orientation
+!  - prepareMCData : preparation des donnes necessaires et remplissage
+!                    des tableaux globaux
+!  - generateTripoliFile : ecriture du fichier de donnees Tripoli4
+!  - generateMCNPFile : ecriture du fichier de donnees MCNP
+!  - generateSerpentFile : ecriture du fichier de donnees Serpent
+!  - putOn80col : formattage sur 80 colonnes style MCNP, d'un buffer a afficher
+!  - findParalleleWithTrans : trouve l'element geometrique a associer avec un
+!                             autre, dans le cas d'une translation
+!
+!-----------------------------------------------------------------------
+!
+module monteCarlo
+    use boundCond
+    use cast
+    use cellulePlaced
+    use constType
+    use constUtiles
+    use segArc
+    use generSAL
+  
+    implicit none
+
+  type t_elemGene
+     logical          :: isPlan !plan -> true      ; cylindre -> false
+     double precision :: x,y    !     -> origine   ;          -> centre
+     double precision :: dx,dy  !     -> extremite ;          -> 0.
+     double precision :: r      !     -> 0.        ;          -> rayon
+     integer          :: limite ! ==Tri_XXXX
+  end type t_elemGene
+
+  type(t_elemGene),dimension(:),allocatable,save :: tabEG
+
+  type t_volume
+     integer                      :: lg      !longueur des tableaux
+     integer                      :: mix     !numero du milieu
+     integer,dimension(:), allocatable :: indElem !indice des elements generiques
+                                             !du contour
+     logical,dimension(:), allocatable :: side    !true -> cote + ; false -> cote -
+     integer,dimension(:), allocatable :: typCL   !condition limite de la face
+  end type t_volume
+
+  type(t_volume),dimension(:),allocatable,save   :: tabVolume
+
+  integer,parameter :: Tri_Not=-1,Tri_Void=0,Tri_Refl=1,Tri_Trans=2,Tri_Cos=3
+
+contains
+
+  subroutine initMC(nbNode)
+    integer,intent(in) :: nbNode
+
+    allocate(tabVolume(nbNode),stat=alloc_ok)
+    if (alloc_ok /= 0) call XABORT("G2S: initMC => allocation pb")
+    tabVolume(1:nbNode)%lg = 0
+  end subroutine initMC
+
+  subroutine destroyMC()
+    integer :: i
+
+    do i = 1,size(tabVolume)
+       deallocate(tabVolume(i)%indElem)
+       deallocate(tabVolume(i)%side)
+       deallocate(tabVolume(i)%typCL)
+    end do
+    deallocate(tabEG,tabVolume)
+  end subroutine destroyMC
+
+  function createElemGene(sa,strDCL)
+    type(t_segArc),intent(in)   :: sa
+    character(len=4),intent(in) :: strDCL
+    type(t_elemGene)            :: createElemGene
+
+    integer :: sunsetCL
+
+    createElemGene%isPlan = (sa%typ==tseg)
+    createElemGene%x = sa%x
+    createElemGene%y = sa%y
+    if (createElemGene%isPlan) then
+       createElemGene%dx = sa%dx
+       createElemGene%dy = sa%dy
+       createElemGene%r  = 0.d0
+    else
+       createElemGene%dx = 0.d0
+       createElemGene%dy = 0.d0
+       createElemGene%r  = sa%r
+    end if
+    createElemGene%limite = Tri_Not
+    sunsetCL = minval((/sa%nodeg,sa%noded/))
+    if (sunsetCL<=0) then
+       sunsetCL = mod(-sunsetCL,100)
+       select case(sunsetCL)
+       case(B_Void,B_Zero)
+          createElemGene%limite = Tri_Void
+       case(B_Refl,B_Ssym,B_Syme,B_Diag)
+          createElemGene%limite = Tri_Refl
+       case(B_Tran)
+          createElemGene%limite = Tri_Trans
+       case(B_Albe)
+          createElemGene%limite = Tri_Cos
+       case(-fooMix,0)
+          !il faut utiliser la CL par defaut
+          if (strDCL=='ALBE') then
+             createElemGene%limite = Tri_Cos
+          else if (strDCL=='REFL') then
+             createElemGene%limite = Tri_Refl
+          else ! => VOID
+             createElemGene%limite = Tri_Void
+          end if
+       case default
+          call XABORT("G2MC: boundary condition not allowed")
+       end select
+    end if
+  end function createElemGene
+
+  function isSameEG(eg1,eg2)
+    type(t_elemGene),intent(in) :: eg1,eg2
+    logical                     :: isSameEG
+
+    if (eg1%isPlan.neqv.eg2%isPlan) then
+       isSameEG = .false.
+    else if (eg1%isPlan) then
+       isSameEG = pointsAlignes(eg1%x,eg1%y,eg1%dx,eg1%dy,eg2%x,eg2%y)  &
+            .and. pointsAlignes(eg1%x,eg1%y,eg1%dx,eg1%dy,eg2%dx,eg2%dy)
+    else
+       isSameEG = isEqualConst(eg1%x,eg2%x) &
+            .and. isEqualConst(eg1%y,eg2%y) &
+            .and. isEqualConst(eg1%r,eg2%r)
+    end if
+    isSameEG = isSameEG .and. (eg1%limite==eg2%limite)
+  end function isSameEG
+
+  function isSameWayEG(eg1,eg2)
+    type(t_elemGene),intent(in) :: eg1,eg2
+    logical                     :: isSameWayEG
+
+    if (eg1%isPlan .and. eg2%isPlan) then
+       isSameWayEG = &
+            (((eg1%dx-eg1%x)*(eg2%dx-eg2%x)+(eg1%dy-eg1%y)*(eg2%dy-eg2%y))>0)
+    else
+       isSameWayEG = .true.
+    end if
+  end function isSameWayEG
+
+  subroutine prepareMCData(szSA,nbNode,szEG)
+    integer,intent(in)    :: szSA,nbNode
+    integer,intent(inout) :: szEG
+
+    integer          :: i,j,volNbr,defautCl,lgMax,indFictive
+    logical          :: found,toPut
+    real             :: albedo
+    character*4      :: strDCL
+    double precision :: ptx,pty,rxx,ryy
+    type(t_elemGene) :: eg
+    integer,dimension(:),allocatable :: indEG,nbPoints
+    logical,dimension(:),allocatable :: isGoodWay,withFictiveEG,inACircle
+    double precision,dimension(:,:),allocatable :: xx,yy
+    type(t_elemGene),dimension(:),allocatable   :: tmpTabEG
+
+    !recuperation des donnees de conditions aux limites
+    call calculDefaultCl(defautCl,albedo,strDCL)
+    !preparation des donnees de remplissage des tableaux globaux
+    allocate(indEG(szSA),isGoodWay(szSA),tmpTabEG(szSA+nbNode),stat=alloc_ok)
+    if (alloc_ok /= 0) call XABORT("G2S: prepareMCData(1) => allocation pb")
+    do i = 1,szSA
+       eg = createElemGene(tabSegArc(i),strDCL)
+       found = .false.
+       do j = 1,szEG
+          found = isSameEG(eg,tmpTabEG(j))
+          if (found) exit
+       end do
+       if (found) then
+          indEG(i) = j
+          isGoodWay(i) = isSameWayEG(eg,tmpTabEG(j))
+       else
+          szEG = szEG + 1
+          tmpTabEG(szEG) = eg
+          indEG(i) = szEG
+          isGoodWay(i) = .true.
+       end if
+    end do
+    !preparation des donnes permettant de savoir si un element fictif
+    !englobant est necessaire
+    allocate(inACircle(nbNode),withFictiveEG(nbNode),stat=alloc_ok)
+    if (alloc_ok /= 0) call XABORT("G2S: prepareMCData(2) => allocation pb")
+    withFictiveEG(:nbNode) = .false. ; inACircle(:nbNode) = .false.
+    !calcul des dimensions des pointeurs du tableau des volumes
+    do i = 1,szSA
+       !cote gauche
+       volNbr = tabSegArc(i)%nodeg
+       if (volNbr>0) then
+          tabVolume(volNbr)%lg = tabVolume(volNbr)%lg + 1
+          if (tabSegArc(i)%typ==tcer .and. .not. inACircle(volNbr)) then
+             !le node est a l'interieur d'un cercle => on n'a pas besoin de
+             !l'englober
+             inACircle(volNbr) = .true.
+             if (withFictiveEG(volNbr)) then
+                !on avait commencer a englober fictivement => on laisse tomber
+                tabVolume(volNbr)%lg = tabVolume(volNbr)%lg - 1
+                withFictiveEG(volNbr) = .false.
+             end if
+          end if
+       end if
+       !cote droit
+       volNbr = tabSegArc(i)%noded
+       if (volNbr>0) then
+          tabVolume(volNbr)%lg = tabVolume(volNbr)%lg + 1
+          if (tabSegArc(i)%typ==tarc .and. .not. withFictiveEG(volNbr) &
+              .and. .not. inACircle(volNbr)) then
+             !le node est a l'exterieur d'un arc => on englobe si pas encore
+             !fait
+             tabVolume(volNbr)%lg = tabVolume(volNbr)%lg + 1
+             withFictiveEG(volNbr) = .true.
+          end if
+       end if
+    end do
+    deallocate(inACircle)
+    lgMax = maxval(tabVolume(:nbNode)%lg)
+    !allocation des pointeurs du tableau des volumes
+    do i = 1,nbNode
+       allocate(tabVolume(i)%indElem(tabVolume(i)%lg),stat=alloc_ok)
+    if (alloc_ok /= 0) call XABORT("G2S: prepareMCData(3) => allocation pb")
+       tabVolume(i)%indElem = 0
+       allocate(tabVolume(i)%side(tabVolume(i)%lg),stat=alloc_ok)
+    if (alloc_ok /= 0) call XABORT("G2S: prepareMCData(4) => allocation pb")
+       allocate(tabVolume(i)%typCL(tabVolume(i)%lg),stat=alloc_ok)
+    if (alloc_ok /= 0) call XABORT("G2S: prepareMCData(5) => allocation pb")
+       tabVolume(i)%typCL = 0
+       !remise a zero du compteur, pour gerer les doublons
+       tabVolume(i)%lg = 0
+    end do
+    !remplissage des pointeurs, avec elimination des doublons
+    allocate(xx(nbNode,lgMax),yy(nbNode,lgMax))
+    allocate(nbPoints(nbNode),stat=alloc_ok)
+    if (alloc_ok /= 0) call XABORT("G2S: prepareMCData(6) => allocation pb")
+    nbPoints(:nbNode) = 0 ; xx(:nbNode,:lgMax) = 0.d0 ; yy(:nbNode,:lgMax) = 0.d0
+    do i = 1,szSA
+       !cote gauche
+       volNbr = tabSegArc(i)%nodeg
+       if (volNbr>0) then
+          toPut = ( count(tabVolume(volNbr)%indElem(:)==indEG(i)) == 0 )
+          if (toPut) then
+             tabVolume(volNbr)%mix = tabSegArc(i)%neutronicMixg
+             tabVolume(volNbr)%lg = tabVolume(volNbr)%lg + 1
+             tabVolume(volNbr)%indElem(tabVolume(volNbr)%lg) = indEG(i)
+             tabVolume(volNbr)%typCL(tabVolume(volNbr)%lg) = &
+                  tmpTabEG(indEG(i))%limite
+             tabVolume(volNbr)%side(tabVolume(volNbr)%lg) = .not.isGoodWay(i)
+          end if
+          !preparation du calcul des coordonnees du cercle englobant
+          if (withFictiveEG(volNbr).and.tabSegArc(i)%typ/=tcer) then
+             call giveOrigine(tabSegArc(i),ptx,pty)
+             nbPoints(volNbr) = nbPoints(volNbr) + 1
+             xx(volNbr,nbPoints(volNbr)) = ptx
+             yy(volNbr,nbPoints(volNbr)) = pty
+          end if
+       end if
+       !cote droit
+       volNbr = tabSegArc(i)%noded
+       if (volNbr>0) then
+          toPut = ( count(tabVolume(volNbr)%indElem(:)==indEG(i)) == 0 )
+          if (toPut) then
+             tabVolume(volNbr)%mix = tabSegArc(i)%neutronicMixd
+             tabVolume(volNbr)%lg = tabVolume(volNbr)%lg + 1
+             tabVolume(volNbr)%indElem(tabVolume(volNbr)%lg) = indEG(i)
+             tabVolume(volNbr)%typCL(tabVolume(volNbr)%lg) = &
+                  tmpTabEG(indEG(i))%limite
+             tabVolume(volNbr)%side(tabVolume(volNbr)%lg) = isGoodWay(i)
+          end if
+          !preparation du calcul des coordonnees du cercle englobant
+          if (withFictiveEG(volNbr).and.tabSegArc(i)%typ/=tcer) then
+             call giveExtremite(tabSegArc(i),ptx,pty)
+             nbPoints(volNbr) = nbPoints(volNbr) + 1
+             xx(volNbr,nbPoints(volNbr)) = ptx
+             yy(volNbr,nbPoints(volNbr)) = pty
+          end if
+       end if
+    end do
+    !calcul des cercles englobants
+    do i = 1,nbNode
+       if (.not.withFictiveEG(i)) cycle
+       !creation du cercle
+       eg%isPlan = .false.
+       eg%x = sum(xx(i,:nbPoints(i)))/nbPoints(i)
+       eg%y = sum(yy(i,:nbPoints(i)))/nbPoints(i)
+       eg%dx = 0.d0 ; eg%dy = 0.d0 ; eg%r = 0.d0 ; eg%limite = Tri_Not
+       do j = 1,nbPoints(i)
+          rxx = xx(i,j)-eg%x ; ryy = yy(i,j)-eg%y
+          eg%r = max(eg%r,longVect(rxx,ryy))
+       end do
+       !integration eventuelle dans le tableau des elements geometriques
+       found = .false.
+       do j = 1,szEG
+          found = isSameEG(eg,tmpTabEG(j))
+          if (found) exit
+       end do
+       if (found) then
+          indFictive = j
+       else
+          szEG = szEG + 1
+          tmpTabEG(szEG) = eg
+          indFictive = szEG
+       end if
+       !ajout eventuel dans la liste des elements geometriques constitutifs
+       !de la zone
+       toPut = ( count(tabVolume(i)%indElem(:)==indFictive) == 0 )
+       if (toPut) then
+          tabVolume(i)%lg = tabVolume(i)%lg + 1
+          tabVolume(i)%indElem(tabVolume(i)%lg) = indFictive
+          tabVolume(i)%typCL(tabVolume(i)%lg) = Tri_Not
+          tabVolume(i)%side(tabVolume(i)%lg) = .false.
+       end if
+    end do
+    deallocate(nbPoints,xx,yy,indEG,isGoodWay,withFictiveEG)
+    !recopie du tableau global des elements geometriques, et liberation
+    !du temporaire
+    allocate(tabEG(szEG),stat=alloc_ok)
+    if (alloc_ok /= 0) call XABORT("G2S: prepareMCData(7) => allocation pb")
+    tabEG(1:szEG) = tmpTabEG(1:szEG)
+    deallocate(tmpTabEG)
+  end subroutine prepareMCData
+
+  subroutine generateTripoliFile(fileTripoli,szSA,nbNode)
+    integer,intent(in) :: fileTripoli,szSA,nbNode
+
+    real    :: a,b,c,d,xx,yy,rr
+    integer :: i,j,szEG,nbPlus,nbMoins,nbMix,nbCLTot,nbCL
+    integer,dimension(:),allocatable :: sp,sm,listMix,volmix,clList
+
+    szEG = 0
+    call initMC(nbNode)
+    call prepareMCData(szSA,nbNode,szEG)
+
+    write(fileTripoli,'(/"GEOMETRY"/)')
+    write(fileTripoli,'("//jeu de donnees geometriques TRIPOLI")')
+    write(fileTripoli,'("//genere pour DRAGON par le module G2MC"/)')
+    write(fileTripoli,'("TITRE geometrie DRAGON pour TRIPOLI"/)')
+    !definition des surfaces
+    write(fileTripoli,'("//definition des surfaces"/)')
+    c = 0.
+    do i = 1,szEG
+       if (tabEG(i)%isPlan) then
+          a = real(tabEG(i)%dy - tabEG(i)%y)
+          b = real(tabEG(i)%x - tabEG(i)%dx)
+          d = real(tabEG(i)%y*tabEG(i)%dx - tabEG(i)%x*tabEG(i)%dy)
+          write(fileTripoli,*) " SURF ",i," PLAN ",a,b,c,d
+       else
+          xx = real(tabEG(i)%x)
+          yy = real(tabEG(i)%y)
+          rr = real(tabEG(i)%r)
+          write(fileTripoli,*) " SURF ",i," CYLZ ",xx,yy,rr
+       end if
+    end do
+    write(fileTripoli,*) " SURF ",szEG+1," PLANZ  0. //plan inferieur" 
+    write(fileTripoli,*) " SURF ",szEG+2," PLANZ 30. //plan superieur"
+    !definition des volumes
+    write(fileTripoli,'(/"//definition des volumes"/)')
+    do i = 1,nbNode
+       nbPlus = count(tabVolume(i)%side(1:tabVolume(i)%lg))
+       nbMoins = tabVolume(i)%lg - nbPlus
+       allocate(sp(nbPlus),sm(nbMoins),stat=alloc_ok)
+       if (alloc_ok /= 0) call XABORT("G2S: generateTripoliFile(1) => allocation pb")
+       sp = pack(tabVolume(i)%indElem(1:tabVolume(i)%lg), &
+                 tabVolume(i)%side(1:tabVolume(i)%lg))
+       sm = pack(tabVolume(i)%indElem(1:tabVolume(i)%lg), &
+                 .not.tabVolume(i)%side(1:tabVolume(i)%lg))
+       write(fileTripoli,*) " VOLU"
+       write(fileTripoli,*) "  ",i
+       write(fileTripoli,*) "   EQUA"
+       write(fileTripoli,*) "     PLUS",nbPlus+1
+       write(fileTripoli,*) "      ",szEG+1,sp(:) !plan inf en premier
+       write(fileTripoli,*) "     MOINS",nbMoins+1
+       write(fileTripoli,*) "      ",sm(:),szEG+2 !plan sup en dernier
+       write(fileTripoli,*) " FINV"
+       write(fileTripoli,*)
+       deallocate(sp,sm)
+    end do
+    write(fileTripoli,'("FINGEOM")')
+    
+    !conditions aux limites
+    write(fileTripoli,'(/"//conditions aux limites")')
+    write(fileTripoli,'("LIMIT")')
+    !calcul du nombre total de lignes de conditions limites particulieres
+    nbCLTot = 2*nbNode
+    do i = 1,nbNode
+       nbCLTot = nbCLTot &
+            + count(tabVolume(i)%typCL(1:tabVolume(i)%lg)==Tri_Refl) &
+            + count(tabVolume(i)%typCL(1:tabVolume(i)%lg)==Tri_Trans) &
+            + count(tabVolume(i)%typCL(1:tabVolume(i)%lg)==Tri_Cos)
+    end do
+    write(fileTripoli,*) nbCLTot
+    !affichage des conditions limites
+    do i = 1,nbNode
+       !Reflexions
+       write(fileTripoli,*) "  ",i," REFLECTION ",szEG+1
+       write(fileTripoli,*) "  ",i," REFLECTION ",szEG+2
+       nbCL = count(tabVolume(i)%typCL(1:tabVolume(i)%lg)==Tri_Refl)
+       if (nbCL/=0) then
+          allocate(clList(nbCL),stat=alloc_ok)
+          if (alloc_ok /= 0) call XABORT("G2S: generateTripoliFile(2) => allocation pb")
+          clList = pack(tabVolume(i)%indElem(1:tabVolume(i)%lg), &
+                        tabVolume(i)%typCL(1:tabVolume(i)%lg)==Tri_Refl)
+          do j = 1,nbCL
+             write(fileTripoli,*) "  ",i," REFLECTION ",clList(j)
+          end do
+          deallocate(clList)
+       end if
+       !Translations
+       nbCL = count(tabVolume(i)%typCL(1:tabVolume(i)%lg)==Tri_Trans)
+       if (nbCL/=0) then
+          allocate(clList(nbCL),stat=alloc_ok)
+          if (alloc_ok /= 0) call XABORT("G2S: generateTripoliFile(3) => allocation pb")
+          clList = pack(tabVolume(i)%indElem(1:tabVolume(i)%lg), &
+                        tabVolume(i)%typCL(1:tabVolume(i)%lg)==Tri_Trans)
+          do j = 1,nbCL
+             write(fileTripoli,*) "  ",i," TRANSLATION ",clList(j)
+          end do
+          deallocate(clList)
+       end if
+       !Cosinus
+       nbCL = count(tabVolume(i)%typCL(1:tabVolume(i)%lg)==Tri_Cos)
+       if (nbCL/=0) then
+          allocate(clList(nbCL),stat=alloc_ok)
+          if (alloc_ok /= 0) call XABORT("G2S: generateTripoliFile(4) => allocation pb")
+          clList = pack(tabVolume(i)%indElem(1:tabVolume(i)%lg), &
+                        tabVolume(i)%typCL(1:tabVolume(i)%lg)==Tri_Cos)
+          do j = 1,nbCL
+             write(fileTripoli,*) "  ",i," COSINUS ",clList(j)
+          end do
+          deallocate(clList)
+       end if
+    end do
+    write(fileTripoli,'("FIN_LIMIT"/)')
+
+    ! insertion des donnees neutroniques des materiaux
+    ! le milieu ne correspond au materiaux nomme "MIX_n"
+    write(fileTripoli,'("//insertion des donnees neutroniques")')
+    write(fileTripoli,'("//des matreriaux")')
+    write(fileTripoli,'("FILE")')
+    write(fileTripoli,'("   mix_data"/)') !le nom est fixe, mais pourra etre
+                                          !change si besoin  ######
+
+    !definition du milieu par volume
+    write(fileTripoli,'("//definition du milieu par volume")')
+    write(fileTripoli,'("GEOMCOMP")')
+    allocate(listMix(nbNode),stat=alloc_ok)
+    if (alloc_ok /= 0) call XABORT("G2S: generateTripoliFile(5) => allocation pb")
+    listMix(:nbNode) = 0
+    nbMix = 0
+    do i = 1,nbNode
+       if (count(listMix(:nbNode)==tabVolume(i)%mix)==0) then
+          nbMix = nbMix + 1
+          listMix(nbMix) = tabVolume(i)%mix
+       end if
+    end do
+    do i = 1,nbMix
+       allocate(volMix(count(tabVolume(:)%mix==listMix(i))),stat=alloc_ok)
+       if (alloc_ok /= 0) call XABORT("G2S: generateTripoliFile(6) => allocation pb")
+       volMix = pack((/(j,j=1,nbNode)/),tabVolume(:)%mix==listMix(i))
+       write(fileTripoli,*) " MIX_"//trim(i2s(listMix(i))),size(volMix)
+       write(fileTripoli,*) "   ",volMix
+       deallocate(volMix)
+    end do
+    deallocate(listMix)
+    write(fileTripoli,'("FIN_GEOMCOMP"/)')
+
+    call destroyMC()
+  end subroutine generateTripoliFile
+
+  subroutine generateMCNPFile(fileMCNP,szSA,nbNode)
+    integer,intent(in) :: fileMCNP,szSA,nbNode
+
+    integer :: i,j,k,szEG,ind,transInd,nbEGCL
+    real    :: a,b,c,d,xx,yy,rr
+    character(len=1)    :: sgn
+    character(len=5)    :: nbr,mix
+    character(len=6)    :: other,signedNbr
+    character(len=12)   :: text12
+    character(len=1000) :: buffer
+    integer,dimension(:),allocatable :: tabEGCL
+
+    szEG = 0
+    call initMC(nbNode)
+    call prepareMCData(szSA,nbNode,szEG)
+    
+    !title card
+    write(fileMCNP,'("Jeu de donnees geometriques MCNP genere pour DRAGON &
+         &par le module G2MC")')
+    !message en commentaire
+    write(fileMCNP,'("C Attention ! Les densites des milieux n''etant pas &
+         &disponibles dans les")')
+    write(fileMCNP,'("C jeux de donnees geometriques DRAGON, il faut passer &
+         &une moulinette pour")')
+    write(fileMCNP,'("C remplacer les ""DENSITE_I"" presents dans ce fichier, &
+         &par leur valeur."/)')
+    !cell card
+    write(fileMCNP,'("C   cell card")')
+    do i = 1,nbNode
+       text12 = i2s(i)
+       nbr = text12(:5)
+       text12 = i2s(tabVolume(i)%mix)
+       mix = text12(:5)
+       buffer =nbr//" "//trim(mix)//" DENSITE_"//trim(mix)
+       ind = len_trim(buffer)+1
+       write(buffer(ind:),*) merge(tabVolume(i)%indElem(1:tabVolume(i)%lg), &
+            -tabVolume(i)%indElem(1:tabVolume(i)%lg), &
+            tabVolume(i)%side(1:tabVolume(i)%lg)) , &
+            szEG+1 , -(szEG+2)
+       call putOn80col(fileMCNP,buffer)
+    end do
+    text12=i2s(nbNode+1)
+    nbr = text12(:5)
+    buffer = nbr//" 0"
+    ind = 9
+    if (geomTyp==HexTyp.or.geomTyp==TriaTyp) then
+       !on a des risques de geometrie non convexe
+       ! => on prend les complementaires
+       do i = 1,nbNode
+          text12 = i2s(i)
+          nbr = text12(:5)
+          buffer(ind:) = "#"//nbr
+          ind = ind + len_trim(nbr) + 2
+       end do
+    else
+       !la geometrie est convexe
+       ! on prend l'union des complementaires des cotes a condition limite
+       nbEGCL = count(tabEG(:)%limite/=Tri_Not)
+       allocate(tabEGCL(nbEGCL),stat=alloc_ok)
+       if (alloc_ok /= 0) call XABORT("G2S: generateTripoliFile(7) => allocation pb")
+       tabEGCL = pack((/(i,i=1,szEG)/),tabEG(:)%limite/=Tri_Not)
+       do i = 1,nbEGCL
+          do j = 1,nbNode
+             do k = 1,tabVolume(j)%lg
+                if (tabVolume(j)%indElem(k)==tabEGCL(i)) then
+                   if (tabVolume(j)%side(k)) tabEGCL(i) = -tabEGCL(i)
+                   exit
+                end if
+             end do
+          end do
+       end do
+       do i = 1,nbEGCL
+          text12 = i2s(tabEGCL(i))
+          signedNbr = text12(:6)
+          buffer(ind:) = trim(signedNbr)//" :"
+          ind = ind + len_trim(signedNbr) + 3
+       end do
+       buffer(ind:) = trim(i2s(-(szEG+1)))//" : "//i2s(szEG+2)
+       deallocate(tabEGCL)
+    end if
+    call putOn80col(fileMCNP,buffer)
+
+    !surface card
+    write(fileMCNP,'(/"C   surface card")')
+    c = 0.
+    do i = 1,szEG
+       text12 = i2s(i)
+       nbr = text12(:5)
+       sgn = ' '
+       other = ' '
+       select case(tabEG(i)%limite)
+       case(Tri_Refl)
+          sgn = '*'
+       case(Tri_Trans)
+          transInd = findParalleleWithTrans(i,szEG)
+          text12=i2s(transInd)
+          if (isSameWayEG(tabEG(i),tabEG(transInd))) then
+             other = '-'//text12(:6)
+          else
+             other = text12(:6)
+          end if
+       case(Tri_Cos)
+          sgn = '+'
+       end select
+       if (tabEG(i)%isPlan) then
+          a = real(tabEG(i)%dy - tabEG(i)%y)
+          b = real(tabEG(i)%x - tabEG(i)%dx)
+          d = real(tabEG(i)%y*tabEG(i)%dx - tabEG(i)%x*tabEG(i)%dy)
+          buffer = adjustl(sgn//nbr//trim(" "//other)//" P")
+          ind = len_trim(buffer) + 1
+          write(buffer(ind:),*) a,b,c,d
+       else
+          xx = real(tabEG(i)%x)
+          yy = real(tabEG(i)%y)
+          rr = real(tabEG(i)%r)
+          buffer = adjustl(sgn//nbr//" C/Z")
+          ind = len_trim(buffer) + 1
+          write(buffer(ind:),*) xx,yy,rr
+       end if
+       call putOn80col(fileMCNP,buffer)
+    end do
+    text12 = i2s(szEG+1)
+    nbr = text12(:5)
+    buffer = "*"//nbr//" PZ  0. $plan inferieur"
+    call putOn80col(fileMCNP,buffer)
+    text12 = i2s(szEG+2)
+    nbr = text12(:5)
+    buffer = "*"//nbr//" PZ 30. $plan superieur"
+    call putOn80col(fileMCNP,buffer)
+
+    !data card
+
+    call destroyMC()
+  end subroutine generateMCNPFile
+
+  subroutine generateSerpentFile(fileSerpent,szSA,nbNode)
+    integer,intent(in) :: fileSerpent,szSA,nbNode
+
+    integer :: i,j,k,szEG,ind,transInd,nbEGCL,nb2,isurf,jjj 
+    real    :: a,b,c,d,xx,yy,rr,cuboid(6)
+    character(len=1)    :: sgn
+    character(len=5)    :: nbr,mix
+    character(len=6)    :: other
+    character(len=12)   :: text12
+    character(len=1000) :: buffer
+    integer,dimension(:),allocatable :: tabEGCL,tmpnod,tmpnod2
+    logical :: lcuboid
+
+    szEG = 0
+    call initMC(nbNode)
+    call prepareMCData(szSA,nbNode,szEG)
+
+    !title card
+    write(fileSerpent,'("% Serpent combinatorial geometry generated by module G2MC:")')
+    !cell card
+    nbEGCL = count(tabEG(:)%limite/=Tri_Not)
+    allocate(tabEGCL(nbEGCL),stat=alloc_ok)
+    if (alloc_ok /= 0) call XABORT("G2S: generateSerpentFile(1) => allocation pb")
+    tabEGCL = pack((/(i,i=1,szEG)/),tabEG(:)%limite/=Tri_Not)
+    do i = 1,nbEGCL
+       do j = 1,nbNode
+          do k = 1,tabVolume(j)%lg
+             if (tabVolume(j)%indElem(k)==tabEGCL(i)) then
+                if (tabVolume(j)%side(k)) tabEGCL(i) = -tabEGCL(i)
+                exit
+             end if
+          end do
+       end do
+    end do
+    write(fileSerpent,'("%   cell cards")')
+    nb2=0
+    cuboid=(/ 0.0, 0.0, 0.0, 0.0, -1.0E10, 1.0E10 /)
+    do i = 1,nbNode
+       allocate(tmpnod(tabVolume(i)%lg),tmpnod2(tabVolume(i)%lg),stat=alloc_ok)
+       if (alloc_ok /= 0) call XABORT("G2S: generateSerpentFile(2) => allocation pb")
+       text12 = i2s(i)
+       nbr = text12(:5)
+       text12 = i2s(tabVolume(i)%mix)
+       mix = text12(:5)
+       buffer ="cell "//nbr//" 0 MIX_"//trim(mix)
+       ind = len_trim(buffer)+1
+       tmpnod(:tabVolume(i)%lg)=merge(-tabVolume(i)%indElem(1:tabVolume(i)%lg), &
+            tabVolume(i)%indElem(1:tabVolume(i)%lg), &
+            tabVolume(i)%side(1:tabVolume(i)%lg))
+       jjj=0
+       lcuboid=.true.
+       do j = 1,tabVolume(i)%lg
+         do k =1,nbEGCL
+           if(tmpnod(j) == tabEGCL(k)) then
+             isurf=abs(tabEGCL(k))
+             if (.not.tabEG(isurf)%isPlan) call XABORT('g2s_generatingMC: unsupported BC')
+             a = real(tabEG(isurf)%dy - tabEG(isurf)%y)
+             b = real(tabEG(isurf)%x - tabEG(isurf)%dx)
+             d = real(tabEG(isurf)%y*tabEG(isurf)%dx - tabEG(isurf)%x*tabEG(isurf)%dy)
+             if ((a == 0.0).and.(tabEGCL(k) < 0)) then
+               cuboid(4)=-d/b
+             else if ((a == 0.0).and.(tabEGCL(k) > 0)) then
+               cuboid(3)=-d/b
+             else if ((b == 0.0).and.(tabEGCL(k) < 0)) then
+               cuboid(1)=-d/a
+             else if ((b == 0.0).and.(tabEGCL(k) > 0)) then
+               cuboid(2)=-d/a
+             endif
+             if(lcuboid) then
+               jjj=jjj+1
+               tmpnod2(jjj)=-szEG-1
+               lcuboid=.false.
+             endif
+             go to 10
+           endif
+         end do
+         jjj=jjj+1
+         tmpnod2(jjj)=-tmpnod(j)
+    10   continue
+       end do
+       if(jjj.gt.0) then
+         write(buffer(ind:),*) tmpnod2(:jjj)
+         write(fileSerpent,'(a)') trim(buffer)
+       endif
+       deallocate(tmpnod2,tmpnod)
+    end do
+    if (nbEGCL > 0) then
+       text12 = i2s(nbNode+1)
+       nbr = text12(:5)
+       text12 = i2s(szEG+1)
+       buffer ="cell "//nbr//" 0 outside         "//trim(text12)
+       write(fileSerpent,'(a)') trim(buffer)
+    endif
+
+    !surface card
+    write(fileSerpent,'(/"%   surface cards")')
+    c = 0.
+    do i = 1,szEG
+       do k = 1,nbEGCL
+         if (i == abs(tabEGCL(k))) go to 20
+       end do
+       text12 = i2s(i)
+       nbr = text12(:5)
+       sgn = ' '
+       other = ' '
+       select case(tabEG(i)%limite)
+       case(Tri_Refl)
+          sgn = '*'
+       case(Tri_Trans)
+          transInd = findParalleleWithTrans(i,szEG)
+          text12=i2s(transInd)
+          if (isSameWayEG(tabEG(i),tabEG(transInd))) then
+             other = '-'//text12(:6)
+          else
+             other = text12(:6)
+          end if
+       case(Tri_Cos)
+          sgn = '+'
+       end select
+       if (tabEG(i)%isPlan) then
+         a = real(tabEG(i)%dy - tabEG(i)%y)
+         b = real(tabEG(i)%x - tabEG(i)%dx)
+         d = -real(tabEG(i)%y*tabEG(i)%dx - tabEG(i)%x*tabEG(i)%dy)
+         buffer = adjustl("surf "//sgn//nbr//trim(" "//other)//" plane")
+         ind = len_trim(buffer) + 1
+         write(buffer(ind:),*) a,b,c,d
+       else
+         xx = real(tabEG(i)%x)
+         yy = real(tabEG(i)%y)
+         rr = real(tabEG(i)%r)
+         buffer = adjustl("surf "//sgn//nbr//" cyl")
+         ind = len_trim(buffer) + 1
+         write(buffer(ind:),*) xx,yy,rr
+       end if
+       write(fileSerpent,'(a)') trim(buffer)
+   20  continue
+    end do
+    if (nbEGCL > 0) then
+       text12 = i2s(szEG+1)
+       nbr = text12(:5)
+       buffer = adjustl("surf "//sgn//nbr//" cuboid ")
+       ind = len_trim(buffer) + 1
+       write(buffer(ind:),*) (cuboid(i),i=1,6)
+       write(fileSerpent,'(a)') trim(buffer)
+    endif
+    deallocate(tabEGCL)
+    !data card
+
+    call destroyMC()
+  end subroutine generateSerpentFile
+
+  subroutine putOn80col(outFile,strIn)
+    integer,intent(in)          :: outFile
+    character(len=*),intent(in) :: strIn
+
+    integer           :: lg,lastBlank
+    character(len=15) :: frm
+    character(len=80) :: truncStr
+    character(len=len(strIn)) :: str
+
+    str = strIn
+    do
+       lg = len_trim(str)
+       if (lg <= 80) then
+          frm = "(a"//i2s(lg)//")"
+          write(outFile,frm) str
+          return
+       end if
+       truncStr = str(1:78)
+       do lastBlank = 78,1,-1
+          if (truncStr(lastBlank:lastBlank)==' ') exit
+       end do
+       truncStr = trim(str(1:lastBlank-1))//' &'
+       frm = "(a"//i2s(len_trim(truncStr))//")"
+       write(outFile,frm) truncStr
+       str = '     '//adjustl(str(lastBlank:))
+    end do
+  end subroutine putOn80col
+
+  function findParalleleWithTrans(ind,szEG)
+    integer,intent(in) :: ind,szEG
+    integer            :: findParalleleWithTrans
+
+    integer :: i
+
+    findParalleleWithTrans = 0
+    do i = 1,szEG
+       if (i==ind) cycle
+       if (.not.tabEG(i)%isPlan) cycle
+       if (tabEG(i)%limite/=Tri_Trans) cycle
+       if (.not.estColi(tabEG(i)%dx-tabEG(i)%x,tabEG(i)%dy-tabEG(i)%y, &
+            tabEG(ind)%dx-tabEG(ind)%x,tabEG(ind)%dy-tabEG(ind)%y)) cycle
+       findParalleleWithTrans = i
+       return
+    end do
+    call XABORT("G2MC: error, no parallel side found for translation &
+         &boundary condition")
+  end function findParalleleWithTrans
+
+end module monteCarlo