V3FIT
primed_grid.f
Go to the documentation of this file.
1 !*******************************************************************************
4 !
5 ! Note separating the Doxygen comment block here so detailed decription is
6 ! found in the Module not the file.
7 !
10 !*******************************************************************************
11  MODULE primed_grid
12  USE stel_kinds, ONLY: rprec
13  USE profiler
15 
16  IMPLICIT NONE
17 
18 !*******************************************************************************
19 ! DERIVED-TYPE DECLARATIONS
20 ! 1) primed grid base class
21 !
22 !*******************************************************************************
23 !-------------------------------------------------------------------------------
26 !-------------------------------------------------------------------------------
29  REAL (rprec) :: dvol
31  REAL (rprec) :: dv
32 
34  REAL (rprec), DIMENSION(:,:,:), POINTER :: x => null()
36  REAL (rprec), DIMENSION(:,:,:), POINTER :: y => null()
38  REAL (rprec), DIMENSION(:,:,:), POINTER :: z => null()
39 
41  REAL (rprec), DIMENSION(:,:,:), POINTER :: j_x => null()
43  REAL (rprec), DIMENSION(:,:,:), POINTER :: j_y => null()
45  REAL (rprec), DIMENSION(:,:,:), POINTER :: j_z => null()
46 
47  REAL (rprec), DIMENSION(:,:,:), POINTER :: vol => null() ! TEMP
48  END TYPE
49 
50  CONTAINS
51 !*******************************************************************************
52 ! CONSTRUCTION SUBROUTINES
53 !*******************************************************************************
54 !-------------------------------------------------------------------------------
65 !-------------------------------------------------------------------------------
66  FUNCTION primed_grid_construct(num_v, flags, siesta_file, &
67  & parallel, io_unit)
69 
70  IMPLICIT NONE
71 
72 ! Declare Arguments
73  TYPE (primed_grid_class), POINTER :: primed_grid_construct
74  INTEGER, INTENT(in) :: num_v
75  INTEGER, INTENT(in) :: flags
76  CHARACTER (len=*), INTENT(in) :: siesta_file
77  TYPE (bmw_parallel_context_class), INTENT(in) :: parallel
78  INTEGER, INTENT(in) :: io_unit
79 
80 ! local variables
81  REAL (rprec) :: start_time
82 
83 ! Start of executable code
84  start_time = profiler_get_start_time()
85 
86  IF (btest(flags, bmw_state_flags_ju)) THEN
88  & primed_grid_construct_ju(num_v, parallel)
89  ELSE IF (btest(flags, bmw_state_flags_jv)) THEN
91  & primed_grid_construct_jv(num_v, parallel)
92  ELSE IF (btest(flags, bmw_state_flags_siesta)) THEN
94  & primed_grid_construct_siesta(num_v, siesta_file, parallel)
95  ELSE
97  & primed_grid_construct_both(num_v, parallel)
98  END IF
99 
100  IF (parallel%offset .eq. 0) THEN
101  WRITE(io_unit,1000)
102  END IF
103 
104  CALL profiler_set_stop_time('primed_grid_construct', start_time)
105 
106 1000 FORMAT('Prime Grid Ready')
107 
108  END FUNCTION
109 
110 !-------------------------------------------------------------------------------
120 !-------------------------------------------------------------------------------
121  FUNCTION primed_grid_construct_both(num_v, parallel)
122  USE stel_constants, ONLY: twopi, mu0
123  USE read_wout_mod, ONLY: mnmax, mnmax_nyq, lasym, isigng, ns, &
124  & xm, xn, xm_nyq, xn_nyq, &
125  & rmnc, zmns, currumnc, currvmnc, &
126  & rmns, zmnc, currumns, currvmns!, gmnc ! TEMP ,gmnc
127 
128  IMPLICIT NONE
129 
130 ! Declare Arguments
132  INTEGER, INTENT(in) :: num_v
133  TYPE (bmw_parallel_context_class), INTENT(in) :: parallel
134 
135 ! local variables
136  REAL (rprec) :: start_time
137  INTEGER :: i
138  REAL (rprec) :: x
139  REAL (rprec) :: ds
140  INTEGER :: si
141  INTEGER :: ui
142  INTEGER :: vi
143  REAL (rprec) :: r
144  REAL (rprec) :: z
145  REAL (rprec) :: ru
146  REAL (rprec) :: zu
147  REAL (rprec) :: rv
148  REAL (rprec) :: zv
149  REAL (rprec) :: ju
150  REAL (rprec) :: jv
151  REAL (rprec) :: jr
152  REAL (rprec) :: jp
153  REAL (rprec), DIMENSION(:), ALLOCATABLE :: cosv
154  REAL (rprec), DIMENSION(:), ALLOCATABLE :: sinv
155  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: cosmu
156  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: sinmu
157  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: cosmu_nyq
158  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: sinmu_nyq
159  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: cosnv
160  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: sinnv
161  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: cosnv_nyq
162  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: sinnv_nyq
163  REAL (rprec), DIMENSION(:), ALLOCATABLE :: cosmn
164  REAL (rprec), DIMENSION(:), ALLOCATABLE :: sinmn
165  REAL (rprec), DIMENSION(:), ALLOCATABLE :: cosmn_nyq
166  REAL (rprec), DIMENSION(:), ALLOCATABLE :: sinmn_nyq
167  INTEGER :: num_v_use
168 
169 ! local parameters
170  INTEGER, PARAMETER :: num_u = 101
171  REAL (rprec), PARAMETER :: du = twopi/num_u
172 
173 ! Start of executable code
174  start_time = profiler_get_start_time()
175 
177 
178  IF (num_v .eq. 1) THEN
179  num_v_use = 360
180  ELSE
181  num_v_use = num_v
182  END IF
183 
184  ALLOCATE(primed_grid_construct_both%x(ns,num_u,num_v_use))
185  ALLOCATE(primed_grid_construct_both%y(ns,num_u,num_v_use))
186  ALLOCATE(primed_grid_construct_both%z(ns,num_u,num_v_use))
187 
188  ALLOCATE(primed_grid_construct_both%j_x(ns,num_u,num_v_use))
189  ALLOCATE(primed_grid_construct_both%j_y(ns,num_u,num_v_use))
190  ALLOCATE(primed_grid_construct_both%j_z(ns,num_u,num_v_use))
191 
192 ! ALLOCATE(primed_grid_construct_both%vol(ns,num_u,num_v_use)) !TEMP
193 
194  ALLOCATE(cosv(num_v_use))
195  ALLOCATE(sinv(num_v_use))
196  ALLOCATE(cosmu(mnmax,num_u))
197  ALLOCATE(sinmu(mnmax,num_u))
198  ALLOCATE(cosnv(mnmax,num_v_use))
199  ALLOCATE(sinnv(mnmax,num_v_use))
200  ALLOCATE(cosmu_nyq(mnmax_nyq,num_u))
201  ALLOCATE(sinmu_nyq(mnmax_nyq,num_u))
202  ALLOCATE(cosnv_nyq(mnmax_nyq,num_v_use))
203  ALLOCATE(sinnv_nyq(mnmax_nyq,num_v_use))
204 
205  ds = 1.0/(ns - 1.0)
206  primed_grid_construct_both%dv = twopi/num_v_use
207 
208  primed_grid_construct_both%dvol = isigng*ds*du &
209  & * primed_grid_construct_both%dv/(2.0*twopi)
210 
211 !$OMP PARALLEL
212 !$OMP& DEFAULT(SHARED)
213 !$OMP& PRIVATE(i, x, si, ui, vi, r, z, ru, rv, zu, zv, ju, jv, jr, jp, &
214 !$OMP& cosmn, sinmn, cosmn_nyq, sinmn_nyq)
215 
216 ! Multi process will do an all reduce so these arrays need to be initalized.
217  IF (parallel%stride .gt. 1) THEN
218 !$OMP WORKSHARE
225  cosv = 0.0
226  sinv = 0.0
227  cosmu = 0.0
228  sinmu = 0.0
229  cosnv = 0.0
230  sinnv = 0.0
231  cosmu_nyq = 0.0
232  sinmu_nyq = 0.0
233  cosnv_nyq = 0.0
234  sinnv_nyq = 0.0
235 !$OMP END WORKSHARE
236  END IF
237 
238 !$OMP DO
239 !$OMP& SCHEDULE(STATIC)
240  DO i = bmw_parallel_context_start(parallel, num_u), &
241  & bmw_parallel_context_end(parallel, num_u)
242  x = (i - 0.5)*du
243 ! WRITE (*,*) x ! TEMP
244  cosmu(:,i) = cos(xm*x)
245  sinmu(:,i) = sin(xm*x)
246  cosmu_nyq(:,i) = cos(xm_nyq*x)
247  sinmu_nyq(:,i) = sin(xm_nyq*x)
248  END DO
249 !$OMP END DO
250 
251 !$OMP DO
252 !$OMP& SCHEDULE(STATIC)
253  DO i = bmw_parallel_context_start(parallel, num_v_use), &
254  & bmw_parallel_context_end(parallel, num_v_use)
255  x = (i - 0.5)*primed_grid_construct_both%dv
256 ! WRITE (*,*) x ! TEMP
257  cosv(i) = cos(x)
258  sinv(i) = sin(x)
259  cosnv(:,i) = cos(xn*x)
260  sinnv(:,i) = sin(xn*x)
261  cosnv_nyq(:,i) = cos(xn_nyq*x)
262  sinnv_nyq(:,i) = sin(xn_nyq*x)
263  END DO
264 !$OMP END DO
265 
266  IF (parallel%stride .gt. 1) THEN
267 !$OMP SINGLE
268  CALL bmw_parallel_context_reduce(parallel, cosv)
269  CALL bmw_parallel_context_reduce(parallel, sinv)
270  CALL bmw_parallel_context_reduce(parallel, cosmu)
271  CALL bmw_parallel_context_reduce(parallel, sinmu)
272  CALL bmw_parallel_context_reduce(parallel, cosmu_nyq)
273  CALL bmw_parallel_context_reduce(parallel, sinmu_nyq)
274  CALL bmw_parallel_context_reduce(parallel, cosnv)
275  CALL bmw_parallel_context_reduce(parallel, sinnv)
276  CALL bmw_parallel_context_reduce(parallel, cosnv_nyq)
277  CALL bmw_parallel_context_reduce(parallel, sinnv_nyq)
278 !$OMP END SINGLE
279  END IF
280 
281  ALLOCATE(cosmn(mnmax))
282  ALLOCATE(sinmn(mnmax))
283  ALLOCATE(cosmn_nyq(mnmax_nyq))
284  IF (lasym) THEN
285  ALLOCATE(sinmn_nyq(mnmax_nyq))
286  END IF
287 
288 !$OMP DO
289 !$OMP& SCHEDULE(STATIC)
290  DO i = bmw_parallel_context_start(parallel, ns*num_u*num_v_use), &
291  & bmw_parallel_context_end(parallel, ns*num_u*num_v_use)
292  si = bmw_parallel_context_i(i, ns)
293  ui = bmw_parallel_context_j(i, ns, num_u)
294  vi = bmw_parallel_context_k(i, ns, num_u)
295 
296  cosmn = cosmu(:,ui)*cosnv(:,vi) + sinmu(:,ui)*sinnv(:,vi)
297  sinmn = sinmu(:,ui)*cosnv(:,vi) - cosmu(:,ui)*sinnv(:,vi)
298  cosmn_nyq = cosmu_nyq(:,ui)*cosnv_nyq(:,vi) &
299  & + sinmu_nyq(:,ui)*sinnv_nyq(:,vi)
300 ! x = (vi - 0.5)*primed_grid_construct_both%dv
301  r = sum(rmnc(:,si)*cosmn(:))
302  z = sum(zmns(:,si)*sinmn(:))
303 
304 ! IF (si .eq. ns) WRITE (20000,*) '{',r,',',x,',',z,'},'
305 
306  ru = -sum(xm*rmnc(:,si)*sinmn(:))
307  rv = sum(xn*rmnc(:,si)*sinmn(:))
308  zu = sum(xm*zmns(:,si)*cosmn(:))
309  zv = -sum(xn*zmns(:,si)*cosmn(:))
310 
311  ju = sum(currumnc(:,si)*cosmn_nyq(:))*mu0
312  jv = sum(currvmnc(:,si)*cosmn_nyq(:))*mu0
313 
314 ! primed_grid_construct_both%vol(si,ui,vi) = &
315 ! & SUM(gmnc(:,si)*cosmn_nyq(:)) ! TEMP
316 
317  IF (lasym) THEN
318  sinmn_nyq = sinmu_nyq(:,ui)*cosnv_nyq(:,vi) &
319  & - cosmu_nyq(:,ui)*sinnv_nyq(:,vi)
320 
321  r = r + sum(rmns(:,si)*sinmn(:))
322  z = z + sum(zmnc(:,si)*cosmn(:))
323 
324  ru = ru + sum(xm*rmns(:,si)*cosmn(:))
325  rv = rv - sum(xn*rmns(:,si)*cosmn(:))
326  zu = zu - sum(xm*zmnc(:,si)*sinmn(:))
327  zv = zv + sum(xn*zmnc(:,si)*sinmn(:))
328 
329  ju = ju + sum(currumns(:,si)*sinmn_nyq(:))*mu0
330  jv = jv + sum(currvmns(:,si)*sinmn_nyq(:))*mu0
331  END IF
332 
333  jr = ju*ru + jv*rv
334  jp = jv*r
335  primed_grid_construct_both%j_z(si,ui,vi) = ju*zu + jv*zv
336 
337  IF (si .eq. 1 .or. si .eq. ns) THEN
338  jr = jr/2.0
339  jp = jp/2.0
340  primed_grid_construct_both%j_z(si,ui,vi) = &
341  & primed_grid_construct_both%j_z(si,ui,vi)/2.0
342 
343 ! primed_grid_construct_both%vol(si,ui,vi) = &
344 ! & primed_grid_construct_both%vol(si,ui,vi)/2.0 ! TEMP
345  END IF
346 
347  primed_grid_construct_both%j_x(si,ui,vi) = jr*cosv(vi) &
348  & - jp*sinv(vi)
349  primed_grid_construct_both%j_y(si,ui,vi) = jr*sinv(vi) &
350  & + jp*cosv(vi)
351 
352  primed_grid_construct_both%x(si,ui,vi) = r*cosv(vi)
353  primed_grid_construct_both%y(si,ui,vi) = r*sinv(vi)
354  primed_grid_construct_both%z(si,ui,vi) = z
355  END DO
356 !$OMP END DO
357 
358  DEALLOCATE(cosmn)
359  DEALLOCATE(sinmn)
360  DEALLOCATE(cosmn_nyq)
361  IF (lasym) THEN
362  DEALLOCATE(sinmn_nyq)
363  END IF
364 !$OMP END PARALLEL
365 
366  DEALLOCATE(cosv)
367  DEALLOCATE(sinv)
368  DEALLOCATE(cosmu)
369  DEALLOCATE(sinmu)
370  DEALLOCATE(cosnv)
371  DEALLOCATE(sinnv)
372  DEALLOCATE(cosmu_nyq)
373  DEALLOCATE(sinmu_nyq)
374  DEALLOCATE(cosnv_nyq)
375  DEALLOCATE(sinnv_nyq)
376 
377 ! Multi process did not fill out the entire array. Get the missing pieces from
378 ! the other processes.
379  IF (parallel%stride .gt. 1) THEN
380  CALL bmw_parallel_context_reduce(parallel, &
382  CALL bmw_parallel_context_reduce(parallel, &
384  CALL bmw_parallel_context_reduce(parallel, &
386  CALL bmw_parallel_context_reduce(parallel, &
388  CALL bmw_parallel_context_reduce(parallel, &
390  CALL bmw_parallel_context_reduce(parallel, &
392  END IF
393 
394  CALL profiler_set_stop_time('primed_grid_construct_both', &
395  & start_time)
396 
397  END FUNCTION
398 
399 !-------------------------------------------------------------------------------
411 !-------------------------------------------------------------------------------
412  FUNCTION primed_grid_construct_ju(num_v, parallel)
413  USE stel_constants, ONLY: twopi, mu0
414  USE read_wout_mod, ONLY: mnmax, mnmax_nyq, lasym, isigng, ns, &
415  & rmnc, rmns, currvmnc, &
416  & zmnc, zmns, currvmns, &
417  & bsupumnc, bsupvmnc, bsupumns, bsupvmns, &
418  & xm, xn, xm_nyq, xn_nyq, presf
419 
420  IMPLICIT NONE
421 
422 ! Declare Arguments
424  INTEGER, INTENT(in) :: num_v
425  TYPE (bmw_parallel_context_class), INTENT(in) :: parallel
426 
427 ! local variables
428  REAL (rprec) :: start_time
429  INTEGER :: i
430  REAL (rprec) :: x
431  REAL (rprec) :: ds
432  INTEGER :: si
433  INTEGER :: ui
434  INTEGER :: vi
435  REAL (rprec) :: r
436  REAL (rprec) :: z
437  REAL (rprec) :: ru
438  REAL (rprec) :: zu
439  REAL (rprec) :: rv
440  REAL (rprec) :: zv
441  REAL (rprec) :: ju
442  REAL (rprec) :: jv
443  REAL (rprec) :: bu
444  REAL (rprec) :: bv
445  REAL (rprec) :: jr
446  REAL (rprec) :: jp
447  REAL (rprec), DIMENSION(:), ALLOCATABLE :: cosv
448  REAL (rprec), DIMENSION(:), ALLOCATABLE :: sinv
449  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: cosmu
450  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: sinmu
451  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: cosmu_nyq
452  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: sinmu_nyq
453  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: cosnv
454  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: sinnv
455  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: cosnv_nyq
456  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: sinnv_nyq
457  REAL (rprec), DIMENSION(:), ALLOCATABLE :: cosmn
458  REAL (rprec), DIMENSION(:), ALLOCATABLE :: sinmn
459  REAL (rprec), DIMENSION(:), ALLOCATABLE :: cosmn_nyq
460  REAL (rprec), DIMENSION(:), ALLOCATABLE :: sinmn_nyq
461  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: rmnch
462  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: rmnsh
463  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: zmnch
464  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: zmnsh
465  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: currvmnch
466  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: currvmnsh
467  REAL (rprec), DIMENSION(:), ALLOCATABLE :: p_prime
468  INTEGER :: num_v_use
469 
470 ! local parameters
471  INTEGER, PARAMETER :: num_u = 101
472  REAL (rprec), PARAMETER :: du = twopi/num_u
473 
474 ! Start of executable code
475  start_time = profiler_get_start_time()
476 
477  ALLOCATE(primed_grid_construct_ju)
478 
479  IF (num_v .eq. 1) THEN
480  num_v_use = 360
481  ELSE
482  num_v_use = num_v
483  END IF
484 
485  ALLOCATE(primed_grid_construct_ju%x(ns - 1,num_u,num_v_use))
486  ALLOCATE(primed_grid_construct_ju%y(ns - 1,num_u,num_v_use))
487  ALLOCATE(primed_grid_construct_ju%z(ns - 1,num_u,num_v_use))
488 
489  ALLOCATE(primed_grid_construct_ju%j_x(ns - 1,num_u,num_v_use))
490  ALLOCATE(primed_grid_construct_ju%j_y(ns - 1,num_u,num_v_use))
491  ALLOCATE(primed_grid_construct_ju%j_z(ns - 1,num_u,num_v_use))
492 
493  ALLOCATE(cosv(num_v_use))
494  ALLOCATE(sinv(num_v_use))
495  ALLOCATE(cosmu(mnmax,num_u))
496  ALLOCATE(sinmu(mnmax,num_u))
497  ALLOCATE(cosnv(mnmax,num_v_use))
498  ALLOCATE(sinnv(mnmax,num_v_use))
499  ALLOCATE(cosmu_nyq(mnmax_nyq,num_u))
500  ALLOCATE(sinmu_nyq(mnmax_nyq,num_u))
501  ALLOCATE(cosnv_nyq(mnmax_nyq,num_v_use))
502  ALLOCATE(sinnv_nyq(mnmax_nyq,num_v_use))
503 
504  ALLOCATE(rmnch(mnmax,ns - 1))
505  ALLOCATE(zmnsh(mnmax,ns - 1))
506  ALLOCATE(currvmnch(mnmax_nyq,ns - 1))
507 
508  IF (lasym) THEN
509  ALLOCATE(rmnsh(mnmax,ns - 1))
510  ALLOCATE(zmnch(mnmax,ns - 1))
511  ALLOCATE(currvmnsh(mnmax_nyq,ns - 1))
512  END IF
513 
514  ALLOCATE(p_prime(ns - 1))
515 
516  ds = 1.0/(ns - 1.0)
517  primed_grid_construct_ju%dv = twopi/num_v_use
518 
519  primed_grid_construct_ju%dvol = isigng*ds*du &
520  & * primed_grid_construct_ju%dv/(2.0*twopi)
521 
522 !$OMP PARALLEL
523 !$OMP& DEFAULT(SHARED)
524 !$OMP& PRIVATE(i, x, si, ui, vi, r, z, ru, rv, zu, zv, ju, jv, bu, bv, &
525 !$OMP& cosmn, sinmn, cosmn_nyq, sinmn_nyq, jr, jp)
526 
527 ! Multi process will do an all reduce so these arrays need to be initalized.
528  IF (parallel%stride .gt. 1) THEN
529 !$OMP WORKSHARE
533  primed_grid_construct_ju%j_x = 0.0
534  primed_grid_construct_ju%j_y = 0.0
535  primed_grid_construct_ju%j_z = 0.0
536  cosv = 0.0
537  sinv = 0.0
538  cosmu = 0.0
539  sinmu = 0.0
540  cosnv = 0.0
541  sinnv = 0.0
542  cosmu_nyq = 0.0
543  sinmu_nyq = 0.0
544  cosnv_nyq = 0.0
545  sinnv_nyq = 0.0
546  rmnch = 0.0
547  zmnsh = 0.0
548  currvmnch = 0.0
549  p_prime = 0.0
550 !$OMP END WORKSHARE
551  IF (lasym) THEN
552 !$OMP WORKSHARE
553  rmnsh = 0.0
554  zmnch = 0.0
555  currvmnsh = 0.0
556 !$OMP END WORKSHARE
557  END IF
558  END IF
559 
560 !$OMP DO
561 !$OMP& SCHEDULE(STATIC)
562  DO i = bmw_parallel_context_start(parallel, ns - 1), &
563  & bmw_parallel_context_end(parallel, ns - 1)
564  rmnch(:,i) = (rmnc(:,i + 1) + rmnc(:,i))/2.0
565  zmnsh(:,i) = (zmns(:,i + 1) + zmns(:,i))/2.0
566 
567  currvmnch(:,i) = (currvmnc(:,i + 1) + currvmnc(:,i))/2.0
568 
569  p_prime(i) = (presf(i + 1) - presf(i))/ds
570 
571  IF (lasym) THEN
572  rmnsh(:,i) = (rmns(:,i + 1) + rmns(:,i))/2.0
573  zmnch(:,i) = (zmnc(:,i + 1) + zmnc(:,i))/2.0
574 
575  currvmnsh(:,i) = (currvmns(:,i + 1) + currvmns(:,i))/2.0
576  END IF
577  END DO
578 !$OMP END DO
579 
580 !$OMP DO
581 !$OMP& SCHEDULE(STATIC)
582  DO i = bmw_parallel_context_start(parallel, num_u), &
583  & bmw_parallel_context_end(parallel, num_u)
584  x = (i - 0.5)*du
585  cosmu(:,i) = cos(xm*x)
586  sinmu(:,i) = sin(xm*x)
587  cosmu_nyq(:,i) = cos(xm_nyq*x)
588  sinmu_nyq(:,i) = sin(xm_nyq*x)
589  END DO
590 !$OMP END DO
591 
592 !$OMP DO
593 !$OMP& SCHEDULE(STATIC)
594  DO i = bmw_parallel_context_start(parallel, num_v_use), &
595  & bmw_parallel_context_end(parallel, num_v_use)
596  x = (i - 0.5)*primed_grid_construct_ju%dv
597  cosv(i) = cos(x)
598  sinv(i) = sin(x)
599  cosnv(:,i) = cos(xn*x)
600  sinnv(:,i) = sin(xn*x)
601  cosnv_nyq(:,i) = cos(xn_nyq*x)
602  sinnv_nyq(:,i) = sin(xn_nyq*x)
603  END DO
604 !$OMP END DO
605 
606  IF (parallel%stride .gt. 1) THEN
607 !$OMP SINGLE
608  CALL bmw_parallel_context_reduce(parallel, cosv)
609  CALL bmw_parallel_context_reduce(parallel, sinv)
610  CALL bmw_parallel_context_reduce(parallel, cosmu)
611  CALL bmw_parallel_context_reduce(parallel, sinmu)
612  CALL bmw_parallel_context_reduce(parallel, cosmu_nyq)
613  CALL bmw_parallel_context_reduce(parallel, sinmu_nyq)
614  CALL bmw_parallel_context_reduce(parallel, cosnv)
615  CALL bmw_parallel_context_reduce(parallel, sinnv)
616  CALL bmw_parallel_context_reduce(parallel, cosnv_nyq)
617  CALL bmw_parallel_context_reduce(parallel, sinnv_nyq)
618  CALL bmw_parallel_context_reduce(parallel, rmnch)
619  CALL bmw_parallel_context_reduce(parallel, zmnsh)
620  CALL bmw_parallel_context_reduce(parallel, currvmnch)
621  CALL bmw_parallel_context_reduce(parallel, p_prime)
622  IF (lasym) THEN
623  CALL bmw_parallel_context_reduce(parallel, rmnsh)
624  CALL bmw_parallel_context_reduce(parallel, zmnch)
625  CALL bmw_parallel_context_reduce(parallel, currvmnsh)
626  END IF
627 !$OMP END SINGLE
628  END IF
629 
630  ALLOCATE(cosmn(mnmax))
631  ALLOCATE(sinmn(mnmax))
632  ALLOCATE(cosmn_nyq(mnmax_nyq))
633  IF (lasym) THEN
634  ALLOCATE(sinmn_nyq(mnmax_nyq))
635  END IF
636 
637 !$OMP DO
638 !$OMP& SCHEDULE(STATIC)
639  DO i = bmw_parallel_context_start(parallel, ns*num_u*num_v_use), &
640  & bmw_parallel_context_end(parallel, ns*num_u*num_v_use)
641  si = bmw_parallel_context_i(i, ns)
642  ui = bmw_parallel_context_j(i, ns, num_u)
643  vi = bmw_parallel_context_k(i, ns, num_u)
644 
645  cosmn = cosmu(:,ui)*cosnv(:,vi) + sinmu(:,ui)*sinnv(:,vi)
646  sinmn = sinmu(:,ui)*cosnv(:,vi) - cosmu(:,ui)*sinnv(:,vi)
647  cosmn_nyq = cosmu_nyq(:,ui)*cosnv_nyq(:,vi) &
648  & + sinmu_nyq(:,ui)*sinnv_nyq(:,vi)
649 
650  r = sum(rmnch(:,si)*cosmn(:))
651  z = sum(zmnsh(:,si)*sinmn(:))
652 
653  ru = -sum(xm*rmnch(:,si)*sinmn(:))
654  rv = sum(xn*rmnch(:,si)*sinmn(:))
655  zu = sum(xm*zmnsh(:,si)*cosmn(:))
656  zv = -sum(xn*zmnsh(:,si)*cosmn(:))
657 
658  jv = sum(currvmnch(:,si)*cosmn_nyq(:))
659 
660  bu = sum(bsupumnc(:,si + 1)*cosmn_nyq(:))
661  bv = sum(bsupvmnc(:,si + 1)*cosmn_nyq(:))
662 
663  IF (lasym) THEN
664  sinmn_nyq = sinmu_nyq(:,ui)*cosnv_nyq(:,vi) &
665  & - cosmu_nyq(:,ui)*sinnv_nyq(:,vi)
666 
667  r = r + sum(rmnsh(:,si)*sinmn(:))
668  z = z + sum(zmnch(:,si)*cosmn(:))
669 
670  ru = ru + sum(xm*rmnsh(:,si)*cosmn(:))
671  rv = rv - sum(xn*rmnsh(:,si)*cosmn(:))
672  zu = zu - sum(xm*zmnch(:,si)*sinmn(:))
673  zv = zv + sum(xn*zmnch(:,si)*sinmn(:))
674 
675  jv = jv + sum(currvmnsh(:,si)*sinmn_nyq(:))
676 
677  bu = bu + sum(bsupumns(:,si + 1)*sinmn_nyq(:))
678  bv = bv + sum(bsupvmns(:,si + 1)*sinmn_nyq(:))
679  END IF
680 
681  ju = (p_prime(si) + jv*bu)/bv*mu0
682  jv = jv*mu0
683 
684  jr = ju*ru + jv*rv
685  jp = jv*r
686 
687  primed_grid_construct_ju%j_z(si,ui,vi) = ju*zu + jv*zv
688  primed_grid_construct_ju%j_x(si,ui,vi) = jr*cosv(vi) &
689  & - jp*sinv(vi)
690  primed_grid_construct_ju%j_y(si,ui,vi) = jr*sinv(vi) &
691  & + jp*cosv(vi)
692 
693  primed_grid_construct_ju%x(si,ui,vi) = r*cosv(vi)
694  primed_grid_construct_ju%y(si,ui,vi) = r*sinv(vi)
695  primed_grid_construct_ju%z(si,ui,vi) = z
696  END DO
697 !$OMP END DO
698 
699  DEALLOCATE(cosmn)
700  DEALLOCATE(sinmn)
701  DEALLOCATE(cosmn_nyq)
702  IF (lasym) THEN
703  DEALLOCATE(sinmn_nyq)
704  END IF
705 !$OMP END PARALLEL
706 
707  DEALLOCATE(cosv)
708  DEALLOCATE(sinv)
709  DEALLOCATE(cosmu)
710  DEALLOCATE(sinmu)
711  DEALLOCATE(cosnv)
712  DEALLOCATE(sinnv)
713  DEALLOCATE(cosmu_nyq)
714  DEALLOCATE(sinmu_nyq)
715  DEALLOCATE(cosnv_nyq)
716  DEALLOCATE(sinnv_nyq)
717 
718  DEALLOCATE(rmnch)
719  DEALLOCATE(zmnsh)
720  DEALLOCATE(currvmnch)
721  DEALLOCATE(p_prime)
722  IF (lasym) THEN
723  DEALLOCATE(rmnsh)
724  DEALLOCATE(zmnch)
725  DEALLOCATE(currvmnsh)
726  END IF
727 
728 ! Multi process did not fill out the entire array. Get the missing pieces from
729 ! the other processes.
730  IF (parallel%stride .gt. 1) THEN
731  CALL bmw_parallel_context_reduce(parallel, &
733  CALL bmw_parallel_context_reduce(parallel, &
735  CALL bmw_parallel_context_reduce(parallel, &
737  CALL bmw_parallel_context_reduce(parallel, &
739  CALL bmw_parallel_context_reduce(parallel, &
741  CALL bmw_parallel_context_reduce(parallel, &
743  END IF
744 
745  CALL profiler_set_stop_time('primed_grid_construct_ju', &
746  & start_time)
747 
748  END FUNCTION
749 
750 !-------------------------------------------------------------------------------
762 !-------------------------------------------------------------------------------
763  FUNCTION primed_grid_construct_jv(num_v, parallel)
764  USE stel_constants, ONLY: twopi, mu0
765  USE read_wout_mod, ONLY: mnmax, mnmax_nyq, lasym, isigng, ns, &
766  & rmnc, rmns, currumnc, &
767  & zmnc, zmns, currumns, &
768  & bsupumnc, bsupvmnc, bsupumns, bsupvmns, &
769  & xm, xn, xm_nyq, xn_nyq, presf
770 
771  IMPLICIT NONE
772 
773 ! Declare Arguments
775  INTEGER, INTENT(in) :: num_v
776  TYPE (bmw_parallel_context_class), INTENT(in) :: parallel
777 
778 ! local variables
779  REAL (rprec) :: start_time
780  INTEGER :: i
781  REAL (rprec) :: x
782  REAL (rprec) :: ds
783  INTEGER :: si
784  INTEGER :: ui
785  INTEGER :: vi
786  REAL (rprec) :: r
787  REAL (rprec) :: z
788  REAL (rprec) :: ru
789  REAL (rprec) :: zu
790  REAL (rprec) :: rv
791  REAL (rprec) :: zv
792  REAL (rprec) :: ju
793  REAL (rprec) :: jv
794  REAL (rprec) :: bu
795  REAL (rprec) :: bv
796  REAL (rprec) :: jr
797  REAL (rprec) :: jp
798  REAL (rprec), DIMENSION(:), ALLOCATABLE :: cosv
799  REAL (rprec), DIMENSION(:), ALLOCATABLE :: sinv
800  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: cosmu
801  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: sinmu
802  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: cosmu_nyq
803  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: sinmu_nyq
804  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: cosnv
805  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: sinnv
806  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: cosnv_nyq
807  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: sinnv_nyq
808  REAL (rprec), DIMENSION(:), ALLOCATABLE :: cosmn
809  REAL (rprec), DIMENSION(:), ALLOCATABLE :: sinmn
810  REAL (rprec), DIMENSION(:), ALLOCATABLE :: cosmn_nyq
811  REAL (rprec), DIMENSION(:), ALLOCATABLE :: sinmn_nyq
812  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: rmnch
813  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: rmnsh
814  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: zmnch
815  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: zmnsh
816  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: currumnch
817  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: currumnsh
818  REAL (rprec), DIMENSION(:), ALLOCATABLE :: p_prime
819  INTEGER :: num_v_use
820 
821 ! local parameters
822  INTEGER, PARAMETER :: num_u = 101
823  REAL (rprec), PARAMETER :: du = twopi/num_u
824 
825 ! Start of executable code
826  start_time = profiler_get_start_time()
827 
828  ALLOCATE(primed_grid_construct_jv)
829 
830  IF (num_v .eq. 1) THEN
831  num_v_use = 360
832  ELSE
833  num_v_use = num_v
834  END IF
835 
836  ALLOCATE(primed_grid_construct_jv%x(ns - 1,num_u,num_v_use))
837  ALLOCATE(primed_grid_construct_jv%y(ns - 1,num_u,num_v_use))
838  ALLOCATE(primed_grid_construct_jv%z(ns - 1,num_u,num_v_use))
839 
840  ALLOCATE(primed_grid_construct_jv%j_x(ns - 1,num_u,num_v_use))
841  ALLOCATE(primed_grid_construct_jv%j_y(ns - 1,num_u,num_v_use))
842  ALLOCATE(primed_grid_construct_jv%j_z(ns - 1,num_u,num_v_use))
843 
844  ALLOCATE(cosv(num_v_use))
845  ALLOCATE(sinv(num_v_use))
846  ALLOCATE(cosmu(mnmax,num_u))
847  ALLOCATE(sinmu(mnmax,num_u))
848  ALLOCATE(cosnv(mnmax,num_v_use))
849  ALLOCATE(sinnv(mnmax,num_v_use))
850  ALLOCATE(cosmu_nyq(mnmax_nyq,num_u))
851  ALLOCATE(sinmu_nyq(mnmax_nyq,num_u))
852  ALLOCATE(cosnv_nyq(mnmax_nyq,num_v_use))
853  ALLOCATE(sinnv_nyq(mnmax_nyq,num_v_use))
854 
855  ALLOCATE(rmnch(mnmax,ns - 1))
856  ALLOCATE(zmnsh(mnmax,ns - 1))
857  ALLOCATE(currumnch(mnmax_nyq,ns - 1))
858 
859  IF (lasym) THEN
860  ALLOCATE(rmnsh(mnmax,ns - 1))
861  ALLOCATE(zmnch(mnmax,ns - 1))
862  ALLOCATE(currumnsh(mnmax_nyq,ns - 1))
863  END IF
864 
865  ALLOCATE(p_prime(ns - 1))
866 
867  ds = 1.0/(ns - 1.0)
868  primed_grid_construct_jv%dv = twopi/num_v_use
869 
870  primed_grid_construct_jv%dvol = isigng*ds*du &
871  & * primed_grid_construct_jv%dv/(2.0*twopi)
872 
873 !$OMP PARALLEL
874 !$OMP& DEFAULT(SHARED)
875 !$OMP& PRIVATE(i, x, si, ui, vi, r, z, ru, rv, zu, zv, ju, jv, bu, bv, &
876 !$OMP& cosmn, sinmn, cosmn_nyq, sinmn_nyq, jr, jp)
877 
878 ! Multi process will do an all reduce so these arrays need to be initalized.
879  IF (parallel%stride .gt. 1) THEN
880 !$OMP WORKSHARE
884  primed_grid_construct_jv%j_x = 0.0
885  primed_grid_construct_jv%j_y = 0.0
886  primed_grid_construct_jv%j_z = 0.0
887  cosv = 0.0
888  sinv = 0.0
889  cosmu = 0.0
890  sinmu = 0.0
891  cosnv = 0.0
892  sinnv = 0.0
893  cosmu_nyq = 0.0
894  sinmu_nyq = 0.0
895  cosnv_nyq = 0.0
896  sinnv_nyq = 0.0
897  rmnch = 0.0
898  zmnsh = 0.0
899  currumnch = 0.0
900  p_prime = 0.0
901 !$OMP END WORKSHARE
902  IF (lasym) THEN
903  rmnsh = 0.0
904  zmnch = 0.0
905  currumnsh = 0.0
906  END IF
907  END IF
908 
909 !$OMP DO
910 !$OMP& SCHEDULE(STATIC)
911  DO i = bmw_parallel_context_start(parallel, ns - 1), &
912  & bmw_parallel_context_end(parallel, ns - 1)
913  si = bmw_parallel_context_i(i, ns)
914 
915  rmnch(:,si) = (rmnc(:,si + 1) + rmnc(:,si))/2.0
916  zmnsh(:,si) = (zmns(:,si + 1) + zmns(:,si))/2.0
917 
918  currumnch(:,si) = (currumnc(:,si + 1) + currumnc(:,si))/2.0
919 
920  p_prime(si) = (presf(si + 1) - presf(si))/ds
921 
922  IF (lasym) THEN
923  rmnsh(:,si) = (rmns(:,si + 1) + rmns(:,si))/2.0
924  zmnch(:,si) = (zmnc(:,si + 1) + zmnc(:,si))/2.0
925 
926  currumnsh(:,si) = (currumns(:,si + 1) + currumns(:,si))/2.0
927  END IF
928  END DO
929 !$OMP END DO
930 
931 !$OMP DO
932 !$OMP& SCHEDULE(STATIC)
933  DO i = bmw_parallel_context_start(parallel, num_u), &
934  & bmw_parallel_context_end(parallel, num_u)
935  x = (i - 0.5)*du
936  cosmu(:,i) = cos(xm*x)
937  sinmu(:,i) = sin(xm*x)
938  cosmu_nyq(:,i) = cos(xm_nyq*x)
939  sinmu_nyq(:,i) = sin(xm_nyq*x)
940  END DO
941 !$OMP END DO
942 
943 !$OMP DO
944 !$OMP& SCHEDULE(STATIC)
945  DO i = bmw_parallel_context_start(parallel, num_v_use), &
946  & bmw_parallel_context_end(parallel, num_v_use)
947  x = (i - 0.5)*primed_grid_construct_jv%dv
948  cosv(i) = cos(x)
949  sinv(i) = sin(x)
950  cosnv(:,i) = cos(xn*x)
951  sinnv(:,i) = sin(xn*x)
952  cosnv_nyq(:,i) = cos(xn_nyq*x)
953  sinnv_nyq(:,i) = sin(xn_nyq*x)
954  END DO
955 !$OMP END DO
956 
957  IF (parallel%stride .gt. 1) THEN
958 !$OMP SINGLE
959  CALL bmw_parallel_context_reduce(parallel, cosv)
960  CALL bmw_parallel_context_reduce(parallel, sinv)
961  CALL bmw_parallel_context_reduce(parallel, cosmu)
962  CALL bmw_parallel_context_reduce(parallel, sinmu)
963  CALL bmw_parallel_context_reduce(parallel, cosmu_nyq)
964  CALL bmw_parallel_context_reduce(parallel, sinmu_nyq)
965  CALL bmw_parallel_context_reduce(parallel, cosnv)
966  CALL bmw_parallel_context_reduce(parallel, sinnv)
967  CALL bmw_parallel_context_reduce(parallel, cosnv_nyq)
968  CALL bmw_parallel_context_reduce(parallel, sinnv_nyq)
969  CALL bmw_parallel_context_reduce(parallel, rmnch)
970  CALL bmw_parallel_context_reduce(parallel, zmnsh)
971  CALL bmw_parallel_context_reduce(parallel, currumnch)
972  CALL bmw_parallel_context_reduce(parallel, p_prime)
973  IF (lasym) THEN
974  CALL bmw_parallel_context_reduce(parallel, rmnsh)
975  CALL bmw_parallel_context_reduce(parallel, zmnch)
976  CALL bmw_parallel_context_reduce(parallel, currumnsh)
977  END IF
978 !$OMP END SINGLE
979  END IF
980 
981  ALLOCATE(cosmn(mnmax))
982  ALLOCATE(sinmn(mnmax))
983  ALLOCATE(cosmn_nyq(mnmax_nyq))
984  IF (lasym) THEN
985  ALLOCATE(sinmn_nyq(mnmax_nyq))
986  END IF
987 
988 !$OMP DO
989 !$OMP& SCHEDULE(STATIC)
990  DO i = bmw_parallel_context_start(parallel, ns*num_u*num_v_use), &
991  & bmw_parallel_context_end(parallel, ns*num_u*num_v_use)
992  si = bmw_parallel_context_i(i, ns)
993  ui = bmw_parallel_context_j(i, ns, num_u)
994  vi = bmw_parallel_context_k(i, ns, num_u)
995 
996  cosmn = cosmu(:,ui)*cosnv(:,vi) + sinmu(:,ui)*sinnv(:,vi)
997  sinmn = sinmu(:,ui)*cosnv(:,vi) - cosmu(:,ui)*sinnv(:,vi)
998  cosmn_nyq = cosmu_nyq(:,ui)*cosnv_nyq(:,vi) &
999  & + sinmu_nyq(:,ui)*sinnv_nyq(:,vi)
1000 
1001  r = sum(rmnch(:,si)*cosmn(:))
1002  z = sum(zmnsh(:,si)*sinmn(:))
1003 
1004  ru = -sum(xm*rmnch(:,si)*sinmn(:))
1005  rv = sum(xn*rmnch(:,si)*sinmn(:))
1006  zu = sum(xm*zmnsh(:,si)*cosmn(:))
1007  zv = -sum(xn*zmnsh(:,si)*cosmn(:))
1008 
1009  ju = sum(currumnch(:,si)*cosmn_nyq(:))
1010 
1011  bu = sum(bsupumnc(:,si + 1)*cosmn_nyq(:))
1012  bv = sum(bsupvmnc(:,si + 1)*cosmn_nyq(:))
1013 
1014  IF (lasym) THEN
1015  sinmn_nyq = sinmu_nyq(:,ui)*cosnv_nyq(:,vi) &
1016  & - cosmu_nyq(:,ui)*sinnv_nyq(:,vi)
1017 
1018  r = r + sum(rmnsh(:,si)*sinmn(:))
1019  z = z + sum(zmnch(:,si)*cosmn(:))
1020 
1021  ru = ru + sum(xm*rmnsh(:,si)*cosmn(:))
1022  rv = rv - sum(xn*rmnsh(:,si)*cosmn(:))
1023  zu = zu - sum(xm*zmnch(:,si)*sinmn(:))
1024  zv = zv + sum(xn*zmnch(:,si)*sinmn(:))
1025 
1026  ju = ju + sum(currumnsh(:,si)*sinmn_nyq(:))
1027 
1028  bu = bu + sum(bsupumns(:,si + 1)*sinmn_nyq(:))
1029  bv = bv + sum(bsupvmns(:,si + 1)*sinmn_nyq(:))
1030  END IF
1031 
1032  jv = (ju*bv - p_prime(si))/bu*mu0
1033  ju = ju*mu0
1034 
1035  jr = ju*ru + jv*rv
1036  jp = jv*r
1037 
1038  primed_grid_construct_jv%j_z(si,ui,vi) = ju*zu + jv*zv
1039  primed_grid_construct_jv%j_x(si,ui,vi) = jr*cosv(vi) &
1040  & - jp*sinv(vi)
1041  primed_grid_construct_jv%j_y(si,ui,vi) = jr*sinv(vi) &
1042  & + jp*cosv(vi)
1043 
1044  primed_grid_construct_jv%x(si,ui,vi) = r*cosv(vi)
1045  primed_grid_construct_jv%y(si,ui,vi) = r*sinv(vi)
1046  primed_grid_construct_jv%z(si,ui,vi) = z
1047  END DO
1048 !$OMP END DO
1049 
1050  DEALLOCATE(cosmn)
1051  DEALLOCATE(sinmn)
1052  DEALLOCATE(cosmn_nyq)
1053  IF (lasym) THEN
1054  DEALLOCATE(sinmn_nyq)
1055  END IF
1056 !$OMP END PARALLEL
1057 
1058  DEALLOCATE(cosv)
1059  DEALLOCATE(sinv)
1060  DEALLOCATE(cosmu)
1061  DEALLOCATE(sinmu)
1062  DEALLOCATE(cosnv)
1063  DEALLOCATE(sinnv)
1064  DEALLOCATE(cosmu_nyq)
1065  DEALLOCATE(sinmu_nyq)
1066  DEALLOCATE(cosnv_nyq)
1067  DEALLOCATE(sinnv_nyq)
1068 
1069  DEALLOCATE(rmnch)
1070  DEALLOCATE(zmnsh)
1071  DEALLOCATE(currumnch)
1072  DEALLOCATE(p_prime)
1073  IF (lasym) THEN
1074  DEALLOCATE(rmnsh)
1075  DEALLOCATE(zmnch)
1076  DEALLOCATE(currumnsh)
1077  END IF
1078 
1079 ! Multi process did not fill out the entire array. Get the missing pieces from
1080 ! the other processes.
1081  IF (parallel%stride .gt. 1) THEN
1082  CALL bmw_parallel_context_reduce(parallel, &
1084  CALL bmw_parallel_context_reduce(parallel, &
1086  CALL bmw_parallel_context_reduce(parallel, &
1088  CALL bmw_parallel_context_reduce(parallel, &
1090  CALL bmw_parallel_context_reduce(parallel, &
1092  CALL bmw_parallel_context_reduce(parallel, &
1094  END IF
1095 
1096  CALL profiler_set_stop_time('primed_grid_construct_jv', &
1097  & start_time)
1098 
1099  END FUNCTION
1100 
1101 !-------------------------------------------------------------------------------
1112 !-------------------------------------------------------------------------------
1113  FUNCTION primed_grid_construct_siesta(num_v, siesta_file_name, &
1114  & parallel)
1115  USE stel_constants, ONLY: twopi, mu0
1116  USE read_wout_mod, ONLY: nfp, ns, rmnc, rmns, zmnc, zmns, isigng, &
1117  & mnmax, xm, xn, lasym
1118  USE siesta_file
1119 
1120  IMPLICIT NONE
1121 
1122 ! Declare Arguments
1124  INTEGER, INTENT(in) :: num_v
1125  CHARACTER (len=*) :: siesta_file_name
1126  TYPE (bmw_parallel_context_class), INTENT(in) :: parallel
1127 
1128 ! local variables
1129  REAL (rprec) :: start_time
1130  INTEGER :: i
1131  REAL (rprec) :: x
1132  REAL (rprec) :: ds
1133  INTEGER :: si
1134  INTEGER :: ui
1135  INTEGER :: vi
1136  REAL (rprec) :: r
1137  REAL (rprec) :: z
1138  REAL (rprec) :: rs
1139  REAL (rprec) :: zs
1140  REAL (rprec) :: ru
1141  REAL (rprec) :: zu
1142  REAL (rprec) :: rv
1143  REAL (rprec) :: zv
1144  REAL (rprec) :: js
1145  REAL (rprec) :: ju
1146  REAL (rprec) :: jv
1147  REAL (rprec) :: jr
1148  REAL (rprec) :: jp
1149  REAL (rprec) :: jz
1150  INTEGER :: m
1151  INTEGER :: n
1152  INTEGER :: ilow
1153  INTEGER :: ihigh
1154  REAL (rprec) :: wlow
1155  REAL (rprec) :: whigh
1156  REAL (rprec), DIMENSION(:), ALLOCATABLE :: cosv
1157  REAL (rprec), DIMENSION(:), ALLOCATABLE :: sinv
1158  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: cosmu
1159  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: sinmu
1160  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: cosnv
1161  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: sinnv
1162  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: cosmu_vmec
1163  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: sinmu_vmec
1164  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: cosnv_vmec
1165  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: sinnv_vmec
1166  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: cosmn
1167  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: sinmn
1168  REAL (rprec), DIMENSION(:), ALLOCATABLE :: cosmn_vmec
1169  REAL (rprec), DIMENSION(:), ALLOCATABLE :: sinmn_vmec
1170  REAL (rprec), DIMENSION(:), ALLOCATABLE :: amnint
1171  INTEGER :: num_v_use
1172  TYPE (siesta_file_class), POINTER :: siesta
1173 
1174 ! local parameters
1175  INTEGER, PARAMETER :: num_u = 101
1176  real (rprec), PARAMETER :: du = twopi/num_u
1177 
1178 ! Start of executable code
1179  start_time = profiler_get_start_time()
1180 
1181  siesta => siesta_file_construct(trim(siesta_file_name))
1182 
1184 
1185  IF (num_v .eq. 1) THEN
1186  num_v_use = 360
1187  ELSE
1188  num_v_use = num_v
1189  END IF
1190 
1191  ALLOCATE(primed_grid_construct_siesta%x(siesta%nrad,num_u, &
1192  & num_v_use))
1193  ALLOCATE(primed_grid_construct_siesta%y(siesta%nrad,num_u, &
1194  & num_v_use))
1195  ALLOCATE(primed_grid_construct_siesta%z(siesta%nrad,num_u, &
1196  & num_v_use))
1197 
1198  ALLOCATE(primed_grid_construct_siesta%j_x(siesta%nrad,num_u, &
1199  & num_v_use))
1200  ALLOCATE(primed_grid_construct_siesta%j_y(siesta%nrad,num_u, &
1201  & num_v_use))
1202  ALLOCATE(primed_grid_construct_siesta%j_z(siesta%nrad,num_u, &
1203  & num_v_use))
1204 
1205  ALLOCATE(cosv(num_v_use))
1206  ALLOCATE(sinv(num_v_use))
1207  ALLOCATE(cosmu_vmec(mnmax,num_u))
1208  ALLOCATE(sinmu_vmec(mnmax,num_u))
1209  ALLOCATE(cosnv_vmec(mnmax,num_v_use))
1210  ALLOCATE(sinnv_vmec(mnmax,num_v_use))
1211  ALLOCATE(cosmu(0:siesta%mpol,num_u))
1212  ALLOCATE(sinmu(0:siesta%mpol,num_u))
1213  ALLOCATE(cosnv(-siesta%ntor:siesta%ntor,num_v_use))
1214  ALLOCATE(sinnv(-siesta%ntor:siesta%ntor,num_v_use))
1215 
1216  ds = 1.0/(ns - 1.0)
1217  primed_grid_construct_siesta%dv = twopi/num_v_use
1218 
1219  primed_grid_construct_siesta%dvol = isigng*ds*du &
1220  & * primed_grid_construct_siesta%dv/(2.0*twopi)
1221 
1222 !$OMP PARALLEL
1223 !$OMP& DEFAULT(SHARED)
1224 !$OMP& PRIVATE(i, x, si, ui, vi, n, m, r, ru, rv, z, zu, zv, rs, zs, &
1225 !$OMP& js, ju, jv, jr, jp, jz, ilow, wlow, ihigh, whigh, &
1226 !$OMP& cosmn, sinmn, cosmn_vmec, sinmn_vmec, amnint)
1227 
1228 ! Multi process will do an all reduce so these arrays need to be initalized.
1229  IF (parallel%stride .gt. 1) THEN
1230 !$OMP WORKSHARE
1237  cosmu_vmec = 0.0
1238  cosnv_vmec = 0.0
1239  sinmu_vmec = 0.0
1240  sinnv_vmec = 0.0
1241  cosmu = 0.0
1242  cosnv = 0.0
1243  sinmu = 0.0
1244  sinnv = 0.0
1245  cosv = 0.0
1246  sinv = 0.0
1247 !$OMP END WORKSHARE
1248  END IF
1249 
1250 !$OMP DO
1251 !$OMP& SCHEDULE(STATIC)
1252  DO i = bmw_parallel_context_start(parallel, num_u), &
1253  & bmw_parallel_context_end(parallel, num_u)
1254  x = (i - 0.5)*du
1255  cosmu_vmec(:,i) = cos(xm*x)
1256  sinmu_vmec(:,i) = sin(xm*x)
1257 
1258  DO m = 0, siesta%mpol
1259  cosmu(m,i) = cos(m*x)
1260  sinmu(m,i) = sin(m*x)
1261  END DO
1262  END DO
1263 !$OMP END DO
1264 
1265 !$OMP DO
1266 !$OMP& SCHEDULE(STATIC)
1267  DO i = bmw_parallel_context_start(parallel, num_v_use), &
1268  & bmw_parallel_context_end(parallel, num_v_use)
1269  x = (i - 0.5)*primed_grid_construct_siesta%dv
1270  cosv(i) = cos(x)
1271  sinv(i) = sin(x)
1272  cosnv_vmec(:,i) = cos(xn*x)
1273  sinnv_vmec(:,i) = sin(xn*x)
1274  DO n = -siesta%ntor, siesta%ntor
1275  cosnv(n,i) = cos(n*x)
1276  sinnv(n,i) = sin(n*x)
1277  END DO
1278  END DO
1279 !$OMP END DO
1280 
1281  IF (parallel%stride .gt. 1) THEN
1282 !$OMP SINGLE
1283  CALL bmw_parallel_context_reduce(parallel, cosv)
1284  CALL bmw_parallel_context_reduce(parallel, sinv)
1285  CALL bmw_parallel_context_reduce(parallel, cosmu_vmec)
1286  CALL bmw_parallel_context_reduce(parallel, sinmu_vmec)
1287  CALL bmw_parallel_context_reduce(parallel, cosnv_vmec)
1288  CALL bmw_parallel_context_reduce(parallel, sinnv_vmec)
1289  CALL bmw_parallel_context_reduce(parallel, cosmu)
1290  CALL bmw_parallel_context_reduce(parallel, sinmu)
1291  CALL bmw_parallel_context_reduce(parallel, cosnv)
1292  CALL bmw_parallel_context_reduce(parallel, sinnv)
1293 !$OMP END SINGLE
1294  END IF
1295 
1296  ALLOCATE(cosmn(0:siesta%mpol,-siesta%ntor:siesta%ntor))
1297  ALLOCATE(sinmn(0:siesta%mpol,-siesta%ntor:siesta%ntor))
1298 
1299  ALLOCATE(cosmn_vmec(mnmax))
1300  ALLOCATE(sinmn_vmec(mnmax))
1301 
1302  ALLOCATE(amnint(mnmax))
1303 
1304 !$OMP DO
1305 !$OMP& SCHEDULE(STATIC)
1306  DO i = bmw_parallel_context_start(parallel, &
1307  & siesta%nrad*num_u*num_v_use), &
1308  & bmw_parallel_context_end(parallel, &
1309  & siesta%nrad*num_u*num_v_use)
1310  si = bmw_parallel_context_i(i, siesta%nrad)
1311  ui = bmw_parallel_context_j(i, siesta%nrad, num_u)
1312  vi = bmw_parallel_context_k(i, siesta%nrad, num_u)
1313 
1314  x = (si - 1.0)/(siesta%nrad - 1.0)
1315 
1316  CALL primed_grid_siesta_interpol(x, ilow, wlow, ihigh, whigh)
1317 
1318  cosmn_vmec = cosmu_vmec(:,ui)*cosnv_vmec(:,vi) &
1319  & + sinmu_vmec(:,ui)*sinnv_vmec(:,vi)
1320  sinmn_vmec = sinmu_vmec(:,ui)*cosnv_vmec(:,vi) &
1321  & - cosmu_vmec(:,ui)*sinnv_vmec(:,vi)
1322 
1323  DO n = -siesta%ntor, siesta%ntor
1324  DO m = 0, siesta%mpol
1325  cosmn(m,n) = cosmu(m, ui)*cosnv(n, vi) &
1326  & - sinmu(m, ui)*sinnv(n, vi)
1327  sinmn(m,n) = sinmu(m, ui)*cosnv(n, vi) &
1328  & + cosmu(m, ui)*sinnv(n, vi)
1329  END DO
1330  END DO
1331 
1332  amnint = wlow*rmnc(:,ilow) + whigh*rmnc(:,ihigh)
1333  r = sum(amnint*cosmn_vmec)
1334  ru = -sum(xm*amnint*sinmn_vmec)
1335  rv = sum(xn*amnint*sinmn_vmec)
1336 
1337  amnint = wlow*zmns(:,ilow) + whigh*zmns(:,ihigh)
1338  z = sum(amnint*sinmn_vmec)
1339  zu = sum(xm*amnint*cosmn_vmec)
1340  zv = -sum(xn*amnint*cosmn_vmec)
1341 
1342  amnint = 2.0*x*(rmnc(:,ihigh) - rmnc(:,ilow))*(ns - 1.0)
1343  rs = sum(amnint*cosmn_vmec)
1344 
1345  amnint = 2.0*x*(zmns(:,ihigh) - zmns(:,ilow))*(ns - 1.0)
1346  zs = sum(amnint*sinmn_vmec)
1347 
1348  js = sum(siesta%jksupsmnsf(:,:,si)*sinmn)
1349  ju = sum(siesta%jksupumncf(:,:,si)*cosmn)
1350  jv = sum(siesta%jksupvmncf(:,:,si)*cosmn)
1351  IF (btest(siesta%flags, siesta_lasym_flag)) THEN
1352  IF (lasym) THEN
1353  amnint = wlow*rmns(:,ilow) + whigh*rmns(:,ihigh)
1354  r = r + sum(amnint*sinmn_vmec)
1355  ru = ru + sum(xm*amnint*cosmn_vmec)
1356  rv = rv - sum(xn*amnint*cosmn_vmec)
1357 
1358  amnint = wlow*zmnc(:,ilow) + whigh*zmnc(:,ihigh)
1359  z = z + sum(amnint*cosmn_vmec)
1360  zu = zu - sum(xm*amnint*sinmn_vmec)
1361  zv = zv + sum(xn*amnint*sinmn_vmec)
1362 
1363  amnint = 2.0*x*(rmns(:,ihigh) - rmns(:,ilow))*(ns - 1.0)
1364  rs = rs + sum(amnint*sinmn_vmec)
1365 
1366  amnint = 2.0*x*(zmnc(:,ihigh) - zmnc(:,ilow))*(ns - 1.0)
1367  zs = zs + sum(amnint*cosmn_vmec)
1368  END IF
1369 
1370  js = js + sum(siesta%jksupsmncf(:,:,si)*cosmn)
1371  ju = ju + sum(siesta%jksupumnsf(:,:,si)*sinmn)
1372  jv = jv + sum(siesta%jksupvmnsf(:,:,si)*sinmn)
1373 
1374  END IF
1375 
1376  js = js/(siesta%b_factor*mu0)
1377  ju = ju/(siesta%b_factor*mu0)
1378  jv = jv/(siesta%b_factor*mu0)
1379 
1380  jr = js*rs + ju*ru + jv*rv
1381  jp = jv*r
1382  jz = js*zs + ju*zu + jv*zv
1383 
1384  IF (si .eq. 1 .or. si .eq. ns) THEN
1385  jr = jr/2.0
1386  jp = jp/2.0
1387  jz = jz/2.0
1388  END IF
1389 
1390  primed_grid_construct_siesta%j_x(si,ui,vi) = jr*cosv(vi) &
1391  & - jp*sinv(vi)
1392  primed_grid_construct_siesta%j_y(si,ui,vi) = jr*sinv(vi) &
1393  & + jp*cosv(vi)
1394  primed_grid_construct_siesta%j_z(si,ui,vi) = jz
1395 
1396  primed_grid_construct_siesta%x(si,ui,vi) = r*cosv(vi)
1397  primed_grid_construct_siesta%y(si,ui,vi) = r*sinv(vi)
1398  primed_grid_construct_siesta%z(si,ui,vi) = z
1399  END DO
1400 !$OMP END DO
1401 
1402  DEALLOCATE(cosmn)
1403  DEALLOCATE(sinmn)
1404 
1405  DEALLOCATE(cosmn_vmec)
1406  DEALLOCATE(sinmn_vmec)
1407 
1408  DEALLOCATE(amnint)
1409 !$OMP END PARALLEL
1410 
1411  DEALLOCATE(cosv)
1412  DEALLOCATE(sinv)
1413  DEALLOCATE(cosmu_vmec)
1414  DEALLOCATE(sinmu_vmec)
1415  DEALLOCATE(cosnv_vmec)
1416  DEALLOCATE(sinnv_vmec)
1417  DEALLOCATE(cosmu)
1418  DEALLOCATE(sinmu)
1419  DEALLOCATE(cosnv)
1420  DEALLOCATE(sinnv)
1421 
1422 ! Multi process did not fill out the entire array. Get the missing pieces from
1423 ! the other processes.
1424  IF (parallel%stride .gt. 1) THEN
1425  CALL bmw_parallel_context_reduce(parallel, &
1427  CALL bmw_parallel_context_reduce(parallel, &
1429  CALL bmw_parallel_context_reduce(parallel, &
1431  CALL bmw_parallel_context_reduce(parallel, &
1433  CALL bmw_parallel_context_reduce(parallel, &
1435  CALL bmw_parallel_context_reduce(parallel, &
1437  END IF
1438 
1439  CALL siesta_file_destruct(siesta)
1440 
1441  CALL profiler_set_stop_time('primed_grid_construct_siesta', &
1442  & start_time)
1443 
1444  END FUNCTION
1445 
1446 !*******************************************************************************
1447 ! DESTRUCTION SUBROUTINES
1448 !*******************************************************************************
1449 !-------------------------------------------------------------------------------
1455 !-------------------------------------------------------------------------------
1456  SUBROUTINE primed_grid_destruct(this)
1458  IMPLICIT NONE
1459 
1460 ! Declare Arguments
1461  TYPE (primed_grid_class), POINTER :: this
1462 
1463 ! Start of executable code
1464  IF (ASSOCIATED(this%x)) THEN
1465  DEALLOCATE(this%x)
1466  this%x => null()
1467  END IF
1468 
1469  IF (ASSOCIATED(this%y)) THEN
1470  DEALLOCATE(this%y)
1471  this%y => null()
1472  END IF
1473 
1474  IF (ASSOCIATED(this%z)) THEN
1475  DEALLOCATE(this%z)
1476  this%z => null()
1477  END IF
1478 
1479  IF (ASSOCIATED(this%j_x)) THEN
1480  DEALLOCATE(this%j_x)
1481  this%j_x => null()
1482  END IF
1483 
1484  IF (ASSOCIATED(this%j_y)) THEN
1485  DEALLOCATE(this%j_y)
1486  this%j_y => null()
1487  END IF
1488 
1489  IF (ASSOCIATED(this%j_z)) THEN
1490  DEALLOCATE(this%j_z)
1491  this%j_z => null()
1492  END IF
1493 
1494  DEALLOCATE(this)
1495 
1496  END SUBROUTINE
1497 
1498 !*******************************************************************************
1499 ! UTILITY SUBROUTINES
1500 !*******************************************************************************
1501  SUBROUTINE primed_grid_siesta_interpol(s, ilow, wlow, &
1502  & ihigh, whigh)
1503  USE read_wout_mod, ONLY: ns
1504 
1505  IMPLICIT NONE
1506 
1507 ! Declare Arguments
1508  REAL (rprec), INTENT(in) :: s
1509  INTEGER, INTENT(out) :: ilow
1510  REAL (rprec), INTENT(out) :: wlow
1511  INTEGER, INTENT(out) :: ihigh
1512  REAL (rprec), INTENT(out) :: whigh
1513 
1514 ! local variables
1515  REAL (rprec) :: start_time
1516  REAL (rprec) :: wlow_r
1517 
1518 ! Start of executable code
1519  start_time = profiler_get_start_time()
1520 
1521  wlow_r = s*s*(ns - 1.0) + 1.0
1522  ilow = floor(wlow_r)
1523  IF (ilow .eq. ns) THEN
1524  ilow = ns - 1
1525  END IF
1526  ihigh = ilow + 1
1527  wlow = -wlow_r + 1.0 + ilow
1528  whigh = 1.0 - wlow
1529 
1530  END SUBROUTINE
1531 
1532  END MODULE
siesta_file::siesta_file_construct
type(siesta_file_class) function, pointer siesta_file_construct(siesta_file_name)
Construct a siesta_file_class object.
Definition: siesta_file.f:73
primed_grid::primed_grid_class
Base class representing a primed grid. This is grid the volume integral will be summed over.
Definition: primed_grid.f:27
primed_grid::primed_grid_construct_jv
type(primed_grid_class) function, pointer primed_grid_construct_jv(num_v, parallel)
Construct a primed_grid_class object.
Definition: primed_grid.f:764
profiler
Defines functions for measuring an tabulating performance of function and subroutine calls....
Definition: profiler.f:13
bmw_state_flags::bmw_state_flags_jv
integer, parameter bmw_state_flags_jv
Bit position for the use curl jv response flag.
Definition: bmw_state_flags.f:33
bmw_parallel_context::bmw_parallel_context_k
pure integer function bmw_parallel_context_k(index, num_i, num_j)
Compute the k index of a flat array.
Definition: bmw_parallel_context.f:566
bmw_state_flags
Contains parameters defining the bit positions for flags that mark different options.
Definition: bmw_state_flags.f:11
primed_grid::primed_grid_construct_both
type(primed_grid_class) function, pointer primed_grid_construct_both(num_v, parallel)
Construct a primed_grid_class object.
Definition: primed_grid.f:122
bmw_parallel_context::bmw_parallel_context_start
pure integer function bmw_parallel_context_start(this, total)
Compute the start index of a flat array.
Definition: bmw_parallel_context.f:429
primed_grid::primed_grid_construct
type(primed_grid_class) function, pointer primed_grid_construct(num_v, flags, siesta_file, parallel, io_unit)
Construct a primed_grid_class object.
Definition: primed_grid.f:68
primed_grid::primed_grid_destruct
subroutine primed_grid_destruct(this)
Deconstruct a primed_grid_class object.
Definition: primed_grid.f:1457
bmw_parallel_context
Defines the base class of the type bmw_parallel_context_class. This contains the state variables need...
Definition: bmw_parallel_context.f:11
siesta_file
Defines the base class of the type siesta_file_class. This contains the output of a siesta equilibriu...
Definition: siesta_file.f:11
bmw_parallel_context::bmw_parallel_context_i
pure integer function bmw_parallel_context_i(index, num_i)
Compute the i index of a flat array.
Definition: bmw_parallel_context.f:515
siesta_file::siesta_lasym_flag
integer, parameter siesta_lasym_flag
Version number.
Definition: siesta_file.f:21
siesta_file::siesta_file_class
Base class representing a siesta output.
Definition: siesta_file.f:31
profiler::profiler_get_start_time
real(rprec) function profiler_get_start_time()
Gets the start time of profiled function.
Definition: profiler.f:194
bmw_state_flags::bmw_state_flags_ju
integer, parameter bmw_state_flags_ju
Bit position for the use curl ju response flag.
Definition: bmw_state_flags.f:31
bmw_parallel_context::bmw_parallel_context_j
pure integer function bmw_parallel_context_j(index, num_i, num_j)
Compute the j index of a flat array.
Definition: bmw_parallel_context.f:540
siesta_file::siesta_file_destruct
subroutine siesta_file_destruct(this)
Deconstruct a siesta_file_class object.
Definition: siesta_file.f:162
primed_grid::primed_grid_construct_siesta
type(primed_grid_class) function, pointer primed_grid_construct_siesta(num_v, siesta_file_name, parallel)
Construct a primed_grid_class object.
Definition: primed_grid.f:1115
bmw_parallel_context::bmw_parallel_context_reduce
Interface for the buffer reduction.
Definition: bmw_parallel_context.f:49
profiler::profiler_set_stop_time
subroutine profiler_set_stop_time(symbol_name, start_time)
Gets the end time of profiled function.
Definition: profiler.f:121
bmw_parallel_context::bmw_parallel_context_class
Base class representing a bmw parallel context. This contains all memory needed parameters needed to ...
Definition: bmw_parallel_context.f:26
bmw_parallel_context::bmw_parallel_context_end
pure integer function bmw_parallel_context_end(this, total)
Compute the end index of a flat array.
Definition: bmw_parallel_context.f:472
primed_grid::primed_grid_construct_ju
type(primed_grid_class) function, pointer primed_grid_construct_ju(num_v, parallel)
Construct a primed_grid_class object.
Definition: primed_grid.f:413
bmw_state_flags::bmw_state_flags_siesta
integer, parameter bmw_state_flags_siesta
Bit position for the use siesta instead of vmec.
Definition: bmw_state_flags.f:35
primed_grid
Defines the base class of the type primed_grid_class. This contains the state variables to define the...
Definition: primed_grid.f:11