V3FIT
diagnostics_mod.f90
1  MODULE diagnostics_mod
2  USE v3_utilities, ONLY: assert
3  USE fourier, ONLY: f_sum, f_cos, f_sin, f_none, f_du, f_dv
4  USE descriptor_mod, ONLY: iam, nprocs, siesta_comm
5  USE nscalingtools, ONLY: startglobrow, endglobrow, mpi_err
6  USE timer_mod
7  USE quantities
8  USE metrics, ONLY: tolowerh, nsh
9  USE utilities, ONLY: to_half_mesh
10 #if defined(MPI_OPT)
11  USE mpi_inc
12 #endif
13  USE island_params, ONLY: fourier_context
14 
15  IMPLICIT NONE
16 
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
22 !
23 ! The next three allocatable variables are kept in memory and used to
24 ! produce output whenever L_DUMP_DIAGNO = .TRUE.
25 !
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. !Set in call to EVOLVE_BGRADP
34  LOGICAL :: lcurr_init = .false.
35 !
36 ! CALL THIS AFTER THE B-FIELD COMPONENTS bi = Jac*B-sup-i HAVE BEEN UPDATED
37 ! IN UPDATE_STATE (CALL FROM UPDATE_STATE IS SAFEST WAY!)
38 
39  CONTAINS
40 
41 !*******************************************************************************
43 !*******************************************************************************
44 !-------------------------------------------------------------------------------
49 !-------------------------------------------------------------------------------
50  SUBROUTINE divb(ns_min, ns_max)
51  USE island_params, ONLY: ohs=>ohs_i
52  USE hessian, ONLY: gather_array
53 
54 ! Declare Arguments
55  INTEGER, INTENT(in) :: ns_min
56  INTEGER, INTENT(in) :: ns_max
57 
58 ! Local Variables
59  INTEGER :: istat
60  INTEGER :: m
61  INTEGER :: n
62  INTEGER :: nl
63  INTEGER :: nh
64  REAL (dp) :: tnorm
65  REAL (dp), DIMENSION(2) :: temp
66  REAL (dp) :: ton
67  REAL (dp) :: toff
68 
69 ! Start of executable code
70  CALL second0(ton)
71 
72  nsmin = ns_min
73  nsmax = ns_max
74  nl = max(2, nsmin)
75  nh = min(ns - 1, nsmax)
76 
77  IF (ALLOCATED(divbmnsf)) THEN
78  DEALLOCATE(divbmnsf)
79  END IF
80  IF (ALLOCATED(divbmncf)) THEN
81  DEALLOCATE(divbmncf)
82  END IF
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')
86 
87  divbmnsf = 0
88  divbmncf = 0
89 
90  DO m = 0, mpol
91  DO n = -ntor, ntor
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))
98  IF (lasym) THEN
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))
105  END IF
106  END DO
107  END DO
108 
109  tnorm = hs_i*sum(jbsupumnch(:,:,nsmin:nsmax)**2 &
110  + jbsupvmnch(:,:,nsmin:nsmax)**2 &
111  + jbsupsmnsh(:,:,nsmin:nsmax)**2)
112 
113  divb_rms = sum(divbmnsf**2 + divbmncf**2)
114 #if defined(MPI_OPT)
115  temp(1) = divb_rms
116  temp(2) = tnorm
117  CALL mpi_allreduce(mpi_in_place,temp,2,mpi_real8, mpi_sum, &
118  siesta_comm,mpi_err)
119  divb_rms = temp(1)
120  tnorm = temp(2)
121 #endif
122  IF (tnorm .NE. zero) THEN
123  divb_rms = sqrt(divb_rms/tnorm)
124  divbmnsf = divbmnsf/sqrt(tnorm)
125  divbmncf = divbmncf/sqrt(tnorm)
126  END IF
127 
128  CALL second0(toff)
129  time_divb = time_divb + (toff-ton)
130 
131  END SUBROUTINE
132 
133  SUBROUTINE write_profiles(fsq_total)
134  USE safe_open_mod
135  USE siesta_init, ONLY: init_state
136 !-----------------------------------------------
137 ! D u m m y A r g u m e n t s
138 !-----------------------------------------------
139  REAL(dp) :: fsq_total
140 !-----------------------------------------------
141 ! L o c a l V a r i a b l e s
142 !-----------------------------------------------
143  INTEGER, PARAMETER :: ifull=0, ihalf=1
144  INTEGER :: istat, js
145  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: pijh, pmnh
146  REAL(dp) :: average(ns)
147 !-----------------------------------------------
148 !
149 ! IT IS ASSUMED COMING IN HERE THAT THE PROFILES HAVE BEEN PREVIOUSLY
150 ! GATHERED (IN UPDATE_STATE) AND PREPARE_QUANTITIES HAS BEEN CALLED
151 !
152  IF (iam .NE. 0) RETURN
153 
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')
157 
158  CALL safe_open(unit35, istat, "siesta_profiles.txt", 'replace', &
159  'formatted')
160  CALL safe_open(unit36, istat, "siesta_profiles_pest.txt", &
161  'replace','formatted')
162 
163  CALL assert(istat.EQ.0,'Error opening siesta profiles data file')
164 
165 !Check that init_state was initialized so start/end globrow = (1,ns)
166  CALL assert(startglobrow.EQ.1 .AND. endglobrow.EQ.ns, &
167  ' INIT_STATE NOT INITIALIZED IN WRITE_PROFILES!')
168 
169 !USE pijh, pmnh scratch arrays to compute pressure and bsupX components on half-mesh
170 
171 !pressure
172  CALL write_spectra(pijf0, pmnh, f_cos, f_sin, ifull, 'p', fsq_total)
173 
174 !BSUPS
175  CALL write_spectra(bsupsijf0, pmnh, f_sin, f_cos, ifull, 'B^s', fsq_total)
176 
177 !BSUPU
178  CALL write_spectra(bsupuijf0, pmnh, f_cos, f_sin, ifull, 'B^u', fsq_total)
179 
180 !BSUPV
181  CALL write_spectra(bsupvijf0, pmnh, f_cos, f_sin, ifull, 'B^v', fsq_total)
182 
183 !KSUPS CRCook
184  pijh = ksupsijf0/jacobf
185  CALL write_spectra(pijh, pmnh, f_sin, f_cos, ifull, 'J^s', fsq_total)
186 
187 !KSUPU CRCook
188  pijh = ksupuijf0/jacobf
189  CALL write_spectra(pijh, pmnh, f_cos, f_sin, ifull, 'J^u', fsq_total)
190 
191 !KSUPV CRCook
192  pijh = ksupvijf0/jacobf
193  CALL write_spectra(pijh, pmnh, f_cos, f_sin, ifull, 'J^v', fsq_total)
194 
195 !K|| = (J*B)/B**2 (SPH:040915)
196  CALL bdotj(bdotjijh)
197  pijh = bdotjijh*jacobh
198  DEALLOCATE(bdotjijh)
199  CALL write_spectra(pijh, pmnh, f_cos, f_sin, ihalf, 'K||', fsq_total)
200 
201 
202 
203 !need to extend this to lasym=T
204  CALL fourier_context%tomnsp(pijh, pmnh, f_cos)
205  pmnh(:,0,:) = 0
206  CALL fourier_context%toijsp(pmnh, pijh, f_none, f_cos)
207  CALL surfaverage(average, abs(pijh), 1, ns)
208  WRITE (unit36, *)
209  WRITE (unit36, *) '<ABS(K||) - nonaxisymmetric>'
210  DO js = 2,ns
211  WRITE (unit36, '(i4,1p,e10.2)') js, average(js)
212  END DO
213 
214 !Displacements X jac (SPH:052114)
215  CALL siesta_profiles(unit35, jvsupsmncf, 'vel^s_cosmn; fsq_total: ', fsq_total)
216  IF (lasym) CALL siesta_profiles(unit35, jvsupsmnsf, 'vel^s_sinmn; fsq_total: ', fsq_total)
217  CALL siesta_profiles(unit35, jvsupumnsf, 'vel^u_sinmn; fsq_total: ', fsq_total)
218  IF (lasym) CALL siesta_profiles(unit35, jvsupumncf, 'vel^u_cosmn; fsq_total: ', fsq_total)
219  CALL siesta_profiles(unit35, jvsupvmnsf, 'vel^v_sinmn; fsq_total: ', fsq_total)
220  IF (lasym) CALL siesta_profiles (unit35, jvsupvmncf, 'vel^v_cosmn; fsq_total: ', fsq_total)
221 
222  DEALLOCATE(pijh, pmnh)
223 
224  CLOSE (unit35)
225  CLOSE (unit36)
226 
227  END SUBROUTINE write_profiles
228 
229 !-------------------------------------------------------------------------------
244  SUBROUTINE write_spectra(xij, xmn, fsym, fasym, igrid, label, fsq_total)
245  USE vmec_info, ONLY: lmns => lmns_i, lmnc => lmnc_i
246  USE siesta_namelist, ONLY: l_vessel
247  USE utilities, ONLY: to_full_mesh
248 
249  IMPLICIT NONE
250 
251 ! Declare Arguments
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
259 
260 ! Local variables
261  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: xijf
262  CHARACTER (LEN=128) :: mylabel
263 
264 ! Start of executable code.
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)
268  xij = xijf/jacobf
269  DEALLOCATE(xijf)
270  END IF
271 
272  IF (fsym .EQ. f_cos) THEN
273  WRITE (mylabel,1000) trim(label)
274  ELSE
275  WRITE (mylabel,1001) trim(label)
276  END IF
277 
278  CALL fourier_context%tomnsp(xij, xmn, fsym)
279  CALL siesta_profiles(unit35, xmn, trim(mylabel), fsq_total)
280  IF (.not.l_vessel) THEN
281  CALL fourier_context%tomnsp(xij, xmn, lmns, lmnc, fsym, lasym)
282  CALL siesta_profiles(unit36, xmn, trim(mylabel), fsq_total)
283  END IF
284  IF (lasym) THEN
285 
286  IF (fasym .EQ. f_cos) THEN
287  WRITE (mylabel,1000) trim(label)
288  ELSE
289  WRITE (mylabel,1001) trim(label)
290  END IF
291 
292  CALL fourier_context%tomnsp(xij, xmn, fasym)
293  CALL siesta_profiles(unit35, xmn, trim(mylabel), fsq_total)
294  IF (.not.l_vessel) THEN
295  CALL fourier_context%tomnsp(xij, xmn, lmns, lmnc, fasym, lasym)
296  CALL siesta_profiles(unit36, xmn, trim(mylabel), fsq_total)
297  END IF
298  END IF
299 
300 1000 FORMAT(a,'_cosmn; fsq_total: ')
301 1001 FORMAT(a,'_sinmn; fsq_total: ')
302 
303  END SUBROUTINE
304 
305  SUBROUTINE siesta_profiles (iunit, arr_value, label, fsq_total)
306 !-----------------------------------------------
307 ! D u m m y A r g u m e n t s
308 !-----------------------------------------------
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
313 !-----------------------------------------------
314 ! L o c a l V a r i a b l e s
315 !-----------------------------------------------
316  INTEGER :: n, m, js
317  INTEGER :: mnmax
318  REAL(dp) :: sj, rms_value
319  CHARACTER (LEN=35) :: p_format
320  INTEGER, DIMENSION(:,:), ALLOCATABLE :: m_array
321  INTEGER, DIMENSION(:,:), ALLOCATABLE :: n_array
322 !-----------------------------------------------
323 
324  WRITE (iunit, *)
325  WRITE (iunit, '(1x,a,1p,e12.4)') trim(label), fsq_total
326 
327  ALLOCATE(m_array(0:mpol,-ntor:ntor))
328  ALLOCATE(n_array(0:mpol,-ntor:ntor))
329 
330  DO n = -ntor, ntor
331  DO m = 0, mpol
332  m_array(m,n) = m
333  n_array(m,n) = n
334  END DO
335  END DO
336 
337  mnmax = (mpol + 1)*(2*ntor + 1)
338 
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
343 
344  WRITE (p_format,1002) mnmax
345  DO js = 2, ns
346  sj = hs_i*(js - 1)
347  rms_value = sqrt(sum(arr_value(:,:,js)**2)/mnmax)
348  WRITE(iunit,p_format) sj, rms_value, arr_value(:,:,js)
349  END DO
350 
351  DEALLOCATE(m_array, n_array)
352 
353 1000 FORMAT('(a,8x,a,6x,a,',i3,'(2x,i12))')
354 1002 FORMAT('(f6.3,2x,es12.5,14x',i4,'(2x,es12.5))')
355 
356  END SUBROUTINE siesta_profiles
357 
358 
359  SUBROUTINE bgradp(ns_min, ns_max)
360 !
361 ! WRITTEN 03-21-13 BY S. HIRSHMAN AS PART OF THE ORNL SIESTA PROJECT (c)
362 !
363 ! PURPOSE: Computes B*GRAD P (normalized to B*P) on the FULL mesh
364 !
365 !-----------------------------------------------
366 ! D u m m y A r g u m e n t s
367 !-----------------------------------------------
368  INTEGER, INTENT(IN) :: ns_min, ns_max
369 !-----------------------------------------------
370 ! L o c a l V a r i a b l e s
371 !-----------------------------------------------
372  INTEGER :: istat
373  REAL(dp) :: s1(2), r1(2), ton, toff
374 !-----------------------------------------------
375 !
376 ! COMPUTE B DOT GRAD P AT FULL-GRID POINTS (sans endpts)
377 !
378  nsmin=ns_min; nsmax=ns_max
379 
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)
383 
384  DEALLOCATE (bgradpf, stat=istat)
385 
386  END SUBROUTINE bgradp
387 
388 !-------------------------------------------------------------------------------
393 !-------------------------------------------------------------------------------
394  SUBROUTINE get_bgradp(ns_min, ns_max)
395 
396  IMPLICIT NONE
397 
398 ! Declare Arguments
399  INTEGER, INTENT(in) :: ns_min
400  INTEGER, INTENT(in) :: ns_max
401 
402 ! local variables
403  INTEGER :: istat
404  INTEGER :: n1
405  INTEGER :: n2
406  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: dpduijf
407  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: dpdvijf
408  REAL (dp), DIMENSION(2) :: temp(2)
409  REAL (dp) :: ton
410  REAL (dp) :: toff
411 
412 ! Start of executable code
413  CALL second0(ton)
414 
415  CALL assert(ALLOCATED(jbsupsmnsh),'jbsupXmn unallocated in get_bgradp')
416 
417  nsmin = ns_min
418  nsmax = min(ns_max + 1, ns)
419 
420  bgradpf = 0
421 
422 ! Compute pressure and angle derivatives
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')
426 
427  CALL to_full_mesh(pijh0_du, dpduijf)
428  CALL to_full_mesh(pijh0_dv, dpdvijf)
429 
430  nsmax = min(endglobrow, ns)
431  n1 = max(2, nsmin)
432  n2 = min(nsh, nsmax)
433 ! Get all values except first and
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)
437 
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))
442 
443 ! Ignore inner s<0.1 for this diagnostic
444  n1 = max(ns/20, nsmin)
445  wbgradp = sum(bgradpf(:,:,n1:n2)**2*wint(:,:,n1:n2))
446  temp(1) = tnorm_bgradp
447  temp(2) = wbgradp
448  IF (parsolver) THEN
449 #if defined(MPI_OPT)
450  CALL mpi_allreduce(mpi_in_place, temp, 2, mpi_real8, mpi_sum, &
451  siesta_comm, mpi_err)
452 #endif
453  END IF
454  tnorm_bgradp = temp(1)
455  wbgradp = temp(2)
456 
457  IF (tnorm_bgradp .GT. zero) THEN
458  bgradp_rms = sqrt(wbgradp/tnorm_bgradp)
459  tnorm_bgradp = sqrt(tnorm_bgradp/ns)
460 
461 ! Ignore first point
462  max_bgradp = maxval(bgradpf(:,:,n1:n2))/tnorm_bgradp
463  min_bgradp = minval(bgradpf(:,:,n1:n2))/tnorm_bgradp
464  temp(1) = max_bgradp
465  temp(2) = -min_bgradp
466  IF (parsolver) THEN
467 #if defined(MPI_OPT)
468  CALL mpi_allreduce(mpi_in_place, temp, 2, mpi_real8, mpi_max, &
469  siesta_comm, mpi_err)
470 #endif
471  END IF
472  max_bgradp = temp(1)
473  min_bgradp = -temp(2)
474  IF (.not.lasym) THEN
475 ! Reflect onto u = pi - 2*pi: bgradp has sin parity.
476  max_bgradp = max(max_bgradp, -min_bgradp)
477  min_bgradp = min(min_bgradp, -max_bgradp)
478  END IF
479  END IF
480 
481  DEALLOCATE (dpdvijf, dpduijf, stat=istat)
482 
483  CALL second0(toff)
484  time_bgradp = time_bgradp + (toff - ton)
485 
486  END SUBROUTINE
487 
488 !-------------------------------------------------------------------------------
490 !-------------------------------------------------------------------------------
491  SUBROUTINE beta
492 
493  avbeta = wp/wb
494 
495  END SUBROUTINE
496 
497 
498  SUBROUTINE tflux
499  USE fourier, ONLY: m0, n0
500 
501 !-----------------------------------------------
502 ! L o c a l V a r i a b l e s
503 !-----------------------------------------------
504  REAL(dp) :: part_sum
505 !-----------------------------------------------
506  nsmin=max(2,startglobrow); nsmax=min(endglobrow,ns)
507 
508 !Averages over all toroidal cross sections (which should be the same)
509  toroidal_flux = sum(jbsupvmnch(m0,n0,nsmin:nsmax))
510  IF (parsolver) THEN
511 #if defined(MPI_OPT)
512  CALL mpi_allreduce(mpi_in_place,toroidal_flux,1,mpi_real8, mpi_sum, &
513  siesta_comm,mpi_err)
514 #endif
515 ! ELSE
516 ! toroidal_flux = part_sum
517  END IF
518  toroidal_flux = signjac*twopi*toroidal_flux*hs_i/b_factor
519 
520  IF (toroidal_flux0 .EQ. zero) toroidal_flux0 = toroidal_flux
521 
522  END SUBROUTINE tflux
523 
524 
525  SUBROUTINE bdotj (bdotjijhA)
526 !
527 ! WRITTEN 08-20-07 BY R. SANCHEZ AS PART OF THE ORNL SIESTA PROJECT (c)
528 !
529 ! PURPOSE: Computes parallel current normalized to the total current (J*B/|JxB|) on the half mesh
530 !-----------------------------------------------
531 ! D u m m y A r g u m e n t
532 !-----------------------------------------------
533  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: bdotjijhA
534 !-----------------------------------------------
535 ! L o c a l V a r i a b l e s
536 !-----------------------------------------------
537  INTEGER :: istat
538  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: &
539 ! bsupsijh, bsupuijh, bsupvijh, &
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
544 !-----------------------------------------------
545  CALL second0(ton)
546  CALL assert(SIZE(bsupsijh0,3).EQ.ns, 'bsupsijh0 wrong size in bdotj')
547 !
548 ! Compute covariant field. Used for normalization.
549 !
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
555 
556  CALL tolowerh(bsupsijh0, bsupuijh0, bsupvijh0, &
557  bsubsijh, bsubuijh, bsubvijh, 1, ns)
558 
559 !
560 ! In initial state, calculate currents (remember they are multiplied by a jacobian!)
561 !
562 !
563 ! Move contravariant currents (*jacob) to the half mesh
564 !
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')
569 
570  CALL to_half_mesh(ksupsijf0, ksupsijh)
571  CALL to_half_mesh(ksupuijf0, ksupuijh)
572  CALL to_half_mesh(ksupvijf0, ksupvijh)
573 !
574 ! Remove jacobian
575 !
576  ksupsijh = ksupsijh/jacobh
577  ksupuijh = ksupuijh/jacobh
578  ksupvijh = ksupvijh/jacobh
579 !
580 ! Get covariant currents
581 !
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
587 
588  CALL tolowerh(ksupsijh, ksupuijh, ksupvijh, &
589  ksubsijh, ksubuijh, ksubvijh, 1, ns)
590 !
591 ! Compute B*J in the half mesh.
592 !
593  IF (ALLOCATED(bdotjmnch) .AND. SIZE(bdotjmnch,3).NE. ns) &
594  DEALLOCATE(bdotjmnch)
595  IF (.NOT. ALLOCATED(bdotjmnch)) THEN ! Should already be allocated in QUANTITIES
596  ALLOCATE(bdotjmnch(0:mpol,-ntor:ntor,ns), stat=istat) ! Otherwise, allocate first time called
597  CALL assert(istat.EQ.0,'Allocation #5 failed in BDOTJ')
598  ENDIF
599  IF (lasym) THEN
600  IF (ALLOCATED(bdotjmnsh) .AND. SIZE(bdotjmnsh,3).NE. ns) &
601  DEALLOCATE(bdotjmnsh)
602  IF (.NOT. ALLOCATED(bdotjmnsh)) THEN ! Should already be allocated in QUANTITIES
603  ALLOCATE(bdotjmnsh(0:mpol,-ntor:ntor,ns), stat=istat) ! Otherwise, allocate first time called
604  CALL assert(istat.EQ.0,'Allocation #5 failed in BDOTJ')
605  ENDIF
606  END IF
607  ALLOCATE(bdotjijha(ntheta,nzeta,ns), stat=istat)
608  CALL assert(istat.EQ.0, 'Allocation #6 failed in BDOTJ')
609  bdotjijha=0
610 
611  bdotjijha(:,:,3:nsh-1) = & ! RS: If I include 2 and nsh, it gets huge
612  bsupsijh0(:,:,3:nsh-1)*ksubsijh(:,:,3:nsh-1) + & ! something funny happens at the boundaries.
613  bsupuijh0(:,:,3:nsh-1)*ksubuijh(:,:,3:nsh-1) + &
614  bsupvijh0(:,:,3:nsh-1)*ksubvijh(:,:,3:nsh-1) ! B*J = (B^s*J_s + B^u*J_u + B^v*J_v)
615 !
616 ! Compute tnorm = |JXB|^2 for normalization; thus we print out J_parallel/|J_perp|.
617 !
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), &
621  stat=istat)
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
625 
626  work1 = (bsubuijh*ksubvijh - bsubvijh*ksubuijh)/jacobh ! (JxB)^s = (1/sqrt(g))*(K_uB_b-K_vB_u)
627  work2 = (bsubvijh*ksubsijh - bsubsijh*ksubvijh)/jacobh ! (JxB)^u = .....
628  work3 = (bsubsijh*ksubuijh - bsubuijh*ksubsijh)/jacobh ! (JxB)^v = .....
629  CALL tolowerh(work1, work2, work3, work4, work5, work6, 1, ns)
630 
631  work1 = work1*work4 + work2*work5 + work3*work6 ! |JxB|**2
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)
637 
638  work1 = bsupsijh0*bsubsijh + bsupuijh0*bsubuijh + & ! |B|**2
639  bsupvijh0*bsubvijh
640  work2 = ksupsijh*ksubsijh + ksupuijh*ksubuijh + & ! |J|**2
641  ksupvijh*ksubvijh
642 
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)) ! |J|*|B|
650  bdotj_rms = sum(bdotjijha*bdotjijha)
651 
652  DEALLOCATE (bsubsijh, bsubuijh, bsubvijh, ksupsijh, ksupuijh, &
653  ksupvijh, ksubsijh, ksubuijh, ksubvijh, &
654  work1, work2)
655 
656  tnorm = max(tnorm, epsilon(tnorm))
657  tnorm2= max(tnorm2,epsilon(tnorm2))
658  bdotj2_rms = sqrt(abs(bdotj_rms))/tnorm2 ! RMS of bdotj/|J|*|B|
659  bdotj_rms = sqrt(abs(bdotj_rms))/tnorm ! RMS of bdotj/|JxB|
660  CALL fourier_context%tomnsp(bdotjijha, bdotjmnch, f_cos) ! Keep harmonics of BDOTJ for output
661  bdotjmnch = bdotjmnch/tnorm
662  IF (lasym) THEN
663  CALL fourier_context%tomnsp(bdotjijha, bdotjmnsh, f_cos) ! Keep harmonics of BDOTJ for output
664  bdotjmnsh = bdotjmnsh/tnorm
665  END IF
666  CALL second0(toff)
667  time_bdotj = time_bdotj+(toff-ton)
668 
669  bdotjijh = bdotjijha
670 
671  END SUBROUTINE bdotj
672 
673  SUBROUTINE bdotj_par
674  USE island_params, ONLY: fourier_context
675 
676  IMPLICIT NONE
677 
678 !
679 ! WRITTEN 08-20-07 BY R. SANCHEZ AS PART OF THE ORNL SIESTA PROJECT (c)
680 ! UPDATED 03-21-13 BY S. HIRSHMAN FOR PARALLEL EXECUTION
681 !
682 ! PURPOSE: Computes parallel current normalized to the total current (J*B/|JxB|) on the half mesh
683 !
684 !-----------------------------------------------
685 ! L o c a l V a r i a b l e s
686 !-----------------------------------------------
687  INTEGER :: istat, n1, n2, ns_span
688  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: bsupsijh, &
689  bsupuijh, bsupvijh, &
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
694 !-----------------------------------------------
695  CALL second0(ton)
696 
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
700 !
701 ! Compute contravariant components of magnetic field
702 !
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')
707 
708  CALL fourier_context%toijsp(jbsupsmnsh(:,:,nsmin:nsmax), bsupsijh, &
709  & f_none, f_sin)
710  CALL fourier_context%toijsp(jbsupumnch(:,:,nsmin:nsmax), bsupuijh, &
711  f_none, f_cos)
712  CALL fourier_context%toijsp(jbsupvmnch(:,:,nsmin:nsmax), bsupvijh, &
713  f_none, f_cos)
714  IF (lasym) THEN
715  CALL fourier_context%toijsp(jbsupsmnch(:,:,nsmin:nsmax), bsupsijh, &
716  f_sum, f_cos)
717  CALL fourier_context%toijsp(jbsupumnsh(:,:,nsmin:nsmax), bsupuijh, &
718  f_sum, f_sin)
719  CALL fourier_context%toijsp(jbsupvmnsh(:,:,nsmin:nsmax), bsupvijh, &
720  f_sum, f_sin)
721  END IF
722 
723  bsupsijh = bsupsijh/jacobh(:,:,nsmin:nsmax)
724  bsupuijh = bsupuijh/jacobh(:,:,nsmin:nsmax)
725  bsupvijh = bsupvijh/jacobh(:,:,nsmin:nsmax)
726 !
727 ! Compute covariant field. Used for normalization.
728 !
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')
733 
734  CALL tolowerh(bsupsijh, bsupuijh, bsupvijh, &
735  bsubsijh, bsubuijh, bsubvijh, nsmin, nsmax)
736 !
737 ! In initial state, calculate currents (remember they are multiplied by a jacobian!)
738 !
739 !
740 ! Move contravariant currents (*jacob) to the half mesh
741 !
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')
747 
748  CALL to_half_mesh(ksupsijf0, ksupsijh)
749  CALL to_half_mesh(ksupuijf0, ksupuijh)
750  CALL to_half_mesh(ksupvijf0, ksupvijh)
751 
752 !
753 ! Remove jacobian
754 !
755  ksupsijh = ksupsijh/jacobh(:,:,nsmin:nsmax)
756  ksupuijh = ksupuijh/jacobh(:,:,nsmin:nsmax)
757  ksupvijh = ksupvijh/jacobh(:,:,nsmin:nsmax)
758 
759 !
760 ! Get covariant currents
761 !
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')
766 
767  CALL tolowerh(ksupsijh, ksupuijh, ksupvijh, &
768  ksubsijh, ksubuijh, ksubvijh, nsmin, nsmax)
769  nsmin=max(1,startglobrow)
770 
771 !
772 ! Compute B*J in the half mesh.
773 !
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')
779  ENDIF
780  IF (lasym) THEN
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')
786  ENDIF
787  ENDIF
788 
789  ALLOCATE(bdotjijh(ntheta,nzeta,nsmin:nsmax), stat=istat)
790  CALL assert(istat.EQ.0,'Allocation #6 failed in BDOTJ_PAR')
791  bdotjijh=0
792 
793  n1 = max(3,nsmin); n2 = min(nsh-1, nsmax)
794  bdotjijh(:,:,n1:n2) = & ! RS: If include 2 and nsh, it gets huge
795  bsupsijh(:,:,n1:n2)*ksubsijh(:,:,n1:n2) + & ! something funny happens at the boundaries.
796  bsupuijh(:,:,n1:n2)*ksubuijh(:,:,n1:n2) + &
797  bsupvijh(:,:,n1:n2)*ksubvijh(:,:,n1:n2) ! B*J = (B^s*J_s + B^u*J_u + B^v*J_v)
798 
799 !
800 ! Compute tnorm = |JXB|^2 for normalization; thus we print out J_parallel/|J_perp|.
801 !
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), &
808  stat=istat)
809  CALL assert(istat.EQ.0,'Allocation #7 failed in BDOTJ_PAR')
810 
811  work1 = bsubuijh*ksubvijh(:,:,nsmin:nsmax) &
812  - bsubvijh*ksubuijh(:,:,nsmin:nsmax) ! (JxB)^s = (1/sqrt(g))*(K_uB_b-K_vB_u)
813  work2 = bsubvijh*ksubsijh(:,:,nsmin:nsmax) &
814  - bsubsijh*ksubvijh(:,:,nsmin:nsmax) ! (JxB)^u = .....
815  work3 = bsubsijh*ksubuijh(:,:,nsmin:nsmax) &
816  - bsubuijh*ksubsijh(:,:,nsmin:nsmax) ! (JxB)^v = .....
817  work1 = work1/jacobh(:,:,nsmin:nsmax)
818  work2 = work2/jacobh(:,:,nsmin:nsmax)
819  work3 = work3/jacobh(:,:,nsmin:nsmax)
820  CALL tolowerh(work1, work2, work3, &
821  work4, work5, work6, nsmin, nsmax)
822 
823  work1 = work1*work4 + work2*work5 + work3*work6 ! |JxB|**2
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)
829 
830  work1 = bsupsijh*bsubsijh + bsupuijh*bsubuijh + & ! |B|**2
831  bsupvijh*bsubvijh
832  work2 = ksupsijh(:,:,nsmin:nsmax)*ksubsijh(:,:,nsmin:nsmax) &
833  + ksupuijh(:,:,nsmin:nsmax)*ksubuijh(:,:,nsmin:nsmax) &
834  + ksupvijh(:,:,nsmin:nsmax)*ksubvijh(:,:,nsmin:nsmax) ! |J|**2
835 
836  IF (nsmin .EQ. 1) THEN
837  work1(:,:,1) = 0; work2(:,:,1) = 0
838  END IF
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))
845 #if defined(MPI_OPT)
846  temp(1)=tnorm; temp(2)=tnorm2; temp(3)=bdotj_rms
847  CALL mpi_allreduce(mpi_in_place,temp,3,mpi_real8, mpi_sum, &
848  siesta_comm,mpi_err)
849  tnorm = sqrt(abs(temp(1)))
850  tnorm2 = sqrt(abs(temp(2))) ! |J|*|B|
851  bdotj_rms = temp(3)
852 #else
853  tnorm = sqrt(abs(tnorm)); tnorm2 = sqrt(abs(tnorm2))
854 #endif
855  DEALLOCATE (bsupsijh, bsupuijh, bsupvijh, bsubsijh, &
856  bsubuijh, bsubvijh, ksupsijh, ksupuijh, &
857  ksupvijh, ksubsijh, ksubuijh, ksubvijh, &
858  work1, work2)
859 
860 
861  tnorm2= max(tnorm2,epsilon(tnorm2))
862  bdotj2_rms = sqrt(abs(bdotj_rms))/tnorm2 ! RMS of bdotj/|J|*|B|
863  bdotj_rms = sqrt(abs(bdotj_rms))/tnorm ! RMS of bdotj/|JxB|
864 
865  bdotjijh = bdotjijh/tnorm
866 
867  CALL fourier_context%tomnsp(bdotjijh, bdotjmnch, f_cos) ! Keep harmonics of BDOTJ for output
868  IF (lasym) THEN
869  CALL fourier_context%tomnsp(bdotjijh, bdotjmnsh, f_sin) ! Keep harmonics of BDOTJ for output
870  END IF
871 
872  DEALLOCATE (bdotjijh)
873 
874  CALL second0(toff)
875  time_bdotj = time_bdotj+(toff-ton)
876 
877  END SUBROUTINE bdotj_par
878 
879 !-------------------------------------------------------------------------------
887 !-------------------------------------------------------------------------------
888  SUBROUTINE divj(ns_min, ns_max)
889  USE island_params, ONLY: ohs=>ohs_i
890 
891  IMPLICIT NONE
892 
893 ! Declare Arguments
894  INTEGER, INTENT(in) :: ns_min
895  INTEGER, INTENT(in) :: ns_max
896 
897 ! Local variables
898  INTEGER :: m
899  INTEGER :: n
900  INTEGER :: istat
901  INTEGER :: n1
902  INTEGER :: n2
903  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: temp1
904  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: temp2
905  REAL(dp) :: tnorm
906  REAL(dp), DIMENSION(2) :: temp
907  REAL(dp) :: ton
908  REAL(dp) :: toff
909 
910 ! Start of executable code.
911  CALL second0(ton)
912 
913  CALL assert(lcurr_init,'MUST CALL init_state(.TRUE.) BEFORE divJ_par')
914 
915  nsmin = ns_min
916  nsmax = ns_max
917 
918  IF (ALLOCATED(divjmnsh)) THEN
919  DEALLOCATE(divjmnsh)
920  END IF
921  IF (ALLOCATED(divjmnch)) THEN
922  DEALLOCATE(divjmnch)
923  END IF
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')
927  divjmnsh = 0
928  divjmnch = 0
929 
930 ! Compute divergence of J. It is a SINE parity quantity.
931  n1 = max(3, nsmin)
932  n2 = min(ns-1, nsmax)
933  DO m = 0, mpol
934  DO n = -ntor,ntor
935  divjmnsh(m,n,n1:n2) = &
936  -m*(ksupumncf(m,n,n1:n2) + ksupumncf(m,n,n1-1:n2-1))/2 &
937  - n*nfp*(ksupvmncf(m,n,n1:n2) + ksupvmncf(m,n,n1-1:n2-1))/2 &
938  + ohs*(ksupsmnsf(m,n,n1:n2) - ksupsmnsf(m,n,n1-1:n2-1))
939  IF (lasym) THEN
940  divjmnch(m,n,n1:n2) = &
941  m*(ksupumnsf(m,n,n1:n2) + ksupumnsf(m,n,n1-1:n2-1))/2 &
942  + n*nfp*(ksupvmnsf(m,n,n1:n2) + ksupvmnsf(m,n,n1-1:n2-1))/2 &
943  + ohs*(ksupsmncf(m,n,n1:n2) - ksupsmncf(m,n,n1-1:n2-1))
944  END IF
945  END DO
946  END DO
947 !
948 ! Compute covariant components of current (times jacobian). Used for normalization.
949 !
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')
955  temp1 = 0
956  temp2 = 0
957 
958 ! Norm: |J|^2
959  temp1 = ksupsijf0(:,:,nsmin:nsmax)*ksubsijf(:,:,nsmin:nsmax) &
960  + ksupuijf0(:,:,nsmin:nsmax)*ksubuijf(:,:,nsmin:nsmax) &
961  + ksupvijf0(:,:,nsmin:nsmax)*ksubvijf(:,:,nsmin:nsmax)
962 
963  CALL to_half_mesh(temp1, temp2)
964 
965  temp2 = temp2/(jacobh(:,:,nsmin:nsmax)*jacobh(:,:,nsmin:nsmax))
966 
967  nsmin = max(1, startglobrow)
968  nsmax = min(endglobrow, ns)
969  tnorm = sum(temp2(:,:,nsmin:nsmax))
970 
971  DEALLOCATE(temp1, temp2, stat=istat)
972 
973  divj_rms = sum(divjmnsh(:,:,n1:n2)**2 &
974  + divjmnch(:,:,n1:n2)**2)
975 #if defined(MPI_OPT)
976  temp(1) = divj_rms
977  temp(2) = tnorm
978  CALL mpi_allreduce(mpi_in_place,temp,2,mpi_real8, mpi_sum, &
979  siesta_comm,mpi_err)
980  divj_rms = temp(1)
981  tnorm = temp(2)
982 #endif
983 
984 ! Compute rms of divergence of J
985  IF (tnorm .ge. zero) THEN
986  divjmnsh = divjmnsh/sqrt(tnorm)
987  divjmnch = divjmnch/sqrt(tnorm)
988  divj_rms = sqrt(divj_rms/tnorm)
989  END IF
990 
991  CALL second0(toff)
992  time_divj = time_divj + (toff - ton)
993 
994  END SUBROUTINE
995 
996 !-------------------------------------------------------------------------------
998 !-------------------------------------------------------------------------------
999  SUBROUTINE dealloc_diagnostics
1000 
1001  IF (ALLOCATED(divbmnsf)) THEN
1002  DEALLOCATE(divbmnsf)
1003  END IF
1004  IF (ALLOCATED(divbmncf)) THEN
1005  DEALLOCATE(divbmncf)
1006  END IF
1007  IF (ALLOCATED(divjmnsh)) THEN
1008  DEALLOCATE(divjmnsh)
1009  END IF
1010  IF (ALLOCATED(divjmnch)) THEN
1011  DEALLOCATE(divjmnch)
1012  END IF
1013  IF (ALLOCATED(bdotjmnch)) THEN
1014  DEALLOCATE(bdotjmnch)
1015  END IF
1016  IF (ALLOCATED(bdotjmnsh)) THEN
1017  DEALLOCATE(bdotjmnsh)
1018  END IF
1019 
1020  END SUBROUTINE
1021 
1022  END MODULE
siesta_namelist::l_vessel
logical l_vessel
If extended grid is to be used using an available vessel file.
Definition: siesta_namelist.f90:140
hessian::gather_array
Definition: hessian.f90:31
quantities
This file contains subroutines for allocating and initializing curvilinear magnetic covariant and pre...
Definition: quantities.f90:11
v3_utilities::assert_eq
Definition: v3_utilities.f:62
island_params::fourier_context
type(fourier_class), pointer fourier_context
Fourier transform object.
Definition: island_params.f90:76
fourier
Module contains subroutines for computing FFT on parallel radial intervals. Converts quantities betwe...
Definition: fourier.f90:13
siesta_init::init_state
subroutine, public init_state(lcurrent_only, lpar_in)
Initialize equilibrium state.
Definition: siesta_init.f90:45
fourier::m0
integer, parameter m0
m = 0 mode.
Definition: fourier.f90:58
v3_utilities::assert
Definition: v3_utilities.f:55
fourier::f_sin
integer, parameter f_sin
Sine parity.
Definition: fourier.f90:27
mpi_inc
Umbrella module avoid multiple inlcudes of the mpif.h header.
Definition: mpi_inc.f:11
island_params
This file contains fix parameters related to the computational grids.
Definition: island_params.f90:10
fourier::n0
integer, parameter n0
n = 0 mode.
Definition: fourier.f90:65
fourier::f_du
integer, parameter f_du
Poloidal derivative flag.
Definition: fourier.f90:36
metrics::tolowerh
subroutine tolowerh(xsupsij, xsupuij, xsupvij,
Converts to covariant component half grid.
Definition: metrics.f90:539
fourier::f_sum
integer, parameter f_sum
Sum fouier real space transformed quantity. This is used when a non-stellarator symmetric parity is b...
Definition: fourier.f90:43
island_params::ohs_i
real(dp) ohs_i
Radial grid derivative factor. d X/ ds = (x_i+1 - x_i)/(s_i+1 - s_i) (1) Where ohs = 1/(s_i+1 - s_i) ...
Definition: island_params.f90:51
metrics
Module for reading in VMEC input and computing metric elements on the SIESTA (sqrt-flux) radial grid....
Definition: metrics.f90:13
fourier::f_cos
integer, parameter f_cos
Cosine parity.
Definition: fourier.f90:25
siesta_namelist
This file contains all the variables and maximum sizes of the inputs for a SIESTA namelist input file...
Definition: siesta_namelist.f90:103
siesta_init
Initializes unperturbed siesta fields and pressure in real space.
Definition: siesta_init.f90:10
fourier::f_dv
integer, parameter f_dv
Toroidal derivative flag.
Definition: fourier.f90:38
fourier::f_none
integer, parameter f_none
Do not sum fouier real space transformed quantity. This is used when a stellarator symmetric parity i...
Definition: fourier.f90:32
Generated on Thu Mar 5 2020 15:49:25 for V3FIT by   doxygen 1.8.17