2 USE v3_utilities,
ONLY:
assert
4 USE descriptor_mod,
ONLY: iam, nprocs, siesta_comm
5 USE nscalingtools,
ONLY: startglobrow, endglobrow, mpi_err
9 USE utilities,
ONLY: to_half_mesh
17 INTEGER,
PRIVATE :: unit35=35, unit36=36, nsmin, nsmax
18 REAL(dp) :: divb_rms, avbeta, bdotj_rms, divj_rms, bdotj2_rms
19 REAL(dp) :: toroidal_flux, toroidal_flux0=0, bgradp_rms, wbgradp
20 REAL(dp) :: max_bgradp, min_bgradp
21 REAL(dp) :: tnorm_bgradp
26 REAL(dp),
ALLOCATABLE,
DIMENSION(:,:,:) :: bdotjmnch, bdotjijh
27 REAL(dp),
ALLOCATABLE,
DIMENSION(:,:,:) :: bdotjmnsh
28 REAL(dp),
ALLOCATABLE,
DIMENSION(:,:,:) :: divjmnsh, divbmnsf
29 REAL(dp),
ALLOCATABLE,
DIMENSION(:,:,:) :: divjmnch, divbmncf
30 REAL(dp),
ALLOCATABLE,
DIMENSION(:,:,:) :: bgradpf, jpmnch_loc
31 REAL(dp),
ALLOCATABLE,
DIMENSION(:,:,:),
TARGET :: fgradpmnch
32 REAL(dp),
DIMENSION(:),
POINTER :: xc, gc
33 LOGICAL :: lextend = .false.
34 LOGICAL :: lcurr_init = .false.
50 SUBROUTINE divb(ns_min, ns_max)
55 INTEGER,
INTENT(in) :: ns_min
56 INTEGER,
INTENT(in) :: ns_max
65 REAL (dp),
DIMENSION(2) :: temp
75 nh = min(ns - 1, nsmax)
77 IF (
ALLOCATED(divbmnsf))
THEN
80 IF (
ALLOCATED(divbmncf))
THEN
83 ALLOCATE(divbmnsf(0:mpol,-ntor:ntor,nsmin:nsmax),
84 divbmncf(0:mpol,-ntor:ntor,nsmin:nsmax), stat=istat)
85 CALL assert(istat.EQ.0,
'Allocation failed in divb')
92 divbmnsf(m,n,nl:nh) = -m*(jbsupumnch(m,n,nl+1:nh+1) +
93 jbsupumnch(m,n,nl:nh))*0.5
94 - n*nfp*(jbsupvmnch(m,n,nl+1:nh+1) +
95 jbsupvmnch(m,n,nl:nh))*0.5
96 + ohs*(jbsupsmnsh(m,n,nl+1:nh+1) -
97 jbsupsmnsh(m,n,nl:nh))
99 divbmncf(m,n,nl:nh) = m*(jbsupumnsh(m,n,nl+1:nh+1) +
100 jbsupumnsh(m,n,nl:nh))*0.5
101 + n*nfp*(jbsupvmnsh(m,n,nl+1:nh+1) +
102 jbsupvmnsh(m,n,nl:nh))*0.5
103 + ohs*(jbsupsmnch(m,n,nl+1:nh+1) -
104 jbsupsmnch(m,n,nl:nh))
109 tnorm = hs_i*sum(jbsupumnch(:,:,nsmin:nsmax)**2
110 + jbsupvmnch(:,:,nsmin:nsmax)**2
111 + jbsupsmnsh(:,:,nsmin:nsmax)**2)
113 divb_rms = sum(divbmnsf**2 + divbmncf**2)
117 CALL mpi_allreduce(mpi_in_place,temp,2,mpi_real8, mpi_sum,
122 IF (tnorm .NE. zero)
THEN
123 divb_rms = sqrt(divb_rms/tnorm)
124 divbmnsf = divbmnsf/sqrt(tnorm)
125 divbmncf = divbmncf/sqrt(tnorm)
129 time_divb = time_divb + (toff-ton)
133 SUBROUTINE write_profiles(fsq_total)
139 REAL(dp) :: fsq_total
143 INTEGER,
PARAMETER :: ifull=0, ihalf=1
145 REAL(dp),
ALLOCATABLE,
DIMENSION(:,:,:) :: pijh, pmnh
146 REAL(dp) :: average(ns)
152 IF (iam .NE. 0)
RETURN
154 ALLOCATE (pijh(ntheta,nzeta,ns),
155 pmnh(0:mpol,-ntor:ntor,ns), stat=istat)
156 CALL assert(istat.EQ.0,
'Allocation failed in WRITE_PROFILES')
158 CALL safe_open(unit35, istat,
"siesta_profiles.txt",
'replace',
160 CALL safe_open(unit36, istat,
"siesta_profiles_pest.txt",
161 'replace',
'formatted')
163 CALL assert(istat.EQ.0,
'Error opening siesta profiles data file'
166 CALL assert(startglobrow.EQ.1 .AND. endglobrow.EQ.ns,
167 ' INIT_STATE NOT INITIALIZED IN WRITE_PROFILES!')
172 CALL write_spectra(pijf0, pmnh,
f_cos,
f_sin, ifull,
'p', fsq_total
175 CALL write_spectra(bsupsijf0, pmnh,
f_sin,
f_cos, ifull,
'B^s', fsq_total
178 CALL write_spectra(bsupuijf0, pmnh,
f_cos,
f_sin, ifull,
'B^u', fsq_total
181 CALL write_spectra(bsupvijf0, pmnh,
f_cos,
f_sin, ifull,
'B^v', fsq_total
184 pijh = ksupsijf0/jacobf
185 CALL write_spectra(pijh, pmnh,
f_sin,
f_cos, ifull,
'J^s', fsq_total
188 pijh = ksupuijf0/jacobf
189 CALL write_spectra(pijh, pmnh,
f_cos,
f_sin, ifull,
'J^u', fsq_total
192 pijh = ksupvijf0/jacobf
193 CALL write_spectra(pijh, pmnh,
f_cos,
f_sin, ifull,
'J^v', fsq_total
197 pijh = bdotjijh*jacobh
199 CALL write_spectra(pijh, pmnh,
f_cos,
f_sin, ihalf,
'K||', fsq_total
207 CALL surfaverage(average, abs(pijh), 1, ns)
209 WRITE (unit36, *)
'<ABS(K||) - nonaxisymmetric>'
211 WRITE (unit36,
'(i4,1p,e10.2)') js, average(js)
215 CALL siesta_profiles(unit35, jvsupsmncf,
'vel^s_cosmn; fsq_total: '
216 IF (lasym)
CALL siesta_profiles(unit35, jvsupsmnsf,
'vel^s_sinmn; fsq_total: '
217 CALL siesta_profiles(unit35, jvsupumnsf,
'vel^u_sinmn; fsq_total: '
218 IF (lasym)
CALL siesta_profiles(unit35, jvsupumncf,
'vel^u_cosmn; fsq_total: '
219 CALL siesta_profiles(unit35, jvsupvmnsf,
'vel^v_sinmn; fsq_total: '
220 IF (lasym)
CALL siesta_profiles (unit35, jvsupvmncf,
'vel^v_cosmn; fsq_total: '
222 DEALLOCATE(pijh, pmnh)
227 END SUBROUTINE write_profiles
244 SUBROUTINE write_spectra(xij, xmn, fsym, fasym, igrid, label, fsq_total)
245 USE vmec_info,
ONLY: lmns => lmns_i, lmnc => lmnc_i
247 USE utilities,
ONLY: to_full_mesh
252 REAL (dp),
DIMENSION(:,:,:),
ALLOCATABLE,
INTENT(inout) :: xij
253 REAL (dp),
DIMENSION(:,:,:),
ALLOCATABLE,
INTENT(inout) :: xmn
254 INTEGER,
INTENT(IN) :: fsym
255 INTEGER,
INTENT(IN) :: fasym
256 INTEGER,
INTENT(IN) :: igrid
257 CHARACTER (len=*),
INTENT(in) :: label
258 REAL (dp),
INTENT(IN) :: fsq_total
261 REAL (dp),
DIMENSION(:,:,:),
ALLOCATABLE :: xijf
262 CHARACTER (LEN=128) :: mylabel
265 IF (igrid .EQ. 1)
THEN
266 ALLOCATE(xijf(
SIZE(xij,1),
SIZE(xij,2),
SIZE(xij,3)))
267 CALL to_full_mesh(xij, xijf)
272 IF (fsym .EQ.
f_cos)
THEN
273 WRITE (mylabel,1000) trim(label)
275 WRITE (mylabel,1001) trim(label)
279 CALL siesta_profiles(unit35, xmn, trim(mylabel), fsq_total)
282 CALL siesta_profiles(unit36, xmn, trim(mylabel), fsq_total)
286 IF (fasym .EQ.
f_cos)
THEN
287 WRITE (mylabel,1000) trim(label)
289 WRITE (mylabel,1001) trim(label)
293 CALL siesta_profiles(unit35, xmn, trim(mylabel), fsq_total)
296 CALL siesta_profiles(unit36, xmn, trim(mylabel), fsq_total)
300 1000
FORMAT(a,
'_cosmn; fsq_total: ')
301 1001
FORMAT(a,
'_sinmn; fsq_total: ')
305 SUBROUTINE siesta_profiles (iunit, arr_value, label, fsq_total)
309 INTEGER,
INTENT(IN) :: iunit
310 REAL(dp),
INTENT(IN) :: arr_value(0:mpol,-ntor:ntor,ns)
311 REAL(dp),
INTENT(IN) :: fsq_total
312 CHARACTER*(*),
INTENT(IN) :: label
318 REAL(dp) :: sj, rms_value
319 CHARACTER (LEN=35) :: p_format
320 INTEGER,
DIMENSION(:,:),
ALLOCATABLE :: m_array
321 INTEGER,
DIMENSION(:,:),
ALLOCATABLE :: n_array
325 WRITE (iunit,
'(1x,a,1p,e12.4)') trim(label), fsq_total
327 ALLOCATE(m_array(0:mpol,-ntor:ntor))
328 ALLOCATE(n_array(0:mpol,-ntor:ntor))
337 mnmax = (mpol + 1)*(2*ntor + 1)
339 WRITE (p_format,1000) mnmax
340 WRITE (iunit,p_format)
' ',
' ',
'MPOL--->', m_array
341 WRITE (p_format,1000) mnmax
342 WRITE (iunit,p_format)
'RADIUS',
'RMS',
'NTOR--->', n_array
344 WRITE (p_format,1002) mnmax
347 rms_value = sqrt(sum(arr_value(:,:,js)**2)/mnmax)
348 WRITE(iunit,p_format) sj, rms_value, arr_value(:,:,js)
351 DEALLOCATE(m_array, n_array)
353 1000
FORMAT(
'(a,8x,a,6x,a,',i3,
'(2x,i12))')
354 1002
FORMAT(
'(f6.3,2x,es12.5,14x',i4,
'(2x,es12.5))')
356 END SUBROUTINE siesta_profiles
359 SUBROUTINE bgradp(ns_min, ns_max)
368 INTEGER,
INTENT(IN) :: ns_min, ns_max
373 REAL(dp) :: s1(2), r1(2), ton, toff
378 nsmin=ns_min; nsmax=ns_max
380 ALLOCATE(bgradpf(ntheta,nzeta,nsmin:nsmax), stat=istat)
381 CALL assert(istat.EQ.0,
'Allocation error in BGRADP')
382 CALL get_bgradp(nsmin, nsmax)
384 DEALLOCATE (bgradpf, stat=istat)
386 END SUBROUTINE bgradp
394 SUBROUTINE get_bgradp(ns_min, ns_max)
399 INTEGER,
INTENT(in) :: ns_min
400 INTEGER,
INTENT(in) :: ns_max
406 REAL (dp),
DIMENSION(:,:,:),
ALLOCATABLE :: dpduijf
407 REAL (dp),
DIMENSION(:,:,:),
ALLOCATABLE :: dpdvijf
408 REAL (dp),
DIMENSION(2) :: temp(2)
415 CALL assert(
ALLOCATED(jbsupsmnsh),
'jbsupXmn unallocated in get_bgradp'
418 nsmax = min(ns_max + 1, ns)
423 ALLOCATE(dpduijf(ntheta,nzeta,nsmin:nsmax),
424 dpdvijf(ntheta,nzeta,nsmin:nsmax), stat=istat)
425 CALL assert_eq(0, istat,
'Allocation failed in get_bgradp')
427 CALL to_full_mesh(pijh0_du, dpduijf)
428 CALL to_full_mesh(pijh0_dv, dpdvijf)
430 nsmax = min(endglobrow, ns)
434 bgradpf(:,:,n1:n2) = bsupsijf0(:,:,n1:n2)*pijf0_ds(:,:,n1:n2)
435 + bsupuijf0(:,:,n1:n2)*dpduijf(:,:,n1:n2)
436 + bsupvijf0(:,:,n1:n2)*dpdvijf(:,:,n1:n2)
438 tnorm_bgradp = sum((bsupsijf0(:,:,n1:nsmax)**2 +
439 bsupuijf0(:,:,n1:nsmax)**2 +
440 bsupvijf0(:,:,n1:nsmax)**2)*
441 pijf0(:,:,n1:nsmax)**2*wint(:,:,n1:nsmax))
444 n1 = max(ns/20, nsmin)
445 wbgradp = sum(bgradpf(:,:,n1:n2)**2*wint(:,:,n1:n2))
446 temp(1) = tnorm_bgradp
450 CALL mpi_allreduce(mpi_in_place, temp, 2, mpi_real8, mpi_sum,
451 siesta_comm, mpi_err)
454 tnorm_bgradp = temp(1)
457 IF (tnorm_bgradp .GT. zero)
THEN
458 bgradp_rms = sqrt(wbgradp/tnorm_bgradp)
459 tnorm_bgradp = sqrt(tnorm_bgradp/ns)
462 max_bgradp = maxval(bgradpf(:,:,n1:n2))/tnorm_bgradp
463 min_bgradp = minval(bgradpf(:,:,n1:n2))/tnorm_bgradp
465 temp(2) = -min_bgradp
468 CALL mpi_allreduce(mpi_in_place, temp, 2, mpi_real8, mpi_max
469 siesta_comm, mpi_err)
473 min_bgradp = -temp(2)
476 max_bgradp = max(max_bgradp, -min_bgradp)
477 min_bgradp = min(min_bgradp, -max_bgradp)
481 DEALLOCATE (dpdvijf, dpduijf, stat=istat)
484 time_bgradp = time_bgradp + (toff - ton)
506 nsmin=max(2,startglobrow); nsmax=min(endglobrow,ns)
509 toroidal_flux = sum(jbsupvmnch(
m0,
n0,nsmin:nsmax))
512 CALL mpi_allreduce(mpi_in_place,toroidal_flux,1,mpi_real8, mpi_sum
518 toroidal_flux = signjac*twopi*toroidal_flux*hs_i/b_factor
520 IF (toroidal_flux0 .EQ. zero) toroidal_flux0 = toroidal_flux
525 SUBROUTINE bdotj (bdotjijhA)
533 REAL(dp),
ALLOCATABLE,
DIMENSION(:,:,:) :: bdotjijhA
538 REAL(dp),
ALLOCATABLE,
DIMENSION(:,:,:) :: &
540 bsubsijh, bsubuijh, bsubvijh, ksubsijh, ksubuijh,
541 ksubvijh, ksupsijh, ksupuijh, ksupvijh, work1, work2,
542 work3, work4, work5, work6
543 REAL(dp) :: tnorm, tnorm2, ton, toff
546 CALL assert(
SIZE(bsupsijh0,3).EQ.ns,
'bsupsijh0 wrong size in bdotj'
550 ALLOCATE(bsubsijh(ntheta,nzeta,ns),
551 bsubuijh(ntheta,nzeta,ns),
552 bsubvijh(ntheta,nzeta,ns), stat=istat)
553 CALL assert(istat.EQ.0,
'Allocation #2 failed in BDOTJ')
554 bsubsijh=0; bsubvijh=0; bsubuijh=0
556 CALL tolowerh(bsupsijh0, bsupuijh0, bsupvijh0,
557 bsubsijh, bsubuijh, bsubvijh, 1, ns)
565 ALLOCATE(ksupsijh(ntheta,nzeta,ns),
566 ksupuijh(ntheta,nzeta,ns),
567 ksupvijh(ntheta,nzeta,ns), stat=istat)
568 CALL assert(istat.EQ.0,
'Allocation #3 failed in BDOTJ')
570 CALL to_half_mesh(ksupsijf0, ksupsijh)
571 CALL to_half_mesh(ksupuijf0, ksupuijh)
572 CALL to_half_mesh(ksupvijf0, ksupvijh)
576 ksupsijh = ksupsijh/jacobh
577 ksupuijh = ksupuijh/jacobh
578 ksupvijh = ksupvijh/jacobh
582 ALLOCATE(ksubsijh(ntheta,nzeta,ns),
583 ksubuijh(ntheta,nzeta,ns),
584 ksubvijh(ntheta,nzeta,ns), stat=istat)
585 CALL assert(istat.EQ.0,
'Allocation #4 failed in BDOTJ')
586 ksubsijh=0; ksubuijh=0; ksubvijh=0
588 CALL tolowerh(ksupsijh, ksupuijh, ksupvijh,
589 ksubsijh, ksubuijh, ksubvijh, 1, ns)
593 IF (
ALLOCATED(bdotjmnch) .AND.
SIZE(bdotjmnch,3).NE. ns)
594 DEALLOCATE(bdotjmnch)
595 IF (.NOT.
ALLOCATED(bdotjmnch))
THEN
596 ALLOCATE(bdotjmnch(0:mpol,-ntor:ntor,ns), stat=istat)
597 CALL assert(istat.EQ.0,
'Allocation #5 failed in BDOTJ')
600 IF (
ALLOCATED(bdotjmnsh) .AND.
SIZE(bdotjmnsh,3).NE. ns)
601 DEALLOCATE(bdotjmnsh)
602 IF (.NOT.
ALLOCATED(bdotjmnsh))
THEN
603 ALLOCATE(bdotjmnsh(0:mpol,-ntor:ntor,ns), stat=istat)
604 CALL assert(istat.EQ.0,
'Allocation #5 failed in BDOTJ')
607 ALLOCATE(bdotjijha(ntheta,nzeta,ns), stat=istat)
608 CALL assert(istat.EQ.0,
'Allocation #6 failed in BDOTJ')
611 bdotjijha(:,:,3:nsh-1) =
612 bsupsijh0(:,:,3:nsh-1)*ksubsijh(:,:,3:nsh-1) +
613 bsupuijh0(:,:,3:nsh-1)*ksubuijh(:,:,3:nsh-1) +
614 bsupvijh0(:,:,3:nsh-1)*ksubvijh(:,:,3:nsh-1)
618 ALLOCATE(work1(ntheta,nzeta,ns), work2(ntheta,nzeta,ns),
619 work3(ntheta,nzeta,ns), work4(ntheta,nzeta,ns),
620 work5(ntheta,nzeta,ns), work6(ntheta,nzeta,ns),
622 CALL assert(istat.EQ.0,
'Allocation #7 failed in BDOTJ')
623 work1=0; work2=0; work3=0
624 work4=0; work5=0; work6=0
626 work1 = (bsubuijh*ksubvijh - bsubvijh*ksubuijh)/jacobh
627 work2 = (bsubvijh*ksubsijh - bsubsijh*ksubvijh)/jacobh
628 work3 = (bsubsijh*ksubuijh - bsubuijh*ksubsijh)/jacobh
629 CALL tolowerh(work1, work2, work3, work4, work5, work6, 1, ns)
631 work1 = work1*work4 + work2*work5 + work3*work6
632 tnorm = sum(work1(:,:,3:nsh-1))
633 IF (iam.EQ.0 .AND. any(work1 .LT. zero))
634 print *,
'ERROR1 IN BDOTJ, JXB^2 < 0!'
635 tnorm = sqrt(abs(tnorm))
636 DEALLOCATE(work3, work4, work5, work6)
638 work1 = bsupsijh0*bsubsijh + bsupuijh0*bsubuijh +
640 work2 = ksupsijh*ksubsijh + ksupuijh*ksubuijh +
643 tnorm2 = sum(work1(:,:,3:nsh-1)*work2(:,:,3:nsh-1))
644 work1(:,:,1) = 0; work2(:,:,1) = 0
645 IF (iam.EQ.0 .AND. any(work1.LT.zero))
646 print *,
'ERROR2 IN BDOTJ: B^2 < 0!'
647 IF (iam.EQ.0 .AND. any(work2.LT.zero))
648 print *,
'ERROR3 IN BDOTJ: J^2 < 0!'
649 tnorm2 = sqrt(abs(tnorm2))
650 bdotj_rms = sum(bdotjijha*bdotjijha)
652 DEALLOCATE (bsubsijh, bsubuijh, bsubvijh, ksupsijh, ksupuijh,
653 ksupvijh, ksubsijh, ksubuijh, ksubvijh,
656 tnorm = max(tnorm, epsilon(tnorm))
657 tnorm2= max(tnorm2,epsilon(tnorm2))
658 bdotj2_rms = sqrt(abs(bdotj_rms))/tnorm2
659 bdotj_rms = sqrt(abs(bdotj_rms))/tnorm
661 bdotjmnch = bdotjmnch/tnorm
664 bdotjmnsh = bdotjmnsh/tnorm
667 time_bdotj = time_bdotj+(toff-ton)
687 INTEGER :: istat, n1, n2, ns_span
688 REAL(dp),
ALLOCATABLE,
DIMENSION(:,:,:) :: bsupsijh, &
690 bsubsijh, bsubuijh, bsubvijh, ksubsijh, ksubuijh,
691 ksubvijh, ksupsijh, ksupuijh, ksupvijh, work1, work2,
692 work3, work4, work5, work6
693 REAL(dp) :: tnorm, tnorm2, temp(3), ton, toff
697 CALL assert(lcurr_init,
'MUST CALL init_state(.TRUE.) BEFORE bdotj_par'
698 nsmin=max(1,startglobrow); nsmax=min(ns,endglobrow)
699 ns_span=nsmax-nsmin+1
703 ALLOCATE(bsupsijh(ntheta,nzeta,nsmin:nsmax),
704 bsupuijh(ntheta,nzeta,nsmin:nsmax),
705 bsupvijh(ntheta,nzeta,nsmin:nsmax), stat=istat)
706 CALL assert(istat.EQ.0,
'Allocation #1 failed in BDOTJ_PAR')
723 bsupsijh = bsupsijh/jacobh(:,:,nsmin:nsmax)
724 bsupuijh = bsupuijh/jacobh(:,:,nsmin:nsmax)
725 bsupvijh = bsupvijh/jacobh(:,:,nsmin:nsmax)
729 ALLOCATE(bsubsijh(ntheta,nzeta,nsmin:nsmax),
730 bsubuijh(ntheta,nzeta,nsmin:nsmax),
731 bsubvijh(ntheta,nzeta,nsmin:nsmax), stat=istat)
732 CALL assert(istat.EQ.0,
'Allocation #2 failed in BDOTJ_PAR')
734 CALL tolowerh(bsupsijh, bsupuijh, bsupvijh,
735 bsubsijh, bsubuijh, bsubvijh, nsmin, nsmax)
742 nsmin=max(1,startglobrow-1)
743 ALLOCATE(ksupsijh(ntheta,nzeta,nsmin:nsmax),
744 ksupuijh(ntheta,nzeta,nsmin:nsmax),
745 ksupvijh(ntheta,nzeta,nsmin:nsmax), stat=istat)
746 CALL assert(istat.EQ.0,
'Allocation #3 failed in BDOTJ_PAR')
748 CALL to_half_mesh(ksupsijf0, ksupsijh)
749 CALL to_half_mesh(ksupuijf0, ksupuijh)
750 CALL to_half_mesh(ksupvijf0, ksupvijh)
755 ksupsijh = ksupsijh/jacobh(:,:,nsmin:nsmax)
756 ksupuijh = ksupuijh/jacobh(:,:,nsmin:nsmax)
757 ksupvijh = ksupvijh/jacobh(:,:,nsmin:nsmax)
762 ALLOCATE(ksubsijh(ntheta,nzeta,nsmin:nsmax),
763 ksubuijh(ntheta,nzeta,nsmin:nsmax),
764 ksubvijh(ntheta,nzeta,nsmin:nsmax), stat=istat)
765 CALL assert(istat.EQ.0,
'Allocation #4 failed in BDOTJ_PAR')
767 CALL tolowerh(ksupsijh, ksupuijh, ksupvijh,
768 ksubsijh, ksubuijh, ksubvijh, nsmin, nsmax)
769 nsmin=max(1,startglobrow)
774 IF (
ALLOCATED(bdotjmnch) .AND.
SIZE(bdotjmnch,3).NE.ns_span)
775 DEALLOCATE(bdotjmnch)
776 IF (.NOT.
ALLOCATED(bdotjmnch))
THEN
777 ALLOCATE(bdotjmnch(0:mpol,-ntor:ntor,nsmin:nsmax), stat=istat)
778 CALL assert(istat.EQ.0,
'Allocation #5A failed in BDOTJ_PAR')
781 IF (
ALLOCATED(bdotjmnsh) .AND.
SIZE(bdotjmnsh,3).NE.ns_span)
782 DEALLOCATE(bdotjmnsh)
783 IF (.NOT.
ALLOCATED(bdotjmnsh))
THEN
784 ALLOCATE(bdotjmnsh(0:mpol,-ntor:ntor,nsmin:nsmax), stat=istat)
785 CALL assert(istat.EQ.0,
'Allocation #5B failed in BDOTJ_PAR')
789 ALLOCATE(bdotjijh(ntheta,nzeta,nsmin:nsmax), stat=istat)
790 CALL assert(istat.EQ.0,
'Allocation #6 failed in BDOTJ_PAR')
793 n1 = max(3,nsmin); n2 = min(nsh-1, nsmax)
794 bdotjijh(:,:,n1:n2) =
795 bsupsijh(:,:,n1:n2)*ksubsijh(:,:,n1:n2) +
796 bsupuijh(:,:,n1:n2)*ksubuijh(:,:,n1:n2) +
797 bsupvijh(:,:,n1:n2)*ksubvijh(:,:,n1:n2)
802 ALLOCATE(work1(ntheta,nzeta,nsmin:nsmax),
803 work2(ntheta,nzeta,nsmin:nsmax),
804 work3(ntheta,nzeta,nsmin:nsmax),
805 work4(ntheta,nzeta,nsmin:nsmax),
806 work5(ntheta,nzeta,nsmin:nsmax),
807 work6(ntheta,nzeta,nsmin:nsmax),
809 CALL assert(istat.EQ.0,
'Allocation #7 failed in BDOTJ_PAR')
811 work1 = bsubuijh*ksubvijh(:,:,nsmin:nsmax)
812 - bsubvijh*ksubuijh(:,:,nsmin:nsmax)
813 work2 = bsubvijh*ksubsijh(:,:,nsmin:nsmax)
814 - bsubsijh*ksubvijh(:,:,nsmin:nsmax)
815 work3 = bsubsijh*ksubuijh(:,:,nsmin:nsmax)
816 - bsubuijh*ksubsijh(:,:,nsmin:nsmax)
817 work1 = work1/jacobh(:,:,nsmin:nsmax)
818 work2 = work2/jacobh(:,:,nsmin:nsmax)
819 work3 = work3/jacobh(:,:,nsmin:nsmax)
821 work4, work5, work6, nsmin, nsmax)
823 work1 = work1*work4 + work2*work5 + work3*work6
824 tnorm = sum(work1(:,:,n1:n2))
825 IF (nsmin .EQ. 1) work1(:,:,1) = 0
826 IF (iam.EQ.0 .AND. any(work1 .LT. zero))
827 print *,
'ERROR1 IN BDOTJ_PAR, JXB^2 < 0!'
828 DEALLOCATE(work3, work4, work5, work6)
830 work1 = bsupsijh*bsubsijh + bsupuijh*bsubuijh +
832 work2 = ksupsijh(:,:,nsmin:nsmax)*ksubsijh(:,:,nsmin:nsmax)
833 + ksupuijh(:,:,nsmin:nsmax)*ksubuijh(:,:,nsmin:nsmax)
834 + ksupvijh(:,:,nsmin:nsmax)*ksubvijh(:,:,nsmin:nsmax)
836 IF (nsmin .EQ. 1)
THEN
837 work1(:,:,1) = 0; work2(:,:,1) = 0
839 IF (iam.EQ.0 .AND. any(work1.LT.zero))
840 print *,
'ERROR2 IN BDOTJ_PAR: B^2 < 0!'
841 IF (iam.EQ.0 .AND. any(work2.LT.zero))
842 print *,
'ERROR3 IN BDOTJ_PAR: J^2 < 0!'
843 tnorm2 = sum(work1(:,:,n1:n2)*work2(:,:,n1:n2))
844 bdotj_rms = sum(bdotjijh(:,:,n1:n2)*bdotjijh(:,:,n1:n2))
846 temp(1)=tnorm; temp(2)=tnorm2; temp(3)=bdotj_rms
847 CALL mpi_allreduce(mpi_in_place,temp,3,mpi_real8, mpi_sum,
849 tnorm = sqrt(abs(temp(1)))
850 tnorm2 = sqrt(abs(temp(2)))
853 tnorm = sqrt(abs(tnorm)); tnorm2 = sqrt(abs(tnorm2))
855 DEALLOCATE (bsupsijh, bsupuijh, bsupvijh, bsubsijh,
856 bsubuijh, bsubvijh, ksupsijh, ksupuijh,
857 ksupvijh, ksubsijh, ksubuijh, ksubvijh,
861 tnorm2= max(tnorm2,epsilon(tnorm2))
862 bdotj2_rms = sqrt(abs(bdotj_rms))/tnorm2
863 bdotj_rms = sqrt(abs(bdotj_rms))/tnorm
865 bdotjijh = bdotjijh/tnorm
872 DEALLOCATE (bdotjijh)
875 time_bdotj = time_bdotj+(toff-ton)
877 END SUBROUTINE bdotj_par
888 SUBROUTINE divj(ns_min, ns_max)
894 INTEGER,
INTENT(in) :: ns_min
895 INTEGER,
INTENT(in) :: ns_max
903 REAL(dp),
DIMENSION(:,:,:),
ALLOCATABLE :: temp1
904 REAL(dp),
DIMENSION(:,:,:),
ALLOCATABLE :: temp2
906 REAL(dp),
DIMENSION(2) :: temp
913 CALL assert(lcurr_init,
'MUST CALL init_state(.TRUE.) BEFORE divJ_par'
918 IF (
ALLOCATED(divjmnsh))
THEN
921 IF (
ALLOCATED(divjmnch))
THEN
924 ALLOCATE(divjmnsh(0:mpol,-ntor:ntor,nsmin:nsmax),
925 divjmnch(0:mpol,-ntor:ntor,nsmin:nsmax), stat=istat)
926 CALL assert_eq(istat, 0,
'Allocation failed in DIVJ')
932 n2 = min(ns-1, nsmax)
935 divjmnsh(m,n,n1:n2) =
936 -m*(ksupumncf(m,n,n1:n2) + ksupumncf(m,n,n1-1:n2-1))
937 - n*nfp*(ksupvmncf(m,n,n1:n2) + ksupvmncf(m,n,n1-1:n2-1))
938 + ohs*(ksupsmnsf(m,n,n1:n2) - ksupsmnsf(m,n,n1-1:n2-1))
940 divjmnch(m,n,n1:n2) =
941 m*(ksupumnsf(m,n,n1:n2) + ksupumnsf(m,n,n1-1:n2
942 + n*nfp*(ksupvmnsf(m,n,n1:n2) + ksupvmnsf(m,n,n1-1:n2
943 + ohs*(ksupsmncf(m,n,n1:n2) - ksupsmncf(m,n,n1-1:n2
950 nsmin = max(1, startglobrow - 1)
951 nsmax = min(endglobrow + 1, ns)
952 ALLOCATE(temp1(ntheta,nzeta,nsmin:nsmax),
953 temp2(ntheta,nzeta,nsmin:nsmax), stat=istat)
954 CALL assert_eq(istat, 0,
'Allocation failed in DIVJ')
959 temp1 = ksupsijf0(:,:,nsmin:nsmax)*ksubsijf(:,:,nsmin:nsmax)
960 + ksupuijf0(:,:,nsmin:nsmax)*ksubuijf(:,:,nsmin:nsmax)
961 + ksupvijf0(:,:,nsmin:nsmax)*ksubvijf(:,:,nsmin:nsmax)
963 CALL to_half_mesh(temp1, temp2)
965 temp2 = temp2/(jacobh(:,:,nsmin:nsmax)*jacobh(:,:,nsmin:nsmax))
967 nsmin = max(1, startglobrow)
968 nsmax = min(endglobrow, ns)
969 tnorm = sum(temp2(:,:,nsmin:nsmax))
971 DEALLOCATE(temp1, temp2, stat=istat)
973 divj_rms = sum(divjmnsh(:,:,n1:n2)**2
974 + divjmnch(:,:,n1:n2)**2)
978 CALL mpi_allreduce(mpi_in_place,temp,2,mpi_real8, mpi_sum,
985 IF (tnorm .ge. zero)
THEN
986 divjmnsh = divjmnsh/sqrt(tnorm)
987 divjmnch = divjmnch/sqrt(tnorm)
988 divj_rms = sqrt(divj_rms/tnorm)
992 time_divj = time_divj + (toff - ton)
999 SUBROUTINE dealloc_diagnostics
1001 IF (
ALLOCATED(divbmnsf))
THEN
1002 DEALLOCATE(divbmnsf)
1004 IF (
ALLOCATED(divbmncf))
THEN
1005 DEALLOCATE(divbmncf)
1007 IF (
ALLOCATED(divjmnsh))
THEN
1008 DEALLOCATE(divjmnsh)
1010 IF (
ALLOCATED(divjmnch))
THEN
1011 DEALLOCATE(divjmnch)
1013 IF (
ALLOCATED(bdotjmnch))
THEN
1014 DEALLOCATE(bdotjmnch)
1016 IF (
ALLOCATED(bdotjmnsh))
THEN
1017 DEALLOCATE(bdotjmnsh)