V3FIT
dumping_mod.f90
1  MODULE dumping_mod
2 !
3 ! WRITTEN: 11-31-06 BY R. SANCHEZ AS PART OF THE ORNL SIESTA PROJECT (c)
4 ! Updated: 12-01-10 by R. Sanchez
5 ! Updated: 04-18-13 by S. Hirshman (parallelized)
6 !
7 ! PURPOSE: COMPUTES QUANTITIES for OUTPUT on R-Z-Phi using a VMEC-uniform mesh.
8 ! Dimensions of the mesh is set by: NSS, NUS and NPHIS in PERTURBATIONS module.
9 ! Requires that L_TRACING or L_SILO_OUTPUT set to true to dump it.
10 ! L_TRACING = T produces a BFIELD ASCII file needed for Poincare plots with XPOINCARE.
11 ! L_SILO_OUTPUT = T requires linking against the SILO library and
12 ! produces 2D and 3D data in SILO format.
13 ! FUTURE work: upgrade XPOINCARE to read in SILO output and get rid of BFIELD ASCII file.
14 !
15 ! Legacy output in PGPLOT format also exists in both R-Z and
16 ! s-Theta uniform meshes for pressure. If R-Z requested, use NRS,
17 ! NZS and NVS to set dimensions of CYL mesh and NSS, NUS
18 ! and NVS in s-Theta. This is the default output. Use NPHIS to set number of
19 ! toroidal planes in pressure output.
20 
21 ! QUANTITIES are provided via their XMN harmonics, together with their parity
22 ! (IPARITY=0, COSINE or =1, SINE) and radial gridding info (IHALF=0, full; =1, half mesh).
23 !
24  USE v3_utilities, ONLY: assert
25  USE stel_kinds
26  USE stel_constants, ONLY: twopi, one, zero
27  USE shared_data, ONLY: lverbose
28  USE vmec_utils
29  USE island_params, ns => ns_i, mpol => mpol_i, ntor => ntor_i, &
30  ohs => ohs_i, nfp => nfp_i, ntheta => nu_i, nzeta => nv_i
31  USE metrics, rmax_v => rmax, rmin_v => rmin, zmax_v => zmax, &
32  zmin_v => zmin
33  USE vmec_info, ONLY: rmnc_i, rmns_i, zmnc_i, zmns_i
34  USE siesta_namelist, ONLY: nphis1=>nphis, nrs1=>nrs, nzs1=>nzs, & ! NPHIS = number of toroidal planes (CYL or VMEC)
35  nus1=>nus, nss1=>nss, nvs1=>nvs, restart_ext, & ! NSS, NUS = resolution(VMEC) in PHI-plane
36  l_tracing1=>l_tracing, & ! NRS, NZS = resolution(CYL) in PHI-plane
37  l_silo_output1=>l_silo_output, l_silo3d1=> l_silo3d
38 #if defined(MPI_OPT)
39  USE nscalingtools, ONLY: mpi_status
40  USE mpi_inc
41 #endif
42  USE descriptor_mod, ONLY: iam, nprocs, siesta_comm
43  USE fourier, ONLY: f_ds, f_du, f_dv, f_cos, f_sin, f_none, f_sum
44  USE utilities, ONLY: to_half_mesh
45  IMPLICIT NONE
46  PRIVATE
47 #if defined(SILO_AVAIL)
48  include "silo.inc"
49 #endif
50 !-----------------------------------------------
51  LOGICAL :: l_tracing
52 !-----------------------------------------------
53 ! Grid flags.
54  INTEGER, PARAMETER :: f_full = 0
55  INTEGER, PARAMETER :: f_half = 1
56 
57  INTEGER, PARAMETER :: unit_trace=40, nunit_pres=50, & ! Unit numbers for output
58  nunit_pres1=51
59  INTEGER, PARAMETER :: nsilo=33 ! Number of vectors dumped to silo
60  INTEGER, PARAMETER :: &
61  unit_pres = 1,unit_bsups= 2,unit_bsupu= 3,unit_bsupv= 4, &
62  unit_br = 5,unit_bz = 6,unit_bphi = 7,unit_ksups= 8, &
63  unit_ksupu= 9,unit_ksupv=10,unit_jr =11,unit_jz =12, &
64  unit_jphi =13,unit_fors =14,unit_foru =15,unit_forv =16, &
65  unit_vsups=17,unit_vsupu=18,unit_vsupv=19,unit_vr =20, &
66  unit_vz =21,unit_vphi =22,unit_bsubs=23,unit_bsubu=24, &
67  unit_bsubv=25,unit_ksubs=26,unit_ksubu=27,unit_ksubv=28, &
68  unit_bgradp=29,unit_bdotj=30,unit_bdot_res=31, &
69  unit_divj =32,unit_divb=33
70  CHARACTER(LEN=10), DIMENSION(nsilo), PARAMETER :: &
71  silo_label = &
72  (/ "Pressure ", "Bsups ", "Bsupu ", "Bsupv ", &
73  "Br ", "Bz ", "Bphi ", "Ksups ", &
74  "Ksupu ", "Ksupv ", "Jr ", "Jz ", &
75  "Jphi ", "Force_s ", "Force_u ", "Force_v ", &
76  "Vsups ", "Vsupu ", "Vsupv ", "Vr ", &
77  "Vz ", "Vphi ", "Bsubs ", "Bsubu ", &
78  "Bsubv ", "Ksubs ", "Ksubu ", "Ksubv ", &
79  "Bgradp ", "Bdotj ", "bdot_res ", "DivJ ", &
80  "DivB " /)
81 
82  REAL(dp):: zmaxx, rmaxx, zminx, rminx ! Plotting box boundaries (CYL)
83  INTEGER:: nrs, nzs, nphis ! Resolution (CYL): PHI, R and Z
84  INTEGER:: nvs, nus, nss ! Resolution (VMEC): V, U, S; NVS set to NPHIS if VMEC coordinates required.
85  REAL(dp):: oss, ous, ovs, ors, ozs
86  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:):: pmnch, &
87  bsubsmnsh, bsubumnch, bsubvmnch, & ! Harmonics (w/o jacobian) of P, J_X and B_X
88  bsupsmnsh, bsupumnch, bsupvmnch, jacmnch ! Also jacobian harmonics
89  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:):: pmnsh, &
90  bsubsmnch, bsubumnsh, bsubvmnsh, & ! Harmonics (w/o jacobian) of P, J_X and B_X
91  bsupsmnch, bsupumnsh, bsupvmnsh, jacmnsh ! Also jacobian harmonics
92  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: &
93  ksubsmnsf, ksubumncf, ksubvmncf ! Contains (jacobian*plasma current) harmonics (only dumped if SILO available)
94  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: &
95  ksubsmncf, ksubumnsf, ksubvmnsf ! Contains (jacobian*plasma current) harmonics (only dumped if SILO available)
96  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: &
97  rmnc_ih, zmns_ih ! Coordinate harmonics of R and Z
98  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: &
99  rmns_ih, zmnc_ih
100  LOGICAL :: l_silo_output, l_silo3D
101  CHARACTER(LEN=70):: tagname
102 !-----------------------------------------------
103 ! SILO specific variables
104 !
105  LOGICAL:: l_second = .true. ! Default: dump second cross-section in SILO 2D file
106  LOGICAL:: l_dodisplacement ! Decide whether to dump displacement vector to SILO files or not depending on whether it is non-zero or not
107  INTEGER, PARAMETER:: ilabel = 15
108  CHARACTER(LEN=5), DIMENSION(ilabel):: label_silo
109  CHARACTER(LEN=100) :: filename_cyl, filename_flx
110  INTEGER:: dbfile2d, dbfile3d
111  INTEGER, PARAMETER:: ndims3d = 3
112  INTEGER, DIMENSION(ndims3d):: dims3d
113  INTEGER, PARAMETER:: ndims2d = 2
114  INTEGER, DIMENSION(ndims2d):: dims2d
115  REAL, ALLOCATABLE, DIMENSION(:) :: tmp_silo
116  INTEGER, ALLOCATABLE, DIMENSION(:) :: idum
117  INTEGER:: nxx, nyy, nzz, ilen2D_1, ilen2D_2, ilen2D_3, ilen2D_4, &
118  ilen3d, nsec, index1, index2
119  CHARACTER(LEN=20):: name2d_1, name2d_2, name2d_3, name2d_4, name3d
120  REAL(dp), DIMENSION(:), ALLOCATABLE :: pprof
121  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: pack_bfld
122  REAL, DIMENSION(:,:), ALLOCATABLE :: silo2d, pack_silo
123  REAL, DIMENSION(:,:,:), ALLOCATABLE :: silo3d
124 
125  PUBLIC :: write_output
126 !-----------------------------------------------
127 
128  CONTAINS
129 
130  SUBROUTINE write_output (wout_file, iter, close_wout)
131  USE stel_kinds
132  USE read_wout_mod, ONLY: lwout_opened, read_wout_file, &
133  read_wout_deallocate
134  USE diagnostics_mod, ONLY: dealloc_diagnostics
135  USE siesta_namelist, ONLY: lrestart
136  CHARACTER(LEN=*), INTENT(IN) :: wout_file
137  INTEGER, INTENT(IN) :: iter
138  LOGICAL, INTENT(IN) :: close_wout
139 !-----------------------------------------------
140 ! L o c a l V a r i a b l e s
141 !-----------------------------------------------
142  INTEGER :: istat
143 !-----------------------------------------------
144  IF (lrestart) RETURN
145  CALL read_wout_file(wout_file, istat)
146  CALL assert(istat.eq.0,'Error in write_output')
147 
148  CALL dealloc_diagnostics
149 
150  CALL dump_driver(wout_file, iter)
151 
152  CALL dealloc_diagnostics
153 
154  IF (lwout_opened .and. close_wout) CALL read_wout_deallocate
155 
156  END SUBROUTINE write_output
157 
158  SUBROUTINE dump_driver(wout_file, iter)
159  USE stel_kinds
160  USE read_wout_mod, ns_w => ns, ntor_w => ntor, mpol_w => mpol, &
161  ntmax_w => ntmax, lthreed_w => lthreed, lasym_w => lasym, &
162  nfp_w => nfp
163 !-----------------------------------------------
164  CHARACTER(LEN=*), INTENT(IN) :: wout_file
165  INTEGER, INTENT(IN):: iter
166 !-----------------------------------------------
167 ! L o c a l V a r i a b l e s
168 !-----------------------------------------------
169  CHARACTER(LEN=10) :: cdum1
170  CHARACTER(LEN=120) :: filename
171  INTEGER :: istat, ik, jk, lk, info, nfe, ndum, nmpi, maxprocs, &
172  npmin, npmax, mpi_err
173  INTEGER :: ndum_s, nmpi_s, maxprocs_s, npmin_s, npmax_s
174  REAL(dp) :: zz, rr, zeta0, fmin, ss, uu, vv
175  REAL(dp), DIMENSION(3) :: r_cyl, c_flx, c_flx_old
176  REAL(dp) :: ton, toff
177  REAL(dp) :: dr, dz
178  LOGICAL :: lbfld, lsilo, lendpt
179 !-----------------------------------------------
180 
181 !MODIFY THIS FOR PARALLELIZED VERSION
182  CALL second0(ton)
183 
184  IF (nphis1.LE.0 .OR. nvs1.LE.0) RETURN ! At least one toroidal X-section to produce output
185  nss = nss1; nus = nus1; nvs = nvs1 ! (NSS,NUS,NVS): resolution in flux coordinates
186  nrs = nrs1; nzs = nzs1; nphis = min(10, nphis1) ! (NRS,NZS,NS): resolution in CYL coordinates - only used for pressure PGPLOT output (MAX = 10)
187 
188  l_tracing = l_tracing1 ! =T, writes on file fields for Poincare section calculation
189  l_silo_output = l_silo_output1 ! =T, produce SILO output
190  l_silo3d = l_silo3d1 ! =F, only 2D silo database is dumped out
191  tagname = wout_file ! Used in output files as TAG
192  index1 = index(tagname,"wout")
193  IF (index1 .EQ. 0) index1 = index(tagname,"WOUT")
194  index1 = max(1,index1)
195  tagname = tagname(index1:)
196 
197  CALL loadrzl ! Load rzl_local from VMEC wout file
198  dr = rmax_v-rmin_v; dz = zmax_v-zmin_v
199  rmaxx = rmax_v+0.05*dr; rminx = rmin_v-0.05*dr ! Set plotting box boundaries
200  zmaxx = zmax_v+0.05*dz; zminx = zmin_v-0.05*dz
201  ors = (rmaxx - rminx)/(nrs-1)
202  ozs = (zmaxx - zminx)/(nzs-1)
203  oss = one/(nss-1) ! OSS and OUS are also used later in SILO and Tracing files
204  ous = (2*pi)/(nus-1)
205  ovs = (2*pi)/(nvs-1) ! Define mesh fineness along PHI: not used for pressure plots
206 
207  nsec = nvs/(2*nfp_w)+1 ! Second cross section in 2D SILO file
208  IF (nsec == 1 .OR. nsec > nvs) l_second = .false. ! Then, dump only one cross-section in SILO
209 
210  CALL prepare_quantities
211  CALL open_dump_files(iter) ! Open files needed for dumping info. Also calculates SILO dimensions.
212 #if defined(SILO_AVAIL)
213  IF (l_silo_output) THEN
214  CALL allocate_silo
215  ENDIF
216 #endif
217  IF (iam.EQ.0 .AND. lverbose) THEN
218  WRITE(6,'(/,1x,a)') '>> Starting dump of graphics data...'
219  WRITE(6,'(6x,a)') 'Creating plot box between:'
220  WRITE(6,49) 'RMIN = ', rminx, 'RMAX = ', rmaxx
221  WRITE(6,49) 'ZMIN = ', zminx, 'ZMAX = ', zmaxx
222  WRITE(6,'(6x,a,i4)') 'Toroidal sections computed: ', nphis
223  END IF
224  49 FORMAT(6x,2(a,1p,e10.3,2x))
225 
226  nmpi=0
227  maxprocs = nphis*nzs*nrs
228  ndum = maxprocs/min(nprocs, maxprocs)
229  npmin=iam*ndum+1
230  npmax=(iam+1)*ndum
231  IF (iam .EQ. min(nprocs-1,maxprocs-1)) npmax = maxprocs
232  ALLOCATE (pprof(ndum+(maxprocs-ndum*nprocs)))
233  index1 = 0
234 
235  DO lk = 1, nphis ! Start producing (PGPLOT) PRESSURE file in cylindrical coords
236  zeta0 = (2*(lk-1)*pi)/(nfp_w*nphis) ! Choose all TOROIDAL X-sections in a period [0, 2*pi/NFP], but avoiding last one.
237  DO jk = 1, nzs
238  zz = zminx + (jk-1)*ozs ! Actual Z
239  DO ik = 1, nrs
240  nmpi = nmpi+1
241  IF (nmpi .LT. npmin) cycle
242  IF (nmpi .GT. npmax) GOTO 99
243  rr = rminx + (ik-1)*ors ! Actual R
244  r_cyl(1) = rr; r_cyl(2) = nfp_w*zeta0; r_cyl(3) = zz ! Update target point
245 !SPH: NOTE THIS LOGIC WILL MAKE THE ANSWERS SLIGHTLY DEPENDENT ON NPROCS,
246 ! SINCE C_FLX_OLD WILL BE SLIGHTLY DIFFERENT IF nmpi=npmin DOESN'T
247 ! START AT IK=1!
248  IF (ik > 1 .AND. nmpi.GT.npmin .AND. (c_flx(1) <= one &
249  .AND. c_flx(1) > zero)) THEN ! In this case, use previous solution
250  c_flx = c_flx_old ! as initial guess
251  ELSE
252  c_flx(1:2) = 0 ! Guess for flux coordinate solution
253  ENDIF
254  c_flx(3) = r_cyl(2)
255  CALL cyl2flx(rzl_local, r_cyl, c_flx, ns_w, ntor_w, & ! Obtain flux coordinates
256  mpol_w, ntmax_w, lthreed_w, lasym_w, info, nfe, fmin) ! FMIN not used in what follows.
257  c_flx_old = c_flx
258  IF (info .ne. 0 .AND. info.ne.-3) THEN ! Note: INFO = -3 => ss > 1, OUTSIDE!
259  info = -3
260  ENDIF
261  ss = min(sqrt(c_flx(1)), one) ! SIESTA radial coordinate is s (not flux: s**2)
262  uu = c_flx(2)
263  vv = c_flx(3)/nfp_w
264  index1 = index1+1
265  CALL dump_pressure(vv, ss, uu, info, index1) ! Now, get quantities and dump them on file
266  ENDDO
267  ENDDO
268  ENDDO
269 
270  99 CONTINUE
271 !
272 ! WRITE filename_cyl
273 !
274  CALL write_pressure_data(nunit_pres, ndum, maxprocs)
275 
276  IF (iam .EQ. 0 .and. lverbose) THEN
277  WRITE(*,'(/,6x,a)') 'Completed.'
278  WRITE(*,'(6x,a)') 'Pressure (PGPLOT) using CYLINDRICAL coords..' ! Use uniformly-spaced grid in (R,v,Z)
279  END IF
280 
281  info = 0
282  nmpi = 0
283  maxprocs = nphis*(nss-1)*nus
284  ndum = maxprocs/min(nprocs, maxprocs)
285  npmin=iam*ndum+1
286  npmax=(iam+1)*ndum
287  IF (iam .EQ. min(nprocs-1,maxprocs-1)) npmax = maxprocs
288  ALLOCATE (pprof(ndum+(maxprocs-ndum*nprocs)))
289  index1 = 0
290 
291  DO lk = 1, nphis ! Start producing (PGPLOT) PRESSURE file in SIESTA flux coords
292  vv = (2*(lk-1)*pi)/(nfp_w*nphis) ! Choose all TOROIDAL X-sections in a period [0, 2*pi/NFP], but avoiding last one.
293  DO ik = 2, nss
294  ss = (ik-1)*oss ! SIESTA coords
295  DO jk = 1, nus
296  nmpi = nmpi+1
297  IF (nmpi .LT. npmin) cycle
298  IF (nmpi .GT. npmax) GOTO 199
299  uu = (jk-1)*ous ! Actual U(=Theta) in [0,2*pi]
300 !! CALL flx_to_cyl(ss, uu, vv, rr, zz)
301  index1 = index1+1
302  CALL dump_pressure(vv, ss, uu, info, index1) ! Now, get quantities and store in pprof
303  ENDDO
304  ENDDO
305  ENDDO
306 
307  199 CONTINUE
308 !
309 ! WRITE filename_flx
310 !
311  CALL write_pressure_data(nunit_pres1, ndum, maxprocs)
312 
313  IF (iam .EQ. 0 .and. lverbose) THEN
314  WRITE(*,'(/,6x,a)') 'Completed.'
315  WRITE(*,'(6x,a)') 'Pressure (PGPLOT) using SIESTA coords..' ! Use uniformly-spaced grid in (s,u,v)
316  END IF
317 
318 !
319 ! Calculation in uniformly-spaced grid for Tracing and SILO output if requested
320 !
321  tracing: IF (l_tracing .OR. l_silo_output) THEN ! Do 3D calculation only if L_TRACING = T or SILO_OUTPUT desired
322  IF (iam .EQ. 0 .and. lverbose) THEN
323  WRITE(*,'(/,6x,a)') 'Creating additional files:'
324  IF (l_tracing) THEN
325  WRITE(*,'(6x,a)') 'Tracing BFIELD file..'
326  WRITE(*,'(6x,a,i4)') 'Toroidal sections computed:', nvs
327  ENDIF
328 #if defined(SILO_AVAIL)
329  IF (l_silo_output) THEN
330  WRITE(*,'(6x,a)') 'SILO 2D database..'
331  WRITE(*,'(6x,2a,1p,e10.3,a)') 'Toroidal sections computed: 0.0', &
332  ' and ', real(360*(nsec-1))/(nvs-1), ' degs.'
333  IF (l_silo3d) THEN
334  WRITE(*,'(6x,a)') 'SILO 3D database..'
335  WRITE(*,'(6x,a,i4)') 'Toroidal sections computed:', nvs
336  ENDIF
337  ENDIF
338 #endif
339  END IF
340 
341  npmin = 1; npmax = -1
342  IF (l_tracing) THEN
343  maxprocs = (nvs-1)*(nss-1)*(nus-1)
344  ndum = maxprocs/min(nprocs, maxprocs)
345  npmin=iam*ndum+1
346  npmax=(iam+1)*ndum
347  IF (iam .EQ. min(nprocs-1,maxprocs-1)) npmax = maxprocs
348  nmpi=ndum+(maxprocs-ndum*nprocs)
349  ALLOCATE (pack_bfld(6, nmpi), stat=istat)
350  CALL assert(istat.eq.0,'tracing storage allocation error')
351  END IF
352 
353  npmin_s = 1; npmax_s = -1
354 #if defined(SILO_AVAIL)
355  IF (l_silo_output) THEN
356  maxprocs_s = nvs*(nss-1)*nus
357  ndum_s = maxprocs_s/min(nprocs, maxprocs_s)
358  npmin_s=iam*ndum_s+1
359  npmax_s=(iam+1)*ndum_s
360  IF (iam .EQ. min(nprocs-1,maxprocs_s-1)) npmax_s = maxprocs_s
361  nmpi_s=ndum_s+(maxprocs_s-ndum_s*nprocs)
362  ALLOCATE (pack_silo(nmpi_s,nsilo), stat=istat)
363  CALL assert(istat.eq.0,'silo storage allocation error')
364  END IF
365 #endif
366 
367  index1=0; index2=0
368  nmpi=0; nmpi_s=0
369  lsilo = .false.
370 !
371 ! VECTOR ORDER: ik, lk, jk
372 !
373  DO jk = 1, nus
374  DO lk = 1, nvs
375  vv = (lk-1)*ovs ! Actual V(=PHI) in [0, 2*pi], including last point. Needed for VISIT, but it will not be written on Poincare BFIELD file.
376  lendpt = (jk.NE.nus .AND. lk.NE.nvs) ! Skip angle end points (0,2p]
377  lbfld = .false.
378  DO ik = 2, nss
379  ss = (ik-1)*oss ! Actual SIESTA S-coordinate
380  IF (lendpt .AND. l_tracing) THEN
381  nmpi = nmpi+1
382  lbfld = (nmpi.GE.npmin .AND. nmpi.LE.npmax)
383  END IF
384  IF (l_silo_output) THEN
385  nmpi_s = nmpi_s+1
386  lsilo = (nmpi_s.GE.npmin_s .AND. nmpi_s.LE.npmax_s)
387  END IF
388  IF (nmpi.GT.npmax .AND. nmpi_s.GT.npmax_s) GOTO 299
389  IF (.NOT.lbfld .AND. .NOT.lsilo) cycle
390  uu = (jk-1)*ous ! Actual U(=Theta) in [0,2*pi], including last point. Needed for VISIT but not written on Poincare BFIELD file.
391  CALL flx_to_cyl(ss, uu, vv, rr, zz)
392  CALL dump_quantities(rr, zz, vv, ss, uu, lbfld, lsilo)
393  ENDDO
394  ENDDO
395  ENDDO
396 
397  299 CONTINUE
398 
399  bfield: IF (l_tracing) THEN
400 !
401 ! WRITE bfield-trace file
402 !
403  CALL write_bfield_data(unit_trace, ndum, maxprocs)
404  IF (iam .EQ. 0 .and. lverbose) WRITE(*,'(/,6x,a)') 'Tracing done'
405  END IF bfield
406 
407  END IF tracing
408 
409 #if defined(SILO_AVAIL)
410  IF (l_silo_output) THEN
411  CALL dump_silo (ndum_s, ndum_s+(maxprocs_s-ndum_s*nprocs)) ! combines 2d, 3d silo dump
412  DEALLOCATE (pack_silo)
413  IF (iam .EQ. 0 .and. lverbose) THEN
414  WRITE(*,'(/,6x,a)') 'SILO 2D database done'
415  IF(l_silo3d) WRITE(*,'(6x,a)') 'SILO 3D database done'
416 ! CALL dealloc_SILO
417  END IF
418  ENDIF
419 #endif
420  CALL second0(toff)
421  IF (iam .EQ. 0 .and. lverbose) WRITE(*,'(/,6x,a,f8.3,a)') &
422  'Graphics output completed in ',toff-ton,' s'
423 
424  CALL close_dump_files
425  CALL dealloc_quantities_dump ! Deallocates quantities if needed.
426 
427  END SUBROUTINE dump_driver
428 
429 
430  SUBROUTINE write_pressure_data(nunit, ndum1, maxprocs)
431  INTEGER, INTENT(IN) :: nunit, ndum1, maxprocs
432  INTEGER :: ndum, to=0, from, tag, ik, lk, mpi_err
433 
434  ndum = ndum1
435  tag = 33+nunit
436  DO lk=1, nprocs
437  IF (lk .EQ. nprocs) ndum=ndum1+(maxprocs-ndum1*nprocs)
438 #if defined(MPI_OPT)
439  from = lk-1
440  IF (lk .GT. 1) THEN
441  IF (iam .EQ. 0) THEN
442  CALL mpi_recv(pprof,ndum,mpi_real8,from,tag, &
443  siesta_comm,mpi_status,mpi_err)
444  ELSE IF (iam .EQ. from) THEN
445  CALL mpi_send(pprof,ndum,mpi_real8,to,tag, &
446  siesta_comm,mpi_err)
447  END IF
448  END IF
449 #endif
450  IF (iam .EQ. to) THEN
451  DO ik=1, ndum
452  WRITE(nunit, '(1e14.6)') pprof(ik)
453  END DO
454  END IF
455  END DO
456 
457  DEALLOCATE (pprof)
458 
459  END SUBROUTINE write_pressure_data
460 
461 
462  SUBROUTINE write_bfield_data(nunit, ndum1, maxprocs)
463  INTEGER, INTENT(IN) :: nunit, ndum1, maxprocs
464  INTEGER :: ndum, to=0, from, tag, ik, lk, mpi_err
465 
466  ndum = ndum1
467  tag = 33+nunit
468 
469  DO lk=1, nprocs
470  IF (lk .EQ. nprocs) ndum=ndum+(maxprocs-ndum*nprocs)
471 #if defined(MPI_OPT)
472  from = lk-1
473  IF (lk .GT. 1) THEN
474  IF (iam .EQ. 0) THEN
475  CALL mpi_recv(pack_bfld,6*ndum,mpi_real8,from, &
476  tag,siesta_comm,mpi_status,mpi_err)
477  ELSE IF (iam .EQ. from) THEN
478  CALL mpi_send(pack_bfld,6*ndum,mpi_real8,to, &
479  tag,siesta_comm,mpi_err)
480  END IF
481  END IF
482 #endif
483 ! WRITEOUT file used by POINCARE code.
484  IF (iam .EQ. 0) THEN
485  DO ik=1,ndum
486  WRITE(unit_trace, 122) pack_bfld(:,ik)
487  END DO
488  END IF
489  END DO
490 
491  DEALLOCATE (pack_bfld)
492 122 FORMAT(3(2x,f12.4),3(2x,1pe16.6)) !Must be consistent with POINCARE program
493 
494  END SUBROUTINE write_bfield_data
495 
496 
497  SUBROUTINE open_dump_files(iter)
498  USE stel_kinds
499  INTEGER, INTENT(IN):: iter
500  INTEGER:: ilen, ierr
501  CHARACTER(LEN=30) :: label, tag
502  CHARACTER(LEN=128) :: siloname3d, siloname2d
503  CHARACTER(LEN=100) :: filename
504  CHARACTER(LEN=10) :: cdum1
505  REAL(dp) :: r1, z1, p1, rs, zs
506  CHARACTER(LEN=20) :: cdum
507  REAL, ALLOCATABLE, DIMENSION(:) :: rr, zz, pp
508  REAL, ALLOCATABLE, DIMENSION(:,:,:) :: xc, yc, zc
509  INTEGER :: istat, err, optlistid, i, j, k, ierr2
510 
511  IF (iam .NE. 0) RETURN
512 
513  tag = trim(restart_ext)
514  WRITE(cdum1,'(i4.4)') min(iter,9999)
515  cdum = '-' // adjustl(cdum1)
516 
517 !
518 ! PGPLOT pressure plots
519 !
520  filename_cyl = trim(tag)//'-pressure-CYL'//trim(cdum)
521  filename=trim(filename_cyl)//'.dat' ! By default, always write out pressure files
522  OPEN(unit=nunit_pres, file=filename, status='replace')
523  WRITE (nunit_pres, '(i6)') nphis ! Headers for CYL PGPLOT-plots
524  WRITE (nunit_pres, '(2i6)') nrs, nzs
525  WRITE (nunit_pres, '(4e14.6)') rminx, rmaxx, zminx, zmaxx
526  label = 'pressure'
527  WRITE (nunit_pres, '(a)') trim(label)
528  WRITE (nunit_pres, '(a)') trim(tag)
529 
530  filename_flx = trim(tag)//'-pressure-FLUX'//trim(cdum)
531  filename=trim(filename_flx)//'.dat' ! By default, always write out pressure files
532  OPEN(unit=nunit_pres1, file=filename, status='replace')
533  WRITE (nunit_pres1, '(i6)') nphis ! Headers for SIESTA PGPLOT-plots
534  WRITE (nunit_pres1, '(2i6)') nus, nss-1
535  WRITE (nunit_pres1, '(4e14.6)') zero, (nus-1)*ous/(2.*pi), oss, one
536  label = 'pressure'
537  WRITE (nunit_pres1, '(a)') trim(label)
538  WRITE (nunit_pres1, '(a)') trim(tag)
539 !
540 ! Tracing BFIELD file for POINCARE
541 !
542  IF (l_tracing) THEN ! Then, l_VMEC_uniform = .TRUE.
543  filename = trim(tag)//'-bfield_tracing'// &
544  trim(cdum)//'.dat' ! Since Poincare plots require flux-uniform mesh
545  OPEN(unit=unit_trace, file = filename, status = 'unknown')
546  WRITE(unit_trace, *) nvs-1, nss-1, nus-1 ! NSS-1, since js = 1 is not done because it corresponds to the magnetic axis
547  ENDIF ! NVS-1, NUS-1 since phi=2*pi and theta=2*pi are not included in POINCARE file
548 
549 #if defined(SILO_AVAIL)
550  IF (.NOT.l_silo_output) RETURN ! Get SILO OUTPUT ready.
551  IF (iam .NE. 0) RETURN
552 
553  siloname2d = trim(tag)//'-2D'//trim(cdum)//'.silo'
554  ilen = len(trim(siloname2d)); dbfile2d = 0
555  ierr = dbcreate(trim(siloname2d), ilen, db_clobber, db_local, &
556  "2Ddata", 6, db_pdb, dbfile2d)
557  IF(l_silo3d) THEN
558  siloname3d = trim(tag)//'-3D'//trim(cdum)//'.silo' ! Output mesh is in cartesian form, but it is uniformly spaced in VMEC 3D grid
559  ilen = len(trim(siloname3d)); dbfile3d = 0
560  ierr = dbcreate(trim(siloname3d), ilen, db_clobber, db_local, &
561  "3Ddata", 6, db_pdb, dbfile3d)
562  ENDIF
563 !
564 ! Update SILO information: compute and dump meshes onto file
565 !
566  ALLOCATE(idum(1)); idum = db_f77null ! Use to dump scalars on database as fake MIXVAR
567  nxx = nss-1; nyy = nvs; nzz =nus
568  name3d = "Cart3D"; ilen3d = len(trim(name3d))
569  name2d_1 = "Flux2D_1"; ilen2d_1 = len(trim(name2d_1))
570  name2d_2 = "Flux2D_2"; ilen2d_2 = len(trim(name2d_2))
571  name2d_3 = "Cart2D_1"; ilen2d_3 = len(trim(name2d_3))
572  name2d_4 = "Cart2D_2"; ilen2d_4 = len(trim(name2d_4))
573  dims3d(1) = nxx; dims3d(2) = nyy; dims3d(3) = nzz
574  dims2d(1) = nxx; dims2d(2) = nzz
575 
576  ALLOCATE(rr(nxx), zz(nyy), pp(nzz))
577  DO i = 1, nxx
578  rr(i) = i*oss ! Radial S
579  ENDDO
580  DO i = 1, nyy
581  zz(i) = (i-1)*ovs/(2.d0*pi) ! Toroidal angle/2pi
582  ENDDO
583  DO i = 1, nzz
584  pp(i) = (i-1)*ous/(2.d0*pi) ! Poloidal angle/2pi
585  ENDDO
586 
587  err = dbmkoptlist(2, optlistid) ! 2D MESH (in SIESTA coordinates)
588  err = dbaddcopt(optlistid, dbopt_xlabel, "s", 1)
589  err = dbaddcopt(optlistid, dbopt_ylabel, "u", 1)
590  err = dbputqm(dbfile2d, trim(name2d_1), ilen2d_1, "s", & ! Use different meshes for different
591  1, "u", 1, "v", 1, rr, pp, zz, dims2d, ndims2d, & ! toroidal cross-sections, since, in
592  db_float, db_collinear, optlistid, ierr) ! stellarator fields, materials can change
593  err = dbputqm(dbfile2d, trim(name2d_2), ilen2d_2, "s", &
594  1, "u", 1, "v", 1, rr, pp, zz, dims2d, ndims2d, &
595  db_float, db_collinear, optlistid, ierr)
596  err = dbfreeoptlist(optlistid)
597  zz = 2.d0*pi*zz; pp = 2.d0*pi*pp
598 
599  ALLOCATE(xc(nxx,nzz,2), zc(nxx,nzz,2)) ! 2D MESH (R-Z in a constant PHI X-section)
600  xc = 0.d0; zc = 0.d0
601  DO i = 1, nxx
602  r1 = rr(i) ! Flux surface label
603  DO j = 1, nzz
604  z1 = zz(1) ! Actual toroidal angle
605  p1 = pp(j) ! Poloidal angle
606  CALL flx_to_cyl(r1, p1, z1, rs, zs) ! Get cylindrical coords (rs, zs, z1) from flux coords (r1, p1, z1) for z1 =0.d0
607  xc(i, j, 1) = rs
608  zc(i, j, 1) = zs
609  z1 = zz(nsec)
610  CALL flx_to_cyl(r1, p1, z1, rs, zs) ! Get cylindrical coords (rs, zs, z1) from flux coords (r1, p1, z1)
611  xc(i, j, 2) = rs
612  zc(i, j, 2) = zs
613  ENDDO
614  ENDDO
615  err = dbmkoptlist(2, optlistid)
616  err = dbaddcopt(optlistid, dbopt_xlabel, "R", 1)
617  err = dbaddcopt(optlistid, dbopt_ylabel, "Z", 1)
618  err = dbputqm(dbfile2d, trim(name2d_3), ilen2d_3, "R", &
619  1, "Z", 1, "P", 1, xc(:,:,1), zc(:,:,1), zz, dims2d, & ! zz is a dummy here
620  ndims2d, db_float, db_noncollinear, optlistid, ierr)
621  err = dbputqm(dbfile2d, trim(name2d_4), ilen2d_4, "R", &
622  1, "Z", 1, "P", 1, xc(:,:,2), zc(:,:,2), zz, dims2d, & ! zz is a dummy here
623  ndims2d, db_float, db_noncollinear, optlistid, ierr)
624  err = dbfreeoptlist(optlistid)
625  DEALLOCATE(xc, zc)
626 
627  IF (l_silo3d) THEN ! Prepare 3D MESH
628  ALLOCATE(xc(nxx, nyy, nzz), yc(nxx, nyy, nzz), & ! Convert to cartesian coordinates to accelerate visit and improve accuracy
629  zc(nxx, nyy, nzz))
630  DO i = 1, nxx
631  r1 = rr(i)
632  DO j = 1, nyy
633  z1 = zz(j) ! Actual toroidal angle
634  DO k = 1, nzz
635  p1 = pp(k) ! Actual poloidal angle
636  CALL flx_to_cyl(r1, p1, z1, rs, zs) ! Get cylindrical coords (rs, zs, z1) from flux coords (r1, z1, p1)
637  xc(i, j, k) = rs*cos(z1)
638  yc(i, j, k) = rs*sin(z1)
639  zc(i, j, k) = zs
640  ENDDO
641  ENDDO
642  ENDDO
643  err = dbmkoptlist(3, optlistid)
644  err = dbaddcopt(optlistid, dbopt_xlabel, "X", 1)
645  err = dbaddcopt(optlistid, dbopt_ylabel, "Y", 1)
646  err = dbaddcopt(optlistid, dbopt_zlabel, "Z", 1)
647  err = dbputqm(dbfile3d, trim(name3d), ilen3d, &
648  "X", 1, "Y", 1, "Z", 1, xc, yc, zc, dims3d, &
649  ndims3d, db_float, db_noncollinear, optlistid, ierr)
650  err = dbfreeoptlist(optlistid)
651  DEALLOCATE(xc, yc, zc)
652  ENDIF
653  DEALLOCATE(rr, zz, pp)
654 #endif
655  END SUBROUTINE open_dump_files
656 
657 
658  SUBROUTINE close_dump_files
659  USE stel_kinds
660  INTEGER:: ierr
661 
662  IF (iam .NE. 0) RETURN
663  CLOSE(unit=nunit_pres)
664  CLOSE(unit=nunit_pres1)
665  IF(l_tracing) CLOSE(unit=unit_trace)
666 
667 #if defined(SILO_AVAIL)
668  IF (l_silo_output) THEN ! SILO Close databases
669  ierr = dbclose(dbfile2d)
670  IF (l_silo3d) ierr = dbclose(dbfile3d)
671  ENDIF
672 #endif
673  END SUBROUTINE close_dump_files
674 
675 
676  SUBROUTINE prepare_quantities
677  USE stel_kinds
678  USE diagnostics_mod, ONLY: divj, divb, bdotj, lcurr_init, &
679  bdotjijh
680  USE siesta_namelist, ONLY: lresistive
681  USE nscalingtools, ONLY: startglobrow, endglobrow
682  USE quantities, ONLY: pijh=>pijh0, jacobh, & ! Remember that the JBSUPs and the KSUBs have a jacobian factor that needs to be removed
683  b_factor, p_factor, jvsupsmncf, &
684  bsupsijh0, bsupuijh0, bsupvijh0, &
685  ksupsmnsf, ksupumncf, ksupvmncf, &
686  ksupsmncf, ksupumnsf, ksupvmnsf, &
687  bsubsijf, bsubuijf, bsubvijf, &
688  ksubsijf, ksubuijf, ksubvijf
689  USE siesta_init, ONLY: init_state
690  USE island_params, ONLY: fourier_context
691 
692  IMPLICIT NONE
693 
694 !-----------------------------------------------
695 ! L o c a l V a r i a b l e s
696 !-----------------------------------------------
697  INTEGER :: istat, start, endr, js
698  LOGICAL :: lresist
699  REAL(dp) :: pp
700  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: &
701  bsubsijh , bsubuijh , bsubvijh
702 !-----------------------------------------------
703 !
704 ! Deal with metrics first
705 !
706  js = SIZE(rmnc_i,3)
707  CALL assert(js.EQ.ns,' rmnc_i NOT allocated on extended mesh!')
708  ALLOCATE(rmnc_ih(0:mpol,-ntor:ntor,ns), &
709  zmns_ih(0:mpol,-ntor:ntor,ns), stat=istat)
710  rmnc_ih = zero
711  zmns_ih = zero
712 
713  CALL to_half_mesh(rmnc_i, rmnc_ih) ! We need RMNC and ZMNS on half mesh.
714  CALL to_half_mesh(zmns_i, zmns_ih) ! RMNC_i, ZMNS_i from METRIC. Deallocated elsewhere.
715  ALLOCATE(jacmnch(0:mpol,-ntor:ntor,ns), stat=istat) ! Get harmonics of jacobian on half mesh
716  CALL assert(istat.eq.0,'Allocation failed in PREPARE_QUANTITIES')
717  CALL fourier_context%tomnsp(jacobh, jacmnch, f_cos)
718  jacmnch(:,:,1) = 0
719 
720  IF (lasym) THEN
721  ALLOCATE(rmns_ih(0:mpol,-ntor:ntor,ns), &
722  zmnc_ih(0:mpol,-ntor:ntor,ns), stat=istat)
723  rmns_ih = zero; zmnc_ih = zero
724 
725  CALL to_half_mesh(rmns_i, rmns_ih) ! We need RMNC and ZMNS on half mesh.
726  CALL to_half_mesh(zmnc_i, zmnc_ih) ! RMNS_i, ZMNC_i from METRIC. Deallocated elsewhere.
727  ALLOCATE(jacmnsh(0:mpol,-ntor:ntor,ns), stat=istat) ! Get harmonics of jacobian on half mesh
728  CALL assert(istat.eq.0,'Allocation failed in PREPARE_QUANTITIES')
729  CALL fourier_context%tomnsp(jacobh, jacmnsh, f_sin)
730  jacmnsh(:,:,1)=0
731  END IF
732 
733 !
734 ! COMPUTE QUANTITIES NEEDED FOR OUTPUT FOR ALL RADII HERE
735 !
736  lresist=lresistive
737  lresistive = .true. ! compute ksubX in cv_currents
738  CALL init_state(.false., .false.) ! get quantities on [1,ns]
739  lresistive = lresist
740 
741 ! SCALE FIELDS TO PHYSICAL UNITS
742  pijh = pijh/p_factor ! Remove normalization (restore at end)
743 
744  ALLOCATE(pmnch(0:mpol,-ntor:ntor,ns), stat=istat) ! ... and then recompute harmonics
745  CALL assert(istat.eq.0,'Allocation failed in DUMP_PRESSURE') ! which are stored in PMNCH
746  CALL fourier_context%tomnsp(pijh, pmnch, f_cos); pmnch(:,:,1) = zero
747 
748  WHERE (pijh .LT. zero) pijh = zero
749  CALL assert(SIZE(pmnch,3).EQ.SIZE(pijh,3), &
750  'SIZE pmnch != pijh in prepare_quantities')
751 
752  IF (lasym) THEN
753  ALLOCATE(pmnsh(0:mpol,-ntor:ntor,ns), stat=istat) ! ... and then recompute harmonics
754  CALL assert(istat.eq.0,'Allocation failed in DUMP_PRESSURE') ! which are stored in PMNCH
755  CALL fourier_context%tomnsp(pijh, pmnsh, f_sin); pmnsh(:,:,1) = zero
756  END IF
757 
758  IF (l_tracing .OR. l_silo_output) THEN
759  bsupsijh0 = bsupsijh0/b_factor
760  bsupuijh0 = bsupuijh0/b_factor
761  bsupvijh0 = bsupvijh0/b_factor
762 
763  ALLOCATE(bsupsmnsh(0:mpol,-ntor:ntor,ns), &
764  bsupumnch(0:mpol,-ntor:ntor,ns), &
765  bsupvmnch(0:mpol,-ntor:ntor,ns), stat=istat)
766  CALL assert(istat.eq.0,'Allocation failed in PREPARE_QUANTITIES')
767  CALL fourier_context%tomnsp(bsupsijh0, bsupsmnsh, f_sin); bsupsmnsh(:,:,1) = zero ! compute harmonics
768  CALL fourier_context%tomnsp(bsupuijh0, bsupumnch, f_cos); bsupumnch(:,:,1) = zero
769  CALL fourier_context%tomnsp(bsupvijh0, bsupvmnch, f_cos); bsupvmnch(:,:,1) = zero
770 
771  IF (lasym) THEN
772  ALLOCATE(bsupsmnch(0:mpol,-ntor:ntor,ns), &
773  bsupumnsh(0:mpol,-ntor:ntor,ns), &
774  bsupvmnsh(0:mpol,-ntor:ntor,ns), stat=istat)
775  CALL assert(istat.eq.0,'Allocation failed in PREPARE_QUANTITIES')
776  CALL fourier_context%tomnsp(bsupsijh0, bsupsmnch, f_cos)
777  bsupsmnch(:,:,1) = zero ! compute harmonics
778  CALL fourier_context%tomnsp(bsupuijh0, bsupumnsh, f_sin)
779  bsupumnsh(:,:,1) = zero
780  CALL fourier_context%tomnsp(bsupvijh0, bsupvmnsh, f_sin)
781  bsupvmnsh(:,:,1) = zero
782  END IF
783  END IF
784 
785 
786 #if defined(SILO_AVAIL)
787  IF (l_silo_output) THEN
788  ALLOCATE (bsubsijh(ntheta,nzeta,ns), &
789  bsubuijh(ntheta,nzeta,ns), &
790  bsubvijh(ntheta,nzeta,ns), stat=istat)
791  CALL assert(istat.eq.0,'Allocation failed in DUMPING PREPARE_QUANTITIES')
792  CALL tolowerh(bsupsijh0, bsupuijh0, bsupvijh0, &
793  bsubsijh , bsubuijh , bsubvijh, 1, ns)
794  ksupsmnsf = ksupsmnsf/b_factor; ksubsijf = ksubsijf/b_factor
795  ksupumncf = ksupumncf/b_factor; ksubuijf = ksubuijf/b_factor
796  ksupvmncf = ksupvmncf/b_factor; ksubvijf = ksubvijf/b_factor
797 
798  ALLOCATE(bsubsmnsh(0:mpol,-ntor:ntor,ns), &
799  bsubumnch(0:mpol,-ntor:ntor,ns), &
800  bsubvmnch(0:mpol,-ntor:ntor,ns), &
801  ksubsmnsf(0:mpol,-ntor:ntor,ns), &
802  ksubumncf(0:mpol,-ntor:ntor,ns), &
803  ksubvmncf(0:mpol,-ntor:ntor,ns), stat=istat)
804  CALL assert(istat.eq.0,'Allocation failed in DUMPING PREPARE_QUANTITIES')
805  CALL fourier_context%tomnsp(bsubsijh, bsubsmnsh, f_sin) ! compute harmonics
806  CALL fourier_context%tomnsp(bsubuijh, bsubumnch, f_cos)
807  CALL fourier_context%tomnsp(bsubvijh, bsubvmnch, f_cos)
808  CALL fourier_context%tomnsp(ksubsijf, ksubsmnsf, f_sin) ! compute harmonics
809  CALL fourier_context%tomnsp(ksubuijf, ksubumncf, f_cos)
810  CALL fourier_context%tomnsp(ksubvijf, ksubvmncf, f_cos)
811 
812  IF (lasym) THEN
813  ksupsmncf = ksupsmncf/b_factor
814  ksupumnsf = ksupumnsf/b_factor
815  ksupvmnsf = ksupvmnsf/b_factor
816 
817  ALLOCATE(bsubsmnch(0:mpol,-ntor:ntor,ns), &
818  bsubumnsh(0:mpol,-ntor:ntor,ns), &
819  bsubvmnsh(0:mpol,-ntor:ntor,ns), &
820  ksubsmncf(0:mpol,-ntor:ntor,ns), &
821  ksubumnsf(0:mpol,-ntor:ntor,ns), &
822  ksubvmnsf(0:mpol,-ntor:ntor,ns), stat=istat)
823  CALL assert(istat.eq.0,'Allocation failed in DUMPING PREPARE_QUANTITIES')
824  CALL fourier_context%tomnsp(bsubsijh, bsubsmnch, f_cos) ! compute harmonics
825  CALL fourier_context%tomnsp(bsubuijh, bsubumnsh, f_sin)
826  CALL fourier_context%tomnsp(bsubvijh, bsubvmnsh, f_sin)
827  CALL fourier_context%tomnsp(ksubsijf, ksubsmncf, f_cos)
828  CALL fourier_context%tomnsp(ksubuijf, ksubumnsf, f_sin)
829  CALL fourier_context%tomnsp(ksubvijf, ksubvmnsf, f_sin)
830  END IF
831 
832  DEALLOCATE (bsubsijh, bsubuijh, bsubvijh)
833 
834  IF (maxval(jvsupsmncf) == zero) THEN ! Decide whether to dump displacement vector or not
835  l_dodisplacement = .false.
836  ELSE
837  l_dodisplacement = .true.
838  ENDIF
839 
840 ! COMPUTE DIAGNOSTIC OUTPUT
841  lcurr_init=.true.
842  CALL divb(1, ns)
843  CALL divj(1, ns)
844  CALL bdotj(bdotjijh)
845  DEALLOCATE(bdotjijh)
846  END IF
847 #endif
848  END SUBROUTINE prepare_quantities
849 
850 
851  SUBROUTINE dealloc_quantities_dump
852  USE stel_kinds
853 
854  DEALLOCATE(rmnc_ih, zmns_ih, jacmnch, pmnch)
855  IF (lasym) DEALLOCATE(rmns_ih, zmnc_ih, jacmnsh, pmnsh)
856  IF (l_silo_output .OR. l_tracing) THEN
857  DEALLOCATE(bsupsmnsh, bsupumnch, bsupvmnch)
858  IF (lasym) DEALLOCATE(bsupsmnch, bsupumnsh, bsupvmnsh)
859  END IF
860 #if defined(SILO_AVAIL)
861  IF (l_silo_output) THEN
862  DEALLOCATE(ksubsmnsf, ksubumncf, ksubvmncf, &
863  bsubsmnsh, bsubumnch, bsubvmnch)
864  IF (lasym) THEN
865  DEALLOCATE(ksubsmncf, ksubumnsf, ksubvmnsf, &
866  bsubsmnch, bsubumnsh, bsubvmnsh)
867  END IF
868  END IF
869 #endif
870  END SUBROUTINE dealloc_quantities_dump
871 
872 
873  SUBROUTINE flx_to_cyl(ss, uu, vv, rr, zz)
874  USE stel_kinds
875  USE quantities, ONLY: pijh=>pijh0
876  REAL(dp), INTENT(IN):: vv, ss, uu ! flux-coordinates: (s, u, v) = (SS, UU, VV)
877  REAL(dp), INTENT(OUT):: rr, zz ! CYL-coordinates: (rr, zz, v) = (R, Z, Phi)
878 
879  CALL toijsp_1p(rmnc_ih, ss, uu, vv, rr, f_none, f_cos, f_half) ! SUM MN-series at target point
880  CALL toijsp_1p(zmns_ih, ss, uu, vv, zz, f_none, f_sin, f_half) ! SUM MN-series at target point
881 
882  IF (lasym) THEN
883  CALL toijsp_1p(rmns_ih, ss, uu, vv, rr, f_sum, f_sin, f_half) ! SUM MN-series at target point
884  CALL toijsp_1p(zmnc_ih, ss, uu, vv, zz, f_sum, f_cos, f_half) ! SUM MN-series at target point
885  END IF
886 
887  END SUBROUTINE flx_to_cyl
888 
889 
890  SUBROUTINE dump_pressure(vv, ss, uu, info, iunit) ! Subroutine used to produce pressure files (plotted with PGPLOT)
891  USE stel_kinds
892  USE quantities, ONLY: pijh=>pijh0
893  INTEGER, INTENT(IN) :: info, iunit ! ... and "s" is radius (i.e., SQRT(VMEC-S))
894  REAL(dp), INTENT(IN) :: vv, ss, uu ! flux-coordinates: (s, u, v) = (SS, UU, VV)
895  REAL(dp):: pp
896 
897  IF (info == -3) THEN
898 ! WRITE(iunit) zero ! s>1 => Outside of the plasma. (ONLY for CYL)
899  pprof(iunit) = 0
900 
901  ELSE IF (info == 0) THEN ! s<1 => DUMP pressure info
902  CALL toijsp_1p(pmnch, ss, uu, vv, pp, f_none, f_cos, f_half) ! SUM MN-series at target point
903  IF (lasym) THEN
904  CALL toijsp_1p(pmnsh, ss, uu, vv, pp, f_sum, f_sin, f_half) ! SUM MN-series at target point
905  END IF
906 
907  IF (pp .LT. zero) pp = 0 !interpolation roundoff
908  pprof(iunit) = pp
909 
910  ELSE
911  CALL assert(.false.,'INFO IS NEITHER 0 NOR -3 IN DUMP_PRESSURE!')
912  ENDIF
913 
914  END SUBROUTINE dump_pressure
915 
916 
917  SUBROUTINE dump_quantities(rr, zz, vv, ss, uu, lbfld, lsilo) ! Subroutine used to produce SILO files and Poincare BFIELD file
918 
919  USE stel_kinds ! i = s; j = phi; k = theta
920  USE quantities, ONLY: jvsupsmncf, jvsupumnsf, jvsupvmnsf, &
921  ksupsmnsf, ksupumncf, ksupvmncf, &
922  jvsupsmnsf, jvsupumncf, jvsupvmncf, &
923  ksupsmncf, ksupumnsf, ksupvmnsf
924  USE diagnostics_mod, ONLY: bdotjmnch, divjmnsh, divbmnsf, &
925  bdotjmnsh, divjmnch, divbmncf
926  LOGICAL, INTENT(IN) :: lbfld, lsilo
927  REAL(dp), INTENT(IN):: rr, zz, vv, ss, uu ! flux-coordinates: (s, u, v) = (SS, UU, VV)
928  REAL(dp):: pp, bs, bu, bv, ru, rv, rs, zu, zv, zs, &
929  br, bphi, bz, aux1, aux2, aux3, aux4, bdp, ps, pu, &
930  pv, bdp0, bdj, divj, divb, ks, ku, kv, jr, jphi, jz, &
931  ws, wu, wv, wr, wphi, wz, jb, dpp, fors, foru, forv
932  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE:: bdotjmn_res
933  REAL(dp):: temp, bf0, curr0
934  INTEGER :: ds, du, dv, js
935 
936  CALL toijsp_1p(jacmnch, ss, uu, vv, jb, f_none, f_cos, f_half) ! Get Jacobian and metric info
937  CALL toijsp_1p(rmnc_ih, ss, uu, vv, rs, f_ds, f_cos, f_half) ! Get R_s
938  CALL toijsp_1p(rmnc_ih, ss, uu, vv, ru, f_du, f_cos, f_half) ! Get R_u
939  CALL toijsp_1p(rmnc_ih, ss, uu, vv, rv, f_dv, f_cos, f_half) ! Get R_v
940  CALL toijsp_1p(zmns_ih, ss, uu, vv, zs, f_ds, f_sin, f_half) ! Get Z_s
941  CALL toijsp_1p(zmns_ih, ss, uu, vv, zu, f_du, f_sin, f_half) ! Get Z_u
942  CALL toijsp_1p(zmns_ih, ss, uu, vv, zv, f_dv, f_sin, f_half) ! Get Z_v
943 
944  CALL toijsp_1p(pmnch, ss, uu, vv, pp, f_none, f_cos, f_half) ! Get pressure first: SUM MN-series at target point
945  CALL toijsp_1p(bsupsmnsh, ss, uu, vv, bs, f_none, f_sin, f_half) ! SUM MN-series at target point to get contravariant Bs
946  CALL toijsp_1p(bsupumnch, ss, uu, vv, bu, f_none, f_cos, f_half) ! SUM MN-series at target point to get contravariant Bu
947  CALL toijsp_1p(bsupvmnch, ss, uu, vv, bv, f_none, f_cos, f_half) ! SUM MN-series at target point to get contravariant By
948 
949  IF (lasym) THEN
950  ds = ior(f_sum, f_ds)
951  du = ior(f_sum, f_du)
952  dv = ior(f_sum, f_dv)
953 
954  CALL toijsp_1p(jacmnsh, ss, uu, vv, jb, f_sum, f_sin, f_half) ! Get Jacobian and metric info
955  CALL toijsp_1p(rmns_ih, ss, uu, vv, rs, ds, f_sin, f_half) ! Get R_s
956  CALL toijsp_1p(rmns_ih, ss, uu, vv, ru, du, f_sin, f_half) ! Get R_u
957  CALL toijsp_1p(rmns_ih, ss, uu, vv, rv, dv, f_sin, f_half) ! Get R_v
958  CALL toijsp_1p(zmnc_ih, ss, uu, vv, zs, ds, f_cos, f_half) ! Get Z_s
959  CALL toijsp_1p(zmnc_ih, ss, uu, vv, zu, du, f_cos, f_half) ! Get Z_u
960  CALL toijsp_1p(zmnc_ih, ss, uu, vv, zv, dv, f_cos, f_half) ! Get Z_v
961 
962  CALL toijsp_1p(pmnsh, ss, uu, vv, pp, f_sum, f_sin, f_half) ! Get pressure first: SUM MN-series at target point
963  CALL toijsp_1p(bsupsmnch, ss, uu, vv, bs, f_sum, f_cos, f_half) ! SUM MN-series at target point to get contravariant Bs
964  CALL toijsp_1p(bsupumnsh, ss, uu, vv, bu, f_sum, f_sin, f_half) ! SUM MN-series at target point to get contravariant Bu
965  CALL toijsp_1p(bsupvmnsh, ss, uu, vv, bv, f_sum, f_sin, f_half) ! SUM MN-series at target point to get contravariant By
966 
967  END IF
968 
969  br = bs*rs + bu*ru + bv*rv ! Form cylindrical fields (B)
970  bz = bs*zs + bu*zu + bv*zv
971  bphi = rr*bv
972  IF (l_tracing .AND. lbfld) THEN ! WRITE file used by POINCARE code. Do not write out if PHI=2*pi or THETA=2*pi
973  index1 = index1+1
974  pack_bfld(1,index1) = rr; pack_bfld(2,index1) = zz
975  pack_bfld(3,index1) = vv; pack_bfld(4,index1) = br
976  pack_bfld(5,index1) = bz; pack_bfld(6,index1) = bphi
977  ENDIF
978 
979 #if defined(SILO_AVAIL)
980  IF (.NOT.l_silo_output .OR. .NOT.lsilo) RETURN
981  index2 = index2+1
982  CALL assert(index2.LE.SIZE(pack_silo,1),'INDEX2 > SIZE(pack_silo)')
983  pack_silo(index2,unit_pres) = pp ! pressure pres_silo
984  pack_silo(index2,unit_bsups) = bs ! field values: bsups_silo
985  pack_silo(index2,unit_bsupu) = bu ! bsupu_silo
986  pack_silo(index2,unit_bsupv) = bv ! bsupv_silo
987  pack_silo(index2,unit_br) = br ! br_silo
988  pack_silo(index2,unit_bz) = bz ! bz_silo
989  pack_silo(index2,unit_bphi) = bphi ! bphi_si
990 
991  CALL toijsp_1p(ksupsmnsf, ss, uu, vv, ks, f_none, f_sin, f_full) ! Get contravariant jacob*current also: SUM MN-series at target point
992  CALL toijsp_1p(ksupumncf, ss, uu, vv, ku, f_none, f_cos, f_full) ! SUM MN-series at target point
993  CALL toijsp_1p(ksupvmncf, ss, uu, vv, kv, f_none, f_cos, f_full) ! SUM MN-series at target point
994 
995  IF (lasym) THEN
996  CALL toijsp_1p(ksupsmncf, ss, uu, vv, ks, f_sum, f_cos, f_full) ! Get contravariant jacob*current also: SUM MN-series at target point
997  CALL toijsp_1p(ksupumnsf, ss, uu, vv, ku, f_sum, f_sin, f_full) ! SUM MN-series at target point
998  CALL toijsp_1p(ksupvmnsf, ss, uu, vv, kv, f_sum, f_sin, f_full) ! SUM MN-series at target point
999  END IF
1000 
1001  pack_silo(index2,unit_ksups) = ks ! ksups_silo
1002  pack_silo(index2,unit_ksupu) = ku ! ksupu_silo
1003  pack_silo(index2,unit_ksupv) = kv ! ksupv_silo
1004  jr = (ks*rs + ku*ru + kv*rv)/jb ! Form cylindrical fields (Plasma Current)
1005  jz = (ks*zs + ku*zu + kv*zv)/jb
1006  jphi = (rr*kv)/jb
1007  pack_silo(index2,unit_jr) = jr ! jr_silo
1008  pack_silo(index2,unit_jz) = jz ! jz_silo
1009  pack_silo(index2,unit_jphi) = jphi ! jphi_silo
1010 
1011 ! Forces (Need contravariant K's and B's: thus, need to be done here!)
1012  CALL toijsp_1p(pmnch, ss, uu, vv, dpp, f_ds, f_cos, f_half) ! Get radial derivative of pressure
1013  IF (lasym) CALL toijsp_1p(pmnsh, ss, uu, vv, dpp, ds, f_sin, f_half)
1014  fors = (ku*bv - kv*bu - dpp) ! forces_silo
1015 
1016  CALL toijsp_1p(pmnch, ss, uu, vv, dpp, f_du, f_cos, f_half) ! Get poloidal derivative of pressure
1017  IF (lasym) CALL toijsp_1p(pmnsh, ss, uu, vv, dpp, du, f_sin, f_half)
1018  foru = (kv*bs - ks*bv - dpp) ! forceu_silo
1019 
1020  CALL toijsp_1p(pmnch, ss, uu, vv, dpp, f_dv, f_cos, f_half) ! Get toroidal derivative of pressure
1021  IF (lasym) CALL toijsp_1p(pmnsh, ss, uu, vv, dpp, dv, f_sin, f_half)
1022  forv = (ks*bu - ku*bs - dpp) ! forcev_silo
1023 
1024  IF (l_dodisplacement) THEN
1025  CALL toijsp_1p(jvsupsmncf, ss, uu, vv, ws, f_none, f_cos, f_full) ! Get last displacement vector: SUM MN-series at target point
1026  CALL toijsp_1p(jvsupumnsf, ss, uu, vv, wu, f_none, f_sin, f_full) ! SUM MN-series at target point
1027  CALL toijsp_1p(jvsupvmnsf, ss, uu, vv, wv, f_none, f_sin, f_full) ! SUM MN-series at target point
1028 
1029  IF (lasym) THEN
1030  CALL toijsp_1p(jvsupsmncf, ss, uu, vv, ws, f_sum, f_sin, f_full) ! Get last displacement vector: SUM MN-series at target point
1031  CALL toijsp_1p(jvsupumnsf, ss, uu, vv, wu, f_sum, f_cos, f_full) ! SUM MN-series at target point
1032  CALL toijsp_1p(jvsupvmnsf, ss, uu, vv, wv, f_sum, f_cos, f_full) ! SUM MN-series at target point
1033  END IF
1034 
1035  pack_silo(index2,unit_vsups) = ws ! wsups_silo
1036  pack_silo(index2,unit_vsupu) = wu ! wsupu_silo
1037  pack_silo(index2,unit_vsupv) = wv ! wsupv_silo
1038  wr = ws*rs + wu*ru + wv*rv ! Cylindrical fields (Displacement)
1039  wz = ws*zs + wu*zu + wv*zv
1040  wphi = rr*wv
1041  pack_silo(index2,unit_vr) = wr ! wr_silo
1042  pack_silo(index2,unit_vz) = wz ! wz_silo
1043  pack_silo(index2,unit_vphi) = wphi ! wphi_silo
1044  ENDIF
1045 
1046 ! Covariant vectors
1047 
1048  temp = bs
1049  CALL toijsp_1p(bsubsmnsh, ss, uu, vv, bs, f_none, f_sin, f_half) ! Get covariant magnetic field also: SUM MN-series at target point
1050  IF (lasym) CALL toijsp_1p(bsubsmnch, ss, uu, vv, bs, f_sum, f_cos, f_half)
1051  bf0 = bs*temp; temp = bu
1052  CALL toijsp_1p(bsubumnch, ss, uu, vv, bu, f_none, f_cos, f_half) ! SUM MN-series at target point
1053  IF (lasym) CALL toijsp_1p(bsubumnsh, ss, uu, vv, bu, f_sum, f_sin, f_half)
1054  bf0 = bf0 + bu*temp; temp = bv
1055  CALL toijsp_1p(bsubvmnch, ss, uu, vv, bv, f_none, f_cos, f_half) ! SUM MN-series at target point
1056  IF (lasym) CALL toijsp_1p(bsubvmnsh, ss, uu, vv, bv, f_sum, f_sin, f_half)
1057  bf0 = bf0 + bv*temp; bf0 = sqrt(bf0) ! Bfield at point
1058  pack_silo(index2,unit_bsubs) = bs ! bsubs_silo
1059  pack_silo(index2,unit_bsubu) = bu ! bsubu_silo
1060  pack_silo(index2,unit_bsubv) = bv ! bsubv_silo
1061 
1062  temp = ks
1063  CALL toijsp_1p(ksubsmnsf, ss, uu, vv, ks, f_none, f_sin, f_full) ! Get covariant jacob*current also: SUM MN-series at target point
1064  IF (lasym) CALL toijsp_1p(ksubsmncf, ss, uu, vv, ks, f_sum, f_cos, f_full)
1065  curr0 = ks*temp; temp = ku
1066  CALL toijsp_1p(ksubumncf, ss, uu, vv, ku, f_none, f_cos, f_full) ! SUM MN-series at target point
1067  IF (lasym) CALL toijsp_1p(ksubumnsf, ss, uu, vv, ku, f_sum, f_sin, f_full)
1068  curr0 = curr0 + ku*temp; temp = kv
1069  CALL toijsp_1p(ksubvmncf, ss, uu, vv, kv, f_none, f_cos, f_full) ! SUM MN-series at target point
1070  IF (lasym) CALL toijsp_1p(ksubvmnsf, ss, uu, vv, kv, f_sum, f_sin, f_full)
1071  curr0 = curr0 + kv*temp; curr0 = sqrt(abs(curr0/jb)) ! Current J0 at point
1072  pack_silo(index2,unit_ksubs) = ks ! ksubs_silo
1073  pack_silo(index2,unit_ksubu) = ku ! ksubu_silo
1074  pack_silo(index2,unit_ksubv) = kv ! ksubv_silo
1075 !
1076 ! DIAGNOSTICS
1077 !
1078  CALL toijsp_1p(pmnch, ss, uu, vv, ps, f_ds, f_cos, f_half) ! BDOT GRAD P: Need that Bs, Bu, Bv contain covariant components of B!!
1079  CALL toijsp_1p(pmnch, ss, uu, vv, pu, f_du, f_cos, f_half)
1080  CALL toijsp_1p(pmnch, ss, uu, vv, pv, f_dv, f_cos, f_half)
1081 
1082  IF (lasym) THEN
1083  CALL toijsp_1p(pmnsh, ss, uu, vv, ps, ds, f_sin, f_half) ! BDOT GRAD P: Need that Bs, Bu, Bv contain covariant components of B!!
1084  CALL toijsp_1p(pmnsh, ss, uu, vv, pu, du, f_sin, f_half)
1085  CALL toijsp_1p(pmnsh, ss, uu, vv, pv, dv, f_sin, f_half)
1086  END IF
1087 
1088  bdp = bs*ps + bu*pu + bv*pv
1089  bdp0 = bf0*pmnch(2,0,0)/rmnc_ih(2,0,0) ! FIXME: Not sure how to handle asym terms here.
1090  bdp = bdp/bdp0 ! Normalize B-DOT-P to B0, p0 and R0.
1091  pack_silo(index2,unit_bgradp) = bdp ! bgradp_silo
1092 
1093  CALL toijsp_1p(bdotjmnch, ss, uu, vv, bdj, f_none, f_cos, f_half) ! BDOTJ
1094  IF (lasym) CALL toijsp_1p(bdotjmnsh, ss, uu, vv, bdj, f_sum, f_sin, f_half) ! BDOTJ
1095  pack_silo(index2,unit_bdotj) = bdj/bdp0 ! bdotj_silo
1096 
1097  ALLOCATE(bdotjmn_res(0:mpol,-ntor:ntor,ns))
1098  bdotjmn_res = zero
1099  bdotjmn_res(:,mres1,nres1) = bdotjmnch(:,mres1,nres1)
1100  bdotjmn_res(:,mres2,nres2) = bdotjmnch(:,mres2,nres2)
1101  CALL toijsp_1p(bdotjmn_res, ss, uu, vv, bdj, f_none, f_cos, f_half) ! BDOTJ_RES
1102  IF (lasym) THEN ! FIXME: Not sure about this.
1103  bdotjmn_res = zero
1104  bdotjmn_res(:,mres1,nres1) = bdotjmnsh(:,mres1,nres1)
1105  bdotjmn_res(:,mres2,nres2) = bdotjmnsh(:,mres2,nres2)
1106  CALL toijsp_1p(bdotjmn_res, ss, uu, vv, bdj, f_sum, f_sin, f_half) ! BDOTJ_RES
1107  END IF
1108  pack_silo(index2,unit_bdot_res) = bdj/bdp0 ! bdot_res_silo
1109  DEALLOCATE(bdotjmn_res)
1110 
1111  CALL toijsp_1p(divjmnsh, ss, uu, vv, divj, f_none, f_sin, f_half) ! DIVJ
1112  IF (lasym) CALL toijsp_1p(divjmnch, ss, uu, vv, divj, f_sum, f_cos, f_half) ! DIVJ
1113  pack_silo(index2,unit_divj) = divj ! divj_silo
1114 
1115  CALL toijsp_1p(divbmnsf, ss, uu, vv, divb, f_none, f_sin, f_full) ! DIVB
1116  IF (lasym) CALL toijsp_1p(divbmncf, ss, uu, vv, divb, f_sum, f_cos, f_full) ! DIVB
1117  pack_silo(index2,unit_divb) = divb ! divb_silo
1118 
1119  IF (curr0 .EQ. zero) curr0 = 1
1120  IF (ss .GT. (one-one/ohs)) THEN
1121  fors = 0; foru = 0; forv = 0
1122  END IF
1123  pack_silo(index2,unit_fors) = fors/(curr0*bf0) ! Normalize ALL forces to J*B at EACH point
1124  pack_silo(index2,unit_foru) = foru/(curr0*bf0)
1125  pack_silo(index2,unit_forv) = forv/(curr0*bf0)
1126 #endif
1127  END SUBROUTINE dump_quantities
1128 
1129 
1130  SUBROUTINE toijsp_1p(XMN, SS, THETA, ZETA, XUV, DFLAG, &
1131  IPAR, IHALF)
1132 !
1133 ! DESCRIPTION: calculates X on ONE POINT in space (SS,THETA,ZETA) from its Fourier harmonics X_MN.
1134 ! IPARITY = 0, means the quantity X (NOT any of its derivatives!) is COSINE (even); = 1, X is SINE (odd)
1135 ! IHALF= 0, means quantity is on the radial full mesh; =1, on the half radial mesh
1136 ! DFLAG= 0b001 NO/YES radial deriv. DFLAG= 0b010 => NO/YES THETA-Deriv. DFLAG= 0b100 => NO/YES ZETA-Deriv.
1137 !
1138  USE stel_kinds
1139  USE fourier, ONLY: b_ds, b_du, b_dv, b_sum
1140  USE island_params, ONLY: fourier_context
1141 
1142 !-----------------------------------------------
1143 ! D u m m y A r g u m e n t s
1144 !-----------------------------------------------
1145  INTEGER, INTENT(IN):: ihalf ! = 0, FULL (radial) MESH); = 1, HALF MESH
1146  REAL(dp), DIMENSION(0:mpol,-ntor:ntor,ns), & !New Style: ns LAST ARG on input
1147  INTENT(IN):: xmn
1148  INTEGER, INTENT(IN):: dflag
1149  REAL(dp), INTENT(IN):: ss, theta, zeta
1150  REAL(dp), INTENT(INOUT):: xuv
1151  INTEGER, INTENT(IN):: ipar ! = 0, cosine (EVEN); = 1, sine (ODD)
1152 !-----------------------------------------------
1153 ! L o c a l V a r i a b l e s
1154 !-----------------------------------------------
1155  INTEGER :: m, n
1156  INTEGER :: ns_in, ns_out, iout, ndum, istat
1157  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:):: work0
1158  REAL(dp), ALLOCATABLE, DIMENSION(:,:):: work1, work2
1159  REAL(dp), ALLOCATABLE, DIMENSION(:):: xuv_i
1160  REAL(dp) :: arg, d1, d2, ss_in, aux1, aux2
1161  REAL(dp):: lcosmu, lsinmu, lcosnv, lsinnv, lcosmum, &
1162  lsinmum, lcosnvn, lsinnvn
1163 !-----------------------------------------------
1164 
1165  iout = 0 ! Assume first that target point not near edge/axis
1166  ss_in = ss*ohs + ihalf*0.5_dp + one
1167  ns_in = int(ss_in) ! Calculate surface closest target from below
1168  IF (ihalf == 1 .AND. ns_in == 1) ns_in = 2 ! Extrapolation needed on the half mesh near axis
1169  IF (ns_in == ns) ns_in = ns - 1 ! Extrapolation needed (almost always on half mesh) near edge
1170  ns_out = ns_in + 1 ! Surface closest to target from above
1171  d1 = (ss_in - ns_in) ! Coefficients for linear interpolation
1172  d2 = one - d1 ! xv(ss_in) = d1*xv(ns_out) + d2*xv(ns_in)
1173 
1174  IF (btest(dflag, b_ds)) THEN
1175  IF (ihalf == 0) THEN
1176  IF (ns_in == 1) iout = -1 ! Marks that point at axis for radial derivative
1177  ELSE
1178  IF (ns_in == 2) iout = -1 ! Marks that point at axis for radial derivative
1179  ENDIF
1180  IF (ns_in == ns - 1) iout = 1 ! Marks that point at edge for radial derivative
1181  ndum = 4
1182  ELSE
1183  ndum = 2
1184  ENDIF
1185 
1186  ALLOCATE(work0(0:mpol,-ntor:ntor,ndum), stat=istat)
1187  CALL assert(istat.eq.0,'Alloc. problem #1 in TOIJSP_1P')
1188  work0 = zero
1189 
1190  work0(:,:,1) = fourier_context%orthonorm(0:mpol,-ntor:ntor)*xmn(:,:,ns_in)
1191  work0(:,:,2) = fourier_context%orthonorm(0:mpol,-ntor:ntor)*xmn(:,:,ns_out)
1192  IF (btest(dflag, b_ds)) THEN
1193  IF (ihalf == 0) THEN
1194  IF (ns_in .GT. 1) work0(:,:,ndum-1) = &
1195  fourier_context%orthonorm(0:mpol,-ntor:ntor)*xmn(:,:,ns_in-1) ! ...compute quantity at previous...
1196  ELSE
1197  IF (ns_in .GT. 2) work0(:,:,ndum-1) = &
1198  fourier_context%orthonorm(0:mpol,-ntor:ntor)*xmn(:,:,ns_in-1)
1199  ENDIF
1200  IF (ns_out .LT. ns) work0(:,:,ndum) = &
1201  fourier_context%orthonorm(0:mpol,-ntor:ntor)*xmn(:,:,ns_out+1) ! ...and next surfaces as well
1202  ENDIF
1203 
1204  ALLOCATE(work1(-ntor:ntor,ndum), &
1205  work2(-ntor:ntor,ndum), stat=istat)
1206  CALL assert(istat.eq.0,'Alloc. problem #2 in TOIJSP_1P')
1207  work1 = 0; work2 = 0
1208 
1209  DO m = 0, mpol ! First DO sum over poloidal modes
1210  arg = m*theta
1211  lcosmu = cos(arg); lsinmu = sin(arg)
1212  lcosmum = m*lcosmu; lsinmum = m*lsinmu
1213  IF (btest(dflag, b_du)) THEN ! First poloidal-derivative requested
1214  work1(:,:) = work1(:,:) - work0(m,:,:)*lsinmum
1215  work2(:,:) = work2(:,:) + work0(m,:,:)*lcosmum
1216  ELSE
1217  work1(:,:) = work1(:,:) + work0(m,:,:)*lcosmu
1218  work2(:,:) = work2(:,:) + work0(m,:,:)*lsinmu
1219  ENDIF
1220  ENDDO
1221  DEALLOCATE(work0)
1222 
1223  ALLOCATE(xuv_i(ndum), stat=istat)
1224  CALL assert(istat.eq.0,'Alloc. problem #3 in TOIJSP_1P')
1225  xuv_i = zero ! Make sure that SPATIAL variable is zeroed
1226 
1227  IF (btest(dflag, b_dv)) THEN
1228  xuv_i = zero
1229  ELSE IF (ipar == f_cos) THEN
1230  xuv_i = work1(0,:)
1231  ELSE
1232  xuv_i = work2(0,:)
1233  END IF
1234 
1235  DO n = 1, ntor ! Then sum over N>0 and N<0 toroidal modes
1236  arg = n*nfp_i*zeta
1237  lcosnv = cos(arg); lsinnv = sin(arg)
1238  lcosnvn = n*nfp_i*lcosnv; lsinnvn = n*nfp_i*lsinnv
1239 
1240  IF (ipar == f_cos) THEN ! COSINE series
1241 
1242  IF (btest(dflag, b_dv)) THEN
1243  xuv_i = xuv_i &
1244  - (work1(n,:) + work1(-n,:))*lsinnvn &
1245  - (work2(n,:) - work2(-n,:))*lcosnvn
1246  ELSE
1247  xuv_i = xuv_i &
1248  + (work1(n,:) + work1(-n,:))*lcosnv &
1249  - (work2(n,:) - work2(-n,:))*lsinnv
1250  ENDIF
1251 
1252  ELSE ! SINE series
1253 
1254  IF (btest(dflag, b_dv)) THEN
1255  xuv_i = xuv_i &
1256  + (work1(n,:) - work1(-n,:))*lcosnvn &
1257  - (work2(n,:) + work2(-n,:))*lsinnvn
1258  ELSE
1259  xuv_i = xuv_i &
1260  + (work1(n,:) - work1(-n,:))*lsinnv &
1261  + (work2(n,:) + work2(-n,:))*lcosnv
1262  ENDIF
1263 
1264  ENDIF
1265  ENDDO
1266 
1267  DEALLOCATE(work1, work2)
1268 
1269  IF(btest(dflag, b_ds)) THEN
1270  IF (iout == 0) THEN
1271  aux1 = 0.5d0*(xuv_i(2) - xuv_i(ndum-1))
1272  aux2 = 0.5d0*(xuv_i(ndum) - xuv_i(1))
1273  ELSE IF (iout == -1) THEN ! near the axis
1274  aux1 = xuv_i(2) - xuv_i(1)
1275  aux2 = 0.5d0*(xuv_i(ndum) - xuv_i(1))
1276  ELSE ! near the edge
1277  aux1 = 0.5d0*(xuv_i(2) - xuv_i(ndum-1))
1278  aux2 = xuv_i(2) - xuv_i(1)
1279  ENDIF
1280  IF (btest(dflag, b_sum)) THEN
1281  xuv = xuv + (d1*aux2 + d2*aux1)*ohs
1282  ELSE
1283  xuv = (d1*aux2 + d2*aux1)*ohs
1284  END IF
1285  ELSE
1286  IF (btest(dflag, b_sum)) THEN
1287  xuv = xuv + d1*xuv_i(2) + d2*xuv_i(1) ! Interpolate/Extrapolate to actual flux(s) value
1288  ELSE
1289  xuv = d1*xuv_i(2) + d2*xuv_i(1) ! Interpolate/Extrapolate to actual flux(s) value
1290  END IF
1291  ENDIF
1292 
1293  DEALLOCATE(xuv_i)
1294 
1295  END SUBROUTINE toijsp_1p
1296 
1297 #if defined(SILO_AVAIL)
1298  SUBROUTINE allocate_silo
1299  USE stel_kinds
1300  INTEGER :: istat = 0 ! NXX, NYY, NZZ defined in OPEN_FILES
1301 
1302  IF (iam .NE. 0) GOTO 100
1303 
1304  ALLOCATE(silo2d(nxx, nzz), stat = istat) ! Use for output
1305  IF (istat .NE. 0) GOTO 100
1306  silo2d = zero
1307 
1308  ALLOCATE(silo3d(nxx, nyy, nzz), stat = istat) ! Use for output
1309  IF (istat .NE. 0) GOTO 100
1310  silo3d = zero
1311 
1312  ALLOCATE(tmp_silo(nxx*nyy*nzz), stat = istat) ! pressure
1313 
1314  100 CONTINUE
1315 
1316  CALL assert(istat.eq.0,'ERROR in silo alloc.!')
1317 
1318  END SUBROUTINE allocate_silo
1319 
1320 
1321  SUBROUTINE dealloc_silo
1322 
1323  IF (iam .EQ. 0) DEALLOCATE(silo2d, silo3d, tmp_silo, idum)
1324 
1325  END SUBROUTINE dealloc_silo
1326 
1327 
1328  SUBROUTINE dump_silo (nchunk, nlast)
1329  USE stel_kinds
1330  INTEGER, INTENT(IN) :: nchunk, nlast
1331  INTEGER :: i, loop
1332  CHARACTER(LEN=10):: label
1333  REAL, DIMENSION(:), ALLOCATABLE :: br_silo, bphi_silo, &
1334  jr_silo, jphi_silo, jz_silo, &
1335  wr_silo, wphi_silo
1336 
1337  IF (iam.EQ.0 .AND. l_silo3d) THEN
1338  ALLOCATE (br_silo(nxx*nyy*nzz), bphi_silo(nxx*nyy*nzz), &
1339  jr_silo(nxx*nyy*nzz), jphi_silo(nxx*nyy*nzz), &
1340  wr_silo(nxx*nyy*nzz), wphi_silo(nxx*nyy*nzz), &
1341  jz_silo(nxx*nyy*nzz), stat=loop)
1342  CALL assert(loop.eq.0,'Allocation error in dump_SILO')
1343  END IF
1344 
1345  DO loop = 1, nsilo
1346 
1347  IF (.NOT.l_dodisplacement .AND. loop.GE.unit_vsups .AND. &
1348  loop.LE.unit_vphi) cycle
1349 
1350  CALL get_tmp_silo (loop, nchunk, nlast)
1351 
1352  IF (iam .EQ. 0) THEN
1353  label = trim(silo_label(loop))
1354  CALL write_to_silodb_2d(label, 1, tmp_silo, 1)
1355  IF (l_second) CALL write_to_silodb_2d(label, nsec, tmp_silo, 2)
1356  IF (l_silo3d) THEN
1357  CALL write_to_silodb_3d(label, tmp_silo)
1358  IF (loop .EQ. unit_br) br_silo = tmp_silo
1359  IF (loop .EQ. unit_bphi) bphi_silo = tmp_silo
1360  IF (loop .EQ. unit_jr) jr_silo = tmp_silo
1361  IF (loop .EQ. unit_jz) jz_silo = tmp_silo
1362  IF (loop .EQ. unit_jphi) jphi_silo = tmp_silo
1363  IF (loop .EQ. unit_vr) wr_silo = tmp_silo
1364  IF (loop .EQ. unit_vphi) wphi_silo = tmp_silo
1365  END IF
1366  END IF
1367 
1368  END DO
1369 
1370 ! COMPUTE CARTESIAN COMPONENT
1371  IF (iam.EQ.0 .AND. l_silo3d) THEN
1372  CALL dump_silo3d(br_silo, bphi_silo, jr_silo, jphi_silo, &
1373  wr_silo, wphi_silo, jz_silo)
1374  DEALLOCATE (br_silo, bphi_silo, jr_silo, jphi_silo, &
1375  wr_silo, wphi_silo, jz_silo)
1376  END IF
1377 
1378  END SUBROUTINE dump_silo
1379 
1380 
1381  SUBROUTINE get_tmp_silo(index0, nchunk_in, nlast_in)
1382 !
1383 ! STORES pack_silo(:,index) chunks into tmp_silo on processor 0
1384 !
1385  INTEGER, INTENT(IN) :: index0, nchunk_in, nlast_in
1386  INTEGER :: nchunk, ioff, lk, mpi_err, tag, to, from, istat
1387 
1388  nchunk = nchunk_in
1389  IF (iam .EQ. 0) tmp_silo(1:nchunk) = pack_silo(1:nchunk,index0)
1390  ioff=1+nchunk
1391 
1392  to = 0
1393  tag = 333+index0
1394 
1395  DO lk=2, nprocs
1396 #if defined(MPI_OPT)
1397  IF (lk .EQ. nprocs) nchunk=nlast_in
1398  from=lk-1
1399  IF (iam .EQ. 0) THEN
1400  CALL mpi_recv(tmp_silo(ioff),nchunk,mpi_real,from, &
1401  tag,siesta_comm,mpi_status,mpi_err)
1402  ioff = ioff+nchunk
1403  ELSE IF (iam .EQ. from) THEN
1404  CALL mpi_send(pack_silo(1,index0),nchunk,mpi_real,to, &
1405  tag,siesta_comm,mpi_err)
1406  END IF
1407 #endif
1408  END DO
1409 
1410  END SUBROUTINE get_tmp_silo
1411 
1412 
1413  SUBROUTINE dump_silo3d(br_silo, bphi_silo, jr_silo, jphi_silo, &
1414  wr_silo, wphi_silo, jz_silo)
1415  USE stel_kinds
1416  REAL, DIMENSION(nxx,nyy,nzz) :: br_silo, bphi_silo, &
1417  jr_silo, jphi_silo, &
1418  wr_silo, wphi_silo, jz_silo
1419  INTEGER :: i
1420  REAL :: pp, max_curr
1421  CHARACTER(LEN=20):: label
1422 
1423 !
1424 ! Cartesian (X,Y) components of the magnetic field
1425 !
1426  DO i = 1, nvs
1427  pp = 2*pi*(i-1)/(nvs-1) ! Cylindrical angle == Phi
1428  silo3d(:,i,:) = br_silo(:,i,:)*cos(pp) - & ! BX
1429  & bphi_silo(:,i,:)*sin(pp)
1430  ENDDO
1431  label = "Bx"
1432  CALL write_to_silodb_3d(label, silo3d)
1433  DO i = 1, nvs
1434  pp = 2*pi*(i-1)/(nvs-1) ! Cylindrical angle == Phi
1435  silo3d(:,i,:) = br_silo(:,i,:)*sin(pp) + & ! BY
1436  & bphi_silo(:,i,:)*cos(pp)
1437  ENDDO
1438  label = "By"
1439  CALL write_to_silodb_3d(label, silo3d)
1440 
1441  max_curr = maxval(jphi_silo*jphi_silo + jz_silo*jz_silo &
1442  + jr_silo*jr_silo)
1443  max_curr = sqrt(max_curr) ! Used to normalized current if too large
1444 !
1445 ! Cartesian (X,Y) components of the current
1446 !
1447  label = "Jx"
1448  DO i = 1, nvs
1449  pp = 2*pi*(i-1)/(nvs-1) ! Cylindrical angle == Phi
1450  silo3d(:,i,:) = jr_silo(:,i,:)*cos(pp) - & ! JX
1451  jphi_silo(:,i,:)*sin(pp)
1452  ENDDO
1453  IF (max_curr > 1.d2 .AND. max_curr <= 1d5) THEN ! Then, write it out in kA
1454  label = "Jx_kA"
1455  silo3d = silo3d/1.d3
1456  ELSE IF (max_curr > 1.d5) THEN !.. or in MA
1457  label = "Jx_MA"
1458  silo3d = silo3d/1.d6
1459  ENDIF
1460  CALL write_to_silodb_3d(label, silo3d)
1461 
1462  label = "Jy"
1463  DO i = 1, nvs
1464  pp = 2*pi*(i-1)/(nvs-1) ! Cylindrical angle == Phi
1465  silo3d(:,i,:) = jr_silo(:,i,:)*sin(pp) + & ! JY
1466  jphi_silo(:,i,:)*cos(pp)
1467  ENDDO
1468  IF (max_curr > 1.d2 .AND. max_curr <= 1d5) THEN ! Then, write it out in kA
1469  label = "Jy_kA"
1470  silo3d = silo3d/1.d3
1471  ELSE IF (max_curr > 1.d5) THEN ! .. or in MA
1472  label = "Jy_MA"
1473  silo3d = silo3d/1.d6
1474  ENDIF
1475  CALL write_to_silodb_3d(label, silo3d)
1476 
1477  IF (max_curr > 1.d2 .AND. max_curr <= 1d5) THEN ! For completeness, need to rewrite out the cartesian Jz in the correct units.
1478  label = "Jz_kA"
1479  jz_silo = jz_silo/1.d3
1480  CALL write_to_silodb_3d(label, jz_silo)
1481  ELSE IF (max_curr > 1.d5) THEN
1482  label = "Jz_MA"
1483  jz_silo = jz_silo/1.d6
1484  CALL write_to_silodb_3d(label, jz_silo)
1485  ENDIF
1486 
1487  IF (l_dodisplacement) THEN
1488 !
1489 ! Cartesian (X,Y) components of the last displacement vector
1490 !
1491  DO i = 1, nvs
1492  pp = 2*pi*(i-1)/(nvs-1) ! Cylindrical angle == Phi
1493  silo3d(:,i,:) = wr_silo(:,i,:)*cos(pp) - & ! VX
1494  wphi_silo(:,i,:)*sin(pp)
1495  ENDDO
1496  label = "Vx"
1497  CALL write_to_silodb_3d(label, silo3d)
1498  DO i = 1, nvs
1499  pp = 2*pi*(i-1)/(nvs-1) ! Cylindrical angle == Phi
1500  silo3d(:,i,:) = wr_silo(:,i,:)*sin(pp) + & ! VY
1501  wphi_silo(:,i,:)*cos(pp)
1502  ENDDO
1503  label = "Vy"
1504  CALL write_to_silodb_3d(label, silo3d)
1505 
1506  ENDIF
1507 
1508  END SUBROUTINE dump_silo3d
1509 
1510 
1511  SUBROUTINE write_to_silodb_3d(label, array)
1512  CHARACTER*(*):: label
1513  REAL, DIMENSION(NXX, NYY, NZZ):: array
1514  INTEGER:: err, ierr, optlistid
1515  INTEGER:: ilen1
1516 
1517  ilen1 = len(trim(label))
1518  err = dbmkoptlist(6, optlistid)
1519  err = dbaddcopt(optlistid, dbopt_label, trim(label), ilen1)
1520  err = dbputqv1(dbfile3d, trim(label), ilen1, name3d, &
1521  ilen3d, array, dims3d, ndims3d, idum, 1, db_float, &
1522  db_nodecent, optlistid, ierr)
1523  err = dbfreeoptlist(optlistid)
1524 
1525  END SUBROUTINE write_to_silodb_3d
1526 
1527 
1528  SUBROUTINE write_to_silodb_2d(label, islice, array, isec)
1529  CHARACTER*(*):: label
1530  REAL, DIMENSION(NXX, NYY,NZZ), INTENT(IN) :: array
1531  INTEGER:: err, ierr, optlistid, isec
1532  INTEGER:: ilen1, ilen2, islice, iangle ! Islice in [1,NZZ]. PHI angle = 2*pi*(islice-1)/(NZZ-1)
1533  REAL:: phiout
1534  CHARACTER(LEN=3):: cdum, cdum1
1535  CHARACTER(LEN=50):: labelp1, labelp2
1536 
1537  CALL assert(islice<=nzz,'Slice does not exist!')
1538  phiout = 360.*real(islice -1)/real(nzz - 1)
1539  iangle = int(phiout)
1540  IF (iangle < 10) THEN
1541  WRITE(cdum1,'(i1)') iangle
1542  cdum = '00'//trim(cdum1)
1543  ELSE IF (iangle > 9 .AND. iangle < 100) THEN
1544  WRITE(cdum1,'(i2)') iangle
1545  cdum = '0'//trim(cdum1)
1546  ELSE IF (iangle > 99 .AND. iangle< 1000) THEN
1547  WRITE(cdum1,'(i3)') iangle
1548  cdum = cdum1
1549  ELSE IF (lverbose) THEN
1550  WRITE(*,*) ' Change format in SILODB_2D!'
1551  ENDIF
1552 
1553  labelp1 = trim(label)//'_flux_'//trim(cdum)//'_deg'
1554  ilen1 = len(trim(labelp1))
1555  labelp2 = trim(label)//'_RCyl_'//trim(cdum)//'_deg'
1556  ilen2 = len(trim(labelp2))
1557  silo2d = array(:,islice,:)
1558 
1559  IF (isec == 1) THEN
1560 
1561  err = dbmkoptlist(6, optlistid)
1562  err = dbaddcopt(optlistid, dbopt_label, trim(labelp1), ilen1)
1563  err = dbputqv1(dbfile2d, trim(labelp1), ilen1, trim(name2d_1),& ! X-section in flux coords
1564  ilen2d_1, silo2d, dims2d, ndims2d, idum, 1, db_float, &
1565  db_nodecent, optlistid, ierr)
1566  err = dbfreeoptlist(optlistid)
1567 
1568  err = dbmkoptlist(6, optlistid)
1569  err = dbaddcopt(optlistid, dbopt_label, trim(labelp2), ilen2)
1570  err = dbputqv1(dbfile2d, trim(labelp2), ilen2, trim(name2d_3),& ! X-section in RCyl coords
1571  ilen2d_3, silo2d, dims2d, ndims2d, idum, 1, db_float, &
1572  db_nodecent, optlistid, ierr)
1573  err = dbfreeoptlist(optlistid)
1574 
1575  ELSE
1576 
1577  err = dbmkoptlist(6, optlistid)
1578  err = dbaddcopt(optlistid, dbopt_label, trim(labelp1), ilen1)
1579  err = dbputqv1(dbfile2d, trim(labelp1), ilen1, trim(name2d_2),& ! X-section in flux coords
1580  ilen2d_2, silo2d, dims2d, ndims2d, idum, 1, db_float, &
1581  db_nodecent, optlistid, ierr)
1582  err = dbfreeoptlist(optlistid)
1583 
1584  err = dbmkoptlist(6, optlistid)
1585  err = dbaddcopt(optlistid, dbopt_label, trim(labelp2), ilen2)
1586  err = dbputqv1(dbfile2d, trim(labelp2), ilen2, trim(name2d_4),& ! X-section in RCyl coords
1587  ilen2d_4, silo2d, dims2d, ndims2d, idum, 1, db_float, &
1588  db_nodecent, optlistid, ierr)
1589  err = dbfreeoptlist(optlistid)
1590 
1591  ENDIF
1592 
1593  END SUBROUTINE write_to_silodb_2d
1594 #endif
1595 
1596  END MODULE dumping_mod
siesta_namelist::lrestart
logical lrestart
Restart SIESTA from pervious run.
Definition: siesta_namelist.f90:128
shared_data::lverbose
logical lverbose
Use verbose screen output.
Definition: shared_data.f90:234
siesta_namelist::nrs
integer nrs
Number of radial grid points.
Definition: siesta_namelist.f90:177
metrics::zmax
real(dp) zmax
Maximum of the grid Z inside the last closed flux surface.
Definition: metrics.f90:81
quantities
This file contains subroutines for allocating and initializing curvilinear magnetic covariant and pre...
Definition: quantities.f90:11
fourier::b_dv
integer, parameter b_dv
Bit position of the f_dv flag.
Definition: fourier.f90:51
siesta_namelist::nzs
integer nzs
Number of vertical grid points.
Definition: siesta_namelist.f90:179
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
v3_utilities::assert
Definition: v3_utilities.f:55
fourier::f_sin
integer, parameter f_sin
Sine parity.
Definition: fourier.f90:27
siesta_namelist::nphis
integer nphis
Number of cylindrical phi planes.
Definition: siesta_namelist.f90:175
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::b_sum
integer, parameter b_sum
Bit position of the f_sum flag.
Definition: fourier.f90:53
fourier::b_du
integer, parameter b_du
Bit position of the f_du flag.
Definition: fourier.f90:49
metrics::rmax
real(dp) rmax
Maximum of the grid R inside the last closed flux surface.
Definition: metrics.f90:77
fourier::b_ds
integer, parameter b_ds
Bit position of the f_ds flag.
Definition: fourier.f90:47
metrics::rmin
real(dp) rmin
Minimum of the grid R inside the last closed flux surface.
Definition: metrics.f90:79
fourier::f_ds
integer, parameter f_ds
Radial derivative flag.
Definition: fourier.f90:34
fourier::f_du
integer, parameter f_du
Poloidal derivative flag.
Definition: fourier.f90:36
metrics::zmin
real(dp) zmin
Minimum of the grid Z inside the last closed flux surface.
Definition: metrics.f90:83
siesta_namelist::lresistive
logical lresistive
Use resistive perturbaton.
Definition: siesta_namelist.f90:126
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
shared_data
This file contains variables and parameters used by many modules in SIESTA.
Definition: shared_data.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