1 | DO igaz=1,ngaz
2 | !! IF ((iin.EQ.0).OR.(iin.EQ.n+1)) THEN
3 | !! tmp4=kbarphi_s(igaz)
4 | !! tmp5=kbar_s(igaz)
5 | !! ELSE
6 | !! tmp4=kbarphi_s(igaz)+
7 | !! & kgbar(igaz,iin)*phig(igaz,iin)*sigma_2
8 | !! tmp5=kbar_s(igaz)+kgbar(igaz,iin)*sigma_2
9 | !! ENDIF
10 | ! print *, 'epoc', igaz, epocg(igaz), tmp5,(tmp5.GT.1.D-10)
11 | ! print *, 'phic', igaz, phicg(igaz), tmp4/tmp5
12 | ! read *
13 |
14 | IF (epocg(igaz).GT.1.D-10) THEN
15 | c-[phi de curtis godson entre les deux points]
16 | ! phicg(igaz) = tmp4 / tmp5
17 | c-[correction du phi initial pour passer à phi curtis godson]
18 | call f_ss(phicg(igaz),1.D+0,0.D+0,kgskgbar(igaz),tmp8)
19 | call f_ss(phig(igaz,in),1.D+0,0.D+0,kgskgbar(igaz),tmp9)
20 | tmp6 = tmp8/tmp9
21 | cor(igaz) = tmp6
22 |
23 | c-[correction curtis godson pour les derivees]
24 | test1=0
25 | test2=0
26 | !igaz=3
27 | IF ((in.EQ.0).OR.(in.EQ.n+1)) test1=1
28 | IF ((iin.EQ.0).OR.(iin.EQ.n+1)) test2=1
29 | IF ((test1.EQ.1).AND.(test2.EQ.1)) THEN !paroi-paroi
30 | correction = 1.d+0
31 | ELSEIF (test1.EQ.1) THEN !paroi-gaz
32 | te = epocg(igaz)
33 | tphi = phicg(igaz)
34 | toto = 1.d+0 + (2*epocg(igaz)/tphi)
35 | tmp7=kgbar(igaz,iin)/epocg(igaz)*(phig(igaz,iin)-phicg(igaz))
36 |
37 | tGa_j(igaz)= - toto**(-0.5) *
38 | & (kgbar(igaz,iin)+tmp7*(1.d+0 + (te/tphi)- toto**0.5))
39 | tGb_j(igaz)= - toto**(-0.5) *
40 | & (kgbar(igaz,iin))
41 |
42 | !! tmp6= tGa_j/ tGb_j
43 | c.rfo.
44 | c.lat. rajoute condition
45 | !! if (hetero.eq.2) correction=tmp6
46 | c.rfo.
47 | ELSEIF (test2.EQ.1) THEN !gaz-paroi
48 | te = epocg(igaz)
49 | tphi = phicg(igaz)
50 | toto = 1.d+0 + (2*epocg(igaz)/tphi)
51 | tmp7=kgbar(igaz,in)/epocg(igaz)*(phig(igaz,in)-phicg(igaz))
52 |
53 | tGa_i(igaz)= - toto**(-0.5) *
54 | & (kgbar(igaz,in)+tmp7*(1.d+0 + (te/tphi)- toto**0.5))
55 | tGb_i(igaz)= - toto**(-0.5) *
56 | & (kgbar(igaz,in))
57 |
58 | !! tmp6= tGa_i/ tGb_i
59 | c.rfo.
60 | c.lat. rajoute condition
61 | !! if (hetero.eq.2) correction=tmp6
62 |
63 | c.rfo.
64 | ELSE !gaz-gaz
65 | c.old. CALL d1fds1(kgbar(igaz,in),kgskgbar(igaz),tmp5,
66 | c.old. & phicg(igaz),phig(igaz,in),tmp6)
67 | c.old. tmp7=1.D+0-tmp6/kgskgbar(igaz)/kgbar(igaz,in)
68 | c.old. CALL d1fds2(kgbar(igaz,iin),kgskgbar(igaz),tmp5,
69 | c.old. & phicg(igaz),phig(igaz,iin),tmp6)
70 | c.old. tmp7=tmp7-tmp6/kgskgbar(igaz)/kgbar(igaz,iin)
71 | c.old. CALL d2fds1ds2(kgbar(igaz,in),kgbar(igaz,iin)
72 | c.old. & ,kgskgbar(igaz)
73 | c.old. & ,tmp5,phicg(igaz)
74 | c.old. & ,phig(igaz,in),phig(igaz,iin),tmp6)
75 | c.old. tmp7=tmp7+tmp6/(kgskgbar(igaz))**2
76 | c.old. & /kgbar(igaz,in)/kgbar(igaz,iin)
77 | tdphi_j=kgbar(igaz,iin)/epocg(igaz)*(phig(igaz,iin)-phicg(igaz))
78 | tdphi_i=kgbar(igaz,in)/epocg(igaz)*(phig(igaz,in)-phicg(igaz))
79 | tphi=phicg(igaz)
80 | tphi_i=phig(igaz,in)
81 | tphi_j=phig(igaz,iin)
82 | tk_i=kgbar(igaz,in)
83 | tk_j=kgbar(igaz,iin)
84 | te=epocg(igaz)
85 | tb=te/tphi
86 | ta=(1.D+0 + 2.D+0*tb)
87 | tdb_i=(tphi*tk_i-te*tdphi_i)/tphi**2
88 | tdb_j=(tphi*tk_j-te*tdphi_j)/tphi**2
89 | tda_i=2.D+0*tdb_i
90 | tda_j=2.D+0*tdb_j
91 | td2phi=-tk_i*tk_j/te**2 * (tphi_i-2.D+0*tphi+tphi_j)
92 | c
93 | tam12=ta**(-0.5D+0)
94 | tap12=ta**(+0.5D+0)
95 | tam32=ta**(-1.5D+0)
96 | c
97 | tGa_i(igaz)=-tam12*(tk_i+tdphi_i*(1.D+0+tb-tap12))
98 | tGa_j(igaz)=-tam12*(tk_j+tdphi_j*(1.D+0+tb-tap12))
99 | c
100 | tdGa_ij(igaz)=0.5D+0*tda_j*tam32*(tk_i+tdphi_i*(1.D+0+tb-tap12))
101 | & -tam12*(td2phi*(1.D+0+tb-tap12)
102 | & +tdphi_i*(tdb_j-0.5D+0*tda_j*tam12))
103 | ! tdGa_ji=0.5D+0*tda_i*tam32*(tk_j+tdphi_j*(1.D+0+tb-tap12))
104 | ! & -tam12*(td2phi*(1.D+0+tb-tap12)
105 | ! & +tdphi_j*(tdb_i-0.5D+0*tda_i*tam12))
106 | ccccccccccccccccccccccc print *,'test',tdGa_ij-tdGa_ji,tdGa_ji
107 | c
108 | !tGc=tGa_i(igaz)*tGa_j(igaz)+tdGa_ij(igaz)
109 | c
110 | tGb_i(igaz)=-tk_i*tam12
111 | tGb_j(igaz)=-tk_j*tam12
112 | c
113 | tdGb_ij(igaz)=tk_i*tk_j/tphi*tam32
114 | c
115 | !tGd=tGb_i(igaz)*tGb_j(igaz)+tdGb_ij(igaz)
116 | c
117 | !! tmp7=tGc/tGd
118 | c.rfo.
119 | c.lat. rajoute condition
120 | !! if (hetero.eq.2) correction=tmp7
121 |
122 | c print *,'test',cor(igaz),tmp7
123 | c.rfo.
124 | ENDIF
125 | ELSE
126 | cor(igaz) = 1.d+0
127 | ! print *, 'cor(igaz) = 1.d+0',cor(igaz)
128 | !les morceaux de dérivée ne servent pas
129 | tGa_i(igaz) = 0.d+0
130 | tGa_j(igaz) = 0.d+0
131 | tGb_i(igaz) = 0.d+0
132 | tGb_j(igaz) = 0.d+0
133 | tdGa_ij(igaz) = 0.d+0
134 | tdGb_ij(igaz) = 0.d+0
135 | ENDIF
136 | ENDDO
137 | c.rfo.
138 | c wcg=cor(3)*correction
139 | c.rfo.
140 | c*****************************************************************************
141 | !IF ((epocg(igaz).GT.1.D-10).and.(hetero.eq.2)) THEN
142 | !IF ((tmp5.GT.1.D-10).and.(hetero.eq.2)) THEN
143 | IF (hetero.eq.2) THEN
144 | !correction = fonction (morceaux de dérivée)
145 | IF ((test1.EQ.1).AND.(test2.EQ.1)) THEN !paroi-paroi
146 | correction = 1.d+0
147 | ELSEIF (test1.EQ.1) THEN !paroi-gaz
148 | !tGa= tGa_j(3)
149 | !tGb= tGb_j(3)
150 | tGa = 0.d+0
151 | tGb = 0.d+0
152 | do igaz = 1, ngaz
153 | tGa = tGa + tGa_j(igaz)
154 | tGb = tGb + tGb_j(igaz)
155 | enddo
156 | correction = tGa / tGb
157 | ELSEIF (test2.EQ.1) THEN !gaz-paroi
158 | !tGa= tGa_i(3)
159 | !tGb= tGb_i(3)
160 | tGa = 0.d+0
161 | tGb = 0.d+0
162 | do igaz = 1, ngaz
163 | tGa = tGa + tGa_i(igaz)
164 | tGb = tGb + tGb_i(igaz)
165 | enddo
166 | correction = tGa / tGb
167 | ELSE !gaz-gaz
168 | !tGc=tGa_i(3)*tGa_j(3)+tdGa_ij(3)
169 | !tGd=tGb_i(3)*tGb_j(3)+tdGb_ij(3)
170 | tGc = 0.d+0
171 | !tmp6 = 0.d+0
172 | tGd = 0.d+0
173 | !tmp8 = 0.d+0
174 | do igaz = 1, ngaz
175 | tmp6 = 0.d+0
176 | tmp8 = 0.d+0
177 | tGc = tGc + tGa_i(igaz)*tGa_j(igaz)+tdGa_ij(igaz)
178 | tGd = tGd + tGb_i(igaz)*tGb_j(igaz)+tdGb_ij(igaz)
179 | do iz= 1, ngaz
180 | if (iz .ne. igaz) then
181 | tmp6 = tmp6 + tGa_i(igaz)*tGa_j(iz)
182 | tmp8 = tmp8 + tGb_i(igaz)*tGb_j(iz)
183 | endif
184 | enddo
185 | tGc = tGc +tmp6
186 | tGd = tGd +tmp8
187 | enddo
188 | correction = tGc / tGd
189 | ENDIF
190 | !finallement s'il n'y a aucun des trois gaz sur le trajet
191 | !on veut faire une correction mais il n'y a pas matiere a...
192 | if ((epocg(1).LT. 1.d-10).and.(epocg(2).LT.1.d-10)
193 | & .and.(epocg(3).LT.1.d-10 )) then
194 | correction = 1.d+0
195 | endif
196 | else
197 | correction = 1.d+0
198 | ! print *, 'correction = 1.d+0',correction
199 | ENDIF
200 |
201 |
202 | c******************************************************************************
203 | ! print*, cor(1), cor(2), cor(3), correction
204 | ! wcg = cor(1) * cor(2) * cor(3) * correction
205 | ! print*, cor(3) , correction
206 | wcg = cor(1)*cor(2)* cor(3) * correction