V3FIT
All Classes Namespaces Files Functions Variables Enumerations Macros Pages
forces.f
1  SUBROUTINE forces_par
2  USE vmec_main, p5 => cp5
3  USE realspace
4  USE vforces, ru12 => pazmn_e, zu12 => parmn_e,
5  & pazmn_e => pazmn_e, parmn_e => parmn_e,
6  & lv_e => pcrmn_e, lu_e => pczmn_e, lu_o => pczmn_o,
7  & pcrmn_e => pcrmn_e, pczmn_e => pczmn_e,
8  & pczmn_o => pczmn_o
9  USE parallel_include_module
10  USE timer_sub
11  IMPLICIT NONE
12 C-----------------------------------------------
13 C L o c a l P a r a m e t e r s
14 C-----------------------------------------------
15  REAL(dp), PARAMETER :: p25 = p5*p5, dshalfds=p25
16 C-----------------------------------------------
17 C L o c a l V a r i a b l e s
18 C-----------------------------------------------
19  INTEGER :: l, ndim
20  INTEGER :: i, j, k, nsmin, nsmax
21  REAL(dp), DIMENSION(:,:), POINTER ::
22  & bsqr, gvvs, guvs, guus
23  REAL(dp), ALLOCATABLE, DIMENSION(:) :: bcastbuf
24 C-----------------------------------------------
25  IF (.NOT.lactive .AND. .NOT.lfreeb) RETURN
26 
27  CALL second0 (tforon)
28 
29  ndim = 1+nrzt
30 
31  nsmin=tlglob; nsmax=t1rglob
32 
33 ! POINTER ALIASES
34  bsqr => pextra1(:,:,1); gvvs => pextra2(:,:,1)
35  guvs => pextra3(:,:,1); guus => pextra4(:,:,1)
36 
37  lu_e(:,1) = 0; lv_e(:,1) = 0
38  pguu(:,1) = 0; pguv(:,1) = 0; pgvv(:,1) = 0
39 
40  DO l = nsmin, nsmax
41  guus(:,l) = pguu(:,l)*pshalf(:,l)
42  guvs(:,l) = pguv(:,l)*pshalf(:,l)
43  gvvs(:,l) = pgvv(:,l)* pshalf(:,l)
44 
45  parmn_e(:,l) = ohs*zu12(:,l)*lu_e(:,l)
46  pazmn_e(:,l) =-ohs*ru12(:,l)*lu_e(:,l)
47  pbrmn_e(:,l) = pbrmn_e(:,l)*lu_e(:,l)
48  pbzmn_e(:,l) =-pbzmn_e(:,l)*lu_e(:,l)
49  bsqr(:,l) = dshalfds*lu_e(:,l)/pshalf(:,l)
50 
51  parmn_o(:,l) = parmn_e(:,l)*pshalf(:,l)
52  pazmn_o(:,l) = pazmn_e(:,l)*pshalf(:,l)
53  pbrmn_o(:,l) = pbrmn_e(:,l)*pshalf(:,l)
54  pbzmn_o(:,l) = pbzmn_e(:,l)*pshalf(:,l)
55  END DO
56 
57 !
58 ! CONSTRUCT CYLINDRICAL FORCE KERNELS
59 ! NOTE: presg(ns+1) == 0, AND WILL BE "FILLED IN" AT EDGE
60 ! FOR FREE-BOUNDARY BY RBSQ
61 !
62 !DIR$ IVDEP
63  nsmin=tlglob; nsmax=min(ns-1,trglob)
64  DO l = nsmin, nsmax
65  pguu(:,l) = p5*(pguu(:,l) + pguu(:,l+1))
66  pgvv(:,l) = p5*(pgvv(:,l) + pgvv(:,l+1))
67  bsqr(:,l) = bsqr(:,l) + bsqr(:,l+1)
68  guus(:,l) = p5*(guus(:,l) + guus(:,l+1))
69  gvvs(:,l) = p5*(gvvs(:,l) + gvvs(:,l+1))
70  END DO
71 
72  IF (trglob .ge. ns) THEN
73  pguu(:,ns) = p5*pguu(:,ns)
74  pgvv(:,ns) = p5*pgvv(:,ns)
75  guus(:,ns) = p5*guus(:,ns)
76  gvvs(:,ns) = p5*gvvs(:,ns)
77  END IF
78 
79 !DIR$ IVDEP
80  nsmin=tlglob; nsmax=min(ns-1,trglob)
81  DO l = nsmin, nsmax
82  parmn_e(:,l) = parmn_e(:,l+1) - parmn_e(:,l)
83  & + p5*(lv_e(:,l) + lv_e(:,l+1))
84  pazmn_e(:,l) = pazmn_e(:,l+1) - pazmn_e(:,l)
85  pbrmn_e(:,l) = p5*(pbrmn_e(:,l) + pbrmn_e(:,l+1))
86  pbzmn_e(:,l) = p5*(pbzmn_e(:,l) + pbzmn_e(:,l+1))
87  END DO
88 
89  parmn_e(:,ns) = - parmn_e(:,ns) + p5*lv_e(:,ns)
90  pazmn_e(:,ns) = - pazmn_e(:,ns)
91  pbrmn_e(:,ns) = p5*pbrmn_e(:,ns)
92  pbzmn_e(:,ns) = p5*pbzmn_e(:,ns)
93 
94  nsmin=tlglob; nsmax=t1rglob
95  DO l = nsmin, nsmax
96  parmn_e(:,l) = parmn_e(:,l)
97  & - (gvvs(:,l)*pr1(:,l,1) + pgvv(:,l)*pr1(:,l,0))
98  pbrmn_e(:,l) = pbrmn_e(:,l) + bsqr(:,l)*pz1(:,l,1)
99  & - (guus(:,l)*pru(:,l,1) + pguu(:,l)*pru(:,l,0))
100  pbzmn_e(:,l) = pbzmn_e(:,l) - (bsqr(:,l)*pr1(:,l,1)
101  & + guus(:,l)*pzu(:,l,1) + pguu(:,l)*pzu(:,l,0))
102  lv_e(:,l) = lv_e(:,l)*pshalf(:,l)
103  lu_o(:,l) = dshalfds*lu_e(:,l)
104  END DO
105 
106  nsmin=tlglob; nsmax=min(ns-1,trglob)
107 !DIR$ IVDEP
108  DO l = nsmin, nsmax
109  parmn_o(:,l) = parmn_o(:,l+1) - parmn_o(:,l)
110  & - pzu(:,l,0)*bsqr(:,l)
111  & + p5*(lv_e(:,l)+lv_e(:,l+1))
112  pazmn_o(:,l) = pazmn_o(:,l+1) - pazmn_o(:,l)
113  & + pru(:,l,0)*bsqr(:,l)
114  pbrmn_o(:,l) = p5*(pbrmn_o(:,l) + pbrmn_o(:,l+1))
115  pbzmn_o(:,l) = p5*(pbzmn_o(:,l) + pbzmn_o(:,l+1))
116  lu_o(:,l) = lu_o(:,l) + lu_o(:,l+1)
117  END DO
118 
119  parmn_o(:,ns) = - parmn_o(:,ns) - pzu(:,ns,0)*bsqr(:,ns)
120  & + p5*lv_e(:,ns)
121  pazmn_o(:,ns) = - pazmn_o(:,ns) + pru(:,ns,0)*bsqr(:,ns)
122  pbrmn_o(:,ns) = p5*pbrmn_o(:,ns)
123  pbzmn_o(:,ns) = p5*pbzmn_o(:,ns)
124  lu_o(:,ns) = lu_o(:,ns)
125 
126  nsmin=tlglob; nsmax=trglob
127  DO l = nsmin, nsmax
128  pguu(:,l) = pguu(:,l) * psqrts(:,l)**2
129  bsqr(:,l) = pgvv(:,l) * psqrts(:,l)**2
130  END DO
131 
132  DO l = nsmin, nsmax
133  parmn_o(:,l) = parmn_o(:,l) - (pzu(:,l,1)*lu_o(:,l)
134  & + bsqr(:,l)*pr1(:,l,1) + gvvs(:,l)*pr1(:,l,0))
135  pazmn_o(:,l) = pazmn_o(:,l) + pru(:,l,1)*lu_o(:,l)
136  pbrmn_o(:,l) = pbrmn_o(:,l) + pz1(:,l,1)*lu_o(:,l)
137  & -(pguu(:,l)*pru(:,l,1) + guus(:,l)*pru(:,l,0))
138  pbzmn_o(:,l) = pbzmn_o(:,l) - (pr1(:,l,1)*lu_o(:,l)
139  & + pguu(:,l)*pzu(:,l,1) + guus(:,l)*pzu(:,l,0))
140  END DO
141 
142  IF (lthreed) THEN
143 !DIR$ IVDEP
144  nsmin=tlglob; nsmax=min(ns-1,trglob)
145  DO l = nsmin, nsmax
146  pguv(:,l) = p5*(pguv(:,l) + pguv(:,l+1))
147  guvs(:,l) = p5*(guvs(:,l) + guvs(:,l+1))
148  END DO
149  pguv(:,ns) = p5*pguv(:,ns)
150  guvs(:,ns) = p5*guvs(:,ns)
151 
152  nsmin=tlglob; nsmax=trglob
153  DO l = nsmin, nsmax
154  pbrmn_e(:,l) = pbrmn_e(:,l)
155  & - (pguv(:,l)*prv(:,l,0) + guvs(:,l)*prv(:,l,1))
156  pbzmn_e(:,l) = pbzmn_e(:,l)
157  & - (pguv(:,l)*pzv(:,l,0) + guvs(:,l)*pzv(:,l,1))
158  pcrmn_e(:,l) = pguv(:,l)*pru(:,l,0) + pgvv(:,l)*prv(:,l,0)
159  & + gvvs(:,l)*prv(:,l,1) + guvs(:,l)*pru(:,l,1)
160  pczmn_e(:,l) = pguv(:,l)*pzu(:,l,0) + pgvv(:,l)*pzv(:,l,0)
161  & + gvvs(:,l)*pzv(:,l,1) + guvs(:,l)*pzu(:,l,1)
162  pguv(:,l) = pguv(:,l)*psqrts(:,l)*psqrts(:,l)
163  pbrmn_o(:,l) = pbrmn_o(:,l)
164  & - (guvs(:,l)*prv(:,l,0) + pguv(:,l)*prv(:,l,1))
165  pbzmn_o(:,l) = pbzmn_o(:,l)
166  & - (guvs(:,l)*pzv(:,l,0) + pguv(:,l)*pzv(:,l,1))
167  pcrmn_o(:,l) = guvs(:,l)*pru(:,l,0) + gvvs(:,l)*prv(:,l,0)
168  & + bsqr(:,l)*prv(:,l,1) + pguv(:,l)*pru(:,l,1)
169  pczmn_o(:,l) = guvs(:,l)*pzu(:,l,0) + gvvs(:,l)*pzv(:,l,0)
170  & + bsqr(:,l)*pzv(:,l,1) + pguv(:,l)*pzu(:,l,1)
171  END DO
172  ENDIF
173 
174 !
175 ! ASSIGN EDGE FORCES (JS = NS) FOR FREE BOUNDARY CALCULATION
176 !
177  IF (ivac .GE. 1) THEN
178 
179  DO k = 1, ntheta3
180  DO j = 1, nzeta
181  l = (k-1)*nzeta + j
182  parmn_e(l,ns) = parmn_e(l,ns) + pzu0(l,ns)*rbsq(l)
183  parmn_o(l,ns) = parmn_o(l,ns) + pzu0(l,ns)*rbsq(l)
184  pazmn_e(l,ns) = pazmn_e(l,ns) - pru0(l,ns)*rbsq(l)
185  pazmn_o(l,ns) = pazmn_o(l,ns) - pru0(l,ns)*rbsq(l)
186  END DO
187  END DO
188 
189  ENDIF
190 
191  100 CONTINUE
192 
193 !
194 ! COMPUTE CONSTRAINT FORCE KERNELS
195 !
196 #ifndef _HBANGLE
197  DO l = nsmin, nsmax
198  prcon(:,l,0) = (prcon(:,l,0)-prcon0(:,l))*pgcon(:,l)
199  pzcon(:,l,0) = (pzcon(:,l,0)-pzcon0(:,l))*pgcon(:,l)
200  pbrmn_e(:,l) = pbrmn_e(:,l) + prcon(:,l,0)
201  pbzmn_e(:,l) = pbzmn_e(:,l) + pzcon(:,l,0)
202  pbrmn_o(:,l) = pbrmn_o(:,l)+ prcon(:,l,0)*psqrts(:,l)
203  pbzmn_o(:,l) = pbzmn_o(:,l)+ pzcon(:,l,0)*psqrts(:,l)
204  prcon(:,l,0) = pru0(:,l) * pgcon(:,l)
205  pzcon(:,l,0) = pzu0(:,l) * pgcon(:,l)
206  prcon(:,l,1) = prcon(:,l,0) * psqrts(:,l)
207  pzcon(:,l,1) = pzcon(:,l,0) * psqrts(:,l)
208  END DO
209 #endif
210  CALL second0 (tforoff)
211  timer(tfor) = timer(tfor) + (tforoff - tforon)
212  forces_time = timer(tfor)
213 
214  END SUBROUTINE forces_par
215 
216  SUBROUTINE forces
217  USE vmec_main, p5 => cp5
218  USE realspace
219  USE vforces, ru12 => azmn_e, zu12 => armn_e,
220  & azmn_e => azmn_e, armn_e => armn_e,
221  & lv_e => crmn_e, lu_e => czmn_e, lu_o => czmn_o,
222  & crmn_e => crmn_e, czmn_e => czmn_e, czmn_o => czmn_o
223  USE timer_sub
224 
225  IMPLICIT NONE
226 C-----------------------------------------------
227 C L o c a l P a r a m e t e r s
228 C-----------------------------------------------
229  REAL(dp), PARAMETER :: p25 = p5*p5, dshalfds=p25
230 C-----------------------------------------------
231 C L o c a l V a r i a b l e s
232 C-----------------------------------------------
233  INTEGER :: l, ndim
234  REAL(dp), DIMENSION(:), POINTER ::
235  & bsqr, gvvs, guvs, guus
236 
237 C-----------------------------------------------
238  ndim = 1+nrzt
239  CALL second0 (tforon)
240 
241 ! POINTER ALIASES
242  bsqr => extra1(:,1); gvvs => extra2(:,1)
243  guvs => extra3(:,1); guus => extra4(:,1)
244 
245 !
246 ! ON ENTRY, ARMN=ZU,BRMN=ZS,AZMN=RU,BZMN=RS,LU=R*BSQ,LV = BSQ*SQRT(G)/R12
247 ! HERE, XS (X=Z,R) DO NOT INCLUDE DERIVATIVE OF EXPLICIT SQRT(S)
248 ! BSQ = |B|**2/2 + p
249 ! GIJ = (BsupI * BsupJ) * SQRT(G) (I,J = U,V)
250 ! IT IS ESSENTIAL THAT LU,LV AT j=1 ARE ZERO INITIALLY
251 !
252 ! SOME OF THE BIGGER LOOPS WERE SPLIT TO FACILITATE CACHE
253 ! HITS, PIPELINING ON RISCS
254 !
255 ! FOR OPTIMIZATION ON CRAY, MUST USE COMPILER DIRECTIVES TO
256 ! GET VECTORIZATION OF LOOPS INVOLVING POINTERS!
257 !
258 !
259 ! ORIGIN OF VARIOUS TERMS
260 !
261 ! LU : VARIATION OF DOMINANT RADIAL-DERIVATIVE TERMS IN JACOBIAN
262 !
263 ! LV : VARIATION OF R-TERM IN JACOBIAN
264 !
265 ! GVV: VARIATION OF R**2-TERM AND Rv**2,Zv**2 IN gvv
266 !
267 ! GUU, GUV: VARIATION OF Ru, Rv, Zu, Zv IN guu, guv
268 !
269 
270  lu_e(1:ndim:ns) = 0; lv_e(1:ndim:ns) = 0
271  guu(1:ndim:ns) = 0; guv(1:ndim:ns) = 0; gvv(1:ndim:ns) = 0
272  guus = guu*shalf; guvs = guv*shalf; gvvs = gvv*shalf
273 
274  armn_e = ohs*zu12 * lu_e
275  azmn_e =-ohs*ru12 * lu_e
276  brmn_e = brmn_e * lu_e
277  bzmn_e =-bzmn_e * lu_e
278  bsqr = dshalfds*lu_e/shalf
279 
280  armn_o(1:ndim) = armn_e(1:ndim) *shalf
281  azmn_o(1:ndim) = azmn_e(1:ndim) *shalf
282  brmn_o(1:ndim) = brmn_e(1:ndim) *shalf
283  bzmn_o(1:ndim) = bzmn_e(1:ndim) *shalf
284 
285 !
286 ! CONSTRUCT CYLINDRICAL FORCE KERNELS
287 ! NOTE: presg(ns+1) == 0, AND WILL BE "FILLED IN" AT EDGE
288 ! FOR FREE-BOUNDARY BY RBSQ
289 !
290 !DIR$ IVDEP
291  DO l = 1, nrzt
292  guu(l) = p5*(guu(l) + guu(l+1))
293  gvv(l) = p5*(gvv(l) + gvv(l+1))
294  bsqr(l) = bsqr(l) + bsqr(l+1)
295  guus(l) = p5*(guus(l) + guus(l+1))
296  gvvs(l) = p5*(gvvs(l) + gvvs(l+1))
297  END DO
298 
299 !DIR$ IVDEP
300  DO l = 1, nrzt
301  armn_e(l) = armn_e(l+1) - armn_e(l) + p5*(lv_e(l) + lv_e(l+1))
302  azmn_e(l) = azmn_e(l+1) - azmn_e(l)
303  brmn_e(l) = p5*(brmn_e(l) + brmn_e(l+1))
304  bzmn_e(l) = p5*(bzmn_e(l) + bzmn_e(l+1))
305  END DO
306 
307  armn_e(:nrzt) = armn_e(:nrzt) - (gvvs(:nrzt)*r1(:nrzt,1)
308  & + gvv(:nrzt)*r1(:nrzt,0))
309  brmn_e(:nrzt) = brmn_e(:nrzt) + bsqr(:nrzt)*z1(:nrzt,1)
310  & -(guus(:nrzt)*ru(:nrzt,1) + guu(:nrzt)*ru(:nrzt,0))
311  bzmn_e(:nrzt) = bzmn_e(:nrzt) - (bsqr(:nrzt)*r1(:nrzt,1)
312  & + guus(:nrzt)*zu(:nrzt,1) + guu(:nrzt)*zu(:nrzt,0))
313  lv_e(1:ndim) = lv_e(1:ndim)*shalf(1:ndim)
314  lu_o(1:ndim) = dshalfds*lu_e(1:ndim)
315 
316 !DIR$ IVDEP
317  DO l = 1, nrzt
318  armn_o(l) = armn_o(l+1) - armn_o(l) - zu(l,0)*bsqr(l)
319  & + p5*(lv_e(l) + lv_e(l+1))
320  azmn_o(l) = azmn_o(l+1) - azmn_o(l) + ru(l,0)*bsqr(l)
321  brmn_o(l) = p5*(brmn_o(l) + brmn_o(l+1))
322  bzmn_o(l) = p5*(bzmn_o(l) + bzmn_o(l+1))
323  lu_o(l) = lu_o(l) + lu_o(l+1)
324  END DO
325 
326  guu(1:nrzt) = guu(1:nrzt) * sqrts(1:nrzt)**2
327  bsqr(1:nrzt) = gvv(1:nrzt) * sqrts(1:nrzt)**2
328 
329  armn_o(:nrzt) = armn_o(:nrzt) - (zu(:nrzt,1)*lu_o(:nrzt)
330  & + bsqr(:nrzt)*r1(:nrzt,1) + gvvs(:nrzt)*r1(:nrzt,0))
331  azmn_o(:nrzt) = azmn_o(:nrzt) + ru(:nrzt,1)*lu_o(:nrzt)
332  brmn_o(:nrzt) = brmn_o(:nrzt) + z1(:nrzt,1)*lu_o(:nrzt)
333  & -(guu(:nrzt)*ru(:nrzt,1) + guus(:nrzt)*ru(:nrzt,0))
334  bzmn_o(:nrzt) = bzmn_o(:nrzt) - (r1(:nrzt,1)*lu_o(:nrzt)
335  & + guu(:nrzt)*zu(:nrzt,1) + guus(:nrzt)*zu(:nrzt,0))
336 
337  IF (lthreed) THEN
338 !DIR$ IVDEP
339  DO l = 1, nrzt
340  guv(l) = p5*(guv(l) + guv(l+1))
341  guvs(l) = p5*(guvs(l) + guvs(l+1))
342  END DO
343 
344  brmn_e(:nrzt) = brmn_e(:nrzt)
345  & - (guv(:nrzt)*rv(:nrzt,0) + guvs(:nrzt)*rv(:nrzt,1))
346  bzmn_e(:nrzt) = bzmn_e(:nrzt)
347  & - (guv(:nrzt)*zv(:nrzt,0) + guvs(:nrzt)*zv(:nrzt,1))
348  crmn_e(:nrzt) = guv(:nrzt) *ru(:nrzt,0)
349  & + gvv(:nrzt) *rv(:nrzt,0)
350  & + gvvs(:nrzt)*rv(:nrzt,1) + guvs(:nrzt)*ru(:nrzt,1)
351  czmn_e(:nrzt) = guv(:nrzt) *zu(:nrzt,0)
352  & + gvv(:nrzt) *zv(:nrzt,0)
353  & + gvvs(:nrzt)*zv(:nrzt,1) + guvs(:nrzt)*zu(:nrzt,1)
354  guv(:nrzt) = guv(:nrzt) *sqrts(:nrzt)*sqrts(:nrzt)
355  brmn_o(:nrzt) = brmn_o(:nrzt)
356  & - (guvs(:nrzt)*rv(:nrzt,0) + guv(:nrzt)*rv(:nrzt,1))
357  bzmn_o(:nrzt) = bzmn_o(:nrzt)
358  & - (guvs(:nrzt)*zv(:nrzt,0) + guv(:nrzt)*zv(:nrzt,1))
359  crmn_o(:nrzt) = guvs(:nrzt)*ru(:nrzt,0)
360  & + gvvs(:nrzt)*rv(:nrzt,0)
361  & + bsqr(:nrzt)*rv(:nrzt,1) + guv(:nrzt) *ru(:nrzt,1)
362  czmn_o(:nrzt) = guvs(:nrzt)*zu(:nrzt,0)
363  & + gvvs(:nrzt)*zv(:nrzt,0)
364  & + bsqr(:nrzt)*zv(:nrzt,1) + guv(:nrzt) *zu(:nrzt,1)
365  ENDIF
366 
367 !
368 ! ASSIGN EDGE FORCES (JS = NS) FOR FREE BOUNDARY CALCULATION
369 !
370  IF (ivac .ge. 1) THEN
371  armn_e(ns:nrzt:ns) = armn_e(ns:nrzt:ns)
372  & + zu0(ns:nrzt:ns)*rbsq(1:nznt)
373  armn_o(ns:nrzt:ns) = armn_o(ns:nrzt:ns)
374  & + zu0(ns:nrzt:ns)*rbsq(1:nznt)
375  azmn_e(ns:nrzt:ns) = azmn_e(ns:nrzt:ns)
376  & - ru0(ns:nrzt:ns)*rbsq(1:nznt)
377  azmn_o(ns:nrzt:ns) = azmn_o(ns:nrzt:ns)
378  & - ru0(ns:nrzt:ns)*rbsq(1:nznt)
379 ! fz00_edge = SUM(wint(ns:nrzt:ns)*ru0(ns:nrzt:ns)*rbsq(1:nznt))
380  ENDIF
381 
382  100 CONTINUE
383 ! CALL Print1DArrayMNSP (gcon,tlglob,trglob,500,.TRUE.,"prcon")
384 
385 !
386 ! COMPUTE CONSTRAINT FORCE KERNELS
387 !
388 #ifndef _HBANGLE
389  rcon(:nrzt,0) = (rcon(:nrzt,0) - rcon0(:nrzt)) * gcon(:nrzt)
390  zcon(:nrzt,0) = (zcon(:nrzt,0) - zcon0(:nrzt)) * gcon(:nrzt)
391  brmn_e(:nrzt) = brmn_e(:nrzt) + rcon(:nrzt,0)
392  bzmn_e(:nrzt) = bzmn_e(:nrzt) + zcon(:nrzt,0)
393  brmn_o(:nrzt) = brmn_o(:nrzt)+ rcon(:nrzt,0)*sqrts(:nrzt)
394  bzmn_o(:nrzt) = bzmn_o(:nrzt)+ zcon(:nrzt,0)*sqrts(:nrzt)
395  rcon(:nrzt,0) = ru0(:nrzt) * gcon(:nrzt)
396  zcon(:nrzt,0) = zu0(:nrzt) * gcon(:nrzt)
397  rcon(:nrzt,1) = rcon(:nrzt,0) * sqrts(:nrzt)
398  zcon(:nrzt,1) = zcon(:nrzt,0) * sqrts(:nrzt)
399 #endif
400 
401  CALL second0 (tforoff)
402  timer(tfor) = timer(tfor) + (tforoff - tforon)
403 
404  END SUBROUTINE forces