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,
6 USE realspace,
ONLY: pextra1, pextra2, pextra3, pextra4,
7 & pguu, pguv, pgvv, pru, pzu,
8 & pr1, prv, pzv, pshalf, pwint, pz1,
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
16 USE precon2d,
ONLY: ictrl_prec2d, lhess_exact,
19 USE vmec_input,
ONLY: nzeta
20 USE vmec_dim,
ONLY: ntheta3
21 USE parallel_include_module
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),
29 INTEGER,
INTENT(inout) :: ier_flag
35 REAL(dp),
PARAMETER :: c1p5 = (one + p5)
39 INTEGER :: l, js, ndim
40 REAL(dp) :: r2, volume, curpol_temp
42 REAL(dp) :: arnorm, aznorm, tcon_mul
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
49 INTEGER :: i, j, k, nsmin, nsmax, istat
50 REAL(dp) :: wblocal(ns), wbtotal
51 REAL(dp) :: wplocal(ns), wptotal
53 REAL(dp) :: fnlocal(ns), fntotal
54 REAL(dp) :: fn1local(ns), fn1total
55 REAL(dp) :: fnLlocal(ns), fnLtotal
57 IF (irst.EQ.2 .AND. iequi.EQ.0)
RETURN
77 nsmin=t1lglob; nsmax=t1rglob
78 pguu(:,nsmin:nsmax) = 0
79 pguv(:,nsmin:nsmax) = 0
80 pgvv(:,nsmin:nsmax) = 0
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)
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
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)
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)*
119 & (phipog(:,l) + phipog(:,l-1)))
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)))
134 nsmin = tlglob; nsmax = t1rglob
136 tau(:,l) = gsqrt(:,l)
137 gsqrt(:,l) = r12(:,l)*tau(:,l)
139 gsqrt(:,1) = gsqrt(:,2)
141 nsmin = max(2,tlglob); nsmax = t1rglob
142 pgvv(:,nsmin:nsmax) = pgvv(:,nsmin:nsmax)
143 & + r12sq(:,nsmin:nsmax)
147 nsmin = max(2,tlglob); nsmax = t1rglob
148 WHERE (gsqrt(:,nsmin:nsmax) .ne. zero)
149 phipog(:,nsmin:nsmax) = one/gsqrt(:,nsmin:nsmax)
156 vp(js) = signgs*sum(gsqrt(:,js)*pwint(:,js))
175 lu(:,l,:) = lu(:,l,:)*lamscale
176 lv(:,l,:) = lv(:,l,:)*lamscale
177 lu(:,l,0) = lu(:,l,0) + phipf(l)
180 nsmin = max(2,t1lglob)
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)))
190 IF (rank .EQ. 0)
THEN
204 #if defined(CHI_FORCE)
205 CALL add_fluxes_par(phipog, bsupu, bsupv, ictrl_prec2d.EQ.0)
207 CALL add_fluxes_par(phipog, bsupu, bsupv, .true.)
216 bsubuh(:,l) = pguu(:,l)*bsupu(:,l) + pguv(:,l)*bsupv(:,l)
217 bsubvh(:,l) = pguv(:,l)*bsupu(:,l) + pgvv(:,l)*bsupv(:,l)
227 bsq(:,l) = p5*(bsupu(:,l)*bsubuh(:,l) + bsupv(:,l)*bsubvh(:,l))
230 nsmin = max(2,tlglob)
231 nsmax = min(ns,t1rglob)
232 pres(nsmin:nsmax) = mass(nsmin:nsmax)/vp(nsmin:nsmax)**gamma
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)
241 CALL gather1xarray(wblocal)
242 wbtotal = sum(wblocal(2:ns))
243 CALL gather1xarray(wplocal)
244 wptotal = sum(wplocal(2:ns))
253 bsq(:,l) = bsq(:,l) + pres(l)
254 lvv(:,l) = phipog(:,l)*pgvv(:,l)
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)
270 lvv(:,l) = lvv(:,l)*pshalf(:,l)
271 bsubu_e(:,l) = pguv(:,l)*bsupu(:,l)
275 nsmax = min(ns - 1, trglob)
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))
281 bsubv_e(:,ns) = bsubv_e(:,ns)
282 & + p5*(lvv(:,ns)*lu(:,ns,1) + bsubu_e(:,ns))
288 CALL calc_fbal_par(bsubuh, bsubvh)
289 rbtor0= c1p5*fpsi(2) - p5*fpsi(3)
290 rbtor = c1p5*fpsi(ns) - p5*fpsi(ns-1)
309 IF (lhess_exact)
THEN
310 lctor = lfreeb .AND. ictrl_prec2d.NE.0
312 lctor = lfreeb .AND. ictrl_prec2d.GT.1
316 IF (ictrl_prec2d .EQ. 2)
THEN
317 ctor_prec2d = p5*(buco(ns) - buco(ns1))
319 ctor = signgs*twopi*(buco(ns)+ctor_prec2d)
321 ctor = signgs*twopi*(c1p5*buco(ns) - p5*buco(ns1))
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
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)
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)
363 IF (iequi .EQ. 1)
THEN
364 lu(:,nsmin:nsmax,0) = bsupv(:,nsmin:nsmax)
365 lv(:,nsmin:nsmax,0) = bsupu(:,nsmin:nsmax)
373 IF ((mod(iter2-iter1,ns4).EQ.0 .AND. iequi.EQ.0) .AND.
374 & ictrl_prec2d.EQ.0)
THEN
377 phipog(:,nsmin:nsmax) = phipog(:,nsmin:nsmax)
378 & * pwint(:,nsmin:nsmax)
380 CALL lamcal_par(phipog, pguu, pguv, pgvv)
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)
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)
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
403 pguu(:,nsmin:nsmax) = pguu(:,nsmin:nsmax)
404 & * r12(:,nsmin:nsmax)**2
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
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))
420 r2 = max(wb,wp)/volume
421 fnorm = one/(fntotal*(r2*r2))
423 fnorml = one/fnltotal
430 tcon0 = min(abs(tcon0), one)
431 tcon_mul = tcon0*(1 + r2*(one/60 + r2/(200*120)))
433 tcon_mul = tcon_mul/((4*r0scale**2)**2)
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)
442 IF (arnorm .eq. zero .or. aznorm .eq. zero)
THEN
443 ier_flag = arz_bad_value_flag
446 tcon(js) = min(abs(ard(js,1)/arnorm),
447 1 abs(azd(js,1)/aznorm))*tcon_mul*(32*hs)**2
449 tcon(ns) = p5*tcon(ns - 1)
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
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)
477 nsmin=max(tlglob,2); nsmax=min(trglob,ns)
479 curpol_temp = fpsi(js) - sum(bsubvh(:,js)*pwint(:,js))
480 bsubvh(:,js) = bsubvh(:,js) + curpol_temp
483 bsubu_e(:,nsmin:nsmax) = bsubuh(:,nsmin:nsmax)
484 bsubv_e(:,nsmin:nsmax) = bsubvh(:,nsmin:nsmax)
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)
494 nsmin = max(tlglob,2)
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)
507 nsmin = max(tlglob,2)
510 lvv(:,l) = gsqrt(:,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)
518 END SUBROUTINE bcovar_par
520 SUBROUTINE bcovar (lu, lv)
521 USE vmec_main, fpsi => bvco, p5 => cp5
522 USE vmec_params,
ONLY: ns4, signgs, pdamp, lamscale
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
530 USE precon2d,
ONLY: ictrl_prec2d, lhess_exact,
533 USE angle_constraints,
ONLY: precondn_rho, ard2, arm2,
534 & azd2, azm2, brd2, brm2, bzd2, bzm2
541 REAL(dp),
DIMENSION(nrzt,0:1),
INTENT(inout) :: lu, lv
547 REAL(dp),
PARAMETER :: c1p5 = (one + p5)
551 INTEGER :: l, js, ndim
552 REAL(dp) :: r2, volume, curpol_temp
554 REAL(dp) :: arnorm, aznorm, tcon_mul
556 REAL(dp),
POINTER,
DIMENSION(:) :: luu, luv, lvv, tau
557 REAL(dp),
DIMENSION(:),
POINTER :: bsupu, bsubuh,
558 & bsupv, bsubvh, r12sq
566 tau => extra1(:,1); luu => extra2(:,1);
567 luv => extra3(:,1); lvv => extra4(:,1)
569 bsupu => luu; bsubuh => bsubu_o
570 bsupv => luv; bsubvh => bsubv_o
578 guu(ndim) = 0; guv = 0; gvv = 0
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))
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)
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
613 r12sq(1:nrzt) = r1(1:nrzt,0)*r1(1:nrzt,0)
614 & + r12sq(1:nrzt)*r1(1:nrzt,1)*r1(1:nrzt,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) +
620 & shalf(l)*(phipog(l) + phipog(l-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)))
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)
637 gvv(2:nrzt) = gvv(2:nrzt) + r12sq(2:nrzt)
640 WHERE (gsqrt(2:ndim) .ne. zero)
641 phipog(2:ndim) = one/gsqrt(2:ndim)
643 phipog(1:ndim:ns) = 0
645 vp(1) = 0; vp(ns+1) = 0
647 vp(js) = signgs*sum(gsqrt(js:nrzt:ns)*wint(js:nrzt:ns))
649 IF (iter2 .eq. 1)
THEN
650 voli = twopi*twopi*hs*sum(vp(2:ns))
668 lu(js:nrzt:ns,0) = lu(js:nrzt:ns,0) + phipf(js)
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)))
676 bsupu(1)=0; bsupu(ndim)=0
677 bsupv(1)=0; bsupv(ndim)=0
692 #if defined(CHI_FORCE)
693 CALL add_fluxes(phipog, bsupu, bsupv, ictrl_prec2d.EQ.0)
695 CALL add_fluxes(phipog, bsupu, bsupv, .true.)
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)
706 bsubuh(ndim) = 0; bsubvh(ndim) = 0
711 bsq(:nrzt) = p5*(bsupu(:nrzt)*bsubuh(:nrzt)
712 & + bsupv(:nrzt)*bsubvh(:nrzt))
714 wb = hs*abs(sum(wint(:nrzt)*gsqrt(:nrzt)*bsq(:nrzt)))
718 CALL an_pressure(bsubu_e, bsubv_e)
721 bsq(2:nrzt) = bsq(2:nrzt) + pperp(2:nrzt)
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)
729 pres(2:ns) = mass(2:ns)/vp(2:ns)**gamma
730 wp = hs*sum(vp(2:ns)*pres(2:ns))
734 bsq(js:nrzt:ns) = bsq(js:nrzt:ns) + pres(js)
740 lvv = phipog(:ndim)*gvv
741 bsubv_e(1:nrzt) = p5*(lvv(1:nrzt)+lvv(2:ndim))*lu(1:nrzt,0)
744 bsubu_e(:nrzt) = guv(:nrzt)*bsupu(:nrzt)*sigma_an(:nrzt)
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))
755 IF (iequi .EQ. 1) papr = pmap*pres/vp
757 CALL calc_fbal(bsubuh, bsubvh)
759 rbtor0= c1p5*fpsi(2) - p5*fpsi(3)
760 rbtor = c1p5*fpsi(ns) - p5*fpsi(ns-1)
779 IF (lhess_exact)
THEN
780 lctor = lfreeb .and. ictrl_prec2d.ne.0
782 lctor = lfreeb .and. ictrl_prec2d.gt.1
785 IF (ictrl_prec2d .eq. 2)
THEN
786 ctor_prec2d = signgs*twopi*p5*(buco(ns) - buco(ns1))
788 ctor = ctor_prec2d + signgs*twopi*buco(ns)
790 ctor = signgs*twopi*(c1p5*buco(ns) - p5*buco(ns1))
798 lvv(l:nrzt:ns) = bdamp(l)
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'
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)
814 IF (iequi .EQ. 1)
THEN
815 lu(:nrzt,0) = bsupv(:nrzt)
816 lv(:nrzt,0) = bsupu(:nrzt)
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,
830 & arm2, ard2, brm2, brd2,
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,
836 & azm2, azd2, bzm2, bzd2,
838 & crd, rru_fac, sin01)
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
848 volume = hs*sum(vp(2:ns))
849 r2 = max(wb,wp)/volume
850 guu(:nrzt) = guu(:nrzt)*r12(:nrzt)**2
851 fnorm = one/(sum(guu(1:nrzt)*wint(1:nrzt))*(r2*r2))
852 fnorm1 = one/sum(xc(1+ns:2*irzloff)**2)
853 fnorml = one/(sum((bsubuh(1:nrzt)**2 + bsubvh(1:nrzt)**2) *
854 & wint(1:nrzt))*lamscale**2)
864 tcon0 = min(abs(tcon0), one)
865 tcon_mul = tcon0*(1 + r2*(one/60 + r2/(200*120)))
867 tcon_mul = tcon_mul/((4*r0scale**2)**2)
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'
877 tcon(js) = min(abs(ard(js,1)/arnorm),
878 & abs(azd(js,1)/aznorm))*tcon_mul*(32*hs)**2
880 tcon(ns) = p5*tcon(ns-1)
881 IF (lasym) tcon = p5*tcon
889 IF (iequi .eq. 1)
THEN
891 DO js = ns - 1, 2, -1
893 bsubvh(l) = 2*bsubv_e(l) - bsubvh(l + 1)
899 curpol_temp = fpsi(js)
900 & - sum(bsubvh(js:nrzt:ns)*wint(js:nrzt:ns))
902 bsubvh(l) = bsubvh(l) + curpol_temp
906 bsubu_e(:nrzt) = bsubuh(:nrzt)
907 bsubv_e(:nrzt) = bsubvh(:nrzt)
909 bsubu_o(:nrzt) = shalf(:nrzt)*bsubu_e(:nrzt)
910 bsubv_o(:nrzt) = shalf(:nrzt)*bsubv_e(:nrzt)
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)
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)
933 END SUBROUTINE bcovar