V3FIT
siesta_equilibrium.f
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
2 ! The @header, @table_section, @table_subsection, @item and @end_table commands
3 ! are custom defined commands in Doxygen.in. They are defined under ALIASES.
4 ! For the page created here, the 80 column limit is exceeded. Arguments of
5 ! aliases are separated by ','. If you intended ',' to be a string you must use
6 ! an escaped comma '\,'.
7 !
61 !*******************************************************************************
64 !
65 ! Note separating the Doxygen comment block here so detailed decription is
66 ! found in the Module not the file.
67 !
72 !*******************************************************************************
73 
74  MODULE siesta_equilibrium
75  USE siesta_context
76  USE vmec_equilibrium
77  USE siesta_run
78 
79  IMPLICIT NONE
80 
81 !*******************************************************************************
82 ! siesta equilibrium module parameters
83 !*******************************************************************************
84 ! Define id's only for values that can change. All others shouldn't be needed
85 ! outside of the vmec interface. Parameter id's start at 14 since 0-13 are
86 ! reserved for non equilibrium model parameters. These numbers will need to be
87 ! updated if any new model parameters are added. Siesta uses the VMEC model
88 ! parameters but overwrites the auxilary parameters.
89 
90 ! Auxiliary Model Parameters
92  INTEGER, PARAMETER :: siesta_pp_ne_b_id = 59
94  INTEGER, PARAMETER :: siesta_pp_ne_as_id = 60
96  INTEGER, PARAMETER :: siesta_pp_ne_af_id = 61
98  INTEGER, PARAMETER :: siesta_pp_sxrem_b_id = 62
100  INTEGER, PARAMETER :: siesta_pp_sxrem_as_id = 63
102  INTEGER, PARAMETER :: siesta_pp_sxrem_af_id = 64
105  INTEGER, PARAMETER :: siesta_pp_te_b_id = 65
108  INTEGER, PARAMETER :: siesta_pp_te_as_id = 66
111  INTEGER, PARAMETER :: siesta_pp_te_af_id = 67
114  INTEGER, PARAMETER :: siesta_pp_ti_b_id = 68
117  INTEGER, PARAMETER :: siesta_pp_ti_as_id = 69
120  INTEGER, PARAMETER :: siesta_pp_ti_af_id = 70
123  INTEGER, PARAMETER :: siesta_pp_ze_b_id = 71
126  INTEGER, PARAMETER :: siesta_pp_ze_as_id = 72
129  INTEGER, PARAMETER :: siesta_pp_ze_af_id = 73
130 
131 ! Skip 74-82 since these parameters have a VMEC value but not a SIESTA analog.
134  INTEGER, PARAMETER :: siesta_phi_offset_id = 83
137  INTEGER, PARAMETER :: siesta_z_offset_id = 84
138 
140  INTEGER, PARAMETER :: siesta_helpert_id = 85
141 
142 !*******************************************************************************
143 ! DERIVED-TYPE DECLARATIONS
144 ! 1) siesta magnetic cache
145 ! 2) siesta base class
146 !
147 !*******************************************************************************
148 !-------------------------------------------------------------------------------
150 !-------------------------------------------------------------------------------
151  TYPE siesta_magnetic_cache
153  REAL (rprec) :: ds
155  REAL (rprec) :: du
157  REAL (rprec) :: dv
158 
160  REAL (rprec) :: du_a
162  REAL (rprec) :: dv_a
163 
165  REAL (rprec) :: du_full
167  REAL (rprec) :: dv_full
168 
170  REAL (rprec), DIMENSION(:,:,:), POINTER :: rsuv => null()
172  REAL (rprec), DIMENSION(:,:,:), POINTER :: zsuv => null()
173 
175  REAL (rprec), DIMENSION(:,:,:), POINTER :: jrsuv => null()
177  REAL (rprec), DIMENSION(:,:,:), POINTER :: jphisuv => null()
179  REAL (rprec), DIMENSION(:,:,:), POINTER :: jzsuv => null()
180 
183  REAL (rprec), DIMENSION(:,:), POINTER :: kruv => null()
186  REAL (rprec), DIMENSION(:,:), POINTER :: kphiuv => null()
189  REAL (rprec), DIMENSION(:,:), POINTER :: kzuv => null()
190 
192  REAL (rprec), DIMENSION(:,:,:), POINTER :: x_prime => null()
195  REAL (rprec), DIMENSION(:,:), POINTER :: kxuv_full => null()
198  REAL (rprec), DIMENSION(:,:), POINTER :: kyuv_full => null()
201  REAL (rprec), DIMENSION(:,:), POINTER :: kzuv_full => null()
202 
204  REAL (rprec), DIMENSION(:,:,:), POINTER :: x_axi => null()
207  REAL (rprec), DIMENSION(:,:), POINTER :: kxuv_axi => null()
210  REAL (rprec), DIMENSION(:,:), POINTER :: kyuv_axi => null()
213  REAL (rprec), DIMENSION(:,:), POINTER :: kzuv_axi => null()
214  END TYPE
215 
216 !-------------------------------------------------------------------------------
220 !-------------------------------------------------------------------------------
221  TYPE siesta_class
223  CHARACTER (len=path_length) :: restart_file_name
225  CHARACTER (len=path_length) :: siesta_file_name
226 
228  TYPE (pprofile_class), POINTER :: ne => null()
231  TYPE (pprofile_class), POINTER :: te => null()
234  TYPE (pprofile_class), POINTER :: ti => null()
237  TYPE (pprofile_pointer), DIMENSION(:), POINTER :: &
238  & sxrem => null()
239 
241  REAL (rprec) :: pol_rad_ratio
243  TYPE (siesta_magnetic_cache), POINTER :: &
244  & magnetic_cache => null()
245 
247  TYPE (siesta_context_class), POINTER :: context => null()
248 
249 ! Extra reconstruction parameters
251  REAL (rprec) :: phi_offset = 0
253  REAL (rprec) :: z_offset = 0
254 
256  TYPE (vmec_class), POINTER :: vmec => null()
257  END TYPE
258 
259 !*******************************************************************************
260 ! INTERFACE BLOCKS
261 !*******************************************************************************
262 !-------------------------------------------------------------------------------
266 !-------------------------------------------------------------------------------
267  INTERFACE siesta_set_magnetic_cache
268  MODULE PROCEDURE siesta_set_magnetic_cache_responce, &
269  & siesta_set_magnetic_cache_point, &
270  & siesta_set_magnetic_cache_calc
271  END INTERFACE
272 
273 !-------------------------------------------------------------------------------
275 !-------------------------------------------------------------------------------
276  INTERFACE siesta_get_ne
277  MODULE PROCEDURE siesta_get_ne_cart, &
278  & siesta_get_ne_radial
279  END INTERFACE
280 
281 !-------------------------------------------------------------------------------
283 !-------------------------------------------------------------------------------
284  INTERFACE siesta_get_gp_ne
285  MODULE PROCEDURE siesta_get_gp_ne_ij, &
286  & siesta_get_gp_ne_pi, &
287  & siesta_get_gp_ne_pp
288  END INTERFACE
289 
290 !-------------------------------------------------------------------------------
292 !-------------------------------------------------------------------------------
293  INTERFACE siesta_get_te
294  MODULE PROCEDURE siesta_get_te_cart, &
295  & siesta_get_te_radial
296  END INTERFACE
297 
298 !-------------------------------------------------------------------------------
301 !-------------------------------------------------------------------------------
302  INTERFACE siesta_get_gp_te
303  MODULE PROCEDURE siesta_get_gp_te_ij, &
304  & siesta_get_gp_te_pi, &
305  & siesta_get_gp_te_pp
306  END INTERFACE
307 
308 
309 !-------------------------------------------------------------------------------
311 !-------------------------------------------------------------------------------
312  INTERFACE siesta_get_ti
313  MODULE PROCEDURE siesta_get_ti_cart, &
314  & siesta_get_ti_radial
315  END INTERFACE
316 
317 !-------------------------------------------------------------------------------
319 !-------------------------------------------------------------------------------
320  INTERFACE siesta_get_gp_ti
321  MODULE PROCEDURE siesta_get_gp_ti_ij, &
322  & siesta_get_gp_ti_pi, &
323  & siesta_get_gp_ti_pp
324  END INTERFACE
325 
326 !-------------------------------------------------------------------------------
328 !-------------------------------------------------------------------------------
329  INTERFACE siesta_get_sxrem
330  MODULE PROCEDURE siesta_get_sxrem_cart, &
331  & siesta_get_sxrem_radial
332  END INTERFACE
333 
334 !-------------------------------------------------------------------------------
337 !-------------------------------------------------------------------------------
338  INTERFACE siesta_get_gp_sxrem
339  MODULE PROCEDURE siesta_get_gp_sxrem_ij, &
340  & siesta_get_gp_sxrem_pi, &
341  & siesta_get_gp_sxrem_pp
342  END INTERFACE
343 
344 !-------------------------------------------------------------------------------
346 !-------------------------------------------------------------------------------
347  INTERFACE siesta_get_p
348  MODULE PROCEDURE siesta_get_p_cart, &
349  & siesta_get_p_radial
350  END INTERFACE
351 
352  CONTAINS
353 
354 !*******************************************************************************
355 ! CONSTRUCTION SUBROUTINES
356 !*******************************************************************************
357 !-------------------------------------------------------------------------------
384 !-------------------------------------------------------------------------------
385  FUNCTION siesta_construct(file_name, restart_file_name, &
386  & ne, te, ti, sxrem, phi_offset, z_offset, &
387  & pol_rad_ratio, iou, eq_comm, recon_comm, &
388  & vmec, state_flags)
389  USE safe_open_mod
390  USE model_state
391  USE file_opts
392 
393  IMPLICIT NONE
394 
395 ! Declare Arguments
396  TYPE (siesta_class), POINTER :: siesta_construct
397  CHARACTER (len=*), INTENT(in) :: file_name
398  CHARACTER (len=*), INTENT(in) :: restart_file_name
399  TYPE (pprofile_class), POINTER :: ne
400  TYPE (pprofile_class), POINTER :: te
401  TYPE (pprofile_class), POINTER :: ti
402  TYPE (pprofile_pointer), DIMENSION(:), POINTER :: sxrem
403  REAL (rprec), INTENT(in) :: phi_offset
404  REAL (rprec), INTENT(in) :: z_offset
405  REAL (rprec), INTENT(in) :: pol_rad_ratio
406  INTEGER, INTENT(in) :: iou
407  INTEGER, INTENT(in) :: eq_comm
408  INTEGER, INTENT(in) :: recon_comm
409  TYPE (vmec_class), POINTER :: vmec
410  INTEGER, INTENT(inout) :: state_flags
411 
412 ! local variables
413  INTEGER :: eq_rank
414  INTEGER :: recon_rank
415  INTEGER :: recon_size
416  INTEGER :: error
417  INTEGER :: iou_nli
418  REAL (rprec) :: start_time
419 
420 ! Start of executable code
421  start_time = profiler_get_start_time()
422 
423  ALLOCATE(siesta_construct)
424 
425  eq_rank = 0
426  recon_rank = 0
427  recon_size = 1
428 #if defined(MPI_OPT)
429  CALL mpi_comm_rank(eq_comm, eq_rank, error)
430 
431  IF (eq_rank .eq. 0) THEN
432  CALL mpi_comm_rank(recon_comm, recon_rank, error)
433  CALL mpi_comm_size(recon_comm, recon_size, error)
434  END IF
435  CALL mpi_bcast(recon_rank, 1, mpi_integer, 0, eq_comm, error)
436 #endif
437 
438  siesta_construct%vmec => vmec
439 
440  IF (eq_rank .eq. 0) THEN
441 ! Save a copy of the jcf orginal jcf file.
442  CALL copy_file(trim(file_name), trim(file_name) // '_save', &
443  & error)
444 
445  WRITE (*,*) ' *** Initializing SIESTA equilibrium from ' // &
446  & 'file ' // trim(file_name)
447  WRITE (iou,*) ' *** Initializing SIESTA equilibrium from ' // &
448  & 'file ' // trim(file_name)
449 
450  siesta_construct%ne => ne
451  siesta_construct%te => te
452  siesta_construct%ti => ti
453  siesta_construct%sxrem => sxrem
454 
455  siesta_construct%phi_offset = phi_offset
456  siesta_construct%z_offset = z_offset
457 
458  siesta_construct%pol_rad_ratio = pol_rad_ratio
459  siesta_construct%siesta_file_name = trim(file_name)
460 
461  iou_nli = 100
462  CALL safe_open(iou_nli, error, trim(file_name), 'old', &
463  & 'formatted')
464  READ (iou_nli, nml=siesta_info)
465  CLOSE (iou_nli, iostat=error)
466  wout_file = vmec%wout_file_name
467  WRITE (siesta_construct%restart_file_name,1000) &
468  & trim(restart_ext)
469 
470  state_flags = ibset(state_flags, model_state_siesta_flag)
471 
472  IF (restart_file_name .ne. '') THEN
473  state_flags = ibclr(state_flags, model_state_siesta_flag)
474  siesta_construct%context => &
475  & siesta_context_construct(restart_file_name)
476  END IF
477  END IF
478 
479  CALL profiler_set_stop_time('siesta_construct', start_time)
480 
481 1000 FORMAT('siesta_',a,'.nc')
482 
483  END FUNCTION
484 
485 !*******************************************************************************
486 ! DESTRUCTION SUBROUTINES
487 !*******************************************************************************
488 !-------------------------------------------------------------------------------
494 !-------------------------------------------------------------------------------
495  SUBROUTINE siesta_destruct(this)
496  USE read_wout_mod, only: read_wout_deallocate
497  USE file_opts
498 
499  IMPLICIT NONE
500 
501 ! Declare Arguments
502  TYPE (siesta_class), POINTER :: this
503 
504 ! local variables
505  INTEGER :: i, error
506 
507 ! Start of executable code
508  IF (ASSOCIATED(this%ne)) THEN
509  CALL pprofile_destruct(this%ne)
510  this%ne => null()
511  END IF
512 
513  IF (ASSOCIATED(this%te)) THEN
514  CALL pprofile_destruct(this%te)
515  this%te => null()
516  END IF
517 
518  IF (ASSOCIATED(this%ti)) THEN
519  CALL pprofile_destruct(this%ti)
520  this%ti => null()
521  END IF
522 
523  IF (ASSOCIATED(this%sxrem)) THEN
524  DO i = 1, SIZE(this%sxrem)
525  IF (ASSOCIATED(this%sxrem(i)%p)) THEN
526  CALL pprofile_destruct(this%sxrem(i)%p)
527  this%sxrem(i)%p => null()
528  END IF
529  END DO
530  DEALLOCATE(this%sxrem)
531  this%sxrem => null()
532  END IF
533 
534  IF (ASSOCIATED(this%magnetic_cache)) THEN
535  IF (ASSOCIATED(this%magnetic_cache%rsuv)) THEN
536  DEALLOCATE(this%magnetic_cache%rsuv)
537  this%magnetic_cache%rsuv => null()
538  END IF
539 
540  IF (ASSOCIATED(this%magnetic_cache%zsuv)) THEN
541  DEALLOCATE(this%magnetic_cache%zsuv)
542  this%magnetic_cache%zsuv => null()
543  END IF
544 
545  IF (ASSOCIATED(this%magnetic_cache%jrsuv)) THEN
546  DEALLOCATE(this%magnetic_cache%jrsuv)
547  this%magnetic_cache%jrsuv => null()
548  END IF
549 
550  IF (ASSOCIATED(this%magnetic_cache%jphisuv)) THEN
551  DEALLOCATE(this%magnetic_cache%jphisuv)
552  this%magnetic_cache%jphisuv => null()
553  END IF
554 
555  IF (ASSOCIATED(this%magnetic_cache%jzsuv)) THEN
556  DEALLOCATE(this%magnetic_cache%jzsuv)
557  this%magnetic_cache%jzsuv => null()
558  END IF
559 
560  IF (ASSOCIATED(this%magnetic_cache%kruv)) THEN
561  DEALLOCATE(this%magnetic_cache%kruv)
562  this%magnetic_cache%kruv => null()
563  END IF
564 
565  IF (ASSOCIATED(this%magnetic_cache%kphiuv)) THEN
566  DEALLOCATE(this%magnetic_cache%kphiuv)
567  this%magnetic_cache%kphiuv => null()
568  END IF
569 
570  IF (ASSOCIATED(this%magnetic_cache%kzuv)) THEN
571  DEALLOCATE(this%magnetic_cache%kzuv)
572  this%magnetic_cache%kzuv => null()
573  END IF
574 
575  IF (ASSOCIATED(this%magnetic_cache%x_prime)) THEN
576  DEALLOCATE(this%magnetic_cache%x_prime)
577  this%magnetic_cache%x_prime => null()
578  END IF
579 
580  IF (ASSOCIATED(this%magnetic_cache%kxuv_full)) THEN
581  DEALLOCATE(this%magnetic_cache%kxuv_full)
582  this%magnetic_cache%kxuv_full => null()
583  END IF
584 
585  IF (ASSOCIATED(this%magnetic_cache%kyuv_full)) THEN
586  DEALLOCATE(this%magnetic_cache%kyuv_full)
587  this%magnetic_cache%kyuv_full => null()
588  END IF
589 
590  IF (ASSOCIATED(this%magnetic_cache%kzuv_full)) THEN
591  DEALLOCATE(this%magnetic_cache%kzuv_full)
592  this%magnetic_cache%kzuv_full => null()
593  END IF
594 
595  IF (ASSOCIATED(this%magnetic_cache%x_axi)) THEN
596  DEALLOCATE(this%magnetic_cache%x_axi)
597  this%magnetic_cache%x_axi => null()
598  END IF
599 
600  IF (ASSOCIATED(this%magnetic_cache%kxuv_axi)) THEN
601  DEALLOCATE(this%magnetic_cache%kxuv_axi)
602  this%magnetic_cache%kxuv_axi => null()
603  END IF
604 
605  IF (ASSOCIATED(this%magnetic_cache%kyuv_axi)) THEN
606  DEALLOCATE(this%magnetic_cache%kyuv_axi)
607  this%magnetic_cache%kyuv_axi => null()
608  END IF
609 
610  IF (ASSOCIATED(this%magnetic_cache%kzuv_axi)) THEN
611  DEALLOCATE(this%magnetic_cache%kzuv_axi)
612  this%magnetic_cache%kzuv_axi => null()
613  END IF
614 
615  DEALLOCATE(this%magnetic_cache)
616  this%magnetic_cache => null()
617  END IF
618 
619 ! Delete the restart file. Errors here can safely be ignored.
620  CALL delete_file(trim(this%restart_file_name) // '_cache', error)
621 
622  IF (ASSOCIATED(this%context)) THEN
623  CALL siesta_context_destruct(this%context)
624  this%context => null()
625  END IF
626 
627  IF (ASSOCIATED(this%vmec)) THEN
628  CALL vmec_destruct(this%vmec)
629  this%vmec => null()
630  END IF
631 
632  DEALLOCATE(this)
633 
634  END SUBROUTINE
635 
636 !*******************************************************************************
637 ! SETTER SUBROUTINES
638 !*******************************************************************************
639 !-------------------------------------------------------------------------------
655 !-------------------------------------------------------------------------------
656  SUBROUTINE siesta_set_param(this, id, i_index, j_index, value, &
657  & eq_comm, state_flags)
658  USE model_state
659 
660  IMPLICIT NONE
661 
662 ! Declare Arguments
663  TYPE (siesta_class), INTENT(inout) :: this
664  INTEGER, INTENT(in) :: id
665  INTEGER, INTENT(in) :: i_index
666  INTEGER, INTENT(in) :: j_index
667  REAL (rprec), INTENT(in) :: value
668  INTEGER, INTENT(in) :: eq_comm
669  INTEGER, INTENT(inout) :: state_flags
670 
671 ! local variables
672  INTEGER :: error
673  REAL (rprec) :: start_time
674 
675 ! Start of executable code
676  start_time = profiler_get_start_time()
677 
678  SELECT CASE (id)
679 
680  CASE (siesta_pp_ne_b_id)
681  state_flags = ibset(state_flags, model_state_ne_flag)
682  this%ne%b(i_index) = value
683 
684  CASE (siesta_pp_ne_as_id)
685  state_flags = ibset(state_flags, model_state_ne_flag)
686  this%ne%as(i_index) = value
687 
688  CASE (siesta_pp_ne_af_id)
689  state_flags = ibset(state_flags, model_state_ne_flag)
690  this%ne%af(i_index) = value
691 
692  CASE (siesta_pp_te_b_id)
693  state_flags = ibset(state_flags, model_state_te_flag)
694  this%te%b(i_index) = value
695 
696  CASE (siesta_pp_te_as_id)
697  state_flags = ibset(state_flags, model_state_te_flag)
698  this%te%as(i_index) = value
699 
700  CASE (siesta_pp_te_af_id)
701  state_flags = ibset(state_flags, model_state_te_flag)
702  this%te%af(i_index) = value
703 
704  CASE (siesta_pp_ti_b_id)
705  state_flags = ibset(state_flags, model_state_ti_flag)
706  this%ti%b(i_index) = value
707 
708  CASE (siesta_pp_ti_as_id)
709  state_flags = ibset(state_flags, model_state_ti_flag)
710  this%ti%as(i_index) = value
711 
712  CASE (siesta_pp_ti_af_id)
713  state_flags = ibset(state_flags, model_state_ti_flag)
714  this%ti%af(i_index) = value
715 
716 ! There are multiple soft x-ray emission profiles. These need to be offset by
717 ! the array index.
718  CASE (siesta_pp_sxrem_b_id)
719  state_flags = ibset(state_flags, model_state_sxrem_flag + &
720  & (i_index - 1))
721  this%sxrem(i_index)%p%b(j_index) = value
722 
723  CASE (siesta_pp_sxrem_as_id)
724  state_flags = ibset(state_flags, model_state_sxrem_flag + &
725  & (i_index - 1))
726  this%sxrem(i_index)%p%as(j_index) = value
727 
728  CASE (siesta_pp_sxrem_af_id)
729  state_flags = ibset(state_flags, model_state_sxrem_flag + &
730  & (i_index - 1))
731  this%sxrem(i_index)%p%af(j_index) = value
732 
733 ! For the phi_offset, the magnetic cache needs to be updated. Since the phi
734 ! offset doesn't alter the equilibrium, the magnetics were not getting updated
735 ! since the equilibrium was not reconverged. Call vmec_set_magnetic_cache_calc
736 ! to update the magnetic cache without reconverging VMEC.
737  CASE (siesta_phi_offset_id)
738  state_flags = ibset(state_flags, model_state_shift_flag)
739  this%phi_offset = value
740  IF (ASSOCIATED(this%magnetic_cache)) THEN
741  CALL siesta_set_magnetic_cache(this)
742  END IF
743 
744  CASE (siesta_z_offset_id)
745  state_flags = ibset(state_flags, model_state_shift_flag)
746  this%z_offset = value
747  IF (ASSOCIATED(this%magnetic_cache)) THEN
748  CALL siesta_set_magnetic_cache(this)
749  END IF
750 
751  CASE (siesta_helpert_id)
752  state_flags = ibset(state_flags, model_state_siesta_flag)
753  helpert(i_index) = value
754 
755  CASE DEFAULT
756  CALL vmec_set_param(this%vmec, id, i_index, j_index, value, &
757  & eq_comm, state_flags)
758 
759  END SELECT
760 
761  CALL profiler_set_stop_time('siesta_set_param', start_time)
762 
763  END SUBROUTINE
764 
765 !-------------------------------------------------------------------------------
775 !-------------------------------------------------------------------------------
776  SUBROUTINE siesta_set_magnetic_cache_responce(this, &
777  & response_object)
778  USE stel_constants, only: twopi
779  USE magnetic_response
780 
781  IMPLICIT NONE
782 
783 ! Declare Arguments
784  TYPE (siesta_class), INTENT(inout) :: this
785  TYPE (magnetic_response_class), INTENT(in) :: response_object
786 
787 ! local variables
788  INTEGER :: numU
789  INTEGER :: ns
790  INTEGER :: mpol
791  REAL (rprec) :: start_time
792 
793 ! Start of executable code.
794  start_time = profiler_get_start_time()
795 
796  IF (.not.ASSOCIATED(this%magnetic_cache)) THEN
797  ALLOCATE(this%magnetic_cache)
798  END IF
799 
800  IF (magnetic_response_use_plasma(response_object)) THEN
801 ! In order to sample the fields on a grid, the sampling frequency of the grid
802 ! must be at least twice the frequency of the highest mode.
803  IF (ASSOCIATED(this%context)) THEN
804  ns = this%context%ns
805  mpol = this%context%mpol
806  ELSE
807  ns = nsin
808  mpol = mpolin
809  END IF
810  numu = max(int(this%pol_rad_ratio*ns), 2*mpol)
811 
812 ! Set the differental volume elements.
813  this%magnetic_cache%ds = 1.0/(ns - 1.0)
814  this%magnetic_cache%du = twopi/numu
815  this%magnetic_cache%dv = twopi/response_object%n_field_periods &
816  & / response_object%num_t
817 
818 ! When stellarator symmetry is used, do not need to store all the toroidal
819 ! planes in a single field period.
820  ALLOCATE(this%magnetic_cache%rsuv(ns,numu, &
821  & SIZE(response_object%a_r)))
822  ALLOCATE(this%magnetic_cache%zsuv(ns,numu, &
823  & SIZE(response_object%a_r)))
824  ALLOCATE(this%magnetic_cache%jrsuv(ns,numu, &
825  & SIZE(response_object%a_r)))
826  ALLOCATE(this%magnetic_cache%jphisuv(ns,numu, &
827  & SIZE(response_object%a_r)))
828  ALLOCATE(this%magnetic_cache%jzsuv(ns,numu, &
829  & SIZE(response_object%a_r)))
830  END IF
831 
832 ! When conducting shell is used allocate the surface vector arrays.
833  IF (magnetic_response_use_shell(response_object)) THEN
834 ! Set the differental area elements.
835  this%magnetic_cache%du_a = twopi/compression_get_dimension1( &
836  & response_object%a_s_r)
837  this%magnetic_cache%dv_a &
838  & = (twopi/response_object%n_field_periods) &
839  & / response_object%num_t_shell
840 
841  ALLOCATE(this%magnetic_cache%kruv( &
842  & compression_get_dimension1(response_object%a_s_r), &
843  & compression_get_dimension2(response_object%a_s_r)))
844  ALLOCATE(this%magnetic_cache%kphiuv( &
845  & compression_get_dimension1(response_object%a_s_r), &
846  & compression_get_dimension2(response_object%a_s_r)))
847  ALLOCATE(this%magnetic_cache%kzuv( &
848  & compression_get_dimension1(response_object%a_s_r), &
849  & compression_get_dimension2(response_object%a_s_r)))
850  END IF
851 
852  CALL profiler_set_stop_time('siesta_set_magnetic_cache_responce', &
853  & start_time)
854 
855  END SUBROUTINE
856 
857 !-------------------------------------------------------------------------------
866 !-------------------------------------------------------------------------------
867  SUBROUTINE siesta_set_magnetic_cache_point(this, use_axi)
868 
869  IMPLICIT NONE
870 
871 ! Declare Arguments
872  TYPE (siesta_class), INTENT(inout) :: this
873  LOGICAL, INTENT(in) :: use_axi
874 
875 ! local variables
876  INTEGER :: u_size, v_size
877  REAL (rprec) :: rbc00, rbc01, zbs01
878  REAL (rprec) :: start_time
879 
880 ! Start of executable code
881  start_time = profiler_get_start_time()
882 
883  IF (.not.ASSOCIATED(this%magnetic_cache)) THEN
884  ALLOCATE(this%magnetic_cache)
885  END IF
886 
887  rbc00 = vmec_get_param_value(this%vmec, vmec_rbc_id, 0, 0)
888  rbc01 = vmec_get_param_value(this%vmec, vmec_rbc_id, 0, 1)
889  zbs01 = vmec_get_param_value(this%vmec, vmec_zbs_id, 0, 1)
890 
891 ! Set the grid size based on the size of the rough size of the plasma.
892  v_size = max(int(twopi*rbc00/magnetic_cache_vc_grid_dim), &
893  & magnetic_cache_vc_min_grid_points)
894  IF (mod(v_size, 2) .eq. 0) THEN
895  v_size = v_size + 1
896  END IF
897  this%magnetic_cache%dv_full = twopi/v_size
898 
899  u_size = max(int(twopi*max(rbc01, zbs01)/ &
900  & magnetic_cache_vc_grid_dim), &
901  & magnetic_cache_vc_min_grid_points)
902  IF (mod(u_size, 2) .eq. 0) THEN
903  u_size = u_size + 1
904  END IF
905  this%magnetic_cache%du_full = twopi/u_size
906 
907  ALLOCATE(this%magnetic_cache%kxuv_full(u_size,v_size))
908  ALLOCATE(this%magnetic_cache%kyuv_full(u_size,v_size))
909  ALLOCATE(this%magnetic_cache%kzuv_full(u_size,v_size))
910 
911  ALLOCATE(this%magnetic_cache%x_prime(u_size,v_size,3))
912 
913  IF (use_axi) THEN
914  ALLOCATE(this%magnetic_cache%kxuv_axi(u_size,v_size))
915  ALLOCATE(this%magnetic_cache%kyuv_axi(u_size,v_size))
916  ALLOCATE(this%magnetic_cache%kzuv_axi(u_size,v_size))
917 
918  ALLOCATE(this%magnetic_cache%x_axi(u_size,v_size,3))
919 
920  END IF
921 
922  CALL profiler_set_stop_time('siesta_set_magnetic_cache_point', &
923  & start_time)
924 
925  END SUBROUTINE
926 
927 !-------------------------------------------------------------------------------
934 !-------------------------------------------------------------------------------
935  SUBROUTINE siesta_set_magnetic_cache_calc(this)
936  USE stel_constants, only: mu0
937  USE read_wout_mod, only: nfp, mnmax, xm, xn, rmnc, rmns, zmnc, &
938  & zmns, lasym, nsvmec => ns
939 
940  IMPLICIT NONE
941 
942 ! Declare Arguments
943  TYPE (siesta_class), INTENT(inout) :: this
944 
945 ! local variables
946  INTEGER :: i, j, k
947  REAL (rprec) :: s, u, v
948  INTEGER :: numS, numU, numV
949  INTEGER :: m, n
950  INTEGER :: mpol, ntor
951  REAL (rprec) :: r, rs, ru, rv
952  REAL (rprec) :: z, zs, zu, zv
953  REAL (rprec) :: js, ju, jv
954  REAL (rprec) :: bs, bu, bv
955  REAL (rprec) :: kr, kp
956  REAL (rprec) :: theta, phi
957  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: bsubuc, bsubus
958  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: bsubvc, bsubvs
959  REAL (rprec), DIMENSION(:), ALLOCATABLE :: rmncint, rmnsint
960  REAL (rprec), DIMENSION(:), ALLOCATABLE :: zmncint, zmnsint
961  REAL (rprec) :: cosphi
962  REAL (rprec) :: sinphi
963  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: cosmn
964  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: sinmn
965  REAL (rprec), DIMENSION(:), ALLOCATABLE :: cosmn_vmec
966  REAL (rprec), DIMENSION(:), ALLOCATABLE :: sinmn_vmec
967  REAL (rprec) :: wlow, whigh
968  INTEGER :: ilow, ihigh
969  REAL (rprec) :: start_time
970 
971 ! Start of executable code.
972  start_time = profiler_get_start_time()
973 
974  IF (ASSOCIATED(this%magnetic_cache%rsuv)) THEN
975 ! Get the size of the SUV grid.
976  nums = SIZE(this%magnetic_cache%rsuv, 1)
977  numu = SIZE(this%magnetic_cache%rsuv, 2)
978  numv = SIZE(this%magnetic_cache%rsuv, 3)
979 
980  IF (ASSOCIATED(this%context)) THEN
981  mpol = this%context%mpol
982  ntor = this%context%ntor
983  ELSE
984  mpol = mpolin
985  ntor = ntorin
986  END IF
987 
988 !$OMP PARALLEL
989 !$OMP& DEFAULT(SHARED)
990 !$OMP& PRIVATE(i,j,k,s,u,v,cosmn,sinmn,cosmn_vmec,sinmn_vmec,m,n,r,z,
991 !$OMP& rs,zs,ru,zu,rv,zv,js,ju,jv,rmncint,rmnsint,zmncint,
992 !$OMP& zmnsint,ilow,ihigh,wlow,whigh)
993 
994  ALLOCATE(cosmn(0:mpol,-ntor:ntor))
995  ALLOCATE(sinmn(0:mpol,-ntor:ntor))
996 
997  ALLOCATE(cosmn_vmec(mnmax))
998  ALLOCATE(sinmn_vmec(mnmax))
999 
1000  ALLOCATE(rmncint(mnmax))
1001  ALLOCATE(zmnsint(mnmax))
1002  IF (lasym) THEN
1003  ALLOCATE(rmnsint(mnmax))
1004  ALLOCATE(zmncint(mnmax))
1005  END IF
1006 
1007 ! Set the suv grids
1008 !
1009 ! ru = dR/du
1010 ! rv = dR/dv Note: dRHat/dv = PhiHat
1011 !
1012 ! zu = dZ/du
1013 ! zu = dZ/dv
1014 !
1015 ! jr = ju*ru + jv*rv
1016 ! jphi = jv*R
1017 ! jz = ju*zu + jv*zv
1018 
1019 !$OMP DO
1020 !$OMP& SCHEDULE(STATIC)
1021  DO k = 1, numv
1022  v = (k - 1)*this%magnetic_cache%dv + this%phi_offset
1023 
1024  DO j = 1, numu
1025  u = (j - 1)*this%magnetic_cache%du
1026 
1027  cosmn_vmec = cos(xm*u - xn*v)
1028  sinmn_vmec = sin(xm*u - xn*v)
1029 
1030  DO n = -ntor, ntor
1031  DO m = 0, mpol
1032  cosmn(m,n) = cos(m*u + nfp*n*v)
1033  sinmn(m,n) = sin(m*u + nfp*n*v)
1034  END DO
1035  END DO
1036 
1037  DO i = 1, nums
1038  s = (i - 1.0)/(nums - 1.0)
1039  ilow = siesta_to_i_vmec_low(s)
1040  IF (ilow .eq. nsvmec) THEN
1041  ilow = nsvmec - 1
1042  END IF
1043  ihigh = ilow + 1
1044  wlow = siesta_to_w_vmec_low(s, ilow)
1045  whigh = 1.0 - wlow
1046 
1047  rmncint = wlow*rmnc(:,ilow) + whigh*rmnc(:,ihigh)
1048  zmnsint = wlow*zmns(:,ilow) + whigh*zmns(:,ihigh)
1049 
1050  this%magnetic_cache%rsuv(i,j,k) = &
1051  & sum(rmncint*cosmn_vmec)
1052  this%magnetic_cache%zsuv(i,j,k) = &
1053  & sum(zmnsint*sinmn_vmec)
1054 
1055  ru = -sum(xm*rmncint*sinmn_vmec)
1056  rv = sum(xn*rmncint*sinmn_vmec)
1057  zu = sum(xm*zmnsint*cosmn_vmec)
1058  zv = -sum(xn*zmnsint*cosmn_vmec)
1059 
1060 ! dR(s_v)/ds_s = dR(s_s)/ds_v*ds_v/ds_s = 2sdR(s_v)/ds_v
1061  rmncint = 2.0*s*(rmnc(:,ihigh) - rmnc(:,ilow)) &
1062  & * (nsvmec - 1.0)
1063 ! dZ(s_v)/ds_s = dZ(s_s)/ds_v*ds_v/ds_s = 2sdZ(s_v)/ds_v
1064  zmnsint = 2.0*s*(zmns(:,ihigh) - zmns(:,ilow)) &
1065  & * (nsvmec - 1.0)
1066 
1067  rs = sum(rmncint*cosmn_vmec)
1068  zs = sum(zmnsint*sinmn_vmec)
1069 
1070  js = sum(this%context%jksupsmnsf(:,:,i)*sinmn)
1071  ju = sum(this%context%jksupumncf(:,:,i)*cosmn)
1072  jv = sum(this%context%jksupvmncf(:,:,i)*cosmn)
1073 
1074 
1075  IF (this%context%l_asym) THEN
1076 
1077  IF (lasym) THEN
1078  rmnsint = wlow*rmns(:,ilow) &
1079  & + whigh*rmns(:,ihigh)
1080  zmncint = wlow*zmnc(:,ilow) &
1081  & + whigh*zmnc(:,ihigh)
1082 
1083  this%magnetic_cache%rsuv(i,j,k) = &
1084  & this%magnetic_cache%rsuv(i,j,k) + &
1085  & sum(rmnsint*sinmn_vmec)
1086  this%magnetic_cache%zsuv(i,j,k) = &
1087  & this%magnetic_cache%zsuv(i,j,k) + &
1088  & sum(zmncint*cosmn_vmec)
1089 
1090  ru = ru + sum(xm*rmnsint*cosmn_vmec)
1091  rv = rv - sum(xn*rmnsint*cosmn_vmec)
1092  zu = zu - sum(xm*zmncint*sinmn_vmec)
1093  zv = zv + sum(xn*zmncint*sinmn_vmec)
1094 
1095 ! dR(s_v)/ds_s = dR(s_s)/ds_v*ds_v/ds_s = 2sdR(s_v)/ds_v
1096  rmnsint = 2.0*s*(rmns(:,ihigh) - rmns(:,ilow)) &
1097  & * (nsvmec - 1.0)
1098 ! dZ(s_v)/ds_s = dZ(s_s)/ds_v*ds_v/ds_s = 2sdZ(s_v)/ds_v
1099  zmncint = 2.0*s*(zmnc(:,ihigh) - zmnc(:,ilow)) &
1100  & * (nsvmec - 1.0)
1101 
1102  rs = rs + sum(rmnsint*sinmn_vmec)
1103  zs = zs + sum(zmncint*cosmn_vmec)
1104 
1105  END IF
1106 
1107  js = js + sum(this%context%jksupsmncf(:,:,i)*cosmn)
1108  ju = ju + sum(this%context%jksupumnsf(:,:,i)*sinmn)
1109  jv = jv + sum(this%context%jksupvmnsf(:,:,i)*sinmn)
1110 
1111  END IF
1112 
1113  js = js/(this%context%b_factor*mu0)
1114  ju = ju/(this%context%b_factor*mu0)
1115  jv = jv/(this%context%b_factor*mu0)
1116 
1117  this%magnetic_cache%jrsuv(i,j,k) = &
1118  & (js*rs + ju*ru + jv*rv)
1119  this%magnetic_cache%jphisuv(i,j,k) = &
1120  & jv*this%magnetic_cache%rsuv(i,j,k)
1121  this%magnetic_cache%jzsuv(i,j,k) = &
1122  & js*zs + ju*zu + jv*zv
1123  END DO
1124  END DO
1125  END DO
1126 !$OMP END DO
1127 
1128  DEALLOCATE(cosmn)
1129  DEALLOCATE(sinmn)
1130 
1131  DEALLOCATE(cosmn_vmec)
1132  DEALLOCATE(sinmn_vmec)
1133 
1134  DEALLOCATE(rmncint)
1135  DEALLOCATE(zmnsint)
1136  IF (lasym) THEN
1137  DEALLOCATE(rmnsint)
1138  DEALLOCATE(zmncint)
1139  END IF
1140 
1141 !$OMP END PARALLEL
1142 
1143  this%magnetic_cache%zsuv = this%magnetic_cache%zsuv &
1144  & + this%z_offset
1145  END IF
1146 
1147 ! Compute the conducting shell if needed. If the kruv array is associated, the
1148 ! the conducting shell is being used. The total conducting shell signal for a
1149 ! magnetic diagnostic is
1150 !
1151 ! S_i = AreaIntegral K(u,v) . R_i(u,v) da (1)
1152 !
1153 ! For now, It will be assumed that the last closed flux surface of the plasma
1154 ! is the conducting shell. This can be refined later. The total magnetic signal
1155 ! in flux coordinates becomes
1156 !
1157 ! S_i = AreaIntegral (-e^s/|e^s| X B/mu0) . R_i * |J||e^s|dudv (2)
1158 !
1159 ! Using the idenity e^s = 1/J*e_u X e_v equation 2 becomes
1160 !
1161 ! S_i = AreaIntegral -e_u X e_v X B/mu0 . R_i |J|/J dudv (3)
1162 !
1163 ! From the triple cross product idenity equation 3 reduces to
1164 !
1165 ! S_i = AreaIntegral -1/mu0(B_u e_v - B_v e_u) . R_i |J|/J dudv (4)
1166 !
1167 ! Computing the conducting shell surface current is simply.
1168 !
1169 ! K(u,v) = -1/mu0(B_u*e_v - B_v*e_u) (5)
1170 !
1171 ! NOTE: The shell current cannot be used with the Z position shift.
1172 !
1173 ! NOTE: Equation 2.5.49a in "Flux Coordinates and Magnetic Field Structure" by
1174 ! W.D.D'haeseleer, W.N.G.Hitchon, J.D.Callen and J.L.Shohet assumes a
1175 ! positive jacobian. The correct equation should be
1176 !
1177 ! dS(i)=|J||e^i| du^j dv^k
1178 !
1179 ! As a result the signal due to the conducting shell contains a |J|/J
1180 ! term. The sign of the jacobian needs to be taken into account. This is
1181 ! taken care of in vmec_get_area_int_element.
1182 !
1183 ! Compute virtual casing points if needed, If the kruv_full array is
1184 ! associated, the point magnetics are being used. The total signal for a point
1185 ! is
1186 !
1187 ! B(x) = mu0/2Pi AreaIntegral K(u,v) x (x - x')/|x - x'|^3 da (6)
1188 !
1189 ! The surface used is the last closed flux surface. Simular to the procedure
1190 ! for the conducting shell, the different area element and surface current in
1191 ! equation 2 reduce the point field to.
1192 !
1193 ! B(x) = mu0/2Pi AreaIntegral -1/mu0(B_u*e_v - B_v*e_u)
1194 ! x (x - x')/|x - x'|^3 |J|/J dudv (7)
1195 !
1196 ! Computing the virtual casing surface current is simply.
1197 !
1198 ! K(u,v) = -1/mu0(B_u*e_v - B_v*e_u)
1199 !
1200 ! NOTE: Unlike the shell current, the virtual casing can be used with the Z
1201 ! position shift.
1202 
1203  IF (ASSOCIATED(this%magnetic_cache%kruv) .or. &
1204  & ASSOCIATED(this%magnetic_cache%kxuv_full)) THEN
1205 
1206  ALLOCATE(bsubuc(0:mpol,-ntor:ntor))
1207  ALLOCATE(bsubvc(0:mpol,-ntor:ntor))
1208 
1209  nums = SIZE(this%context%bsubumnch, 3)
1210 
1211 ! Interpolate the b_u and b_v onto the last closed flux surface.
1212  bsubuc = 1.5*this%context%bsubumnch(:,:,nums) &
1213  & - 0.5*this%context%bsubumnch(:,:,nums - 1)
1214  bsubvc = 1.5*this%context%bsubvmnch(:,:,nums) &
1215  & - 0.5*this%context%bsubvmnch(:,:,nums - 1)
1216 
1217  IF (this%context%l_asym) THEN
1218  ALLOCATE(bsubus(0:mpol,-ntor:ntor))
1219  ALLOCATE(bsubvs(0:mpol,-ntor:ntor))
1220 
1221  bsubus = 1.5*this%context%bsubumnsh(:,:,nums) &
1222  & - 0.5*this%context%bsubumnsh(:,:,nums - 1)
1223  bsubvs = 1.5*this%context%bsubvmnsh(:,:,nums) &
1224  & - 0.5*this%context%bsubvmnsh(:,:,nums - 1)
1225  END IF
1226  END IF
1227 
1228  IF (ASSOCIATED(this%magnetic_cache%kruv)) THEN
1229  numu = SIZE(this%magnetic_cache%kruv, 1)
1230  numv = SIZE(this%magnetic_cache%kruv, 2)
1231 
1232 ! U grid is square with the s grid.
1233 !$OMP PARALLEL
1234 !$OMP& DEFAULT(SHARED)
1235 !$OMP& PRIVATE(u,v,cosmn,sinmn,cosmn_vmec,sinmn_vmec,j,k,m,n,ru,rv,zu,
1236 !$OMP& zv,r,bu,bv)
1237 
1238  ALLOCATE(cosmn(0:mpol,-ntor:ntor))
1239  IF (this%context%l_asym) THEN
1240  ALLOCATE(sinmn(0:mpol,-ntor:ntor))
1241  END IF
1242 
1243  ALLOCATE(cosmn_vmec(mnmax))
1244  ALLOCATE(sinmn_vmec(mnmax))
1245 
1246 !$OMP DO
1247 !$OMP& SCHEDULE(STATIC)
1248  DO k = 1, numv
1249  v = (k - 1)*this%magnetic_cache%dv_a + this%phi_offset
1250 
1251  DO j = 1, numu
1252  u = (j - 1)*this%magnetic_cache%du_a
1253 
1254  cosmn_vmec = cos(xm*u - xn*v)
1255  sinmn_vmec = sin(xm*u - xn*v)
1256 
1257  DO n = -ntor, ntor
1258  DO m = 0, mpol
1259  cosmn(m,n) = cos(m*u + nfp*n*v)
1260  END DO
1261  END DO
1262 
1263  r = sum(rmncint*cosmn_vmec)
1264  ru = -sum(xm*rmnc(:,nsvmec)*sinmn_vmec)
1265  rv = sum(xn*rmnc(:,nsvmec)*sinmn_vmec)
1266  zu = sum(xm*zmns(:,nsvmec)*cosmn_vmec)
1267  zv = -sum(xn*zmns(:,nsvmec)*cosmn_vmec)
1268 
1269  bu = sum(bsubuc*cosmn)
1270  bv = sum(bsubvc*cosmn)
1271 
1272  IF (this%context%l_asym) THEN
1273  IF (lasym) THEN
1274  r = r + sum(rmnsint*sinmn_vmec)
1275  ru = ru + sum(xm*rmns(:,nsvmec)*cosmn_vmec)
1276  rv = rv - sum(xn*rmns(:,nsvmec)*cosmn_vmec)
1277  zu = zu - sum(xm*zmnc(:,nsvmec)*sinmn_vmec)
1278  zv = zv + sum(xn*zmnc(:,nsvmec)*sinmn_vmec)
1279  END IF
1280 
1281  DO n = -ntor, ntor
1282  DO m = 0, mpol
1283  sinmn(m,n) = sin(m*u + nfp*n*v)
1284  END DO
1285  END DO
1286 
1287  bu = bu + sum(bsubus*sinmn)
1288  bv = bv + sum(bsubvs*sinmn)
1289 
1290  END IF
1291 
1292  bu = bu/this%context%b_factor
1293  bv = bv/this%context%b_factor
1294 
1295  this%magnetic_cache%kruv(j,k) = -(bu*rv - bv*ru)/mu0
1296  this%magnetic_cache%kphiuv(j,k) = -bu*r/mu0
1297  this%magnetic_cache%kzuv(j,k) = -(bu*zv - bv*zu)/mu0
1298 
1299  END DO
1300  END DO
1301 !$OMP END DO
1302 
1303  DEALLOCATE(cosmn)
1304  IF (this%context%l_asym) THEN
1305  DEALLOCATE(sinmn)
1306  END IF
1307 
1308  DEALLOCATE(cosmn_vmec)
1309  DEALLOCATE(sinmn_vmec)
1310 !$OMP END PARALLEL
1311 
1312  END IF
1313 
1314  IF (ASSOCIATED(this%magnetic_cache%kxuv_full)) THEN
1315  numu = SIZE(this%magnetic_cache%kxuv_full, 1)
1316  numv = SIZE(this%magnetic_cache%kxuv_full, 2)
1317 !$OMP PARALLEL
1318 !$OMP& DEFAULT(SHARED)
1319 !$OMP& PRIVATE(cosmn,sinmn,cosmn_vmec,sinmn_vmec,k,j,v,u,n,m,r,z,
1320 !$OMP& ru,rv,zu,zv,bu,bv,cosphi,sinphi,kr,kp,i)
1321 
1322  ALLOCATE(cosmn(0:mpol,-ntor:ntor))
1323  IF (this%context%l_asym) THEN
1324  ALLOCATE(sinmn(0:mpol,-ntor:ntor))
1325  END IF
1326 
1327  ALLOCATE(cosmn_vmec(mnmax))
1328  ALLOCATE(sinmn_vmec(mnmax))
1329 
1330 !$OMP DO
1331 !$OMP& SCHEDULE(STATIC)
1332  DO k = 1, numv
1333 ! Rotate the signal with respect to a fixed plasma. The phi and z offsets are
1334 ! applied in siesta_get_ext_b_plasma.
1335  v = (k - 1)*this%magnetic_cache%dv_full
1336 
1337  cosphi = cos(v)
1338  sinphi = sin(v)
1339 
1340  DO j = 1, numu
1341  u = (j - 1)*this%magnetic_cache%du_full
1342 
1343  cosmn_vmec = cos(xm*u - xn*v)
1344  sinmn_vmec = sin(xm*u - xn*v)
1345 
1346  DO n = -ntor, ntor
1347  DO m = 0, mpol
1348  cosmn(m,n) = cos(m*u + nfp*n*v)
1349  END DO
1350  END DO
1351 
1352  r = sum(rmnc(:,nsvmec)*cosmn_vmec)
1353  z = sum(zmns(:,nsvmec)*sinmn_vmec)
1354 
1355  ru = -sum(xm*rmnc(:,nsvmec)*sinmn_vmec)
1356  rv = sum(xn*rmnc(:,nsvmec)*sinmn_vmec)
1357  zu = sum(xm*zmns(:,nsvmec)*cosmn_vmec)
1358  zv = -sum(xn*zmns(:,nsvmec)*cosmn_vmec)
1359 
1360  bu = sum(bsubuc*cosmn)
1361  bv = sum(bsubvc*cosmn)
1362 
1363  IF (this%context%l_asym) THEN
1364 
1365  DO n = -ntor, ntor
1366  DO m = 0, mpol
1367  sinmn(m,n) = sin(m*u + nfp*n*v)
1368  END DO
1369  END DO
1370 
1371  IF (lasym) THEN
1372  r = r + sum(rmns(:,nsvmec)*sinmn_vmec)
1373  z = z + sum(zmnc(:,nsvmec)*cosmn_vmec)
1374 
1375  ru = ru + sum(xm*rmns(:,nsvmec)*cosmn_vmec)
1376  rv = rv - sum(xn*rmns(:,nsvmec)*cosmn_vmec)
1377  zu = zu - sum(xm*zmnc(:,nsvmec)*sinmn_vmec)
1378  zv = zv + sum(xn*zmnc(:,nsvmec)*sinmn_vmec)
1379  END IF
1380 
1381  bu = bu + sum(bsubus*sinmn)
1382  bv = bv + sum(bsubvs*sinmn)
1383  END IF
1384 
1385  bu = bu/this%context%b_factor
1386  bv = bv/this%context%b_factor
1387 
1388  kr = -(bu*rv - bv*ru)/mu0
1389  kp = -bu*r/mu0
1390 
1391  this%magnetic_cache%kxuv_full(j,k) = kr*cosphi &
1392  & - kp*sinphi
1393  this%magnetic_cache%kyuv_full(j,k) = kr*sinphi &
1394  & + kp*cosphi
1395  this%magnetic_cache%kzuv_full(j,k) = -(bu*zv - bv*zu)/mu0
1396 
1397  this%magnetic_cache%x_prime(j,k,1) = r*cosphi
1398  this%magnetic_cache%x_prime(j,k,2) = r*sinphi
1399  this%magnetic_cache%x_prime(j,k,3) = z
1400 
1401  IF (ASSOCIATED(this%magnetic_cache%kxuv_axi)) THEN
1402  r = 0.0
1403  z = 0.0
1404  ru = 0.0
1405  zu = 0.0
1406 
1407  bu = 0.0
1408  bv = 0.0
1409 
1410  DO i = 1, SIZE(xn)
1411  IF (xn(i) .eq. 0.0) THEN
1412  r = r + rmnc(i,nsvmec)*cosmn_vmec(i)
1413  z = z + zmns(i,nsvmec)*sinmn_vmec(i)
1414  ru = ru - xm(i)*rmnc(i,nsvmec)*sinmn_vmec(i)
1415  zu = zu + xm(i)*zmns(i,nsvmec)*cosmn_vmec(i)
1416  END IF
1417  END DO
1418 
1419  bu = sum(bsubuc(:,0)*cosmn(:,0))
1420  bv = sum(bsubvc(:,0)*cosmn(:,0))
1421 
1422  IF (this%context%l_asym) THEN
1423 
1424  DO i = 1, SIZE(xn)
1425  IF (xn(i) .eq. 0.0) THEN
1426  r = r + rmns(i,nsvmec)*sinmn_vmec(i)
1427  z = z + zmnc(i,nsvmec)*cosmn_vmec(i)
1428  ru = ru + xm(i)*rmns(i,nsvmec)*cosmn_vmec(i)
1429  zu = zu - xm(i)*zmnc(i,nsvmec)*sinmn_vmec(i)
1430  END IF
1431  END DO
1432 
1433  bu = bu + sum(bsubus(:,0)*sinmn(:,0))
1434  bv = bv + sum(bsubvs(:,0)*sinmn(:,0))
1435 
1436  END IF
1437 
1438  bu = bu/this%context%b_factor
1439  bv = bv/this%context%b_factor
1440 
1441  kr = bv*ru/mu0
1442  kp = -bu*r/mu0
1443 
1444  this%magnetic_cache%kxuv_axi(j,k) = kr*cosphi &
1445  & - kp*sinphi
1446  this%magnetic_cache%kyuv_axi(j,k) = kr*sinphi &
1447  & + k*cosphi
1448  this%magnetic_cache%kzuv_axi(j,k) = 1.0/mu0*bv*zu
1449 
1450  this%magnetic_cache%x_axi(j,k,1) = r*cosphi
1451  this%magnetic_cache%x_axi(j,k,2) = r*sinphi
1452  this%magnetic_cache%x_axi(j,k,3) = z
1453 
1454  END IF
1455  END DO
1456  END DO
1457 !$OMP END DO
1458 
1459  DEALLOCATE(cosmn)
1460  IF (this%context%l_asym) THEN
1461  DEALLOCATE(sinmn)
1462  END IF
1463 
1464  DEALLOCATE(cosmn_vmec)
1465  DEALLOCATE(sinmn_vmec)
1466 
1467 !$OMP END PARALLEL
1468  END IF
1469 
1470  IF (ASSOCIATED(this%magnetic_cache%kruv) .or. &
1471  & ASSOCIATED(this%magnetic_cache%kxuv_full)) THEN
1472  DEALLOCATE(bsubuc)
1473  DEALLOCATE(bsubvc)
1474 
1475  IF (this%context%l_asym) THEN
1476  DEALLOCATE(bsubus)
1477  DEALLOCATE(bsubvs)
1478  END IF
1479  END IF
1480 
1481  CALL profiler_set_stop_time('siesta_set_magnetic_cache_calc', &
1482  & start_time)
1483 
1484  END SUBROUTINE
1485 
1486 !-------------------------------------------------------------------------------
1494 !-------------------------------------------------------------------------------
1495  SUBROUTINE siesta_set_namelist(this)
1496 
1497  IMPLICIT NONE
1498 
1499 ! Declare Arguments
1500  TYPE (siesta_class), INTENT(in) :: this
1501 
1502 ! local variables
1503  REAL (rprec) :: start_time
1504 
1505 ! Start of executable code
1506  start_time = profiler_get_start_time()
1507 
1508  CALL vmec_set_namelist(this%vmec)
1509 
1510  CALL profiler_set_stop_time('siesta_set_namelist', start_time)
1511 
1512  END SUBROUTINE
1513 
1514 !*******************************************************************************
1515 ! GETTER SUBROUTINES
1516 !*******************************************************************************
1517 !-------------------------------------------------------------------------------
1525 !-------------------------------------------------------------------------------
1526  FUNCTION siesta_get_param_id(this, param_name)
1527 
1528  IMPLICIT NONE
1529 
1530 ! Declare Arguments
1531  INTEGER :: siesta_get_param_id
1532  TYPE (siesta_class), INTENT(in) :: this
1533  CHARACTER (len=*), INTENT(in) :: param_name
1534 
1535 ! local variables
1536  REAL (rprec) :: start_time
1537 
1538 ! Start of executable code
1539  start_time = profiler_get_start_time()
1540 
1541  SELECT CASE (trim(param_name))
1542 
1543  CASE ('pp_ne_b')
1544  siesta_get_param_id = siesta_pp_ne_b_id
1545 
1546  CASE ('pp_ne_as')
1547  siesta_get_param_id = siesta_pp_ne_as_id
1548 
1549  CASE ('pp_ne_af')
1550  siesta_get_param_id = siesta_pp_ne_af_id
1551 
1552  CASE ('pp_sxrem_b_a')
1553  siesta_get_param_id = siesta_pp_sxrem_b_id
1554 
1555  CASE ('pp_sxrem_as_a')
1556  siesta_get_param_id = siesta_pp_sxrem_as_id
1557 
1558  CASE ('pp_sxrem_af_a')
1559  siesta_get_param_id = siesta_pp_sxrem_af_id
1560 
1561  CASE ('pp_te_b')
1562  siesta_get_param_id = siesta_pp_te_b_id
1563 
1564  CASE ('pp_te_as')
1565  siesta_get_param_id = siesta_pp_te_as_id
1566 
1567  CASE ('pp_te_af')
1568  siesta_get_param_id = siesta_pp_te_af_id
1569 
1570  CASE ('pp_ti_b')
1571  siesta_get_param_id = siesta_pp_ti_b_id
1572 
1573  CASE ('pp_ti_as')
1574  siesta_get_param_id = siesta_pp_ti_as_id
1575 
1576  CASE ('pp_ti_af')
1577  siesta_get_param_id = siesta_pp_ti_af_id
1578 
1579  CASE ('phi_offset')
1580  siesta_get_param_id = siesta_phi_offset_id
1581 
1582  CASE ('z_offset')
1583  siesta_get_param_id = siesta_z_offset_id
1584 
1585  CASE ('helpert')
1586  siesta_get_param_id = siesta_helpert_id
1587 
1588  CASE DEFAULT
1589  siesta_get_param_id = vmec_get_param_id(this%vmec, &
1590  & param_name)
1591 
1592  END SELECT
1593 
1594  CALL profiler_set_stop_time('siesta_get_param_id', start_time)
1595 
1596  END FUNCTION
1597 
1598 !-------------------------------------------------------------------------------
1608 !-------------------------------------------------------------------------------
1609  FUNCTION siesta_get_param_value(this, id, i_index, j_index)
1610 
1611  IMPLICIT NONE
1612 
1613 ! Declare Arguments
1614  REAL (rprec) :: siesta_get_param_value
1615  TYPE (siesta_class), INTENT(in) :: this
1616  INTEGER, INTENT(in) :: id
1617  INTEGER, INTENT(in) :: i_index
1618  INTEGER, INTENT(in) :: j_index
1619 
1620 ! local variables
1621  REAL (rprec) :: start_time
1622 
1623 ! Start of executable code
1624  start_time = profiler_get_start_time()
1625 
1626  SELECT CASE (id)
1627 
1628  CASE (siesta_pp_ne_b_id)
1629  siesta_get_param_value = this%ne%b(i_index)
1630 
1631  CASE (siesta_pp_ne_as_id)
1632  siesta_get_param_value = this%ne%as(i_index)
1633 
1634  CASE (siesta_pp_ne_af_id)
1635  siesta_get_param_value = this%ne%af(i_index)
1636 
1637  CASE (siesta_pp_sxrem_b_id)
1638  siesta_get_param_value = this%sxrem(i_index)%p%b(j_index)
1639 
1640  CASE (siesta_pp_sxrem_as_id)
1641  siesta_get_param_value = this%sxrem(i_index)%p%as(j_index)
1642 
1643  CASE (siesta_pp_sxrem_af_id)
1644  siesta_get_param_value = this%sxrem(i_index)%p%af(j_index)
1645 
1646  CASE (siesta_pp_te_b_id)
1647  siesta_get_param_value = this%te%b(i_index)
1648 
1649  CASE (siesta_pp_te_as_id)
1650  siesta_get_param_value = this%te%as(i_index)
1651 
1652  CASE (siesta_pp_te_af_id)
1653  siesta_get_param_value = this%te%af(i_index)
1654 
1655  CASE (siesta_pp_ti_b_id)
1656  siesta_get_param_value = this%ti%b(i_index)
1657 
1658  CASE (siesta_pp_ti_as_id)
1659  siesta_get_param_value = this%ti%as(i_index)
1660 
1661  CASE (siesta_pp_ti_af_id)
1662  siesta_get_param_value = this%ti%af(i_index)
1663 
1664  CASE (siesta_phi_offset_id)
1665  siesta_get_param_value = this%phi_offset
1666 
1667  CASE (siesta_z_offset_id)
1668  siesta_get_param_value = this%z_offset
1669 
1670  CASE (siesta_helpert_id)
1671  siesta_get_param_value = helpert(i_index)
1672 
1673  CASE DEFAULT
1674  siesta_get_param_value = &
1675  & vmec_get_param_value(this%vmec, id, i_index, j_index)
1676 
1677  END SELECT
1678 
1679  CALL profiler_set_stop_time('siesta_get_param_value', start_time)
1680 
1681  END FUNCTION
1682 
1683 !-------------------------------------------------------------------------------
1691 !-------------------------------------------------------------------------------
1692  FUNCTION siesta_get_param_name(this, id)
1693 
1694  IMPLICIT NONE
1695 
1696 ! Declare Arguments
1697  CHARACTER(len=data_name_length) :: siesta_get_param_name
1698  TYPE (siesta_class), INTENT(in) :: this
1699  INTEGER, INTENT(in) :: id
1700 
1701 ! local variables
1702  REAL (rprec) :: start_time
1703 
1704 ! Start of executable code
1705  start_time = profiler_get_start_time()
1706 
1707  SELECT CASE (id)
1708 
1709  CASE (siesta_pp_ne_b_id)
1710  siesta_get_param_name = 'pp_ne_b'
1711 
1712  CASE (siesta_pp_ne_as_id)
1713  siesta_get_param_name = 'pp_ne_as'
1714 
1715  CASE (siesta_pp_ne_af_id)
1716  siesta_get_param_name = 'pp_ne_af'
1717 
1718  CASE (siesta_pp_sxrem_b_id)
1719  siesta_get_param_name = 'pp_sxrem_b_a'
1720 
1721  CASE (siesta_pp_sxrem_as_id)
1722  siesta_get_param_name = 'pp_sxrem_as_a'
1723 
1724  CASE (siesta_pp_sxrem_af_id)
1725  siesta_get_param_name = 'pp_sxrem_af_a'
1726 
1727  CASE (siesta_pp_te_b_id)
1728  siesta_get_param_name = 'pp_te_b'
1729 
1730  CASE (siesta_pp_te_as_id)
1731  siesta_get_param_name = 'pp_te_as'
1732 
1733  CASE (siesta_pp_te_af_id)
1734  siesta_get_param_name = 'pp_te_af'
1735 
1736  CASE (siesta_pp_ti_b_id)
1737  siesta_get_param_name = 'pp_ti_b'
1738 
1739  CASE (siesta_pp_ti_as_id)
1740  siesta_get_param_name = 'pp_ti_as'
1741 
1742  CASE (siesta_pp_ti_af_id)
1743  siesta_get_param_name = 'pp_ti_af'
1744 
1745  CASE (siesta_phi_offset_id)
1746  siesta_get_param_name = 'phi_offset'
1747 
1748  CASE (siesta_z_offset_id)
1749  siesta_get_param_name = 'z_offset'
1750 
1751  CASE (siesta_helpert_id)
1752  siesta_get_param_name = 'helpert'
1753 
1754  CASE DEFAULT
1755  siesta_get_param_name = &
1756  & vmec_get_param_name(this%vmec, id)
1757 
1758  END SELECT
1759 
1760  CALL profiler_set_stop_time('siesta_get_param_name', start_time)
1761 
1762  END FUNCTION
1763 
1764 !-------------------------------------------------------------------------------
1773 !-------------------------------------------------------------------------------
1774  FUNCTION siesta_get_gp_ne_num_hyper_param(this)
1775 
1776  IMPLICIT NONE
1777 
1778 ! Declare Arguments
1779  INTEGER :: siesta_get_gp_ne_num_hyper_param
1780  TYPE (siesta_class), INTENT(in) :: this
1781 
1782 ! local variables
1783  REAL (rprec) :: start_time
1784 
1785 ! Start of executable code
1786  start_time = profiler_get_start_time()
1787 
1788  IF (ASSOCIATED(this%ne)) THEN
1789  siesta_get_gp_ne_num_hyper_param = &
1790  & pprofile_get_gp_num_hyper_param(this%ne)
1791  ELSE
1792  siesta_get_gp_ne_num_hyper_param = 0
1793  END IF
1794 
1795  CALL profiler_set_stop_time('siesta_get_gp_ne_num_hyper_param', &
1796  & start_time)
1797 
1798  END FUNCTION
1799 
1800 !-------------------------------------------------------------------------------
1808 !-------------------------------------------------------------------------------
1809  FUNCTION siesta_get_ne_af(this)
1810 
1811  IMPLICIT NONE
1812 
1813 ! Declare Arguments
1814  REAL (rprec), DIMENSION(:), POINTER :: siesta_get_ne_af
1815  TYPE (siesta_class), INTENT(in) :: this
1816 
1817 ! local variables
1818  REAL (rprec) :: start_time
1819 
1820 ! Start of executable code
1821  start_time = profiler_get_start_time()
1822 
1823  IF (ASSOCIATED(this%ne)) THEN
1824  siesta_get_ne_af => this%ne%af
1825  ELSE
1826  siesta_get_ne_af => null()
1827  END IF
1828 
1829  END FUNCTION
1830 
1831 !-------------------------------------------------------------------------------
1841 !-------------------------------------------------------------------------------
1842  FUNCTION siesta_get_gp_ne_ij(this, i, j)
1843 
1844  IMPLICIT NONE
1845 
1846 ! Declare Arguments
1847  REAL (rprec) :: siesta_get_gp_ne_ij
1848  TYPE (siesta_class), INTENT(in) :: this
1849  INTEGER, INTENT(in) :: i
1850  INTEGER, INTENT(in) :: j
1851 
1852 ! local variables
1853  REAL (rprec) :: start_time
1854 
1855 ! Start of executable code
1856  start_time = profiler_get_start_time()
1857 
1858  IF (ASSOCIATED(this%ne)) THEN
1859  siesta_get_gp_ne_ij = pprofile_get_gp(this%ne, i, j)
1860  ELSE
1861  siesta_get_gp_ne_ij = 0.0
1862  END IF
1863 
1864  CALL profiler_set_stop_time('siesta_get_gp_ne_ij', start_time)
1865 
1866  END FUNCTION
1867 
1868 !-------------------------------------------------------------------------------
1878 !-------------------------------------------------------------------------------
1879  FUNCTION siesta_get_gp_ne_pi(this, x_cart, i)
1880 
1881  IMPLICIT NONE
1882 
1883 ! Declare Arguments
1884  REAL (rprec) :: siesta_get_gp_ne_pi
1885  TYPE (siesta_class), INTENT(in) :: this
1886  REAL (rprec), DIMENSION(3), INTENT(in) :: x_cart
1887  INTEGER, INTENT(in) :: i
1888 
1889 ! local variables
1890  REAL (rprec) :: start_time
1891 
1892 ! Start of executable code
1893  start_time = profiler_get_start_time()
1894 
1895  IF (ASSOCIATED(this%ne)) THEN
1896  siesta_get_gp_ne_pi = pprofile_get_gp(this%ne, &
1897  & siesta_get_p(this, &
1898  & x_cart, &
1899  & .true.), &
1900  & i)
1901  ELSE
1902  siesta_get_gp_ne_pi = 0.0
1903  END IF
1904 
1905  CALL profiler_set_stop_time('siesta_get_gp_ne_pi', start_time)
1906 
1907  END FUNCTION
1908 
1909 !-------------------------------------------------------------------------------
1920 !-------------------------------------------------------------------------------
1921  FUNCTION siesta_get_gp_ne_pp(this, x_cart, y_cart)
1922 
1923  IMPLICIT NONE
1924 
1925 ! Declare Arguments
1926  REAL (rprec) :: siesta_get_gp_ne_pp
1927  TYPE (siesta_class), INTENT(in) :: this
1928  REAL (rprec), DIMENSION(3), INTENT(in) :: x_cart
1929  REAL (rprec), DIMENSION(3), INTENT(in) :: y_cart
1930 
1931 ! local variables
1932  REAL (rprec) :: start_time
1933 
1934 ! Start of executable code
1935  start_time = profiler_get_start_time()
1936 
1937  IF (ASSOCIATED(this%ne)) THEN
1938  siesta_get_gp_ne_pp = &
1939  & pprofile_get_gp(this%ne, siesta_get_p(this, x_cart, .true.), &
1940  & siesta_get_p(this, y_cart, .true.))
1941  ELSE
1942  siesta_get_gp_ne_pp = 0.0
1943  END IF
1944 
1945  CALL profiler_set_stop_time('siesta_get_gp_ne_pp', start_time)
1946 
1947  END FUNCTION
1948 
1949 !-------------------------------------------------------------------------------
1959 !-------------------------------------------------------------------------------
1960  FUNCTION siesta_get_ne_cart(this, x_cart)
1961 
1962  IMPLICIT NONE
1963 
1964 ! Declare Arguments
1965  REAL (rprec) :: siesta_get_ne_cart
1966  TYPE (siesta_class), INTENT(in) :: this
1967  REAL (rprec), DIMENSION(3), INTENT(in) :: x_cart
1968 
1969 ! local variables
1970  REAL (rprec) :: start_time
1971 
1972 ! Start of executable code
1973  start_time = profiler_get_start_time()
1974 
1975  siesta_get_ne_cart = &
1976  & pprofile_get_value(this%ne, siesta_get_p(this, x_cart, .true.))
1977 
1978  CALL profiler_set_stop_time('siesta_get_ne_cart', start_time)
1979 
1980  END FUNCTION
1981 
1982 !-------------------------------------------------------------------------------
1992 !-------------------------------------------------------------------------------
1993  FUNCTION siesta_get_ne_radial(this, s)
1994 
1995  IMPLICIT NONE
1996 
1997 ! Declare Arguments
1998  REAL (rprec) :: siesta_get_ne_radial
1999  TYPE (siesta_class), INTENT(in) :: this
2000  REAL (rprec), INTENT(in) :: s
2001 
2002 ! local variables
2003  REAL (rprec) :: start_time
2004 
2005 ! Start of executable code
2006  start_time = profiler_get_start_time()
2007 
2008  siesta_get_ne_radial = &
2009  & pprofile_get_value(this%ne, siesta_get_p(this, s, .true.))
2010 
2011  CALL profiler_set_stop_time('siesta_get_ne_radial', start_time)
2012 
2013  END FUNCTION
2014 
2015 !-------------------------------------------------------------------------------
2024 !-------------------------------------------------------------------------------
2025  FUNCTION siesta_get_gp_te_num_hyper_param(this)
2026 
2027  IMPLICIT NONE
2028 
2029 ! Declare Arguments
2030  INTEGER :: siesta_get_gp_te_num_hyper_param
2031  TYPE (siesta_class), INTENT(in) :: this
2032 
2033 ! local variables
2034  REAL (rprec) :: start_time
2035 
2036 ! Start of executable code
2037  start_time = profiler_get_start_time()
2038 
2039  IF (ASSOCIATED(this%te)) THEN
2040  siesta_get_gp_te_num_hyper_param = &
2041  & pprofile_get_gp_num_hyper_param(this%te)
2042  ELSE
2043  siesta_get_gp_te_num_hyper_param = 0
2044  END IF
2045 
2046  CALL profiler_set_stop_time('siesta_get_gp_te_num_hyper_param', &
2047  & start_time)
2048 
2049  END FUNCTION
2050 
2051 !-------------------------------------------------------------------------------
2059 !-------------------------------------------------------------------------------
2060  FUNCTION siesta_get_te_af(this)
2061 
2062  IMPLICIT NONE
2063 
2064 ! Declare Arguments
2065  REAL (rprec), DIMENSION(:), POINTER :: siesta_get_te_af
2066  TYPE (siesta_class), INTENT(in) :: this
2067 
2068 ! local variables
2069  REAL (rprec) :: start_time
2070 
2071 ! Start of executable code
2072  start_time = profiler_get_start_time()
2073 
2074  IF (ASSOCIATED(this%te)) THEN
2075  siesta_get_te_af => this%te%af
2076  ELSE
2077  siesta_get_te_af => null()
2078  END IF
2079 
2080  CALL profiler_set_stop_time('siesta_get_te_af', start_time)
2081 
2082  END FUNCTION
2083 
2084 !-------------------------------------------------------------------------------
2094 !-------------------------------------------------------------------------------
2095  FUNCTION siesta_get_gp_te_ij(this, i, j)
2096 
2097  IMPLICIT NONE
2098 
2099 ! Declare Arguments
2100  REAL (rprec) :: siesta_get_gp_te_ij
2101  TYPE (siesta_class), INTENT(in) :: this
2102  INTEGER, INTENT(in) :: i
2103  INTEGER, INTENT(in) :: j
2104 
2105 ! local variables
2106  REAL (rprec) :: start_time
2107 
2108 ! Start of executable code
2109  start_time = profiler_get_start_time()
2110 
2111  IF (ASSOCIATED(this%te)) THEN
2112  siesta_get_gp_te_ij = pprofile_get_gp(this%te, i, j)
2113  ELSE
2114  siesta_get_gp_te_ij = 0.0
2115  END IF
2116 
2117  CALL profiler_set_stop_time('siesta_get_gp_te_ij', start_time)
2118 
2119  END FUNCTION
2120 
2121 !-------------------------------------------------------------------------------
2132 !-------------------------------------------------------------------------------
2133  FUNCTION siesta_get_gp_te_pi(this, x_cart, i)
2134 
2135  IMPLICIT NONE
2136 
2137 ! Declare Arguments
2138  REAL (rprec) :: siesta_get_gp_te_pi
2139  TYPE (siesta_class), INTENT(in) :: this
2140  REAL (rprec), DIMENSION(3), INTENT(in) :: x_cart
2141  INTEGER, INTENT(in) :: i
2142 
2143 ! local variables
2144  REAL (rprec) :: start_time
2145 
2146 ! Start of executable code
2147  start_time = profiler_get_start_time()
2148 
2149  IF (ASSOCIATED(this%te)) THEN
2150  siesta_get_gp_te_pi = pprofile_get_gp(this%te, &
2151  & siesta_get_p(this, &
2152  & x_cart, &
2153  & .true.), &
2154  & i)
2155  ELSE
2156  siesta_get_gp_te_pi = 0.0
2157  END IF
2158 
2159  CALL profiler_set_stop_time('siesta_get_gp_te_pi', start_time)
2160 
2161  END FUNCTION
2162 
2163 !-------------------------------------------------------------------------------
2176 !-------------------------------------------------------------------------------
2177  FUNCTION siesta_get_gp_te_pp(this, x_cart, y_cart)
2178 
2179  IMPLICIT NONE
2180 
2181 ! Declare Arguments
2182  REAL (rprec) :: siesta_get_gp_te_pp
2183  TYPE (siesta_class), INTENT(in) :: this
2184  REAL (rprec), DIMENSION(3), INTENT(in) :: x_cart
2185  REAL (rprec), DIMENSION(3), INTENT(in) :: y_cart
2186 
2187 ! local variables
2188  REAL (rprec) :: start_time
2189 
2190 ! Start of executable code
2191  start_time = profiler_get_start_time()
2192 
2193  IF (ASSOCIATED(this%te)) THEN
2194  siesta_get_gp_te_pp = &
2195  & pprofile_get_gp(this%te, siesta_get_p(this, x_cart, .true.), &
2196  & siesta_get_p(this, y_cart, .true.))
2197  ELSE
2198  siesta_get_gp_te_pp = 0.0
2199  END IF
2200 
2201  CALL profiler_set_stop_time('siesta_get_gp_te_pp', start_time)
2202 
2203  END FUNCTION
2204 
2205 !-------------------------------------------------------------------------------
2215 !-------------------------------------------------------------------------------
2216  FUNCTION siesta_get_te_cart(this, x_cart)
2217 
2218  IMPLICIT NONE
2219 
2220 ! Declare Arguments
2221  REAL (rprec) :: siesta_get_te_cart
2222  TYPE (siesta_class), INTENT(in) :: this
2223  REAL (rprec), DIMENSION(3), INTENT(in) :: x_cart
2224 
2225 ! local variables
2226  REAL (rprec) :: start_time
2227 
2228 ! Start of executable code
2229  start_time = profiler_get_start_time()
2230 
2231  siesta_get_te_cart = &
2232  & pprofile_get_value(this%te, siesta_get_p(this, x_cart, .true.))
2233 
2234  CALL profiler_set_stop_time('siesta_get_te_cart', start_time)
2235 
2236  END FUNCTION
2237 
2238 !-------------------------------------------------------------------------------
2248 !-------------------------------------------------------------------------------
2249  FUNCTION siesta_get_te_radial(this, s)
2250 
2251  IMPLICIT NONE
2252 
2253 ! Declare Arguments
2254  REAL (rprec) :: siesta_get_te_radial
2255  TYPE (siesta_class), INTENT(in) :: this
2256  REAL (rprec), INTENT(in) :: s
2257 
2258 ! local variables
2259  REAL (rprec) :: start_time
2260 
2261 ! Start of executable code
2262  start_time = profiler_get_start_time()
2263 
2264  siesta_get_te_radial = &
2265  & pprofile_get_value(this%te, siesta_get_p(this, s, .true.))
2266 
2267  CALL profiler_set_stop_time('siesta_get_te_radial', start_time)
2268 
2269  END FUNCTION
2270 
2271 !-------------------------------------------------------------------------------
2280 !-------------------------------------------------------------------------------
2281  FUNCTION siesta_get_gp_ti_num_hyper_param(this)
2282 
2283  IMPLICIT NONE
2284 
2285 ! Declare Arguments
2286  INTEGER :: siesta_get_gp_ti_num_hyper_param
2287  TYPE (siesta_class), INTENT(in) :: this
2288 
2289 ! local variables
2290  REAL (rprec) :: start_time
2291 
2292 ! Start of executable code
2293  start_time = profiler_get_start_time()
2294 
2295  IF (ASSOCIATED(this%ti)) THEN
2296  siesta_get_gp_ti_num_hyper_param = &
2297  & pprofile_get_gp_num_hyper_param(this%ti)
2298  ELSE
2299  siesta_get_gp_ti_num_hyper_param = 0
2300  END IF
2301 
2302  CALL profiler_set_stop_time('siesta_get_gp_ti_num_hyper_param', &
2303  & start_time)
2304 
2305  END FUNCTION
2306 
2307 !-------------------------------------------------------------------------------
2315 !-------------------------------------------------------------------------------
2316  FUNCTION siesta_get_ti_af(this)
2317 
2318  IMPLICIT NONE
2319 
2320 ! Declare Arguments
2321  REAL (rprec), DIMENSION(:), POINTER :: siesta_get_ti_af
2322  TYPE (siesta_class), INTENT(in) :: this
2323 
2324 ! local variables
2325  REAL (rprec) :: start_time
2326 
2327 ! Start of executable code
2328  start_time = profiler_get_start_time()
2329 
2330  IF (ASSOCIATED(this%ti)) THEN
2331  siesta_get_ti_af => this%ti%af
2332  ELSE
2333  siesta_get_ti_af => null()
2334  END IF
2335 
2336  CALL profiler_set_stop_time('siesta_get_ti_af', start_time)
2337 
2338  END FUNCTION
2339 
2340 !-------------------------------------------------------------------------------
2350 !-------------------------------------------------------------------------------
2351  FUNCTION siesta_get_gp_ti_ij(this, i, j)
2352 
2353  IMPLICIT NONE
2354 
2355 ! Declare Arguments
2356  REAL (rprec) :: siesta_get_gp_ti_ij
2357  TYPE (siesta_class), INTENT(in) :: this
2358  INTEGER, INTENT(in) :: i
2359  INTEGER, INTENT(in) :: j
2360 
2361 ! local variables
2362  REAL (rprec) :: start_time
2363 
2364 ! Start of executable code
2365  start_time = profiler_get_start_time()
2366 
2367  IF (ASSOCIATED(this%ti)) THEN
2368  siesta_get_gp_ti_ij = pprofile_get_gp(this%ti, i, j)
2369  ELSE
2370  siesta_get_gp_ti_ij = 0.0
2371  END IF
2372 
2373  CALL profiler_set_stop_time('siesta_get_gp_ti_ij', start_time)
2374 
2375  END FUNCTION
2376 
2377 !-------------------------------------------------------------------------------
2387 !-------------------------------------------------------------------------------
2388  FUNCTION siesta_get_gp_ti_pi(this, x_cart, i)
2389 
2390  IMPLICIT NONE
2391 
2392 ! Declare Arguments
2393  REAL (rprec) :: siesta_get_gp_ti_pi
2394  TYPE (siesta_class), INTENT(in) :: this
2395  REAL (rprec), DIMENSION(3), INTENT(in) :: x_cart
2396  INTEGER, INTENT(in) :: i
2397 
2398 ! local variables
2399  REAL (rprec) :: start_time
2400 
2401 ! Start of executable code
2402  start_time = profiler_get_start_time()
2403 
2404  IF (ASSOCIATED(this%ti)) THEN
2405  siesta_get_gp_ti_pi = &
2406  & pprofile_get_gp(this%ti, siesta_get_p(this, x_cart, .true.), &
2407  & i)
2408  ELSE
2409  siesta_get_gp_ti_pi = 0.0
2410  END IF
2411 
2412  CALL profiler_set_stop_time('siesta_get_gp_ti_pi', start_time)
2413 
2414  END FUNCTION
2415 
2416 !-------------------------------------------------------------------------------
2427 !-------------------------------------------------------------------------------
2428  FUNCTION siesta_get_gp_ti_pp(this, x_cart, y_cart)
2429 
2430  IMPLICIT NONE
2431 
2432 ! Declare Arguments
2433  REAL (rprec) :: siesta_get_gp_ti_pp
2434  TYPE (siesta_class), INTENT(in) :: this
2435  REAL (rprec), DIMENSION(3), INTENT(in) :: x_cart
2436  REAL (rprec), DIMENSION(3), INTENT(in) :: y_cart
2437 
2438 ! local variables
2439  REAL (rprec) :: start_time
2440 
2441 ! Start of executable code
2442  start_time = profiler_get_start_time()
2443 
2444  IF (ASSOCIATED(this%ti)) THEN
2445  siesta_get_gp_ti_pp = &
2446  & pprofile_get_gp(this%ti, siesta_get_p(this, x_cart, .true.), &
2447  & siesta_get_p(this, y_cart, .true.))
2448  ELSE
2449  siesta_get_gp_ti_pp = 0.0
2450  END IF
2451 
2452  CALL profiler_set_stop_time('siesta_get_gp_ti_pp', start_time)
2453 
2454  END FUNCTION
2455 
2456 !-------------------------------------------------------------------------------
2466 !-------------------------------------------------------------------------------
2467  FUNCTION siesta_get_ti_cart(this, x_cart)
2468 
2469  IMPLICIT NONE
2470 
2471 ! Declare Arguments
2472  REAL (rprec) :: siesta_get_ti_cart
2473  TYPE (siesta_class), INTENT(in) :: this
2474  REAL (rprec), DIMENSION(3), INTENT(in) :: x_cart
2475 
2476 ! local variables
2477  REAL (rprec) :: start_time
2478 
2479 ! Start of executable code
2480  start_time = profiler_get_start_time()
2481 
2482  siesta_get_ti_cart = &
2483  & pprofile_get_value(this%ti, siesta_get_p(this, x_cart, .true.))
2484 
2485  CALL profiler_set_stop_time('siesta_get_ti_cart', start_time)
2486 
2487  END FUNCTION
2488 
2489 !-------------------------------------------------------------------------------
2499 !-------------------------------------------------------------------------------
2500  FUNCTION siesta_get_ti_radial(this, s)
2501 
2502  IMPLICIT NONE
2503 
2504 ! Declare Arguments
2505  REAL (rprec) :: siesta_get_ti_radial
2506  TYPE (siesta_class), INTENT(in) :: this
2507  REAL (rprec), INTENT(in) :: s
2508 
2509 ! local variables
2510  REAL (rprec) :: start_time
2511 
2512 ! Start of executable code
2513  start_time = profiler_get_start_time()
2514 
2515  siesta_get_ti_radial = &
2516  & pprofile_get_value(this%ti, siesta_get_p(this, s, .true.))
2517 
2518  CALL profiler_set_stop_time('siesta_get_ti_radial', start_time)
2519 
2520  END FUNCTION
2521 
2522 !-------------------------------------------------------------------------------
2532 !-------------------------------------------------------------------------------
2533  FUNCTION siesta_get_gp_sxrem_num_hyper_param(this, index)
2534 
2535  IMPLICIT NONE
2536 
2537 ! Declare Arguments
2538  INTEGER :: siesta_get_gp_sxrem_num_hyper_param
2539  TYPE (siesta_class), INTENT(in) :: this
2540  INTEGER, INTENT(in) :: index
2541 
2542 ! local variables
2543  REAL (rprec) :: start_time
2544 
2545 ! Start of executable code
2546  start_time = profiler_get_start_time()
2547 
2548  IF (ASSOCIATED(this%sxrem) .and. &
2549  & (index .le. SIZE(this%sxrem))) THEN
2550  siesta_get_gp_sxrem_num_hyper_param = &
2551  & pprofile_get_gp_num_hyper_param(this%sxrem(index)%p)
2552  ELSE
2553  siesta_get_gp_sxrem_num_hyper_param = 0
2554  END IF
2555 
2556  CALL profiler_set_stop_time('siesta_get_gp_sxrem_num_hyper_param', &
2557  & start_time)
2558 
2559  END FUNCTION
2560 
2561 !-------------------------------------------------------------------------------
2570 !-------------------------------------------------------------------------------
2571  FUNCTION siesta_get_sxrem_af(this, index)
2572 
2573  IMPLICIT NONE
2574 
2575 ! Declare Arguments
2576  REAL (rprec), DIMENSION(:), POINTER :: siesta_get_sxrem_af
2577  TYPE (siesta_class), INTENT(in) :: this
2578  INTEGER, INTENT(in) :: index
2579 
2580 ! local variables
2581  REAL (rprec) :: start_time
2582 
2583 ! Start of executable code
2584  start_time = profiler_get_start_time()
2585 
2586  IF (ASSOCIATED(this%sxrem) .and. &
2587  & (index .le. SIZE(this%sxrem))) THEN
2588  siesta_get_sxrem_af => this%sxrem(index)%p%af
2589  ELSE
2590  siesta_get_sxrem_af => null()
2591  END IF
2592 
2593  CALL profiler_set_stop_time('siesta_get_sxrem_af', start_time)
2594 
2595  END FUNCTION
2596 
2597 !-------------------------------------------------------------------------------
2608 !-------------------------------------------------------------------------------
2609  FUNCTION siesta_get_gp_sxrem_ij(this, i, j, index)
2610 
2611  IMPLICIT NONE
2612 
2613 ! Declare Arguments
2614  REAL (rprec) :: siesta_get_gp_sxrem_ij
2615  TYPE (siesta_class), INTENT(in) :: this
2616  INTEGER, INTENT(in) :: i
2617  INTEGER, INTENT(in) :: j
2618  INTEGER, INTENT(in) :: index
2619 
2620 ! local variables
2621  REAL (rprec) :: start_time
2622 
2623 ! Start of executable code
2624  start_time = profiler_get_start_time()
2625 
2626  IF (ASSOCIATED(this%sxrem) .and. &
2627  & (index .le. SIZE(this%sxrem))) THEN
2628  siesta_get_gp_sxrem_ij = &
2629  & pprofile_get_gp(this%sxrem(index)%p, i, j)
2630  ELSE
2631  siesta_get_gp_sxrem_ij = 0.0
2632  END IF
2633 
2634  CALL profiler_set_stop_time('siesta_get_gp_sxrem_ij', start_time)
2635 
2636  END FUNCTION
2637 
2638 !-------------------------------------------------------------------------------
2650 !-------------------------------------------------------------------------------
2651  FUNCTION siesta_get_gp_sxrem_pi(this, x_cart, i, index)
2652 
2653  IMPLICIT NONE
2654 
2655 ! Declare Arguments
2656  REAL (rprec) :: siesta_get_gp_sxrem_pi
2657  TYPE (siesta_class), INTENT(in) :: this
2658  REAL (rprec), DIMENSION(3), INTENT(in) :: x_cart
2659  INTEGER, INTENT(in) :: i
2660  INTEGER, INTENT(in) :: index
2661 
2662 ! local variables
2663  REAL (rprec) :: start_time
2664 
2665 ! Start of executable code
2666  start_time = profiler_get_start_time()
2667 
2668  IF (ASSOCIATED(this%sxrem) .and. &
2669  & (index .le. SIZE(this%sxrem))) THEN
2670  siesta_get_gp_sxrem_pi = pprofile_get_gp(this%sxrem(index)%p, &
2671  & siesta_get_p(this, &
2672  & x_cart, &
2673  & .true.), &
2674  & i)
2675  ELSE
2676  siesta_get_gp_sxrem_pi = 0.0
2677  END IF
2678 
2679  CALL profiler_set_stop_time('siesta_get_gp_sxrem_pi', start_time)
2680 
2681  END FUNCTION
2682 
2683 !-------------------------------------------------------------------------------
2697 !-------------------------------------------------------------------------------
2698  FUNCTION siesta_get_gp_sxrem_pp(this, x_cart, y_cart, index)
2699 
2700  IMPLICIT NONE
2701 
2702 ! Declare Arguments
2703  REAL (rprec) :: siesta_get_gp_sxrem_pp
2704  TYPE (siesta_class), INTENT(in) :: this
2705  REAL (rprec), DIMENSION(3), INTENT(in) :: x_cart
2706  REAL (rprec), DIMENSION(3), INTENT(in) :: y_cart
2707  INTEGER, INTENT(in) :: index
2708 
2709 ! local variables
2710  REAL (rprec) :: start_time
2711 
2712 ! Start of executable code
2713  start_time = profiler_get_start_time()
2714 
2715  IF (ASSOCIATED(this%sxrem) .and. &
2716  & (index .le. SIZE(this%sxrem))) THEN
2717  siesta_get_gp_sxrem_pp = &
2718  & pprofile_get_gp(this%sxrem(index)%p, &
2719  & siesta_get_p(this, x_cart, .true.), &
2720  & siesta_get_p(this, y_cart, .true.))
2721  ELSE
2722  siesta_get_gp_sxrem_pp = 0.0
2723  END IF
2724 
2725  CALL profiler_set_stop_time('siesta_get_gp_sxrem_pp', start_time)
2726 
2727  END FUNCTION
2728 
2729 !-------------------------------------------------------------------------------
2740 !-------------------------------------------------------------------------------
2741  FUNCTION siesta_get_sxrem_cart(this, x_cart, index)
2742 
2743  IMPLICIT NONE
2744 
2745 ! Declare Arguments
2746  REAL (rprec) :: siesta_get_sxrem_cart
2747  TYPE (siesta_class), INTENT(in) :: this
2748  REAL (rprec), DIMENSION(3), INTENT(in) :: x_cart
2749  INTEGER, INTENT(in) :: index
2750 
2751 ! local variables
2752  REAL (rprec) :: start_time
2753 
2754 ! Start of executable code
2755  start_time = profiler_get_start_time()
2756 
2757  IF (ASSOCIATED(this%sxrem) .and. &
2758  & (index .le. SIZE(this%sxrem))) THEN
2759  siesta_get_sxrem_cart = &
2760  & pprofile_get_value(this%sxrem(index)%p, &
2761  & siesta_get_p(this, x_cart, .true.))
2762  ELSE
2763  siesta_get_sxrem_cart = 0.0
2764  END IF
2765 
2766  CALL profiler_set_stop_time('siesta_get_sxrem_cart', start_time)
2767 
2768  END FUNCTION
2769 
2770 !-------------------------------------------------------------------------------
2781 !-------------------------------------------------------------------------------
2782  FUNCTION siesta_get_sxrem_radial(this, s, index)
2783 
2784  IMPLICIT NONE
2785 
2786 ! Declare Arguments
2787  REAL (rprec) :: siesta_get_sxrem_radial
2788  TYPE (siesta_class), INTENT(in) :: this
2789  REAL (rprec), INTENT(in) :: s
2790  INTEGER, INTENT(in) :: index
2791 
2792 ! local variables
2793  REAL (rprec) :: start_time
2794 
2795 ! Start of executable code
2796  start_time = profiler_get_start_time()
2797 
2798  IF (ASSOCIATED(this%sxrem) .and. &
2799  & (index .le. SIZE(this%sxrem))) THEN
2800  siesta_get_sxrem_radial = &
2801  & pprofile_get_value(this%sxrem(index)%p, &
2802  & siesta_get_p(this, s, .true.))
2803  ELSE
2804  siesta_get_sxrem_radial = 0.0
2805  END IF
2806 
2807  CALL profiler_set_stop_time('siesta_get_sxrem_radial', start_time)
2808 
2809  END FUNCTION
2810 
2811 !-------------------------------------------------------------------------------
2820 !-------------------------------------------------------------------------------
2821  FUNCTION siesta_get_p_cart(this, x_cart, normalize)
2822 
2823  IMPLICIT NONE
2824 
2825 ! Declare Arguments
2826  REAL (rprec) :: siesta_get_p_cart
2827  TYPE (siesta_class), INTENT(in) :: this
2828  REAL (rprec), DIMENSION(3), INTENT(in) :: x_cart
2829  LOGICAL, INTENT(in) :: normalize
2830 
2831 ! local variables
2832  REAL (rprec) :: start_time
2833 
2834 ! Start of executable code
2835  start_time = profiler_get_start_time()
2836 
2837  siesta_get_p_cart = &
2838  & siesta_get_p_flux(this, siesta_get_suv(this, x_cart), &
2839  & normalize)
2840 
2841  CALL profiler_set_stop_time('siesta_get_p_cart', start_time)
2842 
2843  END FUNCTION
2844 
2845 !-------------------------------------------------------------------------------
2854 !-------------------------------------------------------------------------------
2855  FUNCTION siesta_get_p_radial(this, s, normalize)
2856 
2857  IMPLICIT NONE
2858 
2859 ! Declare Arguments
2860  REAL (rprec) :: siesta_get_p_radial
2861  TYPE (siesta_class), INTENT(in) :: this
2862  REAL (rprec), INTENT(in) :: s
2863  LOGICAL, INTENT(in) :: normalize
2864 
2865 ! local variables
2866  REAL (rprec), DIMENSION(3) :: flux
2867  REAL (rprec) :: start_time
2868 
2869 ! Start of executable code
2870  start_time = profiler_get_start_time()
2871 
2872  flux(1) = min(s, 1.0)
2873  flux(2:3) = 0.0
2874  siesta_get_p_radial = siesta_get_p_flux(this, flux, normalize)
2875 
2876  CALL profiler_set_stop_time('siesta_get_p_cart', start_time)
2877 
2878  END FUNCTION
2879 
2880 !-------------------------------------------------------------------------------
2889 !-------------------------------------------------------------------------------
2890  FUNCTION siesta_get_p_flux(this, flux, normalize)
2891  USE read_wout_mod, only: nfp
2892  USE stel_constants, only: mu0
2893 
2894  IMPLICIT NONE
2895 
2896 ! Declare Arguments
2897  REAL (rprec) :: siesta_get_p_flux
2898  TYPE (siesta_class), INTENT(in) :: this
2899  REAL (rprec), DIMENSION(3), INTENT(in) :: flux
2900  LOGICAL, INTENT(in) :: normalize
2901 
2902 ! local variables
2903  REAL (rprec) :: ds
2904  REAL (rprec) :: wlow
2905  REAL (rprec) :: whigh
2906  INTEGER :: ilow
2907  INTEGER :: ihigh
2908  INTEGER :: mpol
2909  INTEGER :: ntor
2910  INTEGER :: m
2911  INTEGER :: n
2912  REAL (rprec) :: start_time
2913 
2914 ! Start of executable code
2915  start_time = profiler_get_start_time()
2916 
2917  ds = 1.0/(this%context%ns - 1.0)
2918 
2919  IF (flux(1) .lt. ds/2.0) THEN
2920  ilow = 2
2921  ihigh = 3
2922  whigh = flux(1) - ds/2.0
2923  wlow = 1.0 - whigh
2924  ELSE IF (flux(1) .gt. 1.0 - ds/2.0) THEN
2925  ihigh = this%context%ns
2926  ilow = ihigh - 1
2927  wlow = -flux(1)/ds - 0.5 + 1.0/ds
2928  whigh = 1.0 - wlow
2929  ELSE
2930  ilow = floor(flux(1)/ds + 1.5)
2931  ihigh = ilow + 1
2932  whigh = flux(1)/ds + 1.5 - ilow
2933  wlow = 1.0 - whigh
2934  END IF
2935 
2936  mpol = this%context%mpol
2937  ntor = this%context%ntor
2938 
2939  siesta_get_p_flux = 0.0
2940  DO n = -ntor, ntor
2941  DO m = 0, mpol
2942  siesta_get_p_flux = siesta_get_p_flux + &
2943  & (wlow*this%context%pmnch(m,n,ilow) + &
2944  & whigh*this%context%pmnch(m,n,ihigh))* &
2945  & cos(m*flux(2) + nfp*n*flux(3))
2946 
2947  IF (this%context%l_asym) THEN
2948  siesta_get_p_flux = siesta_get_p_flux + &
2949  & (wlow*this%context%pmnsh(m,n,ilow) + &
2950  & whigh*this%context%pmnsh(m,n,ihigh))* &
2951  & sin(m*flux(2) + nfp*n*flux(3))
2952  END IF
2953  END DO
2954  END DO
2955 
2956 ! The normalizes pressure is requested, scale the pressure to the min and max
2957 ! pressure. Otherwise apply constants to comvert to pascals.
2958  IF (normalize) THEN
2959 ! Ensure that the pressure is not outside the bounds of pmin and pmax.
2960  siesta_get_p_flux = min(siesta_get_p_flux, this%context%p_max)
2961  siesta_get_p_flux = max(siesta_get_p_flux, this%context%p_min)
2962 
2963  siesta_get_p_flux = (siesta_get_p_flux - this%context%p_min) &
2964  & / (this%context%p_max - this%context%p_min)
2965  ELSE
2966  siesta_get_p_flux = siesta_get_p_flux &
2967  & / (mu0*this%context%p_factor)
2968  END IF
2969 
2970  CALL profiler_set_stop_time('siesta_get_p_flux', start_time)
2971 
2972  END FUNCTION
2973 
2974 !-------------------------------------------------------------------------------
2984 !-------------------------------------------------------------------------------
2985  FUNCTION siesta_get_b_vec(this, x_cart, cyl)
2986  USE coordinate_utilities, ONLY : cyl_to_cart_vec, cart_to_cyl
2987  USE mgrid_mod
2988 
2989  IMPLICIT NONE
2990 
2991 ! Declare Arguments
2992  REAL (rprec), DIMENSION(3) :: siesta_get_b_vec
2993  TYPE (siesta_class), INTENT(in) :: this
2994  REAL (rprec), DIMENSION(3), INTENT(in) :: x_cart
2995  LOGICAL, INTENT(in) :: cyl
2996 
2997 ! local variables
2998  REAL (rprec), DIMENSION(3) :: r_cyl
2999  REAL (rprec), DIMENSION(3) :: flux
3000  REAL (rprec) :: start_time
3001 
3002 ! Start of executable code
3003  start_time = profiler_get_start_time()
3004 
3005 ! Determine if the point is inside or outside the last closed flux surface.
3006  flux = siesta_get_suv(this, x_cart)
3007 
3008 ! There is no need to shift this.
3009  r_cyl = cart_to_cyl(x_cart)
3010 
3011  siesta_get_b_vec = 0.0
3012  IF (flux(1) .gt. 1.0) THEN
3013  IF (lfreeb) THEN
3014  siesta_get_b_vec = vmec_get_b_vac(this%vmec, r_cyl) &
3015  & + siesta_get_ext_b_plasma(this, r_cyl, &
3016  & .false.)
3017  END IF
3018  ELSE
3019  siesta_get_b_vec = siesta_get_int_b_plasma(this, flux)
3020  END IF
3021 
3022  IF (.not.cyl) THEN
3023  siesta_get_b_vec = cyl_to_cart_vec(r_cyl, siesta_get_b_vec)
3024  END IF
3025 
3026  CALL profiler_set_stop_time('siesta_get_B_vec', start_time)
3027 
3028  END FUNCTION
3029 
3030 !-------------------------------------------------------------------------------
3040 !-------------------------------------------------------------------------------
3041  FUNCTION siesta_get_int_b_dphi(this, r, theta)
3042  USE line_segment, only: line_seg
3043  USE stel_constants, only: twopi
3044  USE read_wout_mod, only: bsubvmnc, ns
3045 
3046  IMPLICIT NONE
3047 
3048 ! Declare Arguments
3049  REAL (rprec) :: siesta_get_int_b_dphi
3050  TYPE (siesta_class), INTENT(in) :: this
3051  REAL (rprec), INTENT(in) :: r
3052  REAL (rprec), INTENT(in) :: theta
3053 
3054 ! local variables
3055  REAL (rprec) :: bsubv00c, ds
3056  INTEGER :: i
3057  REAL (rprec), DIMENSION(2) :: s
3058  REAL (rprec) :: start_time
3059 
3060 ! Start of executable code
3061  start_time = profiler_get_start_time()
3062  ds = 1.0/(ns - 1)
3063 
3064  s(1) = 1.0 - 1.5*ds
3065  s(2) = 1.0 - 0.5*ds
3066 
3067  CALL line_seg(r, bsubv00c, &
3068  & s, this%context%bsubvmnch(0,0,ns - 1:ns), 2)
3069 
3070  siesta_get_int_b_dphi = twopi*bsubv00c/this%context%b_factor
3071 
3072  CALL profiler_set_stop_time('siesta_get_Int_B_dphi', start_time)
3073 
3074  END FUNCTION
3075 
3076 !-------------------------------------------------------------------------------
3086 !-------------------------------------------------------------------------------
3087  FUNCTION siesta_get_suv(this, x_cart)
3088 
3089  IMPLICIT NONE
3090 
3091 ! Declare Arguments
3092  REAL (rprec), DIMENSION(3) :: siesta_get_suv
3093  TYPE (siesta_class), INTENT(in) :: this
3094  REAL (rprec), DIMENSION(3), INTENT(in) :: x_cart
3095 
3096 ! local variables
3097  REAL (rprec), DIMENSION(3) :: flux
3098  REAL (rprec) :: start_time
3099 
3100 ! Start of executable code
3101  start_time = profiler_get_start_time()
3102 
3103  flux = vmec_get_suv(this%vmec, x_cart)
3104 
3105  siesta_get_suv = &
3106  & (/ siesta_to_siesta_s(flux(1)), flux(2), flux(3) /)
3107 
3108  CALL profiler_set_stop_time('siesta_get_suv', start_time)
3109 
3110  END FUNCTION
3111 
3112 !-------------------------------------------------------------------------------
3126 !-------------------------------------------------------------------------------
3127  FUNCTION siesta_get_plasma_edge(this, phi, r, z)
3128 
3129  IMPLICIT NONE
3130 
3131 ! Declare Arguments
3132  INTEGER :: siesta_get_plasma_edge
3133  TYPE (siesta_class), INTENT(in) :: this
3134  REAL (rprec), INTENT (in) :: phi
3135  REAL (rprec), DIMENSION(:), POINTER :: r
3136  REAL (rprec), DIMENSION(:), POINTER :: z
3137 
3138 ! local variables
3139  REAL (rprec) :: start_time
3140 
3141 ! Start of executable code
3142  start_time = profiler_get_start_time()
3143 
3144 ! The inital VMEC defines the boundary.
3145 ! FIXME: Won't be the case for free boundary siesta.
3146  siesta_get_plasma_edge = vmec_get_plasma_edge(this%vmec, phi, r, &
3147  & z)
3148 
3149  CALL profiler_set_stop_time('siesta_get_plasma_edge', start_time)
3150 
3151  END FUNCTION
3152 
3153 !-------------------------------------------------------------------------------
3161 !-------------------------------------------------------------------------------
3162  FUNCTION siesta_get_magnetic_volume_rgrid(this)
3163 
3164  IMPLICIT NONE
3165 
3166 ! Declare Arguments
3167  REAL (rprec), DIMENSION(:,:,:), POINTER :: &
3168  & siesta_get_magnetic_volume_rgrid
3169  TYPE (siesta_class), INTENT(in) :: this
3170 
3171 ! local variables
3172  REAL (rprec) :: start_time
3173 
3174 ! Start of executable code
3175  start_time = profiler_get_start_time()
3176  siesta_get_magnetic_volume_rgrid => this%magnetic_cache%rsuv
3177 
3178  CALL profiler_set_stop_time('siesta_get_magnetic_volume_rgrid', &
3179  & start_time)
3180 
3181  END FUNCTION
3182 
3183 !-------------------------------------------------------------------------------
3191 !-------------------------------------------------------------------------------
3192  FUNCTION siesta_get_magnetic_volume_zgrid(this)
3193 
3194  IMPLICIT NONE
3195 
3196 ! Declare Arguments
3197  REAL (rprec), DIMENSION(:,:,:), POINTER :: &
3198  & siesta_get_magnetic_volume_zgrid
3199  TYPE (siesta_class), INTENT(in) :: this
3200 
3201 ! local variables
3202  REAL (rprec) :: start_time
3203 
3204 ! Start of executable code
3205  start_time = profiler_get_start_time()
3206  siesta_get_magnetic_volume_zgrid => this%magnetic_cache%zsuv
3207 
3208  CALL profiler_set_stop_time('siesta_get_magnetic_volume_zgrid', &
3209  & start_time)
3210 
3211  END FUNCTION
3212 
3213 !-------------------------------------------------------------------------------
3221 !-------------------------------------------------------------------------------
3222  FUNCTION siesta_get_magnetic_volume_jrgrid(this)
3223 
3224  IMPLICIT NONE
3225 
3226 ! Declare Arguments
3227  REAL (rprec), DIMENSION(:,:,:), POINTER :: &
3228  & siesta_get_magnetic_volume_jrgrid
3229  TYPE (siesta_class), INTENT(in) :: this
3230 
3231 ! local variables
3232  REAL (rprec) :: start_time
3233 
3234 ! Start of executable code
3235  start_time = profiler_get_start_time()
3236  siesta_get_magnetic_volume_jrgrid => this%magnetic_cache%jrsuv
3237 
3238  CALL profiler_set_stop_time('siesta_get_magnetic_volume_jrgrid', &
3239  & start_time)
3240 
3241  END FUNCTION
3242 
3243 !-------------------------------------------------------------------------------
3251 !-------------------------------------------------------------------------------
3252  FUNCTION siesta_get_magnetic_volume_jphigrid(this)
3253 
3254  IMPLICIT NONE
3255 
3256 ! Declare Arguments
3257  REAL (rprec), DIMENSION(:,:,:), POINTER :: &
3258  & siesta_get_magnetic_volume_jphigrid
3259  TYPE (siesta_class), INTENT(in) :: this
3260 
3261 ! local variables
3262  REAL (rprec) :: start_time
3263 
3264 ! Start of executable code
3265  start_time = profiler_get_start_time()
3266  siesta_get_magnetic_volume_jphigrid => this%magnetic_cache%jphisuv
3267 
3268  CALL profiler_set_stop_time('siesta_get_magnetic_volume_jphigrid', &
3269  & start_time)
3270 
3271  END FUNCTION
3272 
3273 !-------------------------------------------------------------------------------
3281 !-------------------------------------------------------------------------------
3282  FUNCTION siesta_get_magnetic_volume_jzgrid(this)
3283 
3284  IMPLICIT NONE
3285 
3286 ! Declare Arguments
3287  REAL (rprec), DIMENSION(:,:,:), POINTER :: &
3288  & siesta_get_magnetic_volume_jzgrid
3289  TYPE (siesta_class), INTENT(in) :: this
3290 
3291 ! local variables
3292  REAL (rprec) :: start_time
3293 
3294 ! Start of executable code
3295  start_time = profiler_get_start_time()
3296  siesta_get_magnetic_volume_jzgrid => this%magnetic_cache%jzsuv
3297 
3298  CALL profiler_set_stop_time('siesta_get_magnetic_volume_jzgrid', &
3299  & start_time)
3300 
3301  END FUNCTION
3302 
3303 !-------------------------------------------------------------------------------
3317 !-------------------------------------------------------------------------------
3318  FUNCTION siesta_get_volume_int_element(this)
3319  USE read_wout_mod, only: isigng
3320 
3321  IMPLICIT NONE
3322 
3323 ! Declare Arguments
3324  REAL (rprec) :: siesta_get_volume_int_element
3325  TYPE (siesta_class), INTENT(in) :: this
3326 
3327 ! local variables
3328  REAL (rprec) :: start_time
3329 
3330 ! Start of executable code
3331  start_time = profiler_get_start_time()
3332  siesta_get_volume_int_element = isigng*this%magnetic_cache%ds &
3333  & * this%magnetic_cache%du &
3334  & * this%magnetic_cache%dv
3335 
3336  CALL profiler_set_stop_time('siesta_get_volume_int_element', &
3337  & start_time)
3338 
3339  END FUNCTION
3340 
3341 !-------------------------------------------------------------------------------
3349 !-------------------------------------------------------------------------------
3350  FUNCTION siesta_get_con_surface_krgrid(this)
3351 
3352  IMPLICIT NONE
3353 
3354 ! Declare Arguments
3355  REAL (rprec), DIMENSION(:,:), POINTER :: &
3356  & siesta_get_con_surface_krgrid
3357  TYPE (siesta_class), INTENT(in) :: this
3358 
3359 ! local variables
3360  REAL (rprec) :: start_time
3361 
3362 ! Start of executable code
3363  start_time = profiler_get_start_time()
3364  siesta_get_con_surface_krgrid => this%magnetic_cache%kruv
3365 
3366  CALL profiler_set_stop_time('siesta_get_con_surface_krgrid', &
3367  & start_time)
3368 
3369  END FUNCTION
3370 
3371 !-------------------------------------------------------------------------------
3379 !-------------------------------------------------------------------------------
3380  FUNCTION siesta_get_con_surface_kphigrid(this)
3381 
3382  IMPLICIT NONE
3383 
3384 ! Declare Arguments
3385  REAL (rprec), DIMENSION(:,:), POINTER :: &
3386  & siesta_get_con_surface_kphigrid
3387  TYPE (siesta_class), INTENT(in) :: this
3388 
3389 ! local variables
3390  REAL (rprec) :: start_time
3391 
3392 ! Start of executable code
3393  start_time = profiler_get_start_time()
3394  siesta_get_con_surface_kphigrid => this%magnetic_cache%kphiuv
3395 
3396  CALL profiler_set_stop_time('siesta_get_con_surface_kphigrid', &
3397  & start_time)
3398 
3399  END FUNCTION
3400 
3401 !-------------------------------------------------------------------------------
3409 !-------------------------------------------------------------------------------
3410  FUNCTION siesta_get_con_surface_kzgrid(this)
3411 
3412  IMPLICIT NONE
3413 
3414 ! Declare Arguments
3415  REAL (rprec), DIMENSION(:,:), POINTER :: &
3416  & siesta_get_con_surface_kzgrid
3417  TYPE (siesta_class), INTENT(in) :: this
3418 
3419 ! local variables
3420  REAL (rprec) :: start_time
3421 
3422 ! Start of executable code
3423  start_time = profiler_get_start_time()
3424  siesta_get_con_surface_kzgrid => this%magnetic_cache%kzuv
3425 
3426  CALL profiler_set_stop_time('siesta_get_con_surface_kzgrid', &
3427  & start_time)
3428 
3429  END FUNCTION
3430 
3431 !-------------------------------------------------------------------------------
3445 !-------------------------------------------------------------------------------
3446  FUNCTION siesta_get_area_int_element(this)
3447  USE read_wout_mod, only: isigng
3448 
3449  IMPLICIT NONE
3450 
3451 ! Declare Arguments
3452  REAL (rprec) :: siesta_get_area_int_element
3453  TYPE (siesta_class), INTENT(in) :: this
3454 
3455 ! local variables
3456  REAL (rprec) :: start_time
3457 
3458 ! Start of executable code
3459  start_time = profiler_get_start_time()
3460  siesta_get_area_int_element = this%magnetic_cache%du_a &
3461  & * this%magnetic_cache%dv_a*isigng
3462 
3463  CALL profiler_set_stop_time('siesta_get_area_int_element', &
3464  & start_time)
3465 
3466  END FUNCTION
3467 
3468 !-------------------------------------------------------------------------------
3480 !-------------------------------------------------------------------------------
3481  FUNCTION siesta_get_ext_currents(this, num_currents, &
3482  & scale_currents)
3483 
3484  IMPLICIT NONE
3485 
3486 ! Declare Arguments
3487  REAL (rprec), DIMENSION(:), POINTER :: siesta_get_ext_currents
3488  TYPE (siesta_class), INTENT(in) :: this
3489  INTEGER, INTENT(in) :: num_currents
3490  LOGICAL, INTENT(out) :: scale_currents
3491 
3492 ! local variables
3493  REAL (rprec) :: start_time
3494 
3495 ! Start of executable code
3496  start_time = profiler_get_start_time()
3497 
3498  siesta_get_ext_currents => &
3499  & vmec_get_ext_currents(this%vmec, num_currents, scale_currents)
3500 
3501  CALL profiler_set_stop_time('vmec_get_ext_currents', start_time)
3502 
3503  END FUNCTION
3504 
3505 !-------------------------------------------------------------------------------
3519 !-------------------------------------------------------------------------------
3520  FUNCTION siesta_get_ext_b_plasma(this, position, axi_only)
3521  USE coordinate_utilities
3522  USE read_wout_mod, only: isigng
3523 
3524  IMPLICIT NONE
3525 
3526 ! Declare Arguments
3527  REAL (rprec), DIMENSION(3) :: siesta_get_ext_b_plasma
3528  TYPE (siesta_class), INTENT(in) :: this
3529  REAL (rprec), DIMENSION(3), INTENT(in) :: position
3530  LOGICAL, INTENT(in) :: axi_only
3531 
3532 ! local variables
3533  REAL (rprec), DIMENSION(3) :: x_cart
3534  REAL (rprec), DIMENSION(3) :: r_vec
3535  INTEGER :: u, v
3536  REAL (rprec) :: length
3537  REAL (rprec), DIMENSION(:,:,:), POINTER :: x_prime
3538  REAL (rprec), DIMENSION(:,:), POINTER :: kxuv
3539  REAL (rprec), DIMENSION(:,:), POINTER :: kyuv
3540  REAL (rprec), DIMENSION(:,:), POINTER :: kzuv
3541  REAL (rprec), DIMENSION(3) :: sum_b
3542  REAL (rprec) :: start_time
3543 
3544 ! Start of executable code
3545  start_time = profiler_get_start_time()
3546 
3547 ! Convert the position to cartesian coordinates. Rotate and shift the signal
3548 ! position with respect to the fixed plasma.
3549  x_cart = cyl_to_cart(position - (/ 0.0_dp, &
3550  & this%phi_offset, &
3551  & this%z_offset /))
3552 
3553  IF (axi_only) THEN
3554  CALL assert(ASSOCIATED(this%magnetic_cache%kxuv_axi), &
3555  & 'Cannot get the axisymmtric only ext bfield.')
3556 
3557  x_prime => this%magnetic_cache%x_axi
3558  kxuv => this%magnetic_cache%kxuv_axi
3559  kyuv => this%magnetic_cache%kyuv_axi
3560  kzuv => this%magnetic_cache%kzuv_axi
3561  ELSE
3562  CALL assert(ASSOCIATED(this%magnetic_cache%kxuv_full), &
3563  & 'Cannot get the ext bfield.')
3564 
3565  x_prime => this%magnetic_cache%x_prime
3566  kxuv => this%magnetic_cache%kxuv_full
3567  kyuv => this%magnetic_cache%kyuv_full
3568  kzuv => this%magnetic_cache%kzuv_full
3569  END IF
3570 
3571 ! The magnetic field produced by a sheet of current is
3572 !
3573 ! B(x) = Int K x (x-x')/|x-x'|^3 da (1)
3574  sum_b = 0.0
3575  DO v = 1, SIZE(kxuv,2)
3576  DO u = 1, SIZE(kxuv,1)
3577  r_vec = x_cart - x_prime(u,v,:)
3578  length = sqrt(dot_product(r_vec, r_vec))**3.0
3579 
3580  sum_b(1) = sum_b(1) &
3581  & + (kyuv(u,v)*r_vec(3) - kzuv(u,v)*r_vec(2)) &
3582  & / length
3583  sum_b(2) = sum_b(2) &
3584  & + (kzuv(u,v)*r_vec(1) - kxuv(u,v)*r_vec(3)) &
3585  & / length
3586  sum_b(3) = sum_b(3) &
3587  & + (kxuv(u,v)*r_vec(2) - kyuv(u,v)*r_vec(1)) &
3588  & / length
3589  END DO
3590  END DO
3591 
3592 ! The virtual casing principal acts as if there is a perfectly conducting shell
3593 ! around the last closed flux surface. In this shell, eddy currents induced
3594 ! from the volume currents get canceled out by the shell currents.
3595 !
3596 ! B_v + B_s = 0 (2)
3597 !
3598 ! Thus the shell field is the same magnetude but opposite sign to the volume
3599 ! current.
3600 !
3601 ! B_v = -B_s (3)
3602 !
3603 ! After the surface integral is completed, negate the result to produce a field
3604 ! that would be the same as a volume integral.
3605  siesta_get_ext_b_plasma = -mu0/(4.0*pi)*sum_b*isigng &
3606  & * this%magnetic_cache%du_full &
3607  & * this%magnetic_cache%dv_full
3608 
3609  siesta_get_ext_b_plasma = cart_to_cyl_vec(x_cart, &
3610  & siesta_get_ext_b_plasma)
3611 
3612  CALL profiler_set_stop_time('siesta_get_ext_b_plasma', start_time)
3613 
3614  END FUNCTION
3615 
3616 !-------------------------------------------------------------------------------
3627 !-------------------------------------------------------------------------------
3628  FUNCTION siesta_get_int_b_plasma(this, flux)
3629  USE read_wout_mod, only: nfp, mnmax, xm, xn, rmnc, rmns, zmnc, &
3630  & zmns, lasym, nsvmec => ns
3631 
3632  USE vmec_utils, ONLY : getbcyl_wout
3633 
3634  IMPLICIT NONE
3635 
3636 ! Declare Arguments
3637  REAL (rprec), DIMENSION(3) :: siesta_get_int_b_plasma
3638  TYPE (siesta_class), INTENT(in) :: this
3639  REAL (rprec), DIMENSION(3), INTENT(in) :: flux
3640 
3641 ! local variables
3642  REAL (rprec) :: ds
3643  REAL (rprec) :: wlow_v
3644  REAL (rprec) :: whigh_v
3645  INTEGER :: ilow_v
3646  INTEGER :: ihigh_v
3647  REAL (rprec) :: wlow_s
3648  REAL (rprec) :: whigh_s
3649  INTEGER :: ilow_s
3650  INTEGER :: ihigh_s
3651  INTEGER :: mpol
3652  INTEGER :: ntor
3653  INTEGER :: m
3654  INTEGER :: n
3655  REAL (rprec) :: r
3656  REAL (rprec) :: rs
3657  REAL (rprec) :: ru
3658  REAL (rprec) :: rv
3659  REAL (rprec) :: zu
3660  REAL (rprec) :: zs
3661  REAL (rprec) :: zv
3662  REAL (rprec) :: bs
3663  REAL (rprec) :: bu
3664  REAL (rprec) :: bv
3665  REAL (rprec), DIMENSION(:), ALLOCATABLE :: vmecint
3666  REAL (rprec), DIMENSION(:), ALLOCATABLE :: cosmn_vmec
3667  REAL (rprec), DIMENSION(:), ALLOCATABLE :: sinmn_vmec
3668  REAL (rprec) :: start_time
3669 
3670 ! Start of executable code
3671  start_time = profiler_get_start_time()
3672 
3673  ds = 1.0/(this%context%ns - 1.0)
3674 
3675 ! siesta intepolation parameters.
3676  IF (flux(1) .lt. ds/2.0) THEN
3677  ilow_s = 2
3678  ihigh_s = 3
3679  whigh_s = flux(1) - ds/2.0
3680  wlow_s = 1.0 - whigh_s
3681  ELSE IF (flux(1) .gt. 1.0 - ds/2.0) THEN
3682  ihigh_s = this%context%ns
3683  ilow_s = ihigh_s - 1
3684  wlow_s = -flux(1)/ds - 0.5 + 1.0/ds
3685  whigh_s = 1.0 - wlow_s
3686  ELSE
3687  ilow_s = floor(flux(1)/ds + 1.5)
3688  ihigh_s = ilow_s + 1
3689  whigh_s = flux(1)/ds + 1.5 - ilow_s
3690  wlow_s = 1.0 - whigh_s
3691  END IF
3692 
3693  ilow_v = siesta_to_i_vmec_low(flux(1))
3694  ihigh_v = ilow_v + 1
3695  wlow_v = siesta_to_w_vmec_low(flux(1), ilow_v)
3696  whigh_v = 1.0 - wlow_v
3697 
3698  ALLOCATE(vmecint(mnmax))
3699  ALLOCATE(cosmn_vmec(mnmax))
3700  ALLOCATE(sinmn_vmec(mnmax))
3701 
3702  cosmn_vmec = cos(xm*flux(2) - xn*flux(3))
3703  sinmn_vmec = sin(xm*flux(2) - xn*flux(3))
3704 
3705  vmecint = wlow_v*rmnc(:,ilow_v) + whigh_v*rmnc(:,ihigh_v)
3706  r = sum(vmecint*cosmn_vmec)
3707  ru = -sum(xm*vmecint*sinmn_vmec)
3708  rv = sum(xn*vmecint*sinmn_vmec)
3709 
3710  vmecint = 2.0*flux(1)*(rmnc(:,ihigh_v) - rmnc(:,ilow_v)) &
3711  & * (nsvmec - 1.0)
3712  rs = sum(vmecint*cosmn_vmec)
3713 
3714  vmecint = wlow_v*zmns(:,ilow_v) + whigh_v*zmns(:,ihigh_v)
3715  zu = sum(xm*vmecint*cosmn_vmec)
3716  zv = -sum(xn*vmecint*cosmn_vmec)
3717 
3718  vmecint = 2.0*flux(1)*(zmns(:,ihigh_v) - zmns(:,ilow_v)) &
3719  & * (nsvmec - 1.0)
3720  zs = sum(vmecint*sinmn_vmec)
3721 
3722  mpol = this%context%mpol
3723  ntor = this%context%ntor
3724 
3725  bs = 0.0
3726  bu = 0.0
3727  bv = 0.0
3728 
3729  DO n = -ntor, ntor
3730  DO m = 0, mpol
3731  bs = bs + &
3732  & (wlow_s*this%context%bsupsmnsh(m,n,ilow_s) + &
3733  & whigh_s*this%context%bsupsmnsh(m,n,ihigh_s))* &
3734  & sin(m*flux(2) + nfp*n*flux(3))
3735 
3736  bu = bu + &
3737  & (wlow_s*this%context%bsupumnch(m,n,ilow_s) + &
3738  & whigh_s*this%context%bsupumnch(m,n,ihigh_s))* &
3739  & cos(m*flux(2) + nfp*n*flux(3))
3740 
3741  bv = bv + &
3742  & (wlow_s*this%context%bsupvmnch(m,n,ilow_s) + &
3743  & whigh_s*this%context%bsupvmnch(m,n,ihigh_s))* &
3744  & cos(m*flux(2) + nfp*n*flux(3))
3745  END DO
3746  END DO
3747 
3748  IF (this%context%l_asym) THEN
3749 
3750  IF (lasym) THEN
3751  vmecint = wlow_v*rmns(:,ilow_v) + whigh_v*rmns(:,ihigh_v)
3752  r = r + sum(vmecint*sinmn_vmec)
3753  ru = ru + sum(xm*vmecint*cosmn_vmec)
3754  rv = rv - sum(xn*vmecint*cosmn_vmec)
3755 
3756  vmecint = 2.0*flux(1)*(rmns(:,ihigh_v) - rmns(:,ilow_v)) &
3757  & * (nsvmec - 1.0)
3758  rs = rs + sum(vmecint*sinmn_vmec)
3759 
3760  vmecint = wlow_v*zmnc(:,ilow_v) + whigh_v*zmnc(:,ihigh_v)
3761  zu = zu - sum(xm*vmecint*sinmn_vmec)
3762  zv = zv + sum(xn*vmecint*sinmn_vmec)
3763 
3764  vmecint = 2.0*flux(1)*(zmnc(:,ihigh_v) - zmnc(:,ilow_v)) &
3765  & * (nsvmec - 1.0)
3766  zs = zs + sum(vmecint*cosmn_vmec)
3767  END IF
3768 
3769  DO n = -ntor, ntor
3770  DO m = 0, mpol
3771  bs = bs + &
3772  & (wlow_s*this%context%bsupsmnch(m,n,ilow_s) + &
3773  & whigh_s*this%context%bsupsmnch(m,n,ihigh_s))* &
3774  & cos(m*flux(2) + nfp*n*flux(3))
3775 
3776  bu = bu + &
3777  & (wlow_s*this%context%bsupumnsh(m,n,ilow_s) + &
3778  & whigh_s*this%context%bsupumnsh(m,n,ihigh_s))* &
3779  & sin(m*flux(2) + nfp*n*flux(3))
3780 
3781  bv = bv + &
3782  & (wlow_s*this%context%bsupvmnsh(m,n,ilow_s) + &
3783  & whigh_s*this%context%bsupvmnsh(m,n,ihigh_s))* &
3784  & sin(m*flux(2) + nfp*n*flux(3))
3785  END DO
3786  END DO
3787  END IF
3788 
3789  DEALLOCATE(vmecint)
3790  DEALLOCATE(cosmn_vmec)
3791  DEALLOCATE(sinmn_vmec)
3792 
3793  siesta_get_int_b_plasma(1) = bs*rs + bu*ru + bv*rv
3794  siesta_get_int_b_plasma(2) = bv*r
3795  siesta_get_int_b_plasma(3) = bs*zs + bu*zu + bv*zv
3796 
3797  siesta_get_int_b_plasma = siesta_get_int_b_plasma &
3798  & / this%context%b_factor
3799 
3800  END FUNCTION
3801 
3802 !-------------------------------------------------------------------------------
3809 !-------------------------------------------------------------------------------
3810  FUNCTION siesta_get_grid_size(this)
3811 
3812  IMPLICIT NONE
3813 
3814 ! Declare Arguments
3815  INTEGER :: siesta_get_grid_size
3816  TYPE (siesta_class), INTENT(in) :: this
3817 
3818 ! local variables
3819  REAL (rprec) :: start_time
3820 
3821 ! Start of executable code
3822  start_time = profiler_get_start_time()
3823 
3824  IF (ASSOCIATED(this%context)) THEN
3825  siesta_get_grid_size = this%context%ns
3826  ELSE
3827  siesta_get_grid_size = nsin
3828  END IF
3829 
3830  CALL profiler_set_stop_time('siesta_get_grid_size', start_time)
3831 
3832  END FUNCTION
3833 
3834 !-------------------------------------------------------------------------------
3841 !-------------------------------------------------------------------------------
3842  FUNCTION siesta_get_grid_start(this)
3843 
3844  IMPLICIT NONE
3845 
3846 ! Declare Arguments
3847  REAl (rprec) :: siesta_get_grid_start
3848  TYPE (siesta_class), INTENT(in) :: this
3849 
3850 ! local variables
3851  REAL (rprec) :: start_time
3852 
3853 ! Start of executable code
3854  start_time = profiler_get_start_time()
3855  siesta_get_grid_start = 0.0
3856  CALL profiler_set_stop_time('siesta_get_grid_start', start_time)
3857 
3858  END FUNCTION
3859 
3860 !-------------------------------------------------------------------------------
3867 !-------------------------------------------------------------------------------
3868  FUNCTION siesta_get_grid_stop(this)
3869 
3870  IMPLICIT NONE
3871 
3872 ! Declare Arguments
3873  REAL (rprec) :: siesta_get_grid_stop
3874  TYPE (siesta_class), INTENT(in) :: this
3875 
3876 ! local variables
3877  REAL (rprec) :: start_time
3878 
3879 ! Start of executable code
3880  start_time = profiler_get_start_time()
3881  siesta_get_grid_stop = 1.0
3882  CALL profiler_set_stop_time('siesta_get_grid_stop', start_time)
3883 
3884  END FUNCTION
3885 
3886 !*******************************************************************************
3887 ! QUERY SUBROUTINES
3888 !*******************************************************************************
3889 !-------------------------------------------------------------------------------
3897 !-------------------------------------------------------------------------------
3898  FUNCTION siesta_is_scaler_value(this, id)
3899 
3900  IMPLICIT NONE
3901 
3902 ! Declare Arguments
3903  LOGICAL :: siesta_is_scaler_value
3904  TYPE (siesta_class), INTENT(in) :: this
3905  INTEGER, INTENT(in) :: id
3906 
3907 ! local variables
3908  REAL (rprec) :: start_time
3909 
3910 ! Start of executable code
3911  start_time = profiler_get_start_time()
3912  SELECT CASE (id)
3913 
3914  CASE (siesta_phi_offset_id, siesta_z_offset_id)
3915  siesta_is_scaler_value = .true.
3916 
3917  CASE DEFAULT
3918  siesta_is_scaler_value = vmec_is_scaler_value(this%vmec, id)
3919 
3920  END SELECT
3921 
3922  CALL profiler_set_stop_time('siesta_is_scaler_value', start_time)
3923 
3924  END FUNCTION
3925 
3926 !-------------------------------------------------------------------------------
3934 !-------------------------------------------------------------------------------
3935  FUNCTION siesta_is_1d_array(this, id)
3936 
3937  IMPLICIT NONE
3938 
3939 ! Declare Arguments
3940  LOGICAL :: siesta_is_1d_array
3941  TYPE (siesta_class), INTENT(in) :: this
3942  INTEGER, INTENT(in) :: id
3943 
3944 ! local variables
3945  REAL (rprec) :: start_time
3946 
3947 ! Start of executable code
3948  start_time = profiler_get_start_time()
3949  SELECT CASE (id)
3950 
3951  CASE (siesta_pp_ne_b_id, siesta_pp_ne_as_id, &
3952  & siesta_pp_ne_af_id, siesta_pp_te_b_id, &
3953  & siesta_pp_te_as_id, siesta_pp_te_af_id, &
3954  & siesta_pp_ti_b_id, siesta_pp_ti_as_id, &
3955  & siesta_pp_ti_af_id, siesta_helpert_id)
3956  siesta_is_1d_array = .true.
3957 
3958  CASE DEFAULT
3959  siesta_is_1d_array = vmec_is_1d_array(this%vmec, id)
3960 
3961  END SELECT
3962 
3963  CALL profiler_set_stop_time('siesta_is_1d_array', start_time)
3964 
3965  END FUNCTION
3966 
3967 !-------------------------------------------------------------------------------
3975 !-------------------------------------------------------------------------------
3976  FUNCTION siesta_is_2d_array(this, id)
3977 
3978  IMPLICIT NONE
3979 
3980 ! Declare Arguments
3981  LOGICAL :: siesta_is_2d_array
3982  TYPE (siesta_class), INTENT(in) :: this
3983  INTEGER, INTENT(in) :: id
3984 
3985 ! local variables
3986  REAL (rprec) :: start_time
3987 
3988 ! Start of executable code
3989  start_time = profiler_get_start_time()
3990  SELECT CASE (id)
3991 
3992  CASE (siesta_pp_sxrem_b_id, siesta_pp_sxrem_as_id, &
3993  & siesta_pp_sxrem_af_id)
3994  siesta_is_2d_array = .true.
3995 
3996  CASE DEFAULT
3997  siesta_is_2d_array = vmec_is_2d_array(this%vmec, id)
3998 
3999  END SELECT
4000 
4001  CALL profiler_set_stop_time('siesta_is_2d_array', start_time)
4002 
4003  END FUNCTION
4004 
4005 !-------------------------------------------------------------------------------
4014 !-------------------------------------------------------------------------------
4015  FUNCTION siesta_is_recon_param(this, id)
4016 
4017  IMPLICIT NONE
4018 
4019 ! Declare Arguments
4020  LOGICAL :: siesta_is_recon_param
4021  TYPE (siesta_class), INTENT(in) :: this
4022  INTEGER, INTENT(in) :: id
4023 
4024 ! local variables
4025  REAL (rprec) :: start_time
4026 
4027 ! Start of executable code
4028  start_time = profiler_get_start_time()
4029  SELECT CASE (id)
4030 
4031  CASE (siesta_pp_ne_b_id, siesta_pp_ne_as_id, &
4032  & siesta_pp_ne_af_id, siesta_pp_te_b_id, &
4033  & siesta_pp_te_as_id, siesta_pp_te_af_id, &
4034  & siesta_pp_ti_b_id, siesta_pp_ti_as_id, &
4035  & siesta_pp_ti_af_id, siesta_pp_sxrem_b_id, &
4036  & siesta_pp_sxrem_as_id, siesta_pp_sxrem_af_id, &
4037  & siesta_phi_offset_id, siesta_z_offset_id, &
4038  & siesta_helpert_id)
4039  siesta_is_recon_param = .true.
4040 
4041  CASE DEFAULT
4042  siesta_is_recon_param = &
4043  & vmec_is_recon_param(this%vmec, id)
4044 
4045  END SELECT
4046 
4047  CALL profiler_set_stop_time('siesta_is_recon_param', start_time)
4048 
4049  END FUNCTION
4050 
4051 !-------------------------------------------------------------------------------
4058 !-------------------------------------------------------------------------------
4059  FUNCTION siesta_is_using_point(this)
4060 
4061  IMPLICIT NONE
4062 
4063 ! Declare Arguments
4064  LOGICAL :: siesta_is_using_point
4065  TYPE (siesta_class), INTENT(in) :: this
4066 
4067 ! local variables
4068  REAL (rprec) :: start_time
4069 
4070 ! Start of executable code
4071  start_time = profiler_get_start_time()
4072 
4073 ! Some compilers will generate code that will early terminate when a first
4074 ! operand to a .and. operator evaluates to false. Others are not so clever.
4075 ! Avoid running the second ASSOCIATED check if the magnetic cache is null.
4076  siesta_is_using_point = ASSOCIATED(this%magnetic_cache)
4077  IF (siesta_is_using_point) THEN
4078  siesta_is_using_point = &
4079  & ASSOCIATED(this%magnetic_cache%kxuv_full)
4080  END IF
4081 
4082  CALL profiler_set_stop_time('vmec_is_using_point', start_time)
4083 
4084  END FUNCTION
4085 
4086 !*******************************************************************************
4087 ! UTILITY SUBROUTINES
4088 !*******************************************************************************
4089 !-------------------------------------------------------------------------------
4104 !-------------------------------------------------------------------------------
4105  FUNCTION siesta_converge(this, num_iter, iou, eq_comm, &
4106  & state_flags)
4107  USE model_state
4108  USE data_parameters
4109  USE safe_open_mod
4110  USE siesta_error
4111 
4112  IMPLICIT NONE
4113 
4114 ! Declare Arguments
4115  LOGICAL :: siesta_converge
4116  TYPE (siesta_class), INTENT(inout) :: this
4117  INTEGER, INTENT(inout) :: num_iter
4118  INTEGER, INTENT(in) :: iou
4119  INTEGER, INTENT(in) :: eq_comm
4120  INTEGER, INTENT(in) :: state_flags
4121 
4122 ! local variables
4123  INTEGER :: eq_size
4124  INTEGER :: eq_rank
4125  INTEGER :: iou_nli
4126  INTEGER :: status
4127  INTEGER :: child_comm
4128  TYPE (siesta_run_class), POINTER :: run_context
4129  REAL (rprec) :: start_time
4130 
4131 ! Start of executable code
4132  start_time = profiler_get_start_time()
4133 
4134  siesta_converge = .true.
4135  IF (btest(state_flags, model_state_vmec_flag)) THEN
4136  siesta_converge = vmec_converge(this%vmec, num_iter, iou, &
4137  & eq_comm)
4138  END IF
4139 
4140 ! FIXME: This is currently a hack. Eventuially siesta will be updated to have
4141 ! a run siesta routine like vmec. For now run vmec as an external
4142 ! command.
4143 
4144  eq_size = 1
4145  eq_rank = 0
4146 #if defined(MPI_OPT)
4147  CALL mpi_comm_rank(eq_comm, eq_rank, status)
4148  CALL mpi_comm_size(eq_comm, eq_size, status)
4149 #endif
4150  IF (eq_rank .eq. 0) THEN
4151 
4152  iou_nli = 100
4153  ladd_pert = .true.
4154  lresistive = .false.
4155  lrestart = .false.
4156  CALL safe_open(iou_nli, status, trim(this%siesta_file_name), &
4157  & 'replace', 'formatted', delim_in='quote')
4158  WRITE (iou_nli, nml=siesta_info)
4159  CLOSE (iou_nli, iostat=status)
4160  END IF
4161 
4162 #if defined(MPI_OPT)
4163  CALL mpi_barrier(eq_comm, status)
4164 #endif
4165  run_context => siesta_run_construct(eq_comm, .false.)
4166  CALL siesta_run_set_vmec(run_context)
4167  CALL siesta_run_converge(run_context)
4168  CALL siesta_run_destruct(run_context, .false., .false.)
4169 
4170  IF (siesta_error_state .eq. siesta_error_no_error) THEN
4171 
4172  IF (eq_rank .eq. 0) THEN
4173 
4174  ladd_pert = .false.
4175  lresistive = .false.
4176  lrestart = .true.
4177  CALL safe_open(iou_nli, status, trim(this%siesta_file_name), &
4178  & 'replace', 'formatted', delim_in='quote')
4179  WRITE (iou_nli, nml=siesta_info)
4180  CLOSE (iou_nli, iostat=status)
4181  END IF
4182 
4183 #if defined(MPI_OPT)
4184  CALL mpi_barrier(eq_comm, status)
4185 #endif
4186  run_context => siesta_run_construct(eq_comm, .false.)
4187  CALL siesta_run_set_vmec(run_context)
4188  CALL siesta_run_converge(run_context)
4189  CALL siesta_run_destruct(run_context, .false., .false.)
4190  END IF
4191  siesta_converge = siesta_error_state .eq. siesta_error_no_error
4192 
4193  status = 0
4194  siesta_converge = siesta_converge .and. status .eq. 0
4195 
4196  IF (siesta_converge .and. eq_rank .eq. 0) THEN
4197  IF (ASSOCIATED(this%context)) THEN
4198  CALL siesta_context_read(this%context, &
4199  & this%restart_file_name)
4200  ELSE
4201  this%context => &
4202  & siesta_context_construct(this%restart_file_name)
4203  END IF
4204 
4205  IF (ASSOCIATED(this%magnetic_cache)) THEN
4206  CALL siesta_set_magnetic_cache(this)
4207  END IF
4208  END IF
4209 
4210  CALL profiler_set_stop_time('siesta_converge', start_time)
4211 
4212  END FUNCTION
4213 
4214 !-------------------------------------------------------------------------------
4224 !-------------------------------------------------------------------------------
4225  SUBROUTINE siesta_read_vac_file(this, index, eq_comm)
4226 
4227  IMPLICIT NONE
4228 
4229 ! Declare Arguments
4230  TYPE (siesta_class), INTENT(in) :: this
4231  INTEGER, INTENT(in) :: index
4232  INTEGER, INTENT(in) :: eq_comm
4233 
4234 ! Start of executable code
4235  CALL vmec_read_vac_file(this%vmec, index, eq_comm)
4236 
4237  END SUBROUTINE
4238 
4239 !-------------------------------------------------------------------------------
4245 !-------------------------------------------------------------------------------
4246  SUBROUTINE siesta_save_state(this)
4247  USE file_opts
4248 
4249  IMPLICIT NONE
4250 
4251 ! Declare Arguments
4252  TYPE (siesta_class), INTENT(inout) :: this
4253 
4254 ! local variables
4255  INTEGER :: error
4256  INTEGER :: i
4257  REAL (rprec) :: start_time
4258 
4259 ! Start of executable code
4260  start_time = profiler_get_start_time()
4261 
4262  CALL vmec_save_state(this%vmec)
4263 
4264 ! Cache the restart file by appending _cache to a copy. Use copy here to keep
4265 ! the orginal file intact.
4266  CALL copy_file(trim(this%restart_file_name), &
4267  & trim(this%restart_file_name) // '_cache', &
4268  & error)
4269  CALL assert_eq(error, 0, 'Error copying restart file.')
4270 
4271  IF (ASSOCIATED(this%ne)) THEN
4272  CALL pprofile_save_state(this%ne)
4273  END IF
4274  IF (ASSOCIATED(this%te)) THEN
4275  CALL pprofile_save_state(this%te)
4276  END IF
4277  IF (ASSOCIATED(this%ti)) THEN
4278  CALL pprofile_save_state(this%ti)
4279  END IF
4280  IF (ASSOCIATED(this%sxrem)) THEN
4281  DO i = 1, SIZE(this%sxrem)
4282  IF (ASSOCIATED(this%sxrem(i)%p)) THEN
4283  CALL pprofile_save_state(this%sxrem(i)%p)
4284  END IF
4285  END DO
4286  END IF
4287 
4288  CALL profiler_set_stop_time('siesta_save_state', start_time)
4289 
4290  END SUBROUTINE
4291 
4292 !-------------------------------------------------------------------------------
4298 !-------------------------------------------------------------------------------
4299  SUBROUTINE siesta_reset_state(this)
4300  USE file_opts
4301 
4302  IMPLICIT NONE
4303 
4304 ! Declare Arguments
4305  TYPE (siesta_class), INTENT(inout) :: this
4306 
4307 ! local variables
4308  INTEGER :: error
4309  INTEGER :: i
4310  REAL (rprec) :: start_time
4311 
4312 ! Start of executable code
4313  start_time = profiler_get_start_time()
4314 
4315  CALL vmec_reset_state(this%vmec)
4316 
4317 ! Reset the restart file.
4318  CALL copy_file(trim(this%restart_file_name) // '_cache', &
4319  & trim(this%restart_file_name), error)
4320  CALL assert_eq(error, 0, 'Error moving wout file.')
4321 
4322  CALL siesta_context_read(this%context, this%restart_file_name)
4323 
4324  IF (ASSOCIATED(this%ne)) THEN
4325  CALL pprofile_reset_state(this%ne)
4326  END IF
4327  IF (ASSOCIATED(this%te)) THEN
4328  CALL pprofile_reset_state(this%te)
4329  END IF
4330  IF (ASSOCIATED(this%ti)) THEN
4331  CALL pprofile_reset_state(this%ti)
4332  END IF
4333  IF (ASSOCIATED(this%sxrem)) THEN
4334  DO i = 1, SIZE(this%sxrem)
4335  IF (ASSOCIATED(this%sxrem(i)%p)) THEN
4336  CALL pprofile_reset_state(this%sxrem(i)%p)
4337  END IF
4338  END DO
4339  END IF
4340 
4341  CALL profiler_set_stop_time('siesta_reset_state', start_time)
4342 
4343  END SUBROUTINE
4344 
4345 !-------------------------------------------------------------------------------
4352 !-------------------------------------------------------------------------------
4353  SUBROUTINE siesta_write(this, iou)
4354  USE safe_open_mod
4355  USE siesta_context
4356  USE file_opts
4357 
4358  IMPLICIT NONE
4359 
4360 ! Declare Arguments
4361  TYPE (siesta_class), INTENT(in) :: this
4362  INTEGER, INTENT(in) :: iou
4363 
4364 ! local variables
4365  INTEGER :: i
4366  INTEGER :: iou_nl
4367  INTEGER :: status
4368  REAL (rprec) :: start_time
4369 
4370 ! Start of executable code
4371  start_time = profiler_get_start_time()
4372  WRITE (iou,*)
4373  WRITE (iou,*) 'Equilibrium Type : SIESTA'
4374  IF (ASSOCIATED(this%ne)) THEN
4375  CALL pprofile_write(this%ne, 'pp_ne', iou)
4376  END IF
4377 
4378  IF (ASSOCIATED(this%sxrem)) THEN
4379  DO i = 1, SIZE(this%sxrem)
4380  IF (ASSOCIATED(this%sxrem(i)%p)) THEN
4381  CALL pprofile_write(this%sxrem(i)%p, 'pp_sxrem', iou)
4382  END IF
4383  END DO
4384  END IF
4385 
4386  IF (ASSOCIATED(this%te)) THEN
4387  CALL pprofile_write(this%te, 'pp_te', iou)
4388  END IF
4389 
4390  IF (ASSOCIATED(this%ti)) THEN
4391  CALL pprofile_write(this%ti, 'pp_ti', iou)
4392  END IF
4393 
4394 ! Update the namelist input variables from the equilibrium solution.
4395  CALL siesta_set_namelist(this)
4396 
4397  iou_nl = 0
4398  CALL safe_open(iou_nl, status, &
4399  & trim(this%siesta_file_name) // '_out', &
4400  & 'replace', 'formatted', delim_in='quote')
4401  WRITE (iou_nl, nml=siesta_info)
4402  CLOSE (iou_nl, iostat=status)
4403 
4404  CALL move_file(trim(this%siesta_file_name) // '_save', &
4405  & trim(this%siesta_file_name), status)
4406 
4407  CALL vmec_write(this%vmec, iou)
4408 
4409  CALL profiler_set_stop_time('siesta_write', start_time)
4410 
4411  END SUBROUTINE
4412 
4413 !-------------------------------------------------------------------------------
4422 !-------------------------------------------------------------------------------
4423  SUBROUTINE siesta_write_input(this, current_step)
4424  USE safe_open_mod
4425  USE file_opts
4426  USE siesta_context
4427 
4428  IMPLICIT NONE
4429 
4430 ! Declare Arguments
4431  TYPE (siesta_class), INTENT(in) :: this
4432  INTEGER, INTENT(in) :: current_step
4433 
4434 ! local variables
4435  INTEGER :: iou_nl
4436  INTEGER :: status
4437  CHARACTER (len=path_length) :: filename
4438  REAL (rprec) :: start_time
4439 
4440 ! Start of executable code
4441  start_time = profiler_get_start_time()
4442 
4443 ! Update the namelist input variables from the equilibrium solution.
4444  CALL siesta_set_namelist(this)
4445 
4446 ! Write out the final namelist input file.
4447  iou_nl = 0
4448  status = 0
4449 
4450  WRITE (filename,1000) trim(this%siesta_file_name), current_step
4451 
4452  CALL safe_open(iou_nl, status, trim(filename), 'replace', &
4453  & 'formatted', delim_in='none')
4454  WRITE (iou_nl, nml=siesta_info)
4455  CLOSE (iou_nl, iostat=status)
4456  CALL assert_eq(status, 0, 'Error saving input file.')
4457 
4458 ! Save the restart file by making a copy with a new name.
4459  WRITE (filename,1000) trim(this%restart_file_name), current_step
4460 
4461  CALL copy_file(this%restart_file_name, trim(filename), status)
4462  CALL assert_eq(status, 0, 'Error copying wout file.')
4463 
4464  CALL vmec_write_input(this%vmec, current_step)
4465 
4466  CALL profiler_set_stop_time('vmec_write_input', start_time)
4467 
4468 1000 FORMAT(a,'_',i0.3)
4469 
4470  END SUBROUTINE
4471 
4472 !*******************************************************************************
4473 ! NETCDF SUBROUTINES
4474 !*******************************************************************************
4492 !-------------------------------------------------------------------------------
4504 !-------------------------------------------------------------------------------
4505  SUBROUTINE siesta_def_result(this, result_ncid, maxnsetps_dim_id)
4506 
4507  IMPLICIT NONE
4508 
4509 ! Declare Arguments
4510  TYPE (siesta_class), INTENT(in) :: this
4511  INTEGER, INTENT(in) :: result_ncid
4512  INTEGER, INTENT(in) :: maxnsetps_dim_id
4513 
4514 ! Local variables
4515  REAL (rprec) :: start_time
4516 
4517 ! Start of executable code
4518  start_time = profiler_get_start_time()
4519 
4520  CALL vmec_def_result(this%vmec, result_ncid, maxnsetps_dim_id)
4521 
4522  CALL profiler_set_stop_time('siesta_def_result', start_time)
4523 
4524  END SUBROUTINE
4525 
4526 !-------------------------------------------------------------------------------
4533 !-------------------------------------------------------------------------------
4534  SUBROUTINE siesta_write_init_data(this, result_ncid)
4535 
4536  IMPLICIT NONE
4537 
4538 ! Declare Arguments
4539  TYPE (siesta_class), INTENT(in) :: this
4540  INTEGER, INTENT(in) :: result_ncid
4541 
4542 ! Local variables
4543  REAL (rprec) :: start_time
4544 
4545 ! Start of executable code
4546  start_time = profiler_get_start_time()
4547 
4548 ! vmec_write_init_data will write step data for the first step. Call this
4549 ! before the call to siesta_write_init_data to avoid writing the step twice.
4550  CALL vmec_write_init_data(this%vmec, result_ncid)
4551  CALL siesta_write_step_data(this, result_ncid, 1)
4552 
4553  CALL profiler_set_stop_time('siesta_write_init_data', start_time)
4554 
4555  END SUBROUTINE
4556 
4557 !-------------------------------------------------------------------------------
4565 !-------------------------------------------------------------------------------
4566  SUBROUTINE siesta_write_step_data(this, result_ncid, current_step)
4567 
4568  IMPLICIT NONE
4569 
4570 ! Declare Arguments
4571  TYPE (siesta_class), INTENT(in) :: this
4572  INTEGER, INTENT(in) :: result_ncid
4573  INTEGER, INTENT(in) :: current_step
4574 
4575 ! Local variables
4576  REAL (rprec) :: start_time
4577 
4578 ! Start of executable code
4579  start_time = profiler_get_start_time()
4580 
4581 ! vmec_write_init_data called from siesta_write_init_data already write the
4582 ! first vmec step.
4583  IF (current_step .ne. 1) THEN
4584  CALL vmec_write_step_data(this%vmec, result_ncid, current_step)
4585  END IF
4586 
4587  CALL profiler_set_stop_time('siesta_write_step_data', start_time)
4588 
4589  END SUBROUTINE
4590 
4591 !-------------------------------------------------------------------------------
4599 !-------------------------------------------------------------------------------
4600  SUBROUTINE siesta_restart(this, result_ncid, current_step)
4601 
4602  IMPLICIT NONE
4603 
4604 ! Declare Arguments
4605  TYPE (siesta_class), INTENT(in) :: this
4606  INTEGER, INTENT(in) :: result_ncid
4607  INTEGER, INTENT(in) :: current_step
4608 
4609 ! Local variables
4610  REAL (rprec) :: start_time
4611 
4612 ! Start of executable code
4613  start_time = profiler_get_start_time()
4614 
4615  CALL vmec_restart(this%vmec, result_ncid, current_step)
4616 
4617  CALL profiler_set_stop_time('siesta_write_step_data', start_time)
4618 
4619  END SUBROUTINE
4620 
4621 !*******************************************************************************
4622 ! MPI SUBROUTINES
4623 !*******************************************************************************
4624 !-------------------------------------------------------------------------------
4632 !-------------------------------------------------------------------------------
4633  SUBROUTINE siesta_sync_state(this, recon_comm)
4634  USE file_opts
4635 
4636  IMPLICIT NONE
4637 
4638 ! Declare Arguments
4639  TYPE (siesta_class), INTENT(inout) :: this
4640  INTEGER, INTENT(in) :: recon_comm
4641 
4642 #if defined(MPI_OPT)
4643 ! local variables
4644  INTEGER :: i
4645  INTEGER :: error
4646  INTEGER :: mpi_rank
4647 
4648 ! Start of executable code
4649  CALL mpi_comm_rank(recon_comm, mpi_rank, error)
4650 
4651 ! The barrier ensures that the wout is ready to be copied from the parent
4652 ! directory. Otherwise a child process can try to copy a wout file that does
4653 ! not exist.
4654  CALL mpi_barrier(recon_comm, error)
4655 
4656 ! Copy the vmec state first to ensure the wout file is loaded before seting the
4657 ! siesta magnetic cache.
4658  CALL vmec_sync_state(this%vmec, recon_comm)
4659 
4660 ! If this is the child process, load the wout file.
4661  IF (mpi_rank .gt. 0) THEN
4662  CALL copy_file(build_path('..', this%restart_file_name), &
4663  & trim(this%restart_file_name), error)
4664 
4665  IF (ASSOCIATED(this%context)) THEN
4666  CALL siesta_context_read(this%context, &
4667  & this%restart_file_name)
4668  ELSE
4669  this%context => &
4670  & siesta_context_construct(this%restart_file_name)
4671  END IF
4672 
4673  IF (ASSOCIATED(this%magnetic_cache)) THEN
4674  CALL siesta_set_magnetic_cache(this)
4675  END IF
4676  END IF
4677 #endif
4678  END SUBROUTINE
4679 
4680 !-------------------------------------------------------------------------------
4689 !-------------------------------------------------------------------------------
4690  SUBROUTINE siesta_sync_child(this, index, recon_comm)
4691  USE file_opts
4692 
4693  IMPLICIT NONE
4694 
4695 ! Declare Arguments
4696  TYPE (siesta_class), INTENT(inout) :: this
4697  INTEGER, INTENT(in) :: index
4698  INTEGER, INTENT(in) :: recon_comm
4699 
4700 #if defined(MPI_OPT)
4701 ! local variables
4702  INTEGER :: error
4703  INTEGER :: mpi_rank
4704 
4705 ! Start of executable code
4706  CALL mpi_comm_rank(recon_comm, mpi_rank, error)
4707 
4708 ! Copy the vmec state first to ensure the wout file is loaded before seting the
4709 ! siesta magnetic cache.
4710  CALL vmec_sync_child(this%vmec, index, recon_comm)
4711 
4712 ! If this is the parent process, load the wout file.
4713  IF (mpi_rank .eq. 0) THEN
4714  CALL copy_file(build_path(process_dir(index + 1), &
4715  & this%restart_file_name), &
4716  & trim(this%restart_file_name), error)
4717  CALL assert_eq(error, 0, 'Error reading synced restart file.')
4718 
4719  CALL siesta_context_read(this%context, this%restart_file_name)
4720 
4721  IF (ASSOCIATED(this%magnetic_cache)) THEN
4722  CALL siesta_set_magnetic_cache(this)
4723  END IF
4724  END IF
4725 
4726 #endif
4727  END SUBROUTINE
4728 
4729 !*******************************************************************************
4730 ! PRIVATE
4731 !*******************************************************************************
4732 !-------------------------------------------------------------------------------
4741 !-------------------------------------------------------------------------------
4742  PURE FUNCTION siesta_to_vmec_s(s)
4743 
4744  IMPLICIT NONE
4745 
4746 ! Declare Arguments
4747  REAL (rprec) :: siesta_to_vmec_s
4748  REAL (rprec), INTENT(in) :: s
4749 
4750 ! Start of executable code
4751  siesta_to_vmec_s = s*s
4752 
4753  END FUNCTION
4754 
4755 !-------------------------------------------------------------------------------
4764 !-------------------------------------------------------------------------------
4765  PURE FUNCTION siesta_to_siesta_s(s)
4766 
4767  IMPLICIT NONE
4768 
4769 ! Declare Arguments
4770  REAL (rprec) :: siesta_to_siesta_s
4771  REAL (rprec), INTENT(in) :: s
4772 
4773 ! Start of executable code
4774  siesta_to_siesta_s = sqrt(s)
4775 
4776  END FUNCTION
4777 
4778 !-------------------------------------------------------------------------------
4797 !-------------------------------------------------------------------------------
4798  PURE FUNCTION siesta_to_i_vmec(s)
4799  USE read_wout_mod, only: ns
4800 
4801  IMPLICIT NONE
4802 
4803 ! Declare Arguments
4804  REAL (rprec) :: siesta_to_i_vmec
4805  REAL (rprec), INTENT(in) :: s
4806 
4807 ! Start of executable code
4808  siesta_to_i_vmec = siesta_to_vmec_s(s)*(ns - 1.0) + 1.0
4809 
4810  END FUNCTION
4811 
4812 !-------------------------------------------------------------------------------
4820 !-------------------------------------------------------------------------------
4821  PURE FUNCTION siesta_to_i_vmec_low(s)
4822 
4823  IMPLICIT NONE
4824 
4825 ! Declare Arguments
4826  INTEGER :: siesta_to_i_vmec_low
4827  REAL (rprec), INTENT(in) :: s
4828 
4829 ! Start of executable code
4830  siesta_to_i_vmec_low = floor(siesta_to_i_vmec(s))
4831 
4832  END FUNCTION
4833 
4834 !-------------------------------------------------------------------------------
4860 !-------------------------------------------------------------------------------
4861  PURE FUNCTION siesta_to_w_vmec_low(s, ilow)
4862 
4863  IMPLICIT NONE
4864 
4865 ! Declare Arguments
4866  REAL (rprec) :: siesta_to_w_vmec_low
4867  REAL (rprec), INTENT(in) :: s
4868  INTEGER, INTENT(in) :: ilow
4869 
4870 ! Start of executable code
4871  siesta_to_w_vmec_low = -siesta_to_i_vmec(s) + 1.0 + ilow
4872 
4873  END FUNCTION
4874 
4875  END MODULE
siesta_equilibrium::siesta_get_p_radial
real(rprec) function siesta_get_p_radial(this, s, normalize)
Gets the plasma pressure at a s, u=0, v=0 position.
Definition: siesta_equilibrium.f:2856
model_state::model_state_shift_flag
integer, parameter model_state_shift_flag
Shift parameter changed bit position.
Definition: model_state.f:36
siesta_equilibrium::siesta_get_ti_cart
real(rprec) function siesta_get_ti_cart(this, x_cart)
Gets the ion temperature at a cartesian position.
Definition: siesta_equilibrium.f:2468
siesta_equilibrium::siesta_pp_te_af_id
integer, parameter siesta_pp_te_af_id
1D Array of segment electron temperature profile f values. Defined in pprofile_T.
Definition: siesta_equilibrium.f:111
siesta_equilibrium::siesta_get_con_surface_kphigrid
real(rprec) function, dimension(:,:), pointer siesta_get_con_surface_kphigrid(this)
Get the conducting surface integration kphi grid points.
Definition: siesta_equilibrium.f:3381
siesta_equilibrium::siesta_get_magnetic_volume_rgrid
real(rprec) function, dimension(:,:,:), pointer siesta_get_magnetic_volume_rgrid(this)
Get volume magnetic volume integration radial grid points.
Definition: siesta_equilibrium.f:3163
siesta_equilibrium::siesta_get_gp_ti_ij
real(rprec) function siesta_get_gp_ti_ij(this, i, j)
Get the ion temperature gp kernel value for the two indicies.
Definition: siesta_equilibrium.f:2352
coordinate_utilities
Module is part of the LIBSTELL. This modules containes code to convert from different coordinate syst...
Definition: coordinate_utilities.f:12
vmec_equilibrium::vmec_converge
logical function vmec_converge(this, num_iter, iou, eq_comm)
Solves the VMEC equilibrium.
Definition: vmec_equilibrium.f:5202
siesta_equilibrium::siesta_get_te_cart
real(rprec) function siesta_get_te_cart(this, x_cart)
Gets the electron temperature at a cartesian position.
Definition: siesta_equilibrium.f:2217
siesta_equilibrium::siesta_pp_ne_b_id
integer, parameter siesta_pp_ne_b_id
1D Array of function density profile parameters. Defined in pprofile_T.
Definition: siesta_equilibrium.f:92
vmec_equilibrium::vmec_write
subroutine vmec_write(this, iou)
Write out the equilibrium to an output file.
Definition: vmec_equilibrium.f:5549
siesta_equilibrium::siesta_get_ne_af
real(rprec) function, dimension(:), pointer siesta_get_ne_af(this)
Gets the electron density profile af array.
Definition: siesta_equilibrium.f:1810
siesta_equilibrium::siesta_get_grid_stop
real(rprec) function siesta_get_grid_stop(this)
Get radial grid size.
Definition: siesta_equilibrium.f:3869
siesta_equilibrium::siesta_get_p_flux
real(rprec) function siesta_get_p_flux(this, flux, normalize)
Gets the plasma pressure at a flux position.
Definition: siesta_equilibrium.f:2891
siesta_equilibrium::siesta_get_area_int_element
real(rprec) function siesta_get_area_int_element(this)
Get area integration element.
Definition: siesta_equilibrium.f:3447
vmec_equilibrium::vmec_rbc_id
integer, parameter vmec_rbc_id
2D Array of outer boundary R cosine coeffiecents. Defined in vmec_input.
Definition: vmec_equilibrium.f:173
vmec_equilibrium::vmec_is_1d_array
logical function vmec_is_1d_array(this, id)
Checks if a parameter id is a 1d array.
Definition: vmec_equilibrium.f:5020
siesta_equilibrium::siesta_is_scaler_value
logical function siesta_is_scaler_value(this, id)
Checks if a parameter id is a scaler value.
Definition: siesta_equilibrium.f:3899
vmec_equilibrium::vmec_write_step_data
subroutine vmec_write_step_data(this, result_ncid, current_step)
Write step data to NetCDF result file.
Definition: vmec_equilibrium.f:5915
siesta_equilibrium::siesta_save_state
subroutine siesta_save_state(this)
Save the internal state of the equilibrium.
Definition: siesta_equilibrium.f:4247
vmec_equilibrium::vmec_get_ext_currents
real(rprec) function, dimension(:), pointer vmec_get_ext_currents(this, num_currents, scale_currents)
Get external current.
Definition: vmec_equilibrium.f:4611
vmec_equilibrium::vmec_get_suv
real(rprec) function, dimension(3) vmec_get_suv(this, x_cart)
Gets VMEC suv position.
Definition: vmec_equilibrium.f:4087
file_opts::copy_file
subroutine copy_file(file_source, file_dest, error)
Copies the source file to the destination.
Definition: file_opts.f:78
siesta_equilibrium::siesta_get_suv
real(rprec) function, dimension(3) siesta_get_suv(this, x_cart)
Gets siesta s value at position.
Definition: siesta_equilibrium.f:3088
siesta_equilibrium::siesta_get_te_radial
real(rprec) function siesta_get_te_radial(this, s)
Gets the electron temperature at a s, u=0, v=0 position.
Definition: siesta_equilibrium.f:2250
vmec_equilibrium::vmec_get_param_name
character(len=data_name_length) function vmec_get_param_name(this, id)
Get the name of a reconstruction VMEC parameter.
Definition: vmec_equilibrium.f:2557
magnetic_response::magnetic_response_use_shell
logical function magnetic_response_use_shell(this)
Checks if the conducting shell flag is set.
Definition: magnetic_response.f:831
siesta_equilibrium::siesta_get_sxrem_cart
real(rprec) function siesta_get_sxrem_cart(this, x_cart, index)
Gets the soft x-ray emissivity at a cartesian position.
Definition: siesta_equilibrium.f:2742
siesta_equilibrium::siesta_def_result
subroutine siesta_def_result(this, result_ncid, maxnsetps_dim_id)
Define NetCDF variables for the result file.
Definition: siesta_equilibrium.f:4506
siesta_equilibrium::siesta_get_te
Interface for the siesta electron temperature profile values.
Definition: siesta_equilibrium.f:293
siesta_equilibrium::siesta_magnetic_cache
Object to cache data needed to compute magnetic signals.
Definition: siesta_equilibrium.f:151
siesta_equilibrium::siesta_is_using_point
logical function siesta_is_using_point(this)
Checks if a the point magnetics are being used.
Definition: siesta_equilibrium.f:4060
siesta_equilibrium::siesta_get_gp_sxrem_pp
real(rprec) function siesta_get_gp_sxrem_pp(this, x_cart, y_cart, index)
Get the soft x-ray emissivity gp kernel value for the position and position.
Definition: siesta_equilibrium.f:2699
siesta_equilibrium::siesta_get_con_surface_krgrid
real(rprec) function, dimension(:,:), pointer siesta_get_con_surface_krgrid(this)
Get the conducting surface integration kr grid points.
Definition: siesta_equilibrium.f:3351
siesta_equilibrium::siesta_get_b_vec
real(rprec) function, dimension(3) siesta_get_b_vec(this, x_cart, cyl)
Gets the magnetic field vector at a position.
Definition: siesta_equilibrium.f:2986
siesta_equilibrium::siesta_pp_sxrem_af_id
integer, parameter siesta_pp_sxrem_af_id
1D Array of segment sxrem profile f values. Defined in pprofile_T.
Definition: siesta_equilibrium.f:102
siesta_equilibrium::siesta_get_magnetic_volume_jphigrid
real(rprec) function, dimension(:,:,:), pointer siesta_get_magnetic_volume_jphigrid(this)
Get volume magnetic volume integration jphi grid points.
Definition: siesta_equilibrium.f:3253
vmec_equilibrium::vmec_read_vac_file
subroutine vmec_read_vac_file(this, index, eq_comm)
Loads the vacuum magnetic field file.
Definition: vmec_equilibrium.f:5395
siesta_equilibrium::siesta_get_gp_ne_ij
real(rprec) function siesta_get_gp_ne_ij(this, i, j)
Get the electron density gp kernel value for the two indicies.
Definition: siesta_equilibrium.f:1843
siesta_equilibrium::siesta_get_plasma_edge
integer function siesta_get_plasma_edge(this, phi, r, z)
Gets the r and z positions of the outer surface at a toroidal angle.
Definition: siesta_equilibrium.f:3128
siesta_equilibrium::siesta_pp_te_b_id
integer, parameter siesta_pp_te_b_id
1D Array of function electron temperature profile parameters. Defined in pprofile_T.
Definition: siesta_equilibrium.f:105
vmec_equilibrium::vmec_class
Base class representing a vmec_equilibrium.
Definition: vmec_equilibrium.f:396
magnetic_response::magnetic_response_use_plasma
logical function magnetic_response_use_plasma(this)
Checks if the plasma flag is set.
Definition: magnetic_response.f:800
siesta_context::siesta_context_destruct
subroutine siesta_context_destruct(this)
Deconstruct a siesta_context_class object.
Definition: siesta_context.f:289
siesta_equilibrium::siesta_get_gp_te_num_hyper_param
integer function siesta_get_gp_te_num_hyper_param(this)
Get the number of electron temperature gp kernel hyper parameters.
Definition: siesta_equilibrium.f:2026
siesta_equilibrium::siesta_set_magnetic_cache
Interface for the setting the siesta_magnetic_cache types using siesta_set_magnetic_cache_responce,...
Definition: siesta_equilibrium.f:267
siesta_equilibrium::siesta_converge
logical function siesta_converge(this, num_iter, iou, eq_comm, state_flags)
Solves the siesta equilibrium.
Definition: siesta_equilibrium.f:4107
siesta_equilibrium::siesta_get_grid_size
integer function siesta_get_grid_size(this)
Get radial grid size.
Definition: siesta_equilibrium.f:3811
siesta_equilibrium::siesta_get_gp_ti_pp
real(rprec) function siesta_get_gp_ti_pp(this, x_cart, y_cart)
Get the ion temperature gp kernel value for the position and position.
Definition: siesta_equilibrium.f:2429
siesta_equilibrium::siesta_get_ti_radial
real(rprec) function siesta_get_ti_radial(this, s)
Gets the ion temperature at a s, u=0, v=0 position.
Definition: siesta_equilibrium.f:2501
vmec_equilibrium::vmec_get_plasma_edge
integer function vmec_get_plasma_edge(this, phi, r, z)
Gets the r and z positions of the outer surface at a toroidal angle.
Definition: vmec_equilibrium.f:4227
siesta_equilibrium::siesta_get_param_id
integer function siesta_get_param_id(this, param_name)
Get the id for a reconstruction parameter.
Definition: siesta_equilibrium.f:1527
coordinate_utilities::cart_to_cyl_vec
pure real(rprec) function, dimension(3) cart_to_cyl_vec(cart, vec)
Convert vector from cartesian coordinates to cylindical coordinates.
Definition: coordinate_utilities.f:126
siesta_equilibrium::siesta_pp_ze_as_id
integer, parameter siesta_pp_ze_as_id
1D Array of segment zeff profile s poitions. Defined in pprofile_T.
Definition: siesta_equilibrium.f:126
siesta_equilibrium::siesta_get_ti
Interface for the siesta ion temperature profile values.
Definition: siesta_equilibrium.f:312
siesta_equilibrium::siesta_get_ti_af
real(rprec) function, dimension(:), pointer siesta_get_ti_af(this)
Gets the ion temperature profile af array.
Definition: siesta_equilibrium.f:2317
siesta_equilibrium::siesta_reset_state
subroutine siesta_reset_state(this)
Reset the internal state of the equilibrium.
Definition: siesta_equilibrium.f:4300
siesta_equilibrium::siesta_set_magnetic_cache_calc
subroutine siesta_set_magnetic_cache_calc(this)
Set magnetic cache.
Definition: siesta_equilibrium.f:936
vmec_equilibrium::vmec_get_param_id
integer function vmec_get_param_id(this, param_name)
Get the id for a reconstruction parameter.
Definition: vmec_equilibrium.f:2045
siesta_equilibrium::siesta_get_gp_sxrem_pi
real(rprec) function siesta_get_gp_sxrem_pi(this, x_cart, i, index)
Get the soft x-ray emissivity gp kernel value for the position and index.
Definition: siesta_equilibrium.f:2652
model_state::model_state_te_flag
integer, parameter model_state_te_flag
Temperature profile changed bit position.
Definition: model_state.f:32
siesta_equilibrium::siesta_to_i_vmec_low
pure integer function siesta_to_i_vmec_low(s)
Find the lower vmec s index.
Definition: siesta_equilibrium.f:4822
siesta_equilibrium::siesta_to_i_vmec
pure real(rprec) function siesta_to_i_vmec(s)
Find the vmec s index.
Definition: siesta_equilibrium.f:4799
coordinate_utilities::cyl_to_cart_vec
pure real(rprec) function, dimension(3), public cyl_to_cart_vec(cyl, vec)
Convert vector from cylindical coordinates to cartesian coordinates.
Definition: coordinate_utilities.f:95
siesta_equilibrium::siesta_write_init_data
subroutine siesta_write_init_data(this, result_ncid)
Write inital data to NetCDF result file.
Definition: siesta_equilibrium.f:4535
siesta_equilibrium::siesta_write_input
subroutine siesta_write_input(this, current_step)
Write the current valid input.
Definition: siesta_equilibrium.f:4424
siesta_equilibrium::siesta_get_sxrem_radial
real(rprec) function siesta_get_sxrem_radial(this, s, index)
Gets the soft x-ray emissivity at a s, u=0, v=0 position.
Definition: siesta_equilibrium.f:2783
siesta_equilibrium::siesta_get_ext_currents
real(rprec) function, dimension(:), pointer siesta_get_ext_currents(this, num_currents, scale_currents)
Get external current.
Definition: siesta_equilibrium.f:3483
vmec_equilibrium::vmec_set_namelist
subroutine vmec_set_namelist(this)
Sets namelist variables from the solved VMEC equilibrium.
Definition: vmec_equilibrium.f:1987
siesta_equilibrium::siesta_get_p_cart
real(rprec) function siesta_get_p_cart(this, x_cart, normalize)
Gets the plasma pressure at a cartesian position.
Definition: siesta_equilibrium.f:2822
vmec_equilibrium::vmec_restart
subroutine vmec_restart(this, result_ncid, current_step)
Restart vmec from the result file.
Definition: vmec_equilibrium.f:5970
siesta_equilibrium::siesta_get_gp_te_ij
real(rprec) function siesta_get_gp_te_ij(this, i, j)
Get the electron temperature gp kernel value for the two indicies.
Definition: siesta_equilibrium.f:2096
siesta_equilibrium::siesta_is_1d_array
logical function siesta_is_1d_array(this, id)
Checks if a parameter id is a 1d array.
Definition: siesta_equilibrium.f:3936
model_state
Contains parameters defining the bit positions for flags that mark changes in different parts of the ...
Definition: model_state.f:11
siesta_equilibrium::siesta_get_gp_sxrem_ij
real(rprec) function siesta_get_gp_sxrem_ij(this, i, j, index)
Get the soft x-ray emissivity gp kernel value for the two indicies.
Definition: siesta_equilibrium.f:2610
file_opts::delete_file
subroutine delete_file(file_source, error)
Deletes a file.
Definition: file_opts.f:215
siesta_equilibrium::siesta_get_magnetic_volume_jrgrid
real(rprec) function, dimension(:,:,:), pointer siesta_get_magnetic_volume_jrgrid(this)
Get volume magnetic volume integration jr grid points.
Definition: siesta_equilibrium.f:3223
siesta_equilibrium::siesta_pp_ne_as_id
integer, parameter siesta_pp_ne_as_id
1D Array of segment density profile s poitions. Defined in pprofile_T.
Definition: siesta_equilibrium.f:94
siesta_equilibrium::siesta_set_magnetic_cache_responce
subroutine siesta_set_magnetic_cache_responce(this, response_object)
Set magnetic cache for plasma responce.
Definition: siesta_equilibrium.f:778
siesta_equilibrium::siesta_pp_ti_b_id
integer, parameter siesta_pp_ti_b_id
1D Array of function ion temperature profile parameters. Defined in pprofile_T.
Definition: siesta_equilibrium.f:114
siesta_equilibrium::siesta_get_ext_b_plasma
real(rprec) function, dimension(3) siesta_get_ext_b_plasma(this, position, axi_only)
Get external plasma magnetic field.
Definition: siesta_equilibrium.f:3521
siesta_equilibrium::siesta_get_gp_ti_pi
real(rprec) function siesta_get_gp_ti_pi(this, x_cart, i)
Get the ion temperature gp kernel value for the position and index.
Definition: siesta_equilibrium.f:2389
siesta_equilibrium::siesta_get_te_af
real(rprec) function, dimension(:), pointer siesta_get_te_af(this)
Gets the electron temperature profile af array.
Definition: siesta_equilibrium.f:2061
siesta_equilibrium::siesta_get_sxrem
Interface for the siesta soft x-ray emission profile values.
Definition: siesta_equilibrium.f:329
siesta_equilibrium::siesta_get_p
Interface for the siesta pressure profile values.
Definition: siesta_equilibrium.f:347
siesta_equilibrium::siesta_set_namelist
subroutine siesta_set_namelist(this)
Sets namelist variables from the solved siesta equilibrium.
Definition: siesta_equilibrium.f:1496
siesta_equilibrium::siesta_helpert_id
integer, parameter siesta_helpert_id
Helical perturbation parameter.
Definition: siesta_equilibrium.f:140
siesta_equilibrium::siesta_write_step_data
subroutine siesta_write_step_data(this, result_ncid, current_step)
Write step data to NetCDF result file.
Definition: siesta_equilibrium.f:4567
sxrem
Defines the base class of the type sxrem_class.
Definition: sxrem.f:13
siesta_equilibrium::siesta_destruct
subroutine siesta_destruct(this)
Deconstruct a siesta_class object.
Definition: siesta_equilibrium.f:496
model_state::model_state_ti_flag
integer, parameter model_state_ti_flag
Ion profile changed bit position.
Definition: model_state.f:34
siesta_equilibrium::siesta_get_param_value
real(rprec) function siesta_get_param_value(this, id, i_index, j_index)
Get the value of a reconstruction siesta parameter.
Definition: siesta_equilibrium.f:1610
vmec_equilibrium::vmec_sync_child
subroutine vmec_sync_child(this, index, recon_comm)
Syncronize a child vmec state to the parent.
Definition: vmec_equilibrium.f:6081
model_state::model_state_vmec_flag
integer, parameter model_state_vmec_flag
VMEC Equilibrium changed bit position.
Definition: model_state.f:26
siesta_context::siesta_context_construct
type(siesta_context_class) function, pointer siesta_context_construct(restart_file_name)
Construct a siesta_context_class object.
Definition: siesta_context.f:133
siesta_equilibrium::siesta_is_2d_array
logical function siesta_is_2d_array(this, id)
Checks if a parameter id is a 2d array.
Definition: siesta_equilibrium.f:3977
coordinate_utilities::cyl_to_cart
pure real(rprec) function, dimension(3), public cyl_to_cart(cyl)
Convert a point from cylindical coordinates to cartesian coordinates.
Definition: coordinate_utilities.f:67
file_opts::move_file
subroutine move_file(file_source, file_dest, error)
Moves the source file to the destination.
Definition: file_opts.f:37
file_opts
Contains cross platform routines for manipulating files on the file system. Defines a functions to mo...
Definition: file_opts.f:13
siesta_equilibrium::siesta_pp_ze_af_id
integer, parameter siesta_pp_ze_af_id
1D Array of segment zeff profile f values. Defined in pprofile_T.
Definition: siesta_equilibrium.f:129
siesta_equilibrium::siesta_to_vmec_s
pure real(rprec) function siesta_to_vmec_s(s)
Map the siesta s to the vmec s.
Definition: siesta_equilibrium.f:4743
siesta_equilibrium::siesta_pp_sxrem_as_id
integer, parameter siesta_pp_sxrem_as_id
1D Array of segment sxrem profile s poitions. Defined in pprofile_T.
Definition: siesta_equilibrium.f:100
vmec_equilibrium::vmec_zbs_id
integer, parameter vmec_zbs_id
2D Array of outer boundary Z sine coeffiecents. Defined in vmec_input.
Definition: vmec_equilibrium.f:175
model_state::model_state_sxrem_flag
integer, parameter model_state_sxrem_flag
Soft x-ray emissivity profile changed bit position. This flag needs to always be the last flag an can...
Definition: model_state.f:46
siesta_equilibrium::siesta_to_w_vmec_low
pure real(rprec) function siesta_to_w_vmec_low(s, ilow)
Find the lower interplation weighting.
Definition: siesta_equilibrium.f:4862
siesta_equilibrium::siesta_to_siesta_s
pure real(rprec) function siesta_to_siesta_s(s)
Map the vmec s to the siesta s.
Definition: siesta_equilibrium.f:4766
siesta_equilibrium::siesta_get_gp_ne
Interface for the siesta guassian process density profile values.
Definition: siesta_equilibrium.f:284
siesta_equilibrium::siesta_pp_ne_af_id
integer, parameter siesta_pp_ne_af_id
1D Array of segment density profile f values. Defined in pprofile_T.
Definition: siesta_equilibrium.f:96
siesta_error
This module contains all the code needed to define error.
Definition: siesta_error.f90:10
siesta_equilibrium::siesta_get_ne
Interface for the siesta density profile values.
Definition: siesta_equilibrium.f:276
siesta_equilibrium::siesta_get_gp_sxrem
Interface for the siesta guassian process soft x-ray emissivity profile values.
Definition: siesta_equilibrium.f:338
siesta_equilibrium::siesta_get_gp_ne_pp
real(rprec) function siesta_get_gp_ne_pp(this, x_cart, y_cart)
Get the electron density gp kernel value for the position and position.
Definition: siesta_equilibrium.f:1922
siesta_equilibrium::siesta_get_gp_te_pi
real(rprec) function siesta_get_gp_te_pi(this, x_cart, i)
Get the electron temperature gp kernel value for the position and index.
Definition: siesta_equilibrium.f:2134
file_opts::build_path
character(len=path_length) function build_path(path, name)
Builds a path.
Definition: file_opts.f:382
model_state::model_state_ne_flag
integer, parameter model_state_ne_flag
Denisty profile changed bit position.
Definition: model_state.f:30
vmec_equilibrium::vmec_save_state
subroutine vmec_save_state(this)
Save the internal state of the equilibrium.
Definition: vmec_equilibrium.f:5429
vmec_equilibrium::vmec_is_2d_array
logical function vmec_is_2d_array(this, id)
Checks if a parameter id is a 2d array.
Definition: vmec_equilibrium.f:5067
vmec_equilibrium::vmec_get_b_vac
real(rprec) function, dimension(3) vmec_get_b_vac(this, r_cyl)
Gets the vacuum magnetic field vector at a position.
Definition: vmec_equilibrium.f:4662
siesta_equilibrium::siesta_get_ne_cart
real(rprec) function siesta_get_ne_cart(this, x_cart)
Gets the electron density at a cartesian position.
Definition: siesta_equilibrium.f:1961
data_parameters
This modules contains parameters used by equilibrium models.
Definition: data_parameters.f:10
line_segment
Module is part of the LIBSTELL. This module contains code to create a profile constructed of line sig...
Definition: line_segment.f:13
coordinate_utilities::cart_to_cyl
pure real(rprec) function, dimension(3), public cart_to_cyl(cart)
Convert a point from cartes cartesian coordinates to cylindical coordinates.
Definition: coordinate_utilities.f:40
siesta_equilibrium::siesta_get_int_b_dphi
real(rprec) function siesta_get_int_b_dphi(this, r, theta)
Gets the loop integrated magnetic field at a position.
Definition: siesta_equilibrium.f:3042
siesta_equilibrium::siesta_get_grid_start
real(rprec) function siesta_get_grid_start(this)
Get start of the radial grid.
Definition: siesta_equilibrium.f:3843
siesta_equilibrium::siesta_get_magnetic_volume_zgrid
real(rprec) function, dimension(:,:,:), pointer siesta_get_magnetic_volume_zgrid(this)
Get volume magnetic volume integration z grid points.
Definition: siesta_equilibrium.f:3193
siesta_context
Defines the base class of the type siesta_context_class. This contains the state variables needed by ...
Definition: siesta_context.f:11
vmec_equilibrium::magnetic_cache_vc_grid_dim
real(rprec), parameter magnetic_cache_vc_grid_dim
Virtual casing grid dimension.
Definition: vmec_equilibrium.f:315
vmec_equilibrium::vmec_write_input
subroutine vmec_write_input(this, current_step)
Write the current valid input.
Definition: vmec_equilibrium.f:5617
siesta_equilibrium::siesta_read_vac_file
subroutine siesta_read_vac_file(this, index, eq_comm)
Loads the vacuum magnetic field file.
Definition: siesta_equilibrium.f:4226
vmec_equilibrium::vmec_def_result
subroutine vmec_def_result(this, result_ncid, maxnsetps_dim_id)
Define NetCDF variables for the result file.
Definition: vmec_equilibrium.f:5706
siesta_equilibrium::siesta_get_gp_ti
Interface for the siesta guassian process ion temperature profile values.
Definition: siesta_equilibrium.f:320
siesta_equilibrium::siesta_class
Base class representing a siesta_equilibrium.
Definition: siesta_equilibrium.f:221
vmec_equilibrium::vmec_write_init_data
subroutine vmec_write_init_data(this, result_ncid)
Write inital data to NetCDF result file.
Definition: vmec_equilibrium.f:5819
siesta_equilibrium::siesta_pp_sxrem_b_id
integer, parameter siesta_pp_sxrem_b_id
1D Array of function sxrem profile parameters. Defined in pprofile_T.
Definition: siesta_equilibrium.f:98
siesta_equilibrium::siesta_write
subroutine siesta_write(this, iou)
Write out the equilibrium to an output file.
Definition: siesta_equilibrium.f:4354
siesta_equilibrium::siesta_get_gp_ne_pi
real(rprec) function siesta_get_gp_ne_pi(this, x_cart, i)
Get the electron density gp kernel value for the position and index.
Definition: siesta_equilibrium.f:1880
line_segment::line_seg
subroutine, public line_seg(x, y, xx, yy, n)
Interpolate a point on a line.
Definition: line_segment.f:41
vmec_equilibrium::vmec_destruct
subroutine vmec_destruct(this)
Deconstruct a vmec_class object.
Definition: vmec_equilibrium.f:845
siesta_equilibrium::siesta_pp_ti_as_id
integer, parameter siesta_pp_ti_as_id
1D Array of segment ion temperature profile s poitions. Defined in pprofile_T.
Definition: siesta_equilibrium.f:117
siesta_equilibrium::siesta_get_gp_sxrem_num_hyper_param
integer function siesta_get_gp_sxrem_num_hyper_param(this, index)
Get the number of soft x-ray emission gp kernel hyper parameters.
Definition: siesta_equilibrium.f:2534
siesta_equilibrium::siesta_construct
type(siesta_class) function, pointer siesta_construct(file_name, restart_file_name, ne, te, ti, sxrem, phi_offset, z_offset, pol_rad_ratio, iou, eq_comm, recon_comm, vmec, state_flags)
Construct a siesta_class object.
Definition: siesta_equilibrium.f:389
vmec_equilibrium::vmec_set_param
subroutine vmec_set_param(this, id, i_index, j_index, value, eq_comm, state_flags)
Sets the value of a reconstruction equilibrium parameter.
Definition: vmec_equilibrium.f:1013
vmec_equilibrium::vmec_get_param_value
real(rprec) function vmec_get_param_value(this, id, i_index, j_index)
Get the value of a reconstruction VMEC parameter.
Definition: vmec_equilibrium.f:2295
siesta_equilibrium::siesta_z_offset_id
integer, parameter siesta_z_offset_id
Plasma z offset. This is a parameter to allow changing the vertical shift of the plasma.
Definition: siesta_equilibrium.f:137
siesta_equilibrium::siesta_set_magnetic_cache_point
subroutine siesta_set_magnetic_cache_point(this, use_axi)
Set magnetic cache initialy.
Definition: siesta_equilibrium.f:868
magnetic_response
Defines the base class of the type magnetic_response_class.
Definition: magnetic_response.f:11
siesta_equilibrium::siesta_pp_ze_b_id
integer, parameter siesta_pp_ze_b_id
1D Array of function zeff profile parameters. Defined in pprofile_T.
Definition: siesta_equilibrium.f:123
vmec_equilibrium::vmec_is_recon_param
logical function vmec_is_recon_param(this, id)
Checks if a parameter id is a reconstruction parameter.
Definition: vmec_equilibrium.f:5113
siesta_equilibrium::siesta_pp_ti_af_id
integer, parameter siesta_pp_ti_af_id
1D Array of segment ion temperature profile f values. Defined in pprofile_T.
Definition: siesta_equilibrium.f:120
vmec_equilibrium::vmec_reset_state
subroutine vmec_reset_state(this)
Reset the internal state of the equilibrium.
Definition: vmec_equilibrium.f:5486
siesta_equilibrium::siesta_pp_te_as_id
integer, parameter siesta_pp_te_as_id
1D Array of segment electron temperature profile s poitions. Defined in pprofile_T.
Definition: siesta_equilibrium.f:108
siesta_context::siesta_context_class
Base class representing a siesta_context. This contains a copy of every variable that is needed to de...
Definition: siesta_context.f:54
vmec_equilibrium::magnetic_cache_vc_min_grid_points
integer, parameter magnetic_cache_vc_min_grid_points
Virtual casing grid points for vacuum field measurements.
Definition: vmec_equilibrium.f:313
siesta_equilibrium::siesta_get_sxrem_af
real(rprec) function, dimension(:), pointer siesta_get_sxrem_af(this, index)
Gets the soft x-ray emissivity profile af array.
Definition: siesta_equilibrium.f:2572
siesta_equilibrium::siesta_get_ne_radial
real(rprec) function siesta_get_ne_radial(this, s)
Gets the electron density at a s, u=0, v=0 position.
Definition: siesta_equilibrium.f:1994
vmec_equilibrium
Defines the base class of the type vmec_class. This module contains all the code necessary to interfa...
Definition: vmec_equilibrium.f:135
siesta_equilibrium::siesta_get_gp_ne_num_hyper_param
integer function siesta_get_gp_ne_num_hyper_param(this)
Get the number of electron density gp kernel hyper parameters.
Definition: siesta_equilibrium.f:1775
siesta_equilibrium::siesta_get_gp_te
Interface for the siesta guassian process electron temperature profile values.
Definition: siesta_equilibrium.f:302
siesta_equilibrium::siesta_set_param
subroutine siesta_set_param(this, id, i_index, j_index, value, eq_comm, state_flags)
Sets the value of a reconstruction equilibrium parameter.
Definition: siesta_equilibrium.f:658
siesta_equilibrium
Defines the base class of the type siesta_class. This module contains all the code necessary to inter...
Definition: siesta_equilibrium.f:74
siesta_equilibrium::siesta_sync_state
subroutine siesta_sync_state(this, recon_comm)
Syncronize the siesta state to children.
Definition: siesta_equilibrium.f:4634
vmec_equilibrium::vmec_is_scaler_value
logical function vmec_is_scaler_value(this, id)
Checks if a parameter id is a scaler value.
Definition: vmec_equilibrium.f:4979
siesta_equilibrium::siesta_get_magnetic_volume_jzgrid
real(rprec) function, dimension(:,:,:), pointer siesta_get_magnetic_volume_jzgrid(this)
Get volume magnetic volume integration jz grid points.
Definition: siesta_equilibrium.f:3283
siesta_equilibrium::siesta_get_con_surface_kzgrid
real(rprec) function, dimension(:,:), pointer siesta_get_con_surface_kzgrid(this)
Get the conducting surface integration kz grid points.
Definition: siesta_equilibrium.f:3411
siesta_equilibrium::siesta_get_int_b_plasma
real(rprec) function, dimension(3) siesta_get_int_b_plasma(this, flux)
Get internal plasma magnetic field.
Definition: siesta_equilibrium.f:3629
vmec_equilibrium::vmec_sync_state
subroutine vmec_sync_state(this, recon_comm)
Syncronize the vmec state to children.
Definition: vmec_equilibrium.f:6028
siesta_equilibrium::siesta_restart
subroutine siesta_restart(this, result_ncid, current_step)
Restart from a result file.
Definition: siesta_equilibrium.f:4601
siesta_equilibrium::siesta_get_param_name
character(len=data_name_length) function siesta_get_param_name(this, id)
Get the name of a reconstruction siesta parameter.
Definition: siesta_equilibrium.f:1693
siesta_equilibrium::siesta_sync_child
subroutine siesta_sync_child(this, index, recon_comm)
Syncronize a child siesta state to the parent.
Definition: siesta_equilibrium.f:4691
model_state::model_state_siesta_flag
integer, parameter model_state_siesta_flag
SIESTA Equilibrium changed bit position.
Definition: model_state.f:28
siesta_equilibrium::siesta_phi_offset_id
integer, parameter siesta_phi_offset_id
Plasma Phi offset. This is a parameter to allow changing the phi angle of a quasi helical state in an...
Definition: siesta_equilibrium.f:134
siesta_equilibrium::siesta_get_gp_te_pp
real(rprec) function siesta_get_gp_te_pp(this, x_cart, y_cart)
Get the electron temperature gp kernel value for the position and position.
Definition: siesta_equilibrium.f:2178
siesta_equilibrium::siesta_get_volume_int_element
real(rprec) function siesta_get_volume_int_element(this)
Get volume integration element.
Definition: siesta_equilibrium.f:3319
siesta_context::siesta_context_read
subroutine siesta_context_read(this, restart_file_name)
Read a restart file.
Definition: siesta_context.f:412
siesta_equilibrium::siesta_get_gp_ti_num_hyper_param
integer function siesta_get_gp_ti_num_hyper_param(this)
Get the number of ion temperature gp kernel hyper parameters.
Definition: siesta_equilibrium.f:2282
siesta_equilibrium::siesta_is_recon_param
logical function siesta_is_recon_param(this, id)
Checks if a parameter id is a reconstruction parameter.
Definition: siesta_equilibrium.f:4016
Generated on Thu Mar 5 2020 15:49:26 for V3FIT by   doxygen 1.8.17