V3FIT
bcovar.f
1  SUBROUTINE bcovar_par(lu, lv, tpxc, ier_flag)
2  USE vmec_main, fpsi => bvco, p5 => cp5
3  USE vmec_params, ONLY: ns4, signgs, pdamp, lamscale, ntmax,
4  & bsub_bad_js1_flag, arz_bad_value_flag,
5  & norm_term_flag
6  USE realspace, ONLY: pextra1, pextra2, pextra3, pextra4,
7  & pguu, pguv, pgvv, pru, pzu,
8  & pr1, prv, pzv, pshalf, pwint, pz1,
9  & pru0, pzu0, psqrts
10  USE vforces, r12 => parmn_o, ru12 => pazmn_e, gsqrt => pazmn_o,
11  & rs => pbzmn_e, zs => pbrmn_e, zu12 => parmn_e,
12  & bsubu_e => pclmn_e, bsubv_e => pblmn_e,
13  & bsubu_o => pclmn_o, bsubv_o => pblmn_o,
14  & bsq => pbzmn_o, phipog => pbrmn_o
15  USE xstuff, ONLY: pxc
16  USE precon2d, ONLY: ictrl_prec2d, lhess_exact,
17  & ctor_prec2d
18  USE fbal
19  USE vmec_input, ONLY: nzeta
20  USE vmec_dim, ONLY: ntheta3
21  USE parallel_include_module
22  IMPLICIT NONE
23 !-----------------------------------------------
24 ! D u m m y A r g u m e n t s
25 !-----------------------------------------------
26  REAL(dp), DIMENSION(nznt,ns,0:1), INTENT(INOUT) :: lu, lv
27  REAL(dp), DIMENSION((1+ntor)*(1+mpol1),1:ns,1:2*ntmax),
28  & INTENT(IN) :: tpxc
29  INTEGER, INTENT(inout) :: ier_flag
30 !-----------------------------------------------
31 ! L o c a l P a r a m e t e r s
32 !-----------------------------------------------
33 ! GENERALLY, IF TEMPORAL CONVERGENCE IS POOR, TRY TO INCREASE PDAMP (< 1)
34 ! (STORED IN VMEC_PARAMS)
35  REAL(dp), PARAMETER :: c1p5 = (one + p5)
36 !-----------------------------------------------
37 ! L o c a l V a r i a b l e s
38 !-----------------------------------------------
39  INTEGER :: l, js, ndim
40  REAL(dp) :: r2, volume, curpol_temp
41 #ifndef _HBANGLE
42  REAL(dp) :: arnorm, aznorm, tcon_mul
43 #endif
44  REAL(dp) :: bcastton, bcasttoff
45  REAL(dp), POINTER, DIMENSION(:,:) :: luu, luv, lvv, tau
46  REAL(dp), DIMENSION(:,:), POINTER :: bsupu, bsubuh,
47  & bsupv, bsubvh, r12sq
48  LOGICAL :: lctor
49  INTEGER :: i, j, k, nsmin, nsmax, istat
50  REAL(dp) :: wblocal(ns), wbtotal
51  REAL(dp) :: wplocal(ns), wptotal
52  REAL(dp) :: vptotal
53  REAL(dp) :: fnlocal(ns), fntotal
54  REAL(dp) :: fn1local(ns), fn1total
55  REAL(dp) :: fnLlocal(ns), fnLtotal
56 !-----------------------------------------------
57  IF (irst.EQ.2 .AND. iequi.EQ.0) RETURN
58 
59 !
60 ! POINTER ALIAS ASSIGNMENTS
61 
62  tau => pextra1(:,:,1)
63  luu => pextra2(:,:,1)
64  luv => pextra3(:,:,1)
65  lvv => pextra4(:,:,1)
66 
67  bsupu => luu
68  bsubuh => bsubu_o
69  bsupv => luv
70  bsubvh => bsubv_o
71  r12sq => bsq
72 
73 !
74 ! FOR OPTIMIZATION ON CRAY, MUST USE COMPILER DIRECTIVES TO
75 ! GET VECTORIZATION OF LOOPS INVOLVING (MORE THAN ONE) POINTER!
76 !
77  nsmin=t1lglob; nsmax=t1rglob
78  pguu(:,nsmin:nsmax) = 0
79  pguv(:,nsmin:nsmax) = 0
80  pgvv(:,nsmin:nsmax) = 0
81 
82 !
83 ! COMPUTE METRIC ELEMENTS GIJ ON HALF MESH
84 ! FIRST, GIJ = EVEN PART (ON FULL MESH), LIJ = ODD PART (ON FULL MESH)
85 ! THEN, GIJ(HALF) = < GIJ(even)> + SHALF < GIJ(odd) >
86 
87  DO l = nsmin, nsmax
88  r12sq(:,l) = psqrts(:,l)*psqrts(:,l)
89  pguu(:,l) = pru(:,l,0)*pru(:,l,0) + pzu(:,l,0)*pzu(:,l,0)
90  & + r12sq(:,l)*(pru(:,l,1)*pru(:,l,1)
91  & + pzu(:,l,1)*pzu(:,l,1))
92  luu(:,l) = (pru(:,l,0)*pru(:,l,1) + pzu(:,l,0)*pzu(:,l,1))*2
93  phipog(:,l) = 2*pr1(:,l,0)*pr1(:,l,1)
94  END DO
95 
96  IF (lthreed) THEN
97  DO l = nsmin, nsmax
98  pguv(:,l) = pru(:,l,0) * prv(:,l,0) + pzu(:,l,0)*pzv(:,l,0)
99  & + r12sq(:,l) * (pru(:,l,1)*prv(:,l,1)
100  & + pzu(:,l,1)*pzv(:,l,1))
101  luv(:,l) = pru(:,l,0) * prv(:,l,1) + pru(:,l,1)*prv(:,l,0)
102  & + pzu(:,l,0)*pzv(:,l,1) + pzu(:,l,1)*pzv(:,l,0)
103  pgvv(:,l) = prv(:,l,0) * prv(:,l,0) + pzv(:,l,0)*pzv(:,l,0)
104  & + r12sq(:,l) * (prv(:,l,1)*prv(:,l,1)
105  & + pzv(:,l,1)*pzv(:,l,1) )
106  lvv(:,l) = (prv(:,l,0) * prv(:,l,1) +
107  & pzv(:,l,0)*pzv(:,l,1))*2
108  END DO
109  END IF
110 
111  r12sq(:,nsmin:nsmax) = pr1(:,nsmin:nsmax,0)*pr1(:,nsmin:nsmax,0)
112  & + r12sq(:,nsmin:nsmax)*pr1(:,nsmin:nsmax,1)*
113  & pr1(:,nsmin:nsmax,1)
114 
115  DO l = t1rglob, max(t1lglob,2), -1
116  pguu(:,l) = p5*(pguu(:,l) + pguu(:,l-1) +
117  & + pshalf(:,l)*(luu(:,l) + luu(:,l-1)))
118  r12sq(:,l) = p5*(r12sq(:,l)+r12sq(:,l-1)+pshalf(:,l)* !Comment: r12sq = r12**2
119  & (phipog(:,l) + phipog(:,l-1)))
120  END DO
121 
122  IF (lthreed) THEN
123  DO l = t1rglob, max(t1lglob,2), -1
124  pguv(:,l) = p5*(pguv(:,l) + pguv(:,l-1) +
125  & pshalf(:,l)*(luv(:,l) + luv(:,l-1)))
126  pgvv(:,l) = p5*(pgvv(:,l) + pgvv(:,l-1) +
127  & pshalf(:,l)*(lvv(:,l) + lvv(:,l-1)))
128  END DO
129  END IF
130 
131  pguv(:,1)=0
132  pgvv(:,1)=0
133 
134  nsmin = tlglob; nsmax = t1rglob
135  DO l = nsmin, nsmax
136  tau(:,l) = gsqrt(:,l)
137  gsqrt(:,l) = r12(:,l)*tau(:,l)
138  END DO
139  gsqrt(:,1) = gsqrt(:,2)
140 
141  nsmin = max(2,tlglob); nsmax = t1rglob
142  pgvv(:,nsmin:nsmax) = pgvv(:,nsmin:nsmax)
143  & + r12sq(:,nsmin:nsmax)
144  pgvv(:,1) = 0
145 
146 !CATCH THIS AFTER WHERE LINE BELOW phipog = 0
147  nsmin = max(2,tlglob); nsmax = t1rglob
148  WHERE (gsqrt(:,nsmin:nsmax) .ne. zero)
149  phipog(:,nsmin:nsmax) = one/gsqrt(:,nsmin:nsmax)
150  END WHERE
151  phipog(:,1) = 0
152 
153  vp(1) = 0
154  vp(ns+1) = 0
155  DO js = nsmin, nsmax
156  vp(js) = signgs*sum(gsqrt(:,js)*pwint(:,js))
157  END DO
158 
159 !
160 ! COMPUTE CONTRA-VARIANT COMPONENTS OF B (Bsupu,v) ON RADIAL HALF-MESH
161 ! TO ACCOMODATE LRFP=T CASES, THE OVERALL PHIP FACTOR (PRIOR TO v8.46)
162 ! HAS BEEN REMOVED FROM PHIPOG, SO NOW PHIPOG == 1/GSQRT!
163 !
164 ! NOTE: LU = LAMU == d(LAM)/du, LV = -LAMV == -d(LAM)/dv COMING INTO THIS ROUTINE
165 ! WILL ADD CHIP IN CALL TO ADD_FLUXES. THE NET BSUPU, BSUPV ARE (PHIPOG=1/GSQRT AS NOTED ABOVE):
166 !
167 ! BSUPU = PHIPOG*(chip + LAMV*LAMSCALE),
168 ! BSUPV = PHIPOG*(phip + LAMU*LAMSCALE)
169 !
170 
171  nsmin = t1lglob
172  nsmax = t1rglob
173 
174  DO l = nsmin, nsmax
175  lu(:,l,:) = lu(:,l,:)*lamscale
176  lv(:,l,:) = lv(:,l,:)*lamscale
177  lu(:,l,0) = lu(:,l,0) + phipf(l)
178  END DO
179 
180  nsmin = max(2,t1lglob)
181  nsmax = t1rglob
182  DO l = nsmin, nsmax
183  bsupu(:,l) = p5*phipog(:,l) * (lv(:,l,0) + lv(:,l-1,0)
184  & + pshalf(:,l)*(lv(:,l,1) + lv(:,l-1,1)))
185  bsupv(:,l) = p5*phipog(:,l) * (lu(:,l,0) + lu(:,l-1,0)
186  & + pshalf(:,l)*(lu(:,l,1) + lu(:,l-1,1)))
187  END DO
188 
189 !v8.49: add ndim points
190  IF (rank .EQ. 0) THEN
191  bsupu(:,1) =0; ! bsupu(ndim) = 0
192  bsupv(:,1) =0; ! bsupv(ndim) = 0
193  END IF
194 
195 !
196 ! UPDATE IOTA EITHER OF TWO WAYS:
197 ! 1) FOR ictrl_prec2d = 0, SOLVE THE LINEAR ALGEBRAIC EQUATION <Bsubu> = icurv
198 ! FOR iotas (after testing, this is preferred way)
199 ! 2) FOR ictrl_prec2d > 0, EVOLVE IOTAS IN TIME, USING Force-iota = <Bsubu> - icurv.
200 ! IOTAS IS "STORED" AT LOCATION LAMBDA-SC(0,0) IN XC-ARRAY
201 !
202 
203 ! COMPUTE (IF NEEDED) AND ADD CHIP TO BSUPU
204 #if defined(CHI_FORCE)
205  CALL add_fluxes_par(phipog, bsupu, bsupv, ictrl_prec2d.EQ.0)
206 #else
207  CALL add_fluxes_par(phipog, bsupu, bsupv, .true.)
208 #endif
209 
210 !
211 ! COMPUTE LAMBDA FORCE KERNELS (COVARIANT B COMPONENT bsubu,v) ON RADIAL HALF-MESH
212 !
213  nsmin = t1lglob
214  nsmax = t1rglob
215  DO l = nsmin, nsmax
216  bsubuh(:,l) = pguu(:,l)*bsupu(:,l) + pguv(:,l)*bsupv(:,l)
217  bsubvh(:,l) = pguv(:,l)*bsupu(:,l) + pgvv(:,l)*bsupv(:,l)
218  END DO
219 
220 !v8.49
221 !
222 ! COMPUTE MAGNETIC AND KINETIC PRESSURE ON RADIAL HALF-MESH
223 !
224  nsmin = t1lglob
225  nsmax = t1rglob
226  DO l = nsmin, nsmax
227  bsq(:,l) = p5*(bsupu(:,l)*bsubuh(:,l) + bsupv(:,l)*bsubvh(:,l))
228  END DO
229 
230  nsmin = max(2,tlglob)
231  nsmax = min(ns,t1rglob)
232  pres(nsmin:nsmax) = mass(nsmin:nsmax)/vp(nsmin:nsmax)**gamma
233  pres(1)=0
234 
235  IF (ictrl_prec2d .LE. 1) THEN
236  DO l = tlglob, trglob
237  wblocal(l) = sum(pwint(:,l)*gsqrt(:,l) * bsq(:,l))
238  wplocal(l) = vp(l)*pres(l)
239  END DO
240 
241  CALL gather1xarray(wblocal)
242  wbtotal = sum(wblocal(2:ns))
243  CALL gather1xarray(wplocal)
244  wptotal = sum(wplocal(2:ns))
245  wb = hs*abs(wbtotal)
246  wp = hs*wptotal
247  END IF
248 
249 ! ADD KINETIC PRESSURE TO MAGNETIC PRESSURE
250  nsmin = tlglob
251  nsmax = t1rglob
252  DO l=nsmin, nsmax
253  bsq(:,l) = bsq(:,l) + pres(l)
254  lvv(:,l) = phipog(:,l)*pgvv(:,l)
255  END DO
256 
257 !SPH122407-MOVED HERE: COMPUTE LAMBDA FULL MESH FORCES
258 ! NOTE: bsubu_e is used here ONLY as a temporary array
259 
260  nsmin = tlglob
261  nsmax = min(ns - 1, trglob)
262  bsubv_e(:,nsmin:nsmax) = p5*(lvv(:,nsmin:nsmax)+
263  & lvv(:,nsmin+1:nsmax+1))
264  & * lu(:,nsmin:nsmax,0)
265  bsubv_e(:,ns) = p5*lvv(:,ns)*lu(:,ns,0)
266 
267  nsmin = tlglob
268  nsmax = t1rglob
269  DO l = nsmin, nsmax
270  lvv(:,l) = lvv(:,l)*pshalf(:,l)
271  bsubu_e(:,l) = pguv(:,l)*bsupu(:,l) !*sigma_an(:nrzt) !sigma_an=1 isotropic
272  END DO
273 
274  nsmin = tlglob
275  nsmax = min(ns - 1, trglob)
276  DO l = nsmin, nsmax
277  bsubv_e(:,l) = bsubv_e(:,l)
278  & + p5*((lvv(:,l)+lvv(:,l+1))*lu(:,l,1) +
279  & bsubu_e(:,l) + bsubu_e(:,l+1))
280  END DO
281  bsubv_e(:,ns) = bsubv_e(:,ns)
282  & + p5*(lvv(:,ns)*lu(:,ns,1) + bsubu_e(:,ns))
283 
284 !
285 ! COMPUTE AVERAGE FORCE BALANCE AND TOROIDAL/POLOIDAL CURRENTS
286 !
287 !WAC: UPDATE buco, bvco AFTER pressure called (Gather buco, bvco in calc_fbal_par
288  CALL calc_fbal_par(bsubuh, bsubvh)
289  rbtor0= c1p5*fpsi(2) - p5*fpsi(3)
290  rbtor = c1p5*fpsi(ns) - p5*fpsi(ns-1)
291 !
292 ! (SPH:08/19/04)
293 ! MUST AVOID BREAKING TRI-DIAGONAL RADIAL COUPLING AT EDGE WHEN USING PRECONDITIONER
294 ! CTOR IS PASSED TO VACUUM TO COMPUTE EDGE BSQVAC, SO IT CAN ONLY DEPEND ON NS, NS-1
295 ! THUS, CTOR ~ buco(ns) WORKS, WITH REMAINDER A FIXED CONSTANT.
296 !
297 ! ALSO, IF USING FAST SWEEP IN COMPUTE_BLOCKS, MUST MAKE CTOR CONSTANT
298 ! TO AVOID BREAKING SYMMETRY OF A+(ns-1) AND B-(ns) HESSIAN ELEMENTS
299 !
300 ! TO GET CORRECT HESSIAN, USE THE CTOR=ctor_prec2d +... ASSIGNMENT
301 ! FOR ictrl_prec2d.ne.0 (replace ictrl_prec2d.gt.1 with ictrl_prec2d.ne.0 in IF test below)
302 !
303 !
304 
305 ! NEXT COMPUTE COVARIANT BSUBV COMPONENT ~ lvv ON FULL RADIAL MESH BY AVERAGING HALF-MESH METRICS
306 ! NOTE: EDGE VALUES AT JS=NS DOWN BY 1/2
307 ! THIS IS NEEDED FOR NUMERICAL STABILITY
308 
309  IF (lhess_exact) THEN
310  lctor = lfreeb .AND. ictrl_prec2d.NE.0 !Yields correct hessian near edge
311  ELSE
312  lctor = lfreeb .AND. ictrl_prec2d.GT.1 !Yields better accuracy in solution
313  END IF
314 
315  IF (lctor) THEN
316  IF (ictrl_prec2d .EQ. 2) THEN
317  ctor_prec2d = p5*(buco(ns) - buco(ns1))
318  END IF
319  ctor = signgs*twopi*(buco(ns)+ctor_prec2d)
320  ELSE
321  ctor = signgs*twopi*(c1p5*buco(ns) - p5*buco(ns1))
322  END IF
323 
324 !
325 ! AVERAGE LAMBDA FORCES ONTO FULL RADIAL MESH
326 ! USE BLENDING FOR bsubv_e FOR NUMERICAL STABILITY NEAR AXIS
327 !
328  nsmin = tlglob
329  nsmax = t1rglob
330  DO l = nsmin, nsmax
331  lvv(:,l) = bdamp(l)
332  END DO
333 
334  IF (rank.EQ.0) THEN
335  IF (any(bsubvh(:,1) .ne. zero)) ier_flag = bsub_bad_js1_flag
336  IF (any(bsubuh(:,1) .ne. zero)) ier_flag = bsub_bad_js1_flag
337  END IF
338 
339  nsmin = tlglob
340  nsmax = min(trglob,ns - 1)
341  bsubu_e(:,nsmin:nsmax) = p5*(bsubuh(:,nsmin:nsmax) +
342  & bsubuh(:,nsmin+1:nsmax+1))
343  IF (trglob .EQ. ns) bsubu_e(:,ns) = p5*bsubuh(:,ns)
344 
345  nsmin = tlglob
346  nsmax = min(ns - 1,trglob)
347  bsubv_e(:,nsmin:nsmax) =
348  & bsubv_e(:,nsmin:nsmax)*lvv(:,nsmin:nsmax) +
349  & p5*(1 - lvv(:,nsmin:nsmax))*(bsubvh(:,nsmin:nsmax) +
350  & bsubvh(:,nsmin+1:nsmax+1))
351  IF (trglob .EQ. ns) THEN
352  bsubv_e(:,ns) = bsubv_e(:,ns)*lvv(:,ns)
353  & + p5*(1-lvv(:,ns))*bsubvh(:,ns)
354  END IF
355 
356 !
357 ! COMPUTE R,Z AND LAMBDA PRE-CONDITIONING MATRIX
358 ! ELEMENTS AND FORCE NORMS: NOTE THAT lu=>czmn, lv=>crmn externally
359 ! SO THIS STORES bsupv in czmn_e, bsupu in crmn_e
360 !
361  nsmin = tlglob
362  nsmax = t1rglob
363  IF (iequi .EQ. 1) THEN
364  lu(:,nsmin:nsmax,0) = bsupv(:,nsmin:nsmax)
365  lv(:,nsmin:nsmax,0) = bsupu(:,nsmin:nsmax)
366  END IF
367 
368 !
369 ! COMPUTE PRECONDITIONING (1D) AND SCALING PARAMETERS
370 ! NO NEED TO RECOMPUTE WHEN 2D-PRECONDITIONER ON
371 !
372 
373  IF ((mod(iter2-iter1,ns4).EQ.0 .AND. iequi.EQ.0) .AND.
374  & ictrl_prec2d.EQ.0) THEN
375  nsmin = tlglob
376  nsmax = t1rglob
377  phipog(:,nsmin:nsmax) = phipog(:,nsmin:nsmax)
378  & * pwint(:,nsmin:nsmax)
379 
380  CALL lamcal_par(phipog, pguu, pguv, pgvv)
381 
382  CALL precondn_par(bsupv, bsq, gsqrt, r12, zs, zu12,
383  & pzu(:,:,0), pzu(:,:,1), pz1(:,:,1), arm,
384  & ard, brm, brd, crd, rzu_fac, cos01)
385 
386  CALL precondn_par(bsupv, bsq, gsqrt, r12, rs, ru12,
387  & pru(:,:,0), pru(:,:,1), pr1(:,:,1), azm,
388  & azd, bzm, bzd, crd, rru_fac, sin01)
389 
390  nsmin = max(2,tlglob)
391  nsmax = min(trglob,ns - 1)
392  rzu_fac(nsmin:nsmax) = psqrts(1,nsmin:nsmax)
393  & * rzu_fac(nsmin:nsmax)
394  rru_fac(nsmin:nsmax) = psqrts(1,nsmin:nsmax)
395  & * rru_fac(nsmin:nsmax)
396  frcc_fac(nsmin:nsmax) = one/rzu_fac(nsmin:nsmax)
397  rzu_fac(nsmin:nsmax) = rzu_fac(nsmin:nsmax)/2
398  fzsc_fac(nsmin:nsmax) =-one/rru_fac(nsmin:nsmax)
399  rru_fac(nsmin:nsmax) = rru_fac(nsmin:nsmax)/2
400 
401  nsmin = tlglob
402  nsmax = t1rglob
403  pguu(:,nsmin:nsmax) = pguu(:,nsmin:nsmax)
404  & * r12(:,nsmin:nsmax)**2
405 
406  DO l = max(2,tlglob), trglob
407  fnlocal(l) = sum(pguu(:,l)*pwint(:,l))
408  fn1local(l) = sum(tpxc(2:,l,1:ntmax)**2)
409  & + sum(tpxc(1:,l,ntmax+1:2*ntmax)**2)
410  fnllocal(l) = sum((bsubuh(:,l)**2 + bsubvh(:,l)**2)
411  & * pwint(:,l))*lamscale**2
412  END DO
413 
414  CALL gather1xarray(vp); vptotal = sum(vp(2:ns))
415  CALL gather1xarray(fnlocal); fntotal = sum(fnlocal(2:ns))
416  CALL gather1xarray(fn1local); fn1total= sum(fn1local(2:ns))
417  CALL gather1xarray(fnllocal); fnltotal= sum(fnllocal(2:ns))
418 
419  volume = hs*vptotal
420  r2 = max(wb,wp)/volume
421  fnorm = one/(fntotal*(r2*r2))
422  fnorm1=one/fn1total
423  fnorml = one/fnltotal
424 !
425 ! COMPUTE CONSTRAINT FORCE SCALING FACTOR (TCON)
426 ! OVERRIDE USER INPUT VALUE HERE
427 !
428 #ifndef _HBANGLE
429  r2 = ns
430  tcon0 = min(abs(tcon0), one) !!ignore large tcon0 from old-style files
431  tcon_mul = tcon0*(1 + r2*(one/60 + r2/(200*120)))
432 
433  tcon_mul = tcon_mul/((4*r0scale**2)**2) !!Scaling of ard, azd (2*r0scale**2);
434  !!Scaling of cos**2 in alias (4*r0scale**2)
435  tcon = tcon0
436  DO js = max(2,tlglob), min(ns-1,trglob)
437  arnorm = sum(pwint(:,js)*pru0(:,js)**2)
438  aznorm = sum(pwint(:,js)*pzu0(:,js)**2)
439 ! IF (arnorm .eq. zero .or. aznorm .eq. zero) THEN
440 ! STOP 'arnorm or aznorm=0 in bcovar'
441 ! END IF
442  IF (arnorm .eq. zero .or. aznorm .eq. zero) THEN !SAL 070719
443  ier_flag = arz_bad_value_flag
444  END IF
445 
446  tcon(js) = min(abs(ard(js,1)/arnorm),
447  1 abs(azd(js,1)/aznorm))*tcon_mul*(32*hs)**2
448  END DO
449  tcon(ns) = p5*tcon(ns - 1)
450  IF (lasym) THEN
451  tcon = p5*tcon
452  END IF
453 #endif
454  ENDIF
455 
456  CALL mpi_allreduce(mpi_in_place,ier_flag,1,mpi_integer,
457  1 mpi_max,ns_comm,mpi_err)
458  IF (ier_flag .ne. norm_term_flag) RETURN
459 !
460 ! COMPUTE COVARIANT BSUBU,V (EVEN, ODD) ON HALF RADIAL MESH
461 ! FOR FORCE BALANCE AND RETURN (IEQUI=1)
462 !
463 !
464 ! COMPUTE COVARIANT BSUBU,V (EVEN, ODD) ON HALF RADIAL MESH
465 ! FOR FORCE BALANCE AND RETURN (IEQUI=1)
466 !
467 
468  IF (iequi .EQ. 1) THEN
469  nsmin = max(tlglob,2)
470  nsmax = min(trglob,ns - 1)
471  DO js = nsmax, nsmin, -1
472  bsubvh(:,js) = 2*bsubv_e(:,js) - bsubvh(:,js+1)
473  END DO
474 
475 
476 ! ADJUST <bsubvh> AFTER MESH-BLENDING
477  nsmin=max(tlglob,2); nsmax=min(trglob,ns)
478  DO js = nsmin, nsmax
479  curpol_temp = fpsi(js) - sum(bsubvh(:,js)*pwint(:,js))
480  bsubvh(:,js) = bsubvh(:,js) + curpol_temp
481  END DO
482 
483  bsubu_e(:,nsmin:nsmax) = bsubuh(:,nsmin:nsmax)
484  bsubv_e(:,nsmin:nsmax) = bsubvh(:,nsmin:nsmax)
485 
486  bsubu_o(:,nsmin:nsmax) = pshalf(:,nsmin:nsmax)
487  & * bsubu_e(:,nsmin:nsmax)
488  bsubv_o(:,nsmin:nsmax) = pshalf(:,nsmin:nsmax)
489  & * bsubv_e(:,nsmin:nsmax)
490  RETURN
491  END IF
492 
493 ! MINUS SIGN => HESSIAN DIAGONALS ARE POSITIVE
494  nsmin = max(tlglob,2)
495  nsmax = trglob
496  bsubu_e(:,nsmin:nsmax) = -lamscale*bsubu_e(:,nsmin:nsmax)
497  bsubv_e(:,nsmin:nsmax) = -lamscale*bsubv_e(:,nsmin:nsmax)
498  bsubu_o(:,nsmin:nsmax) = psqrts(:,nsmin:nsmax)
499  & * bsubu_e(:,nsmin:nsmax)
500  bsubv_o(:,nsmin:nsmax) = psqrts(:,nsmin:nsmax)
501  & * bsubv_e(:,nsmin:nsmax)
502 
503 !
504 ! STORE LU * LV COMBINATIONS USED IN FORCES
505 !
506 
507  nsmin = max(tlglob,2)
508  nsmax = t1rglob
509  DO l = nsmin, nsmax
510  lvv(:,l) = gsqrt(:,l) !*sigma_an(:,l)
511  pguu(:,l) = bsupu(:,l)*bsupu(:,l)*lvv(:,l)
512  pguv(:,l) = bsupu(:,l)*bsupv(:,l)*lvv(:,l)
513  pgvv(:,l) = bsupv(:,l)*bsupv(:,l)*lvv(:,l)
514  lv(:,l,0) = bsq(:,l)*tau(:,l)
515  lu(:,l,0) = bsq(:,l)*r12(:,l)
516  END DO
517 
518  END SUBROUTINE bcovar_par
519 
520  SUBROUTINE bcovar (lu, lv)
521  USE vmec_main, fpsi => bvco, p5 => cp5
522  USE vmec_params, ONLY: ns4, signgs, pdamp, lamscale
523  USE realspace
524  USE vforces, r12 => armn_o, ru12 => azmn_e, gsqrt => azmn_o,
525  & rs => bzmn_e, zs => brmn_e, zu12 => armn_e,
526  & bsubu_e => clmn_e, bsubv_e => blmn_e,
527  & bsubu_o => clmn_o, bsubv_o => blmn_o,
528  & bsq => bzmn_o, phipog => brmn_o
529  USE xstuff, ONLY: xc
530  USE precon2d, ONLY: ictrl_prec2d, lhess_exact,
531  & ctor_prec2d
532 #ifdef _HBANGLE
533  USE angle_constraints, ONLY: precondn_rho, ard2, arm2,
534  & azd2, azm2, brd2, brm2, bzd2, bzm2
535 #endif
536  USE fbal
537  IMPLICIT NONE
538 !-----------------------------------------------
539 ! D u m m y A r g u m e n t s
540 !-----------------------------------------------
541  REAL(dp), DIMENSION(nrzt,0:1), INTENT(inout) :: lu, lv
542 !-----------------------------------------------
543 ! L o c a l P a r a m e t e r s
544 !-----------------------------------------------
545 ! GENERALLY, IF TEMPORAL CONVERGENCE IS POOR, TRY TO INCREASE PDAMP (< 1)
546 ! (STORED IN VMEC_PARAMS)
547  REAL(dp), PARAMETER :: c1p5 = (one + p5)
548 !-----------------------------------------------
549 ! L o c a l V a r i a b l e s
550 !-----------------------------------------------
551  INTEGER :: l, js, ndim
552  REAL(dp) :: r2, volume, curpol_temp
553 #ifndef _HBANGLE
554  REAL(dp) :: arnorm, aznorm, tcon_mul
555 #endif
556  REAL(dp), POINTER, DIMENSION(:) :: luu, luv, lvv, tau
557  REAL(dp), DIMENSION(:), POINTER :: bsupu, bsubuh,
558  & bsupv, bsubvh, r12sq
559  LOGICAL :: lctor
560 !-----------------------------------------------
561  ndim = 1 + nrzt
562 
563 !
564 ! POINTER ALIAS ASSIGNMENTS
565 !
566  tau => extra1(:,1); luu => extra2(:,1);
567  luv => extra3(:,1); lvv => extra4(:,1)
568 
569  bsupu => luu; bsubuh => bsubu_o
570  bsupv => luv; bsubvh => bsubv_o
571  r12sq => bsq
572 
573 
574 !
575 ! FOR OPTIMIZATION ON CRAY, MUST USE COMPILER DIRECTIVES TO
576 ! GET VECTORIZATION OF LOOPS INVOLVING (MORE THAN ONE) POINTER!
577 !
578  guu(ndim) = 0; guv = 0; gvv = 0
579 
580 !
581 ! COMPUTE METRIC ELEMENTS GIJ ON HALF MESH
582 ! FIRST, GIJ = EVEN PART (ON FULL MESH), LIJ = ODD PART (ON FULL MESH)
583 ! THEN, GIJ(HALF) = < GIJ(even)> + SHALF < GIJ(odd) >
584 !
585 
586  r12sq(1:nrzt) = sqrts(1:nrzt)*sqrts(1:nrzt)
587  guu(1:nrzt) = ru(1:nrzt,0)*ru(1:nrzt,0)
588  & + zu(1:nrzt,0)*zu(1:nrzt,0) + r12sq(1:nrzt)*
589  & ( ru(1:nrzt,1)*ru(1:nrzt,1)
590  & + zu(1:nrzt,1)*zu(1:nrzt,1))
591 
592  luu(1:nrzt) = (ru(1:nrzt,0)*ru(1:nrzt,1)
593  & + zu(1:nrzt,0)*zu(1:nrzt,1))*2
594  phipog(1:nrzt) = 2*r1(1:nrzt,0)*r1(1:nrzt,1)
595 
596  IF (lthreed) THEN
597  guv(1:nrzt) = ru(1:nrzt,0)*rv(1:nrzt,0)
598  & + zu(1:nrzt,0)*zv(1:nrzt,0)
599  & + r12sq(1:nrzt)*(ru(1:nrzt,1)*rv(1:nrzt,1) +
600  & zu(1:nrzt,1)*zv(1:nrzt,1))
601  luv(1:nrzt) = ru(1:nrzt,0)*rv(1:nrzt,1)
602  & + ru(1:nrzt,1)*rv(1:nrzt,0)
603  & + zu(1:nrzt,0)*zv(1:nrzt,1)
604  & + zu(1:nrzt,1)*zv(1:nrzt,0)
605  gvv(1:nrzt) = rv(1:nrzt,0)*rv(1:nrzt,0)
606  & + zv(1:nrzt,0)*zv(1:nrzt,0)
607  & + r12sq(1:nrzt)*(rv(1:nrzt,1)*rv(1:nrzt,1) +
608  & zv(1:nrzt,1)*zv(1:nrzt,1))
609  lvv(1:nrzt) =(rv(1:nrzt,0)*rv(1:nrzt,1) +
610  & zv(1:nrzt,0)*zv(1:nrzt,1))*2
611  END IF
612 
613  r12sq(1:nrzt) = r1(1:nrzt,0)*r1(1:nrzt,0)
614  & + r12sq(1:nrzt)*r1(1:nrzt,1)*r1(1:nrzt,1)
615 
616 !DIR$ IVDEP
617  DO l = nrzt, 2, -1
618  guu(l) = p5*(guu(l) + guu(l-1) + shalf(l)*(luu(l) + luu(l-1)))
619  r12sq(l) = p5*(r12sq(l) + r12sq(l-1) + !Comment: r12sq = r12**2
620  & shalf(l)*(phipog(l) + phipog(l-1)))
621  END DO
622 
623  IF (lthreed) THEN
624 !DIR$ IVDEP
625  DO l = nrzt, 2, -1
626  guv(l) = p5*(guv(l) + guv(l-1) +
627  & shalf(l)*(luv(l) + luv(l-1)))
628  gvv(l) = p5*(gvv(l) + gvv(l-1) +
629  & shalf(l)*(lvv(l) + lvv(l-1)))
630  END DO
631  END IF
632 
633  tau(1:nrzt) = gsqrt(1:nrzt)
634  gsqrt(1:nrzt) = r12(1:nrzt)*tau(1:nrzt)
635  gsqrt(1:nrzt:ns) = gsqrt(2:nrzt:ns)
636 
637  gvv(2:nrzt) = gvv(2:nrzt) + r12sq(2:nrzt)
638 
639 !CATCH THIS AFTER WHERE LINE BELOW phipog = 0
640  WHERE (gsqrt(2:ndim) .ne. zero)
641  phipog(2:ndim) = one/gsqrt(2:ndim)
642  END WHERE
643  phipog(1:ndim:ns) = 0
644 
645  vp(1) = 0; vp(ns+1) = 0
646  DO js = 2, ns
647  vp(js) = signgs*sum(gsqrt(js:nrzt:ns)*wint(js:nrzt:ns))
648  END DO
649  IF (iter2 .eq. 1) THEN
650  voli = twopi*twopi*hs*sum(vp(2:ns))
651  END IF
652 
653 !
654 ! COMPUTE CONTRA-VARIANT COMPONENTS OF B (Bsupu,v) ON RADIAL HALF-MESH
655 ! TO ACCOMODATE LRFP=T CASES, THE OVERALL PHIP FACTOR (PRIOR TO v8.46)
656 ! HAS BEEN REMOVED FROM PHIPOG, SO NOW PHIPOG == 1/GSQRT!
657 !
658 ! NOTE: LU = LAMU == d(LAM)/du, LV = -LAMV == -d(LAM)/dv COMING INTO THIS ROUTINE
659 ! WILL ADD CHIP IN CALL TO ADD_FLUXES. THE NET BSUPU, BSUPV ARE (PHIPOG=1/GSQRT AS NOTED ABOVE):
660 !
661 ! BSUPU = PHIPOG*(chip - LAMV*LAMSCALE),
662 ! BSUPV = PHIPOG*(phip + LAMU*LAMSCALE)
663 !
664  lu = lu*lamscale
665  lv = lv*lamscale
666 
667  DO js = 1, ns
668  lu(js:nrzt:ns,0) = lu(js:nrzt:ns,0) + phipf(js)
669  END DO
670 
671  bsupu(2:nrzt) = p5*phipog(2:nrzt)*(lv(2:nrzt,0) + lv(1:nrzt-1,0)
672  & + shalf(2:nrzt)*(lv(2:nrzt,1) + lv(1:nrzt-1,1)))
673  bsupv(2:nrzt) = p5*phipog(2:nrzt)*(lu(2:nrzt,0) + lu(1:nrzt-1,0)
674  & + shalf(2:nrzt)*(lu(2:nrzt,1) + lu(1:nrzt-1,1)))
675 !v8.49: add ndim points
676  bsupu(1)=0; bsupu(ndim)=0
677  bsupv(1)=0; bsupv(ndim)=0
678 
679 !
680 ! UPDATE IOTA EITHER OF TWO WAYS:
681 ! 1) FOR ictrl_prec2d = 0, SOLVE THE LINEAR ALGEBRAIC EQUATION <Bsubu> = icurv
682 ! FOR iotas
683 ! 2) FOR ictrl_prec2d > 0, EVOLVE IOTAS IN TIME, USING Force-iota = <Bsubu> - icurv IN TOMNSP.
684 !
685 ! NEED TO DO IT WAY (#2) TO EASILY COMPUTE THE HESSIAN ELEMENTS DUE TO LAMBDA-VARIATIONS.
686 ! IOTAS IS "STORED" AT LOCATION LAMBDA-SC(0,0) IN XC-ARRAY [USE THIS COMPONENT SO IT
687 ! WILL WORK EVEN FOR 2D PLASMA], ALTHOUGH ITS VARIATION IS LIKE THAT OF LV-CS(0,0),
688 ! WITH N -> 1 IN THE HESSIAN CALCULATION ROUTINES (Compute_Hessian_Flam_lam, etc.)
689 !
690 
691 ! COMPUTE (IF NEEDED) AND ADD CHIP TO BSUPU
692 #if defined(CHI_FORCE)
693  CALL add_fluxes(phipog, bsupu, bsupv, ictrl_prec2d.EQ.0)
694 #else
695  CALL add_fluxes(phipog, bsupu, bsupv, .true.)
696 #endif
697 
698 !
699 ! COMPUTE LAMBDA FORCE KERNELS (COVARIANT B COMPONENT bsubu,v) ON RADIAL HALF-MESH
700 !
701  bsubuh(1:nrzt) = guu(1:nrzt)*bsupu(1:nrzt)
702  & + guv(1:nrzt)*bsupv(1:nrzt)
703  bsubvh(1:nrzt) = guv(1:nrzt)*bsupu(1:nrzt)
704  & + gvv(1:nrzt)*bsupv(1:nrzt)
705 !v8.49
706  bsubuh(ndim) = 0; bsubvh(ndim) = 0
707 
708 !
709 ! COMPUTE MAGNETIC AND KINETIC PRESSURE ON RADIAL HALF-MESH
710 !
711  bsq(:nrzt) = p5*(bsupu(:nrzt)*bsubuh(:nrzt)
712  & + bsupv(:nrzt)*bsubvh(:nrzt))
713 
714  wb = hs*abs(sum(wint(:nrzt)*gsqrt(:nrzt)*bsq(:nrzt)))
715 
716 #ifdef _ANIMEC
717 !SPH: MAKE CALL HERE (bsubX_e are used for scratch arrays)
718  CALL an_pressure(bsubu_e, bsubv_e)
719 
720 ! ADD KINETIC PRESSURE TO MAGNETIC PRESSURE
721  bsq(2:nrzt) = bsq(2:nrzt) + pperp(2:nrzt)
722 
723 !WAC-SPH: MODIFY EFFECTIVE CURRENT K = curl(sigma_an*B)
724  phipog(1:nrzt) = phipog(1:nrzt)*sigma_an(1:nrzt)
725  bsubuh(1:nrzt) = bsubuh(1:nrzt)*sigma_an(1:nrzt)
726  bsubvh(1:nrzt) = bsubvh(1:nrzt)*sigma_an(1:nrzt)
727 
728 #else
729  pres(2:ns) = mass(2:ns)/vp(2:ns)**gamma
730  wp = hs*sum(vp(2:ns)*pres(2:ns))
731 
732 ! ADD KINETIC PRESSURE TO MAGNETIC PRESSURE
733  DO js = 2, ns
734  bsq(js:nrzt:ns) = bsq(js:nrzt:ns) + pres(js)
735  END DO
736 #endif
737 
738 !SPH122407-MOVED HERE: COMPUTE LAMBDA FULL MESH FORCES
739 ! NOTE: bsubu_e is used here ONLY as a temporary array
740  lvv = phipog(:ndim)*gvv
741  bsubv_e(1:nrzt) = p5*(lvv(1:nrzt)+lvv(2:ndim))*lu(1:nrzt,0)
742 
743  lvv = lvv*shalf
744  bsubu_e(:nrzt) = guv(:nrzt)*bsupu(:nrzt)*sigma_an(:nrzt) !sigma_an=1 isotropic
745  bsubu_e(ndim) = 0
746  bsubv_e(1:nrzt) = bsubv_e(1:nrzt)
747  & + p5*((lvv(1:nrzt) + lvv(2:ndim))*lu(1:nrzt,1) +
748  & bsubu_e(1:nrzt) + bsubu_e(2:ndim))
749 
750 !
751 ! COMPUTE AVERAGE FORCE BALANCE AND TOROIDAL/POLOIDAL CURRENTS
752 !
753 !WAC: UPDATE buco, bvco AFTER pressure called
754 #ifdef _ANIMEC
755  IF (iequi .EQ. 1) papr = pmap*pres/vp
756 #endif
757  CALL calc_fbal(bsubuh, bsubvh)
758 
759  rbtor0= c1p5*fpsi(2) - p5*fpsi(3)
760  rbtor = c1p5*fpsi(ns) - p5*fpsi(ns-1)
761 !
762 ! (SPH:08/19/04)
763 ! MUST AVOID BREAKING TRI-DIAGONAL RADIAL COUPLING AT EDGE WHEN USING PRECONDITIONER
764 ! CTOR IS PASSED TO VACUUM TO COMPUTE EDGE BSQVAC, SO IT CAN ONLY DEPEND ON NS, NS-1
765 ! THUS, CTOR ~ buco(ns) WORKS, WITH REMAINDER A FIXED CONSTANT.
766 !
767 ! ALSO, IF USING FAST SWEEP IN COMPUTE_BLOCKS, MUST MAKE CTOR CONSTANT
768 ! TO AVOID BREAKING SYMMETRY OF A+(ns-1) AND B-(ns) HESSIAN ELEMENTS
769 !
770 ! TO GET CORRECT HESSIAN, USE THE CTOR=ctor_prec2d +... ASSIGNMENT
771 ! FOR ictrl_prec2d.ne.0 (replace ictrl_prec2d.gt.1 with ictrl_prec2d.ne.0 in IF test below)
772 !
773 !
774 
775 ! NEXT COMPUTE COVARIANT BSUBV COMPONENT ~ lvv ON FULL RADIAL MESH BY AVERAGING HALF-MESH METRICS
776 ! NOTE: EDGE VALUES AT JS=NS DOWN BY 1/2
777 ! THIS IS NEEDED FOR NUMERICAL STABILITY
778 
779  IF (lhess_exact) THEN
780  lctor = lfreeb .and. ictrl_prec2d.ne.0 !Yields correct hessian near edge
781  ELSE
782  lctor = lfreeb .and. ictrl_prec2d.gt.1 !Yields better accuracy in solution
783  END IF
784  IF (lctor) THEN
785  IF (ictrl_prec2d .eq. 2) THEN
786  ctor_prec2d = signgs*twopi*p5*(buco(ns) - buco(ns1))
787  END IF
788  ctor = ctor_prec2d + signgs*twopi*buco(ns)
789  ELSE
790  ctor = signgs*twopi*(c1p5*buco(ns) - p5*buco(ns1))
791  END IF
792 
793 !
794 ! AVERAGE LAMBDA FORCES ONTO FULL RADIAL MESH
795 ! USE BLENDING FOR bsubv_e FOR NUMERICAL STABILITY NEAR AXIS
796 !
797  DO l = 1, ns
798  lvv(l:nrzt:ns) = bdamp(l)
799  END DO
800 
801  IF (any(bsubuh(1:ndim:ns) .ne. zero)) stop 'BSUBUH != 0 AT JS=1'
802  IF (any(bsubvh(1:ndim:ns) .ne. zero)) stop 'BSUBVH != 0 AT JS=1'
803 
804  bsubu_e(1:nrzt) = p5*(bsubuh(1:nrzt) + bsubuh(2:ndim))
805  bsubv_e(1:nrzt) = bsubv_e(1:nrzt)*lvv(1:nrzt)
806  & + p5*(1 - lvv(1:nrzt))*(bsubvh(1:nrzt)
807  & + bsubvh(2:ndim))
808 
809 !
810 ! COMPUTE R,Z AND LAMBDA PRE-CONDITIONING MATRIX
811 ! ELEMENTS AND FORCE NORMS: NOTE THAT lu=>czmn, lv=>crmn externally
812 ! SO THIS STORES bsupv in czmn_e, bsupu in crmn_e
813 !
814  IF (iequi .EQ. 1) THEN
815  lu(:nrzt,0) = bsupv(:nrzt)
816  lv(:nrzt,0) = bsupu(:nrzt)
817  END IF
818 
819 !
820 ! COMPUTE PRECONDITIONING (1D) AND SCALING PARAMETERS
821 ! NO NEED TO RECOMPUTE WHEN 2D-PRECONDITIONER ON
822 !
823  IF ((mod(iter2-iter1,ns4) .eq. 0 .and. iequi .eq. 0) .and.
824  & ictrl_prec2d.eq.0) THEN
825  phipog(:nrzt) = phipog(:nrzt)*wint(:nrzt)
826  CALL lamcal(phipog, guu, guv, gvv)
827  CALL precondn(bsupv, bsq, gsqrt, r12, zs, zu12, zu,z u(1,1),
828  & z1(1,1), arm, ard, brm, brd,
829 #ifdef _HBANGLE
830  & arm2, ard2, brm2, brd2,
831 #endif
832  & crd, rzu_fac, cos01)
833  CALL precondn(bsupv, bsq, gsqrt, r12, rs, ru12, ru, ru(1,1),
834  & r1(1,1), azm, azd, bzm, bzd,
835 #ifdef _HBANGLE
836  & azm2, azd2, bzm2, bzd2,
837 #endif
838  & crd, rru_fac, sin01)
839 
840  rzu_fac(2:ns-1) = sqrts(2:ns-1)*rzu_fac(2:ns-1)
841  rru_fac(2:ns-1) = sqrts(2:ns-1)*rru_fac(2:ns-1)
842  frcc_fac(2:ns-1) = one/rzu_fac(2:ns-1); rzu_fac = rzu_fac/2
843  fzsc_fac(2:ns-1) =-one/rru_fac(2:ns-1); rru_fac = rru_fac/2
844 #ifdef _HBANGLE
845  CALL precondn_rho
846 #endif
847 
848  volume = hs*sum(vp(2:ns))
849  r2 = max(wb,wp)/volume
850  guu(:nrzt) = guu(:nrzt)*r12(:nrzt)**2 !R12 from RP in force
851  fnorm = one/(sum(guu(1:nrzt)*wint(1:nrzt))*(r2*r2)) !Norm, unpreconditioned R,Z forces
852  fnorm1 = one/sum(xc(1+ns:2*irzloff)**2) !Norm for preconditioned R,Z forces
853  fnorml = one/(sum((bsubuh(1:nrzt)**2 + bsubvh(1:nrzt)**2) *
854  & wint(1:nrzt))*lamscale**2) !Norm for unpreconditioned Lambda force
855 ! r3 = one/(2*r0scale)**2
856 ! fnorm2 = one/MAX(SUM(xc(2*irzloff+1:3*irzloff)**2),r3/4) !Norm for preconditioned Lambda force
857 
858 !
859 ! COMPUTE CONSTRAINT FORCE SCALING FACTOR (TCON)
860 ! OVERRIDE USER INPUT VALUE HERE
861 !
862 #ifndef _HBANGLE
863  r2 = ns
864  tcon0 = min(abs(tcon0), one) !!ignore large tcon0 from old-style files
865  tcon_mul = tcon0*(1 + r2*(one/60 + r2/(200*120)))
866 
867  tcon_mul = tcon_mul/((4*r0scale**2)**2) !!Scaling of ard, azd (2*r0scale**2);
868  !!Scaling of cos**2 in alias (4*r0scale**2)
869 
870  DO js = 2, ns - 1
871  arnorm = sum(wint(js:nrzt:ns)*ru0(js:nrzt:ns)**2)
872  aznorm = sum(wint(js:nrzt:ns)*zu0(js:nrzt:ns)**2)
873  IF (arnorm.eq.zero .or. aznorm.eq.zero) THEN
874  stop 'arnorm or aznorm=0 in bcovar'
875  END IF
876 
877  tcon(js) = min(abs(ard(js,1)/arnorm),
878  & abs(azd(js,1)/aznorm))*tcon_mul*(32*hs)**2
879  END DO
880  tcon(ns) = p5*tcon(ns-1)
881  IF (lasym) tcon = p5*tcon
882 #endif
883  ENDIF
884 
885 !
886 ! COMPUTE COVARIANT BSUBU,V (EVEN, ODD) ON HALF RADIAL MESH
887 ! FOR FORCE BALANCE AND RETURN (IEQUI=1)
888 !
889  IF (iequi .eq. 1) THEN
890 
891  DO js = ns - 1, 2, -1
892  DO l = js, nrzt, ns
893  bsubvh(l) = 2*bsubv_e(l) - bsubvh(l + 1)
894  END DO
895  END DO
896 
897 ! ADJUST <bsubvh> AFTER MESH-BLENDING
898  DO js = 2, ns
899  curpol_temp = fpsi(js)
900  & - sum(bsubvh(js:nrzt:ns)*wint(js:nrzt:ns))
901  DO l = js, nrzt, ns
902  bsubvh(l) = bsubvh(l) + curpol_temp
903  END DO
904  END DO
905 
906  bsubu_e(:nrzt) = bsubuh(:nrzt)
907  bsubv_e(:nrzt) = bsubvh(:nrzt)
908 
909  bsubu_o(:nrzt) = shalf(:nrzt)*bsubu_e(:nrzt)
910  bsubv_o(:nrzt) = shalf(:nrzt)*bsubv_e(:nrzt)
911 
912  RETURN
913 
914  END IF
915 
916 ! MINUS SIGN => HESSIAN DIAGONALS ARE POSITIVE
917  bsubu_e = -lamscale*bsubu_e
918  bsubv_e = -lamscale*bsubv_e
919  bsubu_o(:nrzt) = sqrts(:nrzt)*bsubu_e(:nrzt)
920  bsubv_o(:nrzt) = sqrts(:nrzt)*bsubv_e(:nrzt)
921 
922 !
923 ! STORE LU * LV COMBINATIONS USED IN FORCES
924 !
925 !WAC, SPH122407: sigma_an (=1 for isotropic case)
926  lvv(2:nrzt) = gsqrt(2:nrzt)*sigma_an(2:nrzt)
927  guu(2:nrzt) = bsupu(2:nrzt)*bsupu(2:nrzt)*lvv(2:nrzt)
928  guv(2:nrzt) = bsupu(2:nrzt)*bsupv(2:nrzt)*lvv(2:nrzt)
929  gvv(2:nrzt) = bsupv(2:nrzt)*bsupv(2:nrzt)*lvv(2:nrzt)
930  lv(2:nrzt,0) = bsq(2:nrzt)*tau(2:nrzt)
931  lu(2:nrzt,0) = bsq(2:nrzt)*r12(2:nrzt)
932 
933  END SUBROUTINE bcovar
Generated on Thu Mar 5 2020 15:49:26 for V3FIT by   doxygen 1.8.17