~maddevelopers/mg5amcnlo/2.5.3_lep

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
      subroutine jiotxx(fi,fo,tc,gc,gt,vmass,vwidth , jiot)
c
c This subroutine computes an off-shell vector boson wavefunction from a
c flowing-out fermion, a flowing-in fermion and a tensor boson.
c
c input:
c       complex fi(6)          : flow-in  fermion                   |fi>
c       complex fo(6)          : flow-out fermion                   <fo|
c       complex tc(18)         : input    tensor                       T
c       complex gc(2)          : coupling constants                  gvf
c       complex gt             : coupling constant      gtfv=-1/Lambda/2
c       real    vmass          : mass  of output vector v
c       real    vwidth         : width of output vector v
c
c output:
c       complex jiot(6)        : vector boson          j^mu(<fo|v,T|fi>)
c
c- by Q.Li - OCT. 2006
c
      implicit none
      double complex fi(6), fo(6), tc(18), gc(2), gt, jiot(6)
      double precision vmass, vwidth

      double complex ft(6,4)
      double complex d, T00, T12, T13, T14, T23, T24, T34
      double precision pv(4), pv2
      integer i
      
      double precision rZero, rTwo
      parameter( rZero = 0.0d0, rTwo = 2.0d0 )
      double complex cone
      parameter( cone = ( 0.0d0, 1.0d0 ))


      ft(1,1) = tc(1)
      ft(1,2) = tc(2)
      ft(1,3) = tc(3)
      ft(1,4) = tc(4)
      ft(2,1) = tc(5)
      ft(2,2) = tc(6)
      ft(2,3) = tc(7)
      ft(2,4) = tc(8)
      ft(3,1) = tc(9)
      ft(3,2) = tc(10)
      ft(3,3) = tc(11)
      ft(3,4) = tc(12)
      ft(4,1) = tc(13)
      ft(4,2) = tc(14)
      ft(4,3) = tc(15)
      ft(4,4) = tc(16)
      ft(5,1) = tc(17)
      ft(6,1) = tc(18)

      jiot(5) = fo(5) + ft(5,1) -fi(5)
      jiot(6) = fo(6) + ft(6,1) -fi(6)

      pv(1) = dreal(jiot(5))
      pv(2) = dreal(jiot(6))
      pv(3) = dimag(jiot(6))
      pv(4) = dimag(jiot(5))
 
      pv2=pv(1)**2-pv(2)**2-pv(3)**2-pv(4)**2
     

      T00 = ft(1,1)-ft(2,2)-ft(3,3)-ft(4,4)
      T12 = ft(1,2) + ft(2,1)
      T13 = ft(1,3) + ft(3,1)
      T14 = ft(1,4) + ft(4,1)
      T23 = ft(2,3) + ft(3,2)
      T24 = ft(2,4) + ft(4,2)
      T34 = ft(3,4) + ft(4,3)
      
      if ( vmass.gt.rZero ) then
         d =  1.0d0/dcmplx( pv2-vmass**2, vmass*vwidth )
      else
         d =  1.0d0/dcmplx( pv2, rZero )
      end if

      if (vmass.gt.rzero) then
         jiot(1) = fi(3)*gc(2)*(-((rtwo*T00*fo(2)*pv(1)*(-pv(2)
     &	 - cone*pv(3)))/vmass**2) + 
     &       fo(1)*((rtwo*T00) 
     &- (rtwo*T00*pv(1)*(pv(1) - pv(4)))/vmass**2)) + 
     &    fi(2)*gc(1)*(-((rtwo*T00*fo(3)*pv(1)
     &*(pv(2) - cone*pv(3)))/vmass**2) + 
     &       fo(4)*((rtwo*T00) - (rtwo*T00*pv(1)
     &*(pv(1) - pv(4)))/vmass**2)) + 
     &    fi(4)*gc(2)*(-((rtwo*T00*fo(1)*pv(1)
     &*(-pv(2) + cone*pv(3)))/vmass**2) + 
     &       fo(2)*((rtwo*T00) 
     &- (rtwo*T00*pv(1)*(pv(1) + pv(4)))/vmass**2)) + 
     &    fi(1)*gc(1)*(-((rtwo*T00*fo(4)*pv(1)*(pv(2)
     & + cone*pv(3)))/vmass**2) + 
     &       fo(3)*((rtwo*T00) 
     &- (rtwo*T00*pv(1)*(pv(1) + pv(4)))/vmass**2)) + 
     &    fi(4)*gc(2)*(fo(1)*(ft(1,2) - cone*ft(1,3) 
     &+ ft(2,1) - cone*ft(3,1))
     & + fo(2)*(-2*ft(1,1) - ft(1,4) - ft(4,1))) + 
     &    fi(1)*gc(1)*(fo(4)*(-ft(1,2) - cone*ft(1,3) - ft(2,1) 
     &- cone*ft(3,1)) + fo(3)*(-2*ft(1,1) - ft(1,4) - ft(4,1))) + 
     &    fi(3)*gc(2)*(fo(2)*(ft(1,2) + cone*ft(1,3) + ft(2,1) 
     &+ cone*ft(3,1)) + fo(1)*(-2*ft(1,1) + ft(1,4) + ft(4,1))) + 
     &    fi(2)*gc(1)*(fo(3)*(-ft(1,2) + cone*ft(1,3) - ft(2,1)
     & + cone*ft(3,1)) + fo(4)*(-2*ft(1,1) + ft(1,4) + ft(4,1))) + 
     &    (pv(1)*(fi(4)*gc(2)*(fo(1)*(-(T12*pv(1)) + cone*T13*pv(1)
     & + T23*(-(cone*pv(2)) + pv(3)) + T24*pv(4) - cone*T34*pv(4) + 
     &               rtwo*pv(2)*ft(2,2) - cone*rtwo*pv(3)*ft(3,3)) + 
     &            fo(2)*(-(T12*pv(2)) - T24*pv(2) - T13*pv(3) - 
     &T34*pv(3) - T14*(-pv(1) + pv(4)) + rtwo*pv(1)*ft(1,1) - 
     &               rtwo*pv(4)*ft(4,4))) + fi(1)*gc(1)*
     &          (fo(4)*(T12*pv(1) + cone*T13*pv(1) 
     &+ T23*(-(cone*pv(2)) - pv(3)) - T24*pv(4) - cone*T34*pv(4) - 
     &               rtwo*pv(2)*ft(2,2) - cone*rtwo*pv(3)*ft(3,3)) + 
     &            fo(3)*(-(T12*pv(2)) - T24*pv(2) - T13*pv(3) 
     &- T34*pv(3) - T14*(-pv(1) + pv(4)) + rtwo*pv(1)*ft(1,1) - 
     &               rtwo*pv(4)*ft(4,4))) + fi(3)*gc(2)*
     &          (fo(2)*(-(T12*pv(1)) - cone*T13*pv(1) 
     &+ T23*(cone*pv(2) + pv(3)) + T24*pv(4) + cone*T34*pv(4) + 
     &               rtwo*pv(2)*ft(2,2) + cone*rtwo*pv(3)*ft(3,3)) + 
     &            fo(1)*(-(T12*pv(2)) + T24*pv(2) - T13*pv(3)
     & + T34*pv(3) - T14*(pv(1) + pv(4)) + rtwo*pv(1)*ft(1,1) + 
     &               rtwo*pv(4)*ft(4,4))) + fi(2)*gc(1)*
     &          (fo(3)*(T12*pv(1) - cone*T13*pv(1) 
     &+ T23*(cone*pv(2) - pv(3)) - T24*pv(4) + cone*T34*pv(4) - 
     &               rtwo*pv(2)*ft(2,2) + cone*rtwo*pv(3)*ft(3,3)) + 
     &            fo(4)*(-(T12*pv(2)) + T24*pv(2) - T13*pv(3) 
     &+ T34*pv(3) - T14*(pv(1) + pv(4)) + rtwo*pv(1)*ft(1,1) + 
     &               rtwo*pv(4)*ft(4,4)))))/vmass**2
      
         jiot(2) = fi(3)*gc(2)*(fo(2)*((rtwo*T00) 
     &- (rtwo*T00*pv(2)*(-pv(2) - cone*pv(3)))/vmass**2) - 
     &       (rtwo*T00*fo(1)*pv(2)*(pv(1) - pv(4)))/vmass**2) + 
     &    fi(2)*gc(1)*(fo(3)*(-rtwo*T00 - (rtwo*T00*pv(2)*(pv(2)
     & - cone*pv(3)))/vmass**2) - 
     &       (rtwo*T00*fo(4)*pv(2)*(pv(1) - pv(4)))/vmass**2) + 
     &    fi(4)*gc(2)*(fo(1)*((rtwo*T00)
     & - (rtwo*T00*pv(2)*(-pv(2) + cone*pv(3)))/vmass**2) - 
     &       (rtwo*T00*fo(2)*pv(2)*(pv(1) + pv(4)))/vmass**2) + 
     &    fi(1)*gc(1)*(fo(4)*(-rtwo*T00 - (rtwo*T00*pv(2)*(pv(2) 
     &+ cone*pv(3)))/vmass**2) - 
     &       (rtwo*T00*fo(3)*pv(2)*(pv(1) + pv(4)))/vmass**2) + 
     &    fi(4)*gc(2)*(fo(1)*(2*ft(2,2) - cone*ft(2,3) - cone*ft(3,2))
     & + fo(2)*(-ft(1,2) - ft(2,1) - ft(2,4) - ft(4,2))) + 
     &    fi(1)*gc(1)*(fo(4)*(-2*ft(2,2) - cone*ft(2,3) - cone*ft(3,2)) 
     &+ fo(3)*(-ft(1,2) - ft(2,1) - ft(2,4) - ft(4,2))) + 
     &    fi(3)*gc(2)*(fo(2)*(2*ft(2,2) + cone*ft(2,3) + cone*ft(3,2)) 
     &+ fo(1)*(-ft(1,2) - ft(2,1) + ft(2,4) + ft(4,2))) + 
     &    fi(2)*gc(1)*(fo(3)*(-2*ft(2,2) + cone*ft(2,3) + cone*ft(3,2))
     & + fo(4)*(-ft(1,2) - ft(2,1) + ft(2,4) + ft(4,2))) + 
     &    (pv(2)*(fi(4)*gc(2)*(fo(1)*(-(T12*pv(1)) + cone*T13*pv(1)
     & + T23*(-(cone*pv(2)) + pv(3)) + T24*pv(4) - cone*T34*pv(4) + 
     &               rtwo*pv(2)*ft(2,2) - cone*rtwo*pv(3)*ft(3,3)) + 
     &            fo(2)*(-(T12*pv(2)) - T24*pv(2) - T13*pv(3)
     & - T34*pv(3) - T14*(-pv(1) + pv(4)) + rtwo*pv(1)*ft(1,1) - 
     &               rtwo*pv(4)*ft(4,4))) + fi(1)*gc(1)*
     &          (fo(4)*(T12*pv(1) + cone*T13*pv(1) + T23*(-(cone*pv(2))
     & - pv(3)) - T24*pv(4) - cone*T34*pv(4) - 
     &               rtwo*pv(2)*ft(2,2) - cone*rtwo*pv(3)*ft(3,3)) + 
     &            fo(3)*(-(T12*pv(2)) - T24*pv(2) - T13*pv(3)
     & - T34*pv(3) - T14*(-pv(1) + pv(4)) + rtwo*pv(1)*ft(1,1) - 
     &               rtwo*pv(4)*ft(4,4))) + fi(3)*gc(2)*
     &          (fo(2)*(-(T12*pv(1)) - cone*T13*pv(1) 
     &+ T23*(cone*pv(2) + pv(3)) + T24*pv(4) + cone*T34*pv(4) + 
     &               rtwo*pv(2)*ft(2,2) + cone*rtwo*pv(3)*ft(3,3)) + 
     &            fo(1)*(-(T12*pv(2)) + T24*pv(2) - T13*pv(3)
     & + T34*pv(3) - T14*(pv(1) + pv(4)) + rtwo*pv(1)*ft(1,1) + 
     &               rtwo*pv(4)*ft(4,4))) + fi(2)*gc(1)*
     &          (fo(3)*(T12*pv(1) - cone*T13*pv(1) 
     &+ T23*(cone*pv(2) - pv(3)) - T24*pv(4) + cone*T34*pv(4) - 
     &               rtwo*pv(2)*ft(2,2) + cone*rtwo*pv(3)*ft(3,3)) + 
     &            fo(4)*(-(T12*pv(2)) + T24*pv(2) - T13*pv(3) 
     &+ T34*pv(3) - T14*(pv(1) + pv(4)) + rtwo*pv(1)*ft(1,1) + 
     &               rtwo*pv(4)*ft(4,4)))))/vmass**2

         jiot(3) =fi(3)*gc(2)*(fo(2)*((cone*rtwo*T00)
     &	 - (rtwo*T00*pv(3)*(-pv(2) - cone*pv(3)))/vmass**2) - 
     &       (rtwo*T00*fo(1)*pv(3)*(pv(1) - pv(4)))/vmass**2) + 
     &    fi(2)*gc(1)*(fo(3)*((cone*rtwo*T00) 
     &- (rtwo*T00*pv(3)*(pv(2) - cone*pv(3)))/vmass**2) - 
     &       (rtwo*T00*fo(4)*pv(3)*(pv(1) - pv(4)))/vmass**2) + 
     &    fi(4)*gc(2)*(fo(1)*(-cone*rtwo*T00 
     &- (rtwo*T00*pv(3)*(-pv(2) + cone*pv(3)))/vmass**2) - 
     &       (rtwo*T00*fo(2)*pv(3)*(pv(1) + pv(4)))/vmass**2) + 
     &    fi(1)*gc(1)*(fo(4)*(-cone*rtwo*T00
     & - (rtwo*T00*pv(3)*(pv(2) + cone*pv(3)))/vmass**2) - 
     &       (rtwo*T00*fo(3)*pv(3)*(pv(1) + pv(4)))/vmass**2) + 
     &    fi(4)*gc(2)*(fo(1)*(ft(2,3) + ft(3,2) 
     &- 2*cone*ft(3,3)) + fo(2)*(-ft(1,3) - ft(3,1) 
     &- ft(3,4) - ft(4,3))) + 
     &    fi(1)*gc(1)*(fo(4)*(-ft(2,3) - ft(3,2) 
     &- 2*cone*ft(3,3)) + fo(3)*(-ft(1,3) - ft(3,1) 
     &- ft(3,4) - ft(4,3))) + 
     &    fi(3)*gc(2)*(fo(2)*(ft(2,3) + ft(3,2) + 2*cone*ft(3,3))
     & + fo(1)*(-ft(1,3) - ft(3,1) + ft(3,4) + ft(4,3))) + 
     &    fi(2)*gc(1)*(fo(3)*(-ft(2,3) - ft(3,2) + 2*cone*ft(3,3)) 
     &+ fo(4)*(-ft(1,3) - ft(3,1) + ft(3,4) + ft(4,3))) + 
     &    (pv(3)*(fi(4)*gc(2)*(fo(1)*(-(T12*pv(1)) + cone*T13*pv(1) 
     &+ T23*(-(cone*pv(2)) + pv(3)) + T24*pv(4) - cone*T34*pv(4) + 
     &               rtwo*pv(2)*ft(2,2) - cone*rtwo*pv(3)*ft(3,3)) + 
     &            fo(2)*(-(T12*pv(2)) - T24*pv(2) - T13*pv(3)
     &- T34*pv(3) - T14*(-pv(1) + pv(4)) + rtwo*pv(1)*ft(1,1) - 
     &               rtwo*pv(4)*ft(4,4))) + fi(1)*gc(1)*
     &          (fo(4)*(T12*pv(1) + cone*T13*pv(1) 
     &+ T23*(-(cone*pv(2)) - pv(3)) - T24*pv(4) - cone*T34*pv(4) - 
     &               rtwo*pv(2)*ft(2,2) - cone*rtwo*pv(3)*ft(3,3)) + 
     &            fo(3)*(-(T12*pv(2)) - T24*pv(2) - T13*pv(3)
     & - T34*pv(3) - T14*(-pv(1) + pv(4)) + rtwo*pv(1)*ft(1,1) - 
     &               rtwo*pv(4)*ft(4,4))) + fi(3)*gc(2)*
     &          (fo(2)*(-(T12*pv(1)) - cone*T13*pv(1) 
     &+ T23*(cone*pv(2) + pv(3)) + T24*pv(4) + cone*T34*pv(4) + 
     &               rtwo*pv(2)*ft(2,2) + cone*rtwo*pv(3)*ft(3,3)) + 
     &            fo(1)*(-(T12*pv(2)) + T24*pv(2) - T13*pv(3)
     & + T34*pv(3) - T14*(pv(1) + pv(4)) + rtwo*pv(1)*ft(1,1) + 
     &               rtwo*pv(4)*ft(4,4))) + fi(2)*gc(1)*
     &          (fo(3)*(T12*pv(1) - cone*T13*pv(1) 
     &+ T23*(cone*pv(2) - pv(3)) - T24*pv(4) + cone*T34*pv(4) - 
     &               rtwo*pv(2)*ft(2,2) + cone*rtwo*pv(3)*ft(3,3)) + 
     &            fo(4)*(-(T12*pv(2)) + T24*pv(2) - T13*pv(3)
     & + T34*pv(3) - T14*(pv(1) + pv(4)) + rtwo*pv(1)*ft(1,1) + 
     &               rtwo*pv(4)*ft(4,4)))))/vmass**2

         jiot(4) =fi(3)*gc(2)*(-((rtwo*T00*fo(2)*(-pv(2) 
     &	- cone*pv(3))*pv(4))/vmass**2) + 
     &       fo(1)*((rtwo*T00) - (rtwo*T00*(pv(1)
     & - pv(4))*pv(4))/vmass**2)) + 
     &    fi(2)*gc(1)*(-((rtwo*T00*fo(3)*(pv(2)
     & - cone*pv(3))*pv(4))/vmass**2) + 
     &       fo(4)*((rtwo*T00) - (rtwo*T00*(pv(1)
     & - pv(4))*pv(4))/vmass**2)) + 
     &    fi(4)*gc(2)*(-((rtwo*T00*fo(1)*(-pv(2) 
     &+ cone*pv(3))*pv(4))/vmass**2) + 
     &       fo(2)*(-rtwo*T00
     & - (rtwo*T00*pv(4)*(pv(1) + pv(4)))/vmass**2)) + 
     &    fi(1)*gc(1)*(-((rtwo*T00*fo(4)*(pv(2)
     & + cone*pv(3))*pv(4))/vmass**2) + 
     &       fo(3)*(-rtwo*T00 - (rtwo*T00*pv(4)*(pv(1)
     & + pv(4)))/vmass**2)) + 
     &    fi(4)*gc(2)*(fo(1)*(ft(2,4) - cone*ft(3,4)
     & + ft(4,2) - cone*ft(4,3)) + fo(2)*(-ft(1,4) 
     &- ft(4,1) - 2*ft(4,4))) + 
     &    fi(1)*gc(1)*(fo(4)*(-ft(2,4) - cone*ft(3,4)
     & - ft(4,2) - cone*ft(4,3)) + fo(3)*(-ft(1,4)
     & - ft(4,1) - 2*ft(4,4))) + 
     &    fi(3)*gc(2)*(fo(2)*(ft(2,4) + cone*ft(3,4) 
     &+ ft(4,2) + cone*ft(4,3)) + fo(1)*(-ft(1,4) 
     &- ft(4,1) + 2*ft(4,4))) + 
     &    fi(2)*gc(1)*(fo(3)*(-ft(2,4) + cone*ft(3,4)
     & - ft(4,2) + cone*ft(4,3)) + fo(4)*(-ft(1,4)
     & - ft(4,1) + 2*ft(4,4))) + 
     &    (pv(4)*(fi(4)*gc(2)*(fo(1)*(-(T12*pv(1)) 
     &+ cone*T13*pv(1) + T23*(-(cone*pv(2)) + pv(3))
     & + T24*pv(4) - cone*T34*pv(4) + 
     &               rtwo*pv(2)*ft(2,2) - cone*rtwo*pv(3)*ft(3,3)) + 
     &            fo(2)*(-(T12*pv(2)) - T24*pv(2) - T13*pv(3)
     & - T34*pv(3) - T14*(-pv(1) + pv(4)) + rtwo*pv(1)*ft(1,1) - 
     &               rtwo*pv(4)*ft(4,4))) + fi(1)*gc(1)*
     &          (fo(4)*(T12*pv(1) + cone*T13*pv(1) 
     &+ T23*(-(cone*pv(2)) - pv(3)) - T24*pv(4) - cone*T34*pv(4) - 
     &               rtwo*pv(2)*ft(2,2) - cone*rtwo*pv(3)*ft(3,3)) + 
     &            fo(3)*(-(T12*pv(2)) - T24*pv(2) - T13*pv(3)
     & - T34*pv(3) - T14*(-pv(1) + pv(4)) + rtwo*pv(1)*ft(1,1) - 
     &               rtwo*pv(4)*ft(4,4))) + fi(3)*gc(2)*
     &          (fo(2)*(-(T12*pv(1)) - cone*T13*pv(1) 
     &+ T23*(cone*pv(2) + pv(3)) + T24*pv(4) + cone*T34*pv(4) + 
     &               rtwo*pv(2)*ft(2,2) + cone*rtwo*pv(3)*ft(3,3)) + 
     &            fo(1)*(-(T12*pv(2)) + T24*pv(2) - T13*pv(3) 
     &+ T34*pv(3) - T14*(pv(1) + pv(4)) + rtwo*pv(1)*ft(1,1) + 
     &               rtwo*pv(4)*ft(4,4))) + fi(2)*gc(1)*
     &          (fo(3)*(T12*pv(1) - cone*T13*pv(1) + T23*(cone*pv(2) 
     &- pv(3)) - T24*pv(4) + cone*T34*pv(4) - 
     &               rtwo*pv(2)*ft(2,2) + cone*rtwo*pv(3)*ft(3,3)) + 
     &            fo(4)*(-(T12*pv(2)) + T24*pv(2) - T13*pv(3)
     & + T34*pv(3) - T14*(pv(1) + pv(4)) + rtwo*pv(1)*ft(1,1) + 
     &               rtwo*pv(4)*ft(4,4)))))/vmass**2

      else

         jiot(1) = rtwo*T00*fi(1)*fo(3)*gc(1)
     & + rtwo*T00*fi(2)*fo(4)*gc(1) + rtwo*T00*fi(3)*fo(1)*gc(2) + 
     &    rtwo*T00*fi(4)*fo(2)*gc(2) + fi(4)*gc(2)*(fo(1)*(ft(1,2) 
     &- cone*ft(1,3) + ft(2,1) - cone*ft(3,1)) + 
     &       fo(2)*(-2*ft(1,1) - ft(1,4) - ft(4,1))) + 
     &    fi(1)*gc(1)*(fo(4)*(-ft(1,2) - cone*ft(1,3) - ft(2,1) 
     &- cone*ft(3,1)) + fo(3)*(-2*ft(1,1) - ft(1,4) - ft(4,1))) + 
     &    fi(3)*gc(2)*(fo(2)*(ft(1,2) + cone*ft(1,3) + ft(2,1) 
     &+ cone*ft(3,1)) + fo(1)*(-2*ft(1,1) + ft(1,4) + ft(4,1))) + 
     &    fi(2)*gc(1)*(fo(3)*(-ft(1,2) + cone*ft(1,3) - ft(2,1)
     & + cone*ft(3,1)) + fo(4)*(-2*ft(1,1) + ft(1,4) + ft(4,1)))

         jiot(2) = -(rtwo*T00*fi(2)*fo(3)*gc(1)) 
     &  - rtwo*T00*fi(1)*fo(4)*gc(1)
     &	 + rtwo*T00*fi(4)*fo(1)*gc(2) + 
     &    rtwo*T00*fi(3)*fo(2)*gc(2) + fi(4)*gc(2)*(fo(1)*(2*ft(2,2) 
     &- cone*ft(2,3) - cone*ft(3,2)) + 
     &       fo(2)*(-ft(1,2) - ft(2,1) - ft(2,4) - ft(4,2))) + 
     &    fi(1)*gc(1)*(fo(4)*(-2*ft(2,2) - cone*ft(2,3) - cone*ft(3,2))
     & + fo(3)*(-ft(1,2) - ft(2,1) - ft(2,4) - ft(4,2))) + 
     &    fi(3)*gc(2)*(fo(2)*(2*ft(2,2) + cone*ft(2,3) + cone*ft(3,2))
     & + fo(1)*(-ft(1,2) - ft(2,1) + ft(2,4) + ft(4,2))) + 
     &    fi(2)*gc(1)*(fo(3)*(-2*ft(2,2) + cone*ft(2,3) + cone*ft(3,2)) 
     &+ fo(4)*(-ft(1,2) - ft(2,1) + ft(2,4) + ft(4,2)))

         jiot(3) = cone*rtwo*T00*fi(2)*fo(3)*gc(1) 
     &	- cone*rtwo*T00*fi(1)*fo(4)*gc(1)
     & - cone*rtwo*T00*fi(4)*fo(1)*gc(2) + 
     &    cone*rtwo*T00*fi(3)*fo(2)*gc(2) + fi(4)*gc(2)*
     &     (fo(1)*(ft(2,3) + ft(3,2) - 2*cone*ft(3,3))
     & + fo(2)*(-ft(1,3) - ft(3,1) - ft(3,4) - ft(4,3))) + 
     &    fi(1)*gc(1)*(fo(4)*(-ft(2,3) - ft(3,2) - 2*cone*ft(3,3))
     & + fo(3)*(-ft(1,3) - ft(3,1) - ft(3,4) - ft(4,3))) + 
     &    fi(3)*gc(2)*(fo(2)*(ft(2,3) + ft(3,2) + 2*cone*ft(3,3)) 
     &+ fo(1)*(-ft(1,3) - ft(3,1) + ft(3,4) + ft(4,3))) + 
     &    fi(2)*gc(1)*(fo(3)*(-ft(2,3) - ft(3,2) + 2*cone*ft(3,3)) 
     &+ fo(4)*(-ft(1,3) - ft(3,1) + ft(3,4) + ft(4,3)))
       
         jiot(4) = -(rtwo*T00*fi(1)*fo(3)*gc(1))
     &	 + rtwo*T00*fi(2)*fo(4)*gc(1) + rtwo*T00*fi(3)*fo(1)*gc(2) - 
     &    rtwo*T00*fi(4)*fo(2)*gc(2) + fi(4)*gc(2)*(fo(1)*(ft(2,4) 
     &- cone*ft(3,4) + ft(4,2) - cone*ft(4,3)) + 
     &       fo(2)*(-ft(1,4) - ft(4,1) - 2*ft(4,4))) + 
     &    fi(1)*gc(1)*(fo(4)*(-ft(2,4) - cone*ft(3,4) - ft(4,2)
     & - cone*ft(4,3)) + fo(3)*(-ft(1,4) - ft(4,1) - 2*ft(4,4))) + 
     &    fi(3)*gc(2)*(fo(2)*(ft(2,4) + cone*ft(3,4) + ft(4,2)
     & + cone*ft(4,3)) + fo(1)*(-ft(1,4) - ft(4,1) + 2*ft(4,4))) + 
     &    fi(2)*gc(1)*(fo(3)*(-ft(2,4) + cone*ft(3,4) - ft(4,2) 
     &+ cone*ft(4,3)) + fo(4)*(-ft(1,4) - ft(4,1) + 2*ft(4,4)))  

      endif

      do i = 1,4
         jiot(i) = -jiot(i)*d*gt
      end do

      return
      end