1 |       SUBROUTINE genere_kg(cas_genere_kg,ep_genere_kg,numgaz,sumk)
   2 |       IMPLICIT NONE
   3 | c
   4 | c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
   5 | c generation de "kg" pour "cecile.f"
   6 | c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
   7 | c
   8 | c.......................................................................
   9 | c
  10 |       INCLUDE 'cecile.inc'
  11 | 
  12 | c     ngaz_mx est defini dans propradia
  13 | c     et utilise dans radiatif
  14 |       include 'propradia.inc'
  15 |       include 'radiatif.inc'
  16 | c
  17 | c.......................................................................
  18 | c
  19 | c cas : si cas=1 c'est une emission de paroi
  20 | c          si cas=2 c'est une emission de volume
  21 | c ep : epaisseur de la maille d'emission (seulement pour cas=2)
  22 | c
  23 | c numgaz: numero du gaz pour lequel on tire kg
  24 | c sumk: est la somme des k sauf celui de numgaz 
  25 | c       ils ont etaient tires ou calcule grossierement (suie inclus)
  26 | c.......................................................................
  27 | c
  28 |       INTEGER cas_genere_kg, numgaz
  29 |       DOUBLE PRECISION ep_genere_kg
  30 |       DOUBLE PRECISION kgc_genere_kg
  31 |       double precision sumk
  32 | c
  33 | c.......................................................................
  34 | c
  35 |       DOUBLE PRECISION tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8,tmp9
  36 |       INTEGER itmp1,itmp2,itmp3,itmp4,itmp5
  37 | c
  38 | c.......................................................................
  39 | c
  40 | c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  41 | c
  42 | c.......................................................................
  43 | c emission de paroi
  44 | c.......................................................................
  45 | c
  46 | c l'optimisation de l'echange paroi / reste du systeme etait 
  47 | c optimise pour l'emission de la paroi avec une paroi;  fss
  48 | c finallement si on optimise "l'echange" surface/ volume on peut
  49 | c faire comme si on emettait du volume donc
  50 | c cas_genere_kg prend tjrs la valeur 2
  51 | c donc cas_genere_kg =1 doit tjours etre faux
  52 | 
  53 |       IF ((cas_genere_kg.EQ.1).and.(2.eq.1)) THEN
  54 | cc      IF ((cas_genere_kg.EQ.1).and.(2.eq.1)) THEN
  55 | c
  56 | c -----{tirage selon pdf(kg)=fss(kg,0)}
  57 | c
  58 |         CALL rass(phig(numgaz,in),kgbar(numgaz,in),0.D+0,kg(numgaz))
  59 | c
  60 | c.......................................................................
  61 | c emission de volume
  62 | c.......................................................................
  63 | c
  64 |       ELSE
  65 | c
  66 | c -----{calcul de kgc}
  67 | c
  68 | c        kgc_genere_kg=1.D+0/ep_genere_kg - ks(in)
  69 |         kgc_genere_kg=1.D+0/ep_genere_kg - sumk
  70 | c        print *, ' kgc_genere_kg:  ',  kgc_genere_kg
  71 |         IF (kgc_genere_kg.LE.0.D+0) THEN
  72 | c
  73 | c -------{si le coefficient kgc est negatif : "epais"}
  74 | c
  75 |           CALL rass(phig(numgaz,in),kgbar(numgaz,in),0.D+0,kg(numgaz))
  76 |         ELSE
  77 | c
  78 | c -------{sinon calcul de la ponderation "mince" / "epais"}
  79 | c
  80 |           CALL cdss(phig(numgaz,in),kgbar(numgaz,in),0.D+0,
  81 |      &              kgc_genere_kg,tmp3)
  82 | c          print *, 'phig(numgaz,in),kgbar(numgaz,in),0.D+0,
  83 | c     &              kgc_genere_kg,tmp3)'
  84 | c          print *, phig(numgaz,in),kgbar(numgaz,in),0.D+0,
  85 | c     &              kgc_genere_kg,tmp3
  86 |           CALL cdss(phig(numgaz,in),kgbar(numgaz,in),ep_genere_kg,
  87 |      &              kgc_genere_kg,tmp4)
  88 |           CALL trss(phig(numgaz,in),kgbar(numgaz,in),ep_genere_kg,tmp5)
  89 |           tmp8 = dexp(-(sumk)* ep_genere_kg)
  90 |           tmp6=(tmp3 - tmp8*tmp5*tmp4)/(1.D+0 - tmp8*tmp5)
  91 |           
  92 | c          if (in.eq.1) then
  93 | c           if(itir .eq. 5161) then
  94 | c             print *, 'tmp6, tmp8,(1-f)', tmp6,tmp8, (1.D+0 - tmp8*tmp5)
  95 | c             print *, 'tmp3,tmp4', tmp3,tmp4
  96 | c             read *
  97 | c          endif
  98 | 
  99 | ccc debug
 100 | c si tmp5 = 1 toutes les valeurs font 1 = "0/0"
 101 |           if (tmp5 .eq. 1.d+0) then
 102 |              tmp6 = 1.d+0
 103 |           endif
 104 | 
 105 | c
 106 | c -------{test de bernoulli}
 107 | c
 108 |           call rand_uniforme(tmp7)
 109 |           IF (tmp7.LE.tmp6) THEN
 110 | c
 111 | c ---------{mince}
 112 | c
 113 |             CALL rags(phig(numgaz,in),kgbar(numgaz,in),0.D+0,kg(numgaz))
 114 |           ELSE
 115 | c
 116 | c ---------{epais}
 117 | c
 118 |             CALL rass(phig(numgaz,in),kgbar(numgaz,in),0.D+0,kg(numgaz))
 119 |           ENDIF
 120 | c
 121 | c -------{calcul du poids en fonc tion de la proba choisie}
 122 | c
 123 | c          if(itir .eq. 5161) print *, 'GENERE_KG:poids avt', poids
 124 | c           if(itir .eq. 5161) print *, 'GENERE_KG:poids avt', 
 125 | c     & kg(numgaz)/kgbar(numgaz,in), tmp6
 126 |           poids=poids/(tmp6*kg(numgaz)/kgbar(numgaz,in)+(1.D+0-tmp6))
 127 | c          if(itir .eq. 5161) print *, 'GENERE_KG:poids ape', poids
 128 | c
 129 |         ENDIF
 130 | c
 131 |       ENDIF
 132 |       kgskgbar(numgaz)=kg(numgaz)/kgbar(numgaz,in)
 133 | !      PRINT *, 'GENE:',kg(numgaz),kgbar(numgaz,in),kgskgbar(numgaz)
 134 | c
 135 | c.......................................................................
 136 | c
 137 |       RETURN
 138 |       END
 139 | 
 140 | 
 141 | 
 142 | 
 143 | 
 144 | 
 145 | 
 146 | 
 147 | 
 148 | 
 149 | 
 150 | 
 151 | 
 152 |