save_simulation_parameters Subroutine

public subroutine save_simulation_parameters(params, spp, F, P)

This subroutine saves to the HDF5 file "simulation_parameters.h5" all the relevant simulation parameters of KORC, most of them being part of the input file, but also including some derived quantities from the input parameters. This file is intended to facilitate the post-processing of KORC data using any software that supports the HDF5 software.

Arguments

Type IntentOptional AttributesName
type(KORC_PARAMS), intent(in) :: params

Core KORC simulation parameters.

type(SPECIES), intent(in), DIMENSION(:), ALLOCATABLE:: spp

An instance of KORC's derived type SPECIES containing all the information of different electron species. See korc_types.

type(FIELDS), intent(in) :: F

An instance of KORC's derived type FIELDS containing all the information about the fields used in the simulation. See korc_types and korc_fields.

type(PROFILES), intent(in) :: P

An instance of KORC's derived type PROFILES containing all the information about the plasma profiles used in the simulation. See korc_types and korc_profiles.


Contents


Source Code

  subroutine save_simulation_parameters(params,spp,F,P)
    !! @note Subroutine to save to a HDF5 file all the relevant simulation
    !! parameters. @endnote
    !! This subroutine saves to the HDF5 file "<a>simulation_parameters.h5</a>"
    !! all the relevant simulation parameters of KORC, most of them being part
    !! of the input file, but also including some derived quantities from the
    !! input parameters. This file is intended to facilitate the
    !! post-processing of KORC data using any software that supports
    !! the HDF5 software.
    TYPE(KORC_PARAMS), INTENT(IN) 				:: params
    !!Core KORC simulation parameters.
    TYPE(SPECIES), DIMENSION(:), ALLOCATABLE, INTENT(IN) 	:: spp
    !! An instance of KORC's derived type SPECIES containing all
    !! the information of different electron species. See [[korc_types]].
    TYPE(FIELDS), INTENT(IN) 					:: F
    !! An instance of KORC's derived type FIELDS containing all the information
    !! about the fields used in the simulation. See [[korc_types]]
    !! and [[korc_fields]].
    TYPE(PROFILES), INTENT(IN) 					:: P
    !! An instance of KORC's derived type PROFILES containing all the
    !! information about the plasma profiles used in the simulation.
    !! See [[korc_types]] and [[korc_profiles]].
    CHARACTER(MAX_STRING_LENGTH) 				:: filename
    !! String containing the name of the HDF5 file.
    CHARACTER(MAX_STRING_LENGTH) 				:: gname
    !! String containing the group name of a set of KORC parameters.
    CHARACTER(MAX_STRING_LENGTH) 				:: dset
    !! Name of data set to be saved to file.
    INTEGER(HID_T) 						:: h5file_id
    !!  HDF5 file identifier.
    INTEGER(HID_T) 						:: group_id
    !! HDF5 group identifier.
    INTEGER(HSIZE_T), DIMENSION(:), ALLOCATABLE 		:: dims
    !! Dimensions of data saved to HDF5 file.
    REAL(rp), DIMENSION(:), ALLOCATABLE 			:: rdata
    !! 1-D array of real data to be saved to HDF5 file.
    INTEGER, DIMENSION(:), ALLOCATABLE 				:: idata
    !! 1-D array of integer data to be saved to HDF5 file.
    CHARACTER(MAX_STRING_LENGTH), DIMENSION(:), ALLOCATABLE :: attr_array
    !! An 1-D array with attributes of 1-D real or integer arrays that are
    !! passed to KORC interfaces of HDF5 I/O subroutines.
    CHARACTER(MAX_STRING_LENGTH) 				:: attr
    !!  A single attributes of real or integer data that is passed to KORC
    !! interfaces of HDF5 I/O subroutines.
    INTEGER 							:: h5error
    !! HDF5 error status.
    CHARACTER(19) 						:: tmp_str
    !! Temporary string used to manipulate various strings.
    REAL(rp) 							:: units
    !! Temporary variable used to add physical units to KORC parameters.

    ! * * * Error handling * * * !
    call h5eset_auto_f(params%HDF5_error_handling, h5error)
    ! Turn off: 0_idef. Turn on: 1_idef

    if (.NOT.(params%restart)) then

    if (params%mpi_params%rank .EQ. 0) then
       write(6,'("Saving simulations parameters")')
    end if

       
       if (SIZE(params%outputs_list).GT.1_idef) then
          write(tmp_str,'(I18)') params%mpi_params%rank
          filename = TRIM(params%path_to_outputs) // "file_"  &
               // TRIM(ADJUSTL(tmp_str)) // ".h5"
          call h5fcreate_f(TRIM(filename), H5F_ACC_TRUNC_F, h5file_id, h5error)
          call h5fclose_f(h5file_id, h5error)
       end if

       if (params%mpi_params%rank .EQ. 0) then
          filename = TRIM(params%path_to_outputs) // "simulation_parameters.h5"

          call h5fcreate_f(TRIM(filename), H5F_ACC_TRUNC_F, h5file_id, h5error)

          ! Simulation parameters group
          gname = "simulation"
          call h5gcreate_f(h5file_id, TRIM(gname), group_id, h5error)

          ALLOCATE(attr_array(1))
          ALLOCATE(idata(1))

          dset = TRIM(gname) // "/field_model"
          call save_string_parameter(h5file_id,dset,(/params%field_model/))

          dset = TRIM(gname) // "/profile_model"
          call save_string_parameter(h5file_id,dset,(/params%profile_model/))

          dset = TRIM(gname) // "/simulation_time"
          attr = "Total aimed simulation time in seconds"
          call save_to_hdf5(h5file_id,dset,params%simulation_time* &
               params%cpp%time,attr)

          dset = TRIM(gname) // "/snapshot_frequency"
          attr = "Time between snapshots in seconds"
          call save_to_hdf5(h5file_id,dset,params%snapshot_frequency* &
               params%cpp%time,attr)

          dset = TRIM(gname) // "/dt"
          attr = "Time step in secs"
          call save_to_hdf5(h5file_id,dset,params%dt*params%cpp%time,attr)

          dset = TRIM(gname) // "/t_steps"
          attr_array(1) = "Number of time steps"
          idata = params%t_steps
          call save_1d_array_to_hdf5(h5file_id,dset,idata,attr_array)

          dset = TRIM(gname) // "/num_omp_threads"
          attr = "Number of omp threads"
          call save_to_hdf5(h5file_id,dset, params%num_omp_threads,attr)

          dset = TRIM(gname) // "/output_cadence"
          attr_array(1) = "Cadence of output files"
          idata = params%output_cadence
          call save_1d_array_to_hdf5(h5file_id,dset,idata,attr_array)

          dset = TRIM(gname) // "/HDF5_error_handling"
          attr_array(1) = "Error handling option: 0=OFF, 1=ON"
          idata = params%HDF5_error_handling
          call save_1d_array_to_hdf5(h5file_id,dset,idata,attr_array)

          dset = TRIM(gname) // "/restart_output_cadence"
          attr_array(1) = "Cadence of output files"
          idata = params%restart_output_cadence
          call save_1d_array_to_hdf5(h5file_id,dset,idata,attr_array)

          dset = TRIM(gname) // "/num_snapshots"
          attr_array(1) = "Number of outputs for each variable"
          idata = params%num_snapshots
          call save_1d_array_to_hdf5(h5file_id,dset,idata,attr_array)

          dset = TRIM(gname) // "/num_species"
          attr = "Number of particle species"
          call save_to_hdf5(h5file_id,dset,params%num_species,attr)

          dset = TRIM(gname) // "/nmpi"
          attr = "Number of mpi processes"
          call save_to_hdf5(h5file_id,dset,params%mpi_params%nmpi,attr)

          dset = TRIM(gname) // "/minimum_particle_energy"
          attr = "Minimum energy of simulated particles in eV"
          call save_to_hdf5(h5file_id,dset,params%minimum_particle_energy* &
               params%cpp%energy/C_E,attr)

          dset = TRIM(gname) // "/minimum_particle_g"
          attr = "Minimum relativistic factor gamma of simulated particles"
          call save_to_hdf5(h5file_id,dset,params%minimum_particle_g,attr)

          dset = TRIM(gname) // "/radiation"
          attr = "Radiation losses included in simulation"
          if(params%radiation) then
             call save_to_hdf5(h5file_id,dset,1_idef,attr)
          else
             call save_to_hdf5(h5file_id,dset,0_idef,attr)
          end if

          dset = TRIM(gname) // "/collisions"
          attr = "Collisions included in simulation"
          if(params%collisions) then
             call save_to_hdf5(h5file_id,dset,1_idef,attr)
          else
             call save_to_hdf5(h5file_id,dset,0_idef,attr)
          end if

          dset = TRIM(gname) // "/outputs_list"
          call save_string_parameter(h5file_id,dset,params%outputs_list)

          dset = TRIM(gname) // "/orbit_model"
          call save_string_parameter(h5file_id,dset,(/params%orbit_model/))

          dset = TRIM(gname) // "/field_eval"
          call save_string_parameter(h5file_id,dset,(/params%field_eval/))

          DEALLOCATE(idata)
          DEALLOCATE(attr_array)

          call h5gclose_f(group_id, h5error)


          ! Plasma species group
          gname = "species"
          call h5gcreate_f(h5file_id, TRIM(gname), group_id, h5error)

          ALLOCATE(attr_array(params%num_species))

          dset = TRIM(gname) // "/spatial_distribution"
          call save_string_parameter(h5file_id,dset,spp%spatial_distribution)

          dset = TRIM(gname) // "/energy_distribution"
          call save_string_parameter(h5file_id,dset,spp%energy_distribution)

          dset = TRIM(gname) // "/pitch_distribution"
          call save_string_parameter(h5file_id,dset,spp%pitch_distribution)

          dset = TRIM(gname) // "/ppp"
          attr_array(1) = "Particles per (mpi) process"
          call save_1d_array_to_hdf5(h5file_id,dset,spp%ppp,attr_array)

          dset = TRIM(gname) // "/q"
          attr_array(1) = "Electric charge"
          call save_1d_array_to_hdf5(h5file_id,dset,spp%q* &
               params%cpp%charge,attr_array)

          dset = TRIM(gname) // "/m"
          attr_array(1) = "Species mass in kg"
          call save_1d_array_to_hdf5(h5file_id,dset,spp%m* &
               params%cpp%mass,attr_array)

          dset = TRIM(gname) // "/Eo"
          attr_array(1) = "Initial (average) energy in eV"
          call save_1d_array_to_hdf5(h5file_id,dset,spp%Eo* &
               params%cpp%energy/C_E,attr_array)

          dset = TRIM(gname) // "/go"
          attr_array(1) = "Initial relativistic g factor."
          call save_1d_array_to_hdf5(h5file_id,dset,spp%go,attr_array)

          dset = TRIM(gname) // "/etao"
          attr_array(1) = "Initial pitch angle in degrees"
          call save_1d_array_to_hdf5(h5file_id,dset,spp%etao,attr_array)

          dset = TRIM(gname) // "/wc"
          attr_array(1) = "Average relativistic cyclotron frequency in Hz"
          call save_1d_array_to_hdf5(h5file_id,dset,spp%wc_r/ &
               params%cpp%time,attr_array)

          dset = TRIM(gname) // "/Ro"
          attr_array(1) = "Initial radial position of population"
          call save_1d_array_to_hdf5(h5file_id,dset,spp%Ro* &
               params%cpp%length,attr_array)

          dset = TRIM(gname) // "/PHIo"
          attr_array(1) = "Azimuthal angle in degrees."
          call save_1d_array_to_hdf5(h5file_id,dset,spp%PHIo* &
               180.0_rp/C_PI,attr_array)

          dset = TRIM(gname) // "/Zo"
          attr_array(1) = "Initial Z position of population"
          call save_1d_array_to_hdf5(h5file_id,dset,spp%Zo* &
               params%cpp%length,attr_array)

          dset = TRIM(gname) // "/ri"
          attr_array(1) = "Inner radius of initial spatial distribution"
          call save_1d_array_to_hdf5(h5file_id,dset,spp%r_inner* &
               params%cpp%length,attr_array)

          dset = TRIM(gname) // "/ro"
          attr_array(1) = "Outter radius of initial spatial distribution"
          call save_1d_array_to_hdf5(h5file_id,dset,spp%r_outter* &
               params%cpp%length,attr_array)

          dset = TRIM(gname) // "/falloff_rate"
          attr_array(1) = "Falloff of gaussian or exponential radial &
               profile in m"
          call save_1d_array_to_hdf5(h5file_id,dset,spp%falloff_rate/ &
               params%cpp%length,attr_array)

          dset = TRIM(gname) // "/shear_factor"
          attr_array(1) = "Shear factor (in case ELLIPTIC-TORUS  &
               spatial distribution is used."
          call save_1d_array_to_hdf5(h5file_id,dset,spp%shear_factor, &
               attr_array)

          dset = TRIM(gname) // "/sigmaR"
          attr_array(1) = "Variance of first dimension of 2D spatial & 
               distribution."
          call save_1d_array_to_hdf5(h5file_id,dset, &
               spp%sigmaR*params%cpp%length,attr_array)

          dset = TRIM(gname) // "/sigmaZ"
          attr_array(1) = "Variance of second dimension of 2D spatial &
               distribution."
          call save_1d_array_to_hdf5(h5file_id,dset, &
               spp%sigmaZ*params%cpp%length,attr_array)

          dset = TRIM(gname) // "/theta_gauss"
          attr_array(1) = "Angle of rotation of 2D spatial distribution."
          call save_1d_array_to_hdf5(h5file_id,dset,spp%theta_gauss,attr_array)

          dset = TRIM(gname) // "/psi_max"
          attr_array(1) = "Indicator function level of the argument of &
               the 2D gaussian exponential."
          call save_1d_array_to_hdf5(h5file_id,dset,spp%psi_max,attr_array)

          dset = TRIM(gname) // "/dth"
          attr_array(1) = "Variance of sampling normal variate for pitch angle."
          call save_1d_array_to_hdf5(h5file_id,dset,spp%dth,attr_array)

          dset = TRIM(gname) // "/dgam"
          attr_array(1) = "Variance of sampling normal variate for gamma."
          call save_1d_array_to_hdf5(h5file_id,dset,spp%dgam,attr_array)

          dset = TRIM(gname) // "/dR"
          attr_array(1) = "Variance of sampling normal variate for R."
          call save_1d_array_to_hdf5(h5file_id,dset,spp%dR*params%cpp%length, &
               attr_array)

          dset = TRIM(gname) // "/dZ"
          attr_array(1) = "Variance of sampling normal variate for Z."
          call save_1d_array_to_hdf5(h5file_id,dset,spp%dZ*params%cpp%length, &
               attr_array)
          
          call h5gclose_f(group_id, h5error)

          DEALLOCATE(attr_array)


          ! Plasma profiles group
!          if (params%collisions) then
             gname = "profiles"
             call h5gcreate_f(h5file_id, TRIM(gname), group_id, h5error)

             dset = TRIM(gname) // "/density_profile"
             call save_string_parameter(h5file_id,dset,(/P%ne_profile/))

             dset = TRIM(gname) // "/temperature_profile"
             call save_string_parameter(h5file_id,dset,(/P%Te_profile/))

             dset = TRIM(gname) // "/Zeff_profile"
             call save_string_parameter(h5file_id,dset,(/P%Zeff_profile/))

             dset = TRIM(gname) // "/neo"
             attr = "Density at the magnetic axis (m^-3)"
             call save_to_hdf5(h5file_id,dset,P%neo*params%cpp%density,attr)

             dset = TRIM(gname) // "/Teo"
             attr = "Temperature at the magnetic axis (eV)"
             call save_to_hdf5(h5file_id,dset,P%Teo* &
                  params%cpp%temperature/C_E,attr)

             dset = TRIM(gname) // "/Zeffo"
             attr = "Zeff at the magnetic axis"
             call save_to_hdf5(h5file_id,dset,P%Zeffo,attr)

             if (TRIM(params%profile_model) .EQ. 'ANALYTICAL') then
                dset = TRIM(gname) // "/n_ne"
                attr = "Exponent of tanh(x)^n for density profile"
                call save_to_hdf5(h5file_id,dset,P%n_ne,attr)

                dset = TRIM(gname) // "/a_ne"
                attr = "Coefficients f=ao+a1*r+a2*r^2+a3*r^3.  &
                     a_ne=[a0,a1,a2,a3]"
                call save_1d_array_to_hdf5(h5file_id,dset,P%a_ne)

                dset = TRIM(gname) // "/n_Te"
                attr = "Exponent of tanh(x)^n for density profile"
                call save_to_hdf5(h5file_id,dset,P%n_Te,attr)

                dset = TRIM(gname) // "/a_Te"
                attr = "Coefficients f=ao+a1*r+a2*r^2+a3*r^3.  &
                     a_Te=[a0,a1,a2,a3]"
                call save_1d_array_to_hdf5(h5file_id,dset,P%a_Te)

                dset = TRIM(gname) // "/n_Zeff"
                attr = "Exponent of tanh(x)^n for Zeff profile"
                call save_to_hdf5(h5file_id,dset,P%n_Zeff,attr)

                dset = TRIM(gname) // "/a_Zeff"
                attr = "Coefficients f=ao+a1*r+a2*r^2+a3*r^3.  &
                     a_Zeff=[a0,a1,a2,a3]"
                call save_1d_array_to_hdf5(h5file_id,dset,P%a_Zeff)

                if  (params%field_eval.EQ.'interp') then

                   ALLOCATE(attr_array(1))
                   dset = TRIM(gname) // "/dims"
                   attr_array(1) = "Mesh dimension of the profile (NR,NPHI,NZ)"
                   call save_1d_array_to_hdf5(h5file_id,dset,F%dims,attr_array)

                   dset = TRIM(gname) // "/R"
                   attr_array(1) = "Radial position of the magnetic field grid nodes"
                   call save_1d_array_to_hdf5(h5file_id,dset, &
                        F%X%R*params%cpp%length,attr_array)

                   if (ALLOCATED(F%X%PHI)) then
                      dset = TRIM(gname) // "/PHI"
                      attr_array(1) = "Azimuthal angle of the magnetic &
                           field grid nodes"
                      call save_1d_array_to_hdf5(h5file_id,dset,F%X%PHI,attr_array)
                   end if
                   
                   dset = TRIM(gname) // "/Z"
                   attr_array(1) = "Z position of the magnetic field grid nodes"
                   call save_1d_array_to_hdf5(h5file_id,dset,F%X%Z* &
                        params%cpp%length,attr_array)

                   dset = TRIM(gname) // "/ne"
                   units = params%cpp%density
                   call rsave_2d_array_to_hdf5(h5file_id, dset, &
                        units*P%ne_2D)

                   dset = TRIM(gname) // "/Te"
                   units = params%cpp%temperature
                   call rsave_2d_array_to_hdf5(h5file_id, dset, &
                        units*P%Te_2D)

                   dset = TRIM(gname) // "/Zeff"
                   call rsave_2d_array_to_hdf5(h5file_id, dset, &
                        P%Zeff_2D)

                   DEALLOCATE(attr_array)
                end if
                
             else if (params%profile_model(1:8) .EQ. 'EXTERNAL') then
                ALLOCATE(attr_array(1))
                dset = TRIM(gname) // "/dims"
                attr_array(1) = "Mesh dimension of the profiles (NR,NPHI,NZ)"
                call save_1d_array_to_hdf5(h5file_id,dset,P%dims,attr_array)

                dset = TRIM(gname) // "/R"
                attr_array(1) = "Grid nodes of profiles along the &
                     radial position"
                call save_1d_array_to_hdf5(h5file_id,dset,P%X%R* &
                     params%cpp%length,attr_array)

                if (ALLOCATED(F%X%PHI)) then
                   dset = TRIM(gname) // "/PHI"
                   attr_array(1) = "Grid nodes of profiles along the &
                        azimuthal position"
                   call save_1d_array_to_hdf5(h5file_id,dset, &
                        P%X%PHI,attr_array)
                end if

                dset = TRIM(gname) // "/Z"
                attr_array(1) = "Grid nodes of profiles along the Z position"
                call save_1d_array_to_hdf5(h5file_id,dset, &
                     P%X%Z*params%cpp%length,attr_array)

                dset = TRIM(gname) // "/ne"
                units = params%cpp%density
                call rsave_2d_array_to_hdf5(h5file_id, dset, &
                     units*P%ne_2D)

                dset = TRIM(gname) // "/Te"
                units = params%cpp%temperature
                call rsave_2d_array_to_hdf5(h5file_id, dset, &
                     units*P%Te_2D)

                dset = TRIM(gname) // "/Zeff"
                call rsave_2d_array_to_hdf5(h5file_id, dset, &
                     P%Zeff_2D)

                DEALLOCATE(attr_array)
             else if (params%profile_model .EQ. 'UNIFORM') then
                ! Something
             end if

             call h5gclose_f(group_id, h5error)
          !end if


          ! Electromagnetic fields group

          gname = "fields"
          call h5gcreate_f(h5file_id, TRIM(gname), group_id, h5error)

          if (TRIM(params%field_model(1:10)) .EQ. 'ANALYTICAL') then
             dset = TRIM(gname) // "/Bo"
             attr = "Toroidal field at the magnetic axis in T"
             call save_to_hdf5(h5file_id,dset,F%Bo*params%cpp%Bo,attr)

             dset = TRIM(gname) // "/current_direction"
             call save_string_parameter(h5file_id,dset, &
                  (/F%AB%current_direction/))

             dset = TRIM(gname) // "/a"
             attr = "Minor radius in m"
             call save_to_hdf5(h5file_id,dset,F%AB%a*params%cpp%length,attr)

             dset = TRIM(gname) // "/Ro"
             attr = "Magnetic axis radial position"
             call save_to_hdf5(h5file_id,dset,F%Ro*params%cpp%length,attr)

             dset = TRIM(gname) // "/Zo"
             attr = "Magnetic axis vertical position"
             call save_to_hdf5(h5file_id,dset,F%Zo*params%cpp%length,attr)
             
             dset = TRIM(gname) // "/qa"
             attr = "Safety factor at minor radius"
             call save_to_hdf5(h5file_id,dset,F%AB%qa,attr)

             dset = TRIM(gname) // "/qo"
             attr = "Safety factor at the magnetic axis"
             call save_to_hdf5(h5file_id,dset,F%AB%qo,attr)

             dset = TRIM(gname) // "/lambda"
             attr = "Parameter lamda in m"
             call save_to_hdf5(h5file_id,dset,F%AB%lambda* &
                  params%cpp%length,attr)

             dset = TRIM(gname) // "/Bpo"
             attr = "Poloidal magnetic field in T"
             call save_to_hdf5(h5file_id,dset,F%AB%Bpo*params%cpp%Bo,attr)
             
             dset = TRIM(gname) // "/Eo"
             attr = "Electric field at the magnetic axis in V/m"
             call save_to_hdf5(h5file_id,dset,F%Eo*params%cpp%Eo,attr)

             if (params%SC_E) then
                dset = TRIM(gname) // "/dt_E_SC"
                attr = "Time step for self-consistent E calculation"
                call save_to_hdf5(h5file_id,dset,F%dt_E_SC,attr)

                dset = TRIM(gname) // "/Ip_exp"
                attr = "Scaling for self-consistent current density"
                call save_to_hdf5(h5file_id,dset,F%Ip_exp,attr)

                ALLOCATE(attr_array(1))
                dset = TRIM(gname) // "/r_1D"
                attr_array(1) = "1D minor radial mesh for &
                     self-consistent fields"
                call save_1d_array_to_hdf5(h5file_id,dset,F%r_1D,attr_array)
                DEALLOCATE(attr_array)
             end if
             
             if  (params%field_eval.EQ.'interp') then

                ALLOCATE(attr_array(1))
                dset = TRIM(gname) // "/dims"
                attr_array(1) = "Mesh dimension of the magnetic  &
                     field (NR,NPHI,NZ)"
                call save_1d_array_to_hdf5(h5file_id,dset,F%dims,attr_array)

                dset = TRIM(gname) // "/R"
                attr_array(1) = "Radial position of the magnetic field grid nodes"
                call save_1d_array_to_hdf5(h5file_id,dset, &
                     F%X%R*params%cpp%length,attr_array)

                if (ALLOCATED(F%X%PHI)) then
                   dset = TRIM(gname) // "/PHI"
                   attr_array(1) = "Azimuthal angle of the magnetic &
                        field grid nodes"
                   call save_1d_array_to_hdf5(h5file_id,dset,F%X%PHI,attr_array)
                end if

                dset = TRIM(gname) // "/Z"
                attr_array(1) = "Z position of the magnetic field grid nodes"
                call save_1d_array_to_hdf5(h5file_id,dset,F%X%Z* &
                     params%cpp%length,attr_array)

                if (ALLOCATED(F%PSIp)) then
                   dset = TRIM(gname) // "/psi_p"
                   units = params%cpp%Bo*params%cpp%length**2
                   call rsave_2d_array_to_hdf5(h5file_id, dset, &
                        units*F%PSIp)
                end if
                
                if(params%field_model(12:13).eq.'2D') then
                
                   dset = TRIM(gname) // "/BR"
                   units = params%cpp%Bo
                   call rsave_2d_array_to_hdf5(h5file_id, dset, &
                        units*F%B_2D%R)

                   dset = TRIM(gname) // "/BPHI"
                   units = params%cpp%Bo
                   call rsave_2d_array_to_hdf5(h5file_id, dset, &
                        units*F%B_2D%PHI)

                   dset = TRIM(gname) // "/BZ"
                   units = params%cpp%Bo
                   call rsave_2d_array_to_hdf5(h5file_id, dset, &
                        units*F%B_2D%Z)
                   
                   if (ALLOCATED(F%FLAG2D)) then
                      dset = TRIM(gname) // "/Flag"
                      call rsave_2d_array_to_hdf5(h5file_id, dset, &
                           F%FLAG2D)
                   end if
                
                   if  (params%orbit_model(3:5).EQ.'pre') then

                      dset = TRIM(gname) // "/gradBR"
                      units = params%cpp%Bo/params%cpp%length
                      call rsave_2d_array_to_hdf5(h5file_id, dset, &
                           units*F%gradB_2D%R)

                      dset = TRIM(gname) // "/gradBPHI"
                      units = params%cpp%Bo/params%cpp%length
                      call rsave_2d_array_to_hdf5(h5file_id, dset, &
                           units*F%gradB_2D%PHI)

                      dset = TRIM(gname) // "/gradBZ"
                      units = params%cpp%Bo/params%cpp%length
                      call rsave_2d_array_to_hdf5(h5file_id, dset, &
                           units*F%gradB_2D%Z)

                      dset = TRIM(gname) // "/curlbR"
                      units = 1./params%cpp%length
                      call rsave_2d_array_to_hdf5(h5file_id, dset, &
                           units*F%curlb_2D%R)

                      dset = TRIM(gname) // "/curlbPHI"
                      units = 1./params%cpp%length
                      call rsave_2d_array_to_hdf5(h5file_id, dset, &
                           units*F%curlb_2D%PHI)

                      dset = TRIM(gname) // "/curlbZ"
                      units =1./params%cpp%length
                      call rsave_2d_array_to_hdf5(h5file_id, dset, &
                           units*F%curlb_2D%Z)

                   end if

                else if(params%field_model(12:13).eq.'3D') then
                
                   dset = TRIM(gname) // "/BR"
                   units = params%cpp%Bo
                   call rsave_3d_array_to_hdf5(h5file_id, dset, &
                        units*F%B_3D%R)

                   dset = TRIM(gname) // "/BPHI"
                   units = params%cpp%Bo
                   call rsave_3d_array_to_hdf5(h5file_id, dset, &
                        units*F%B_3D%PHI)

                   dset = TRIM(gname) // "/BZ"
                   units = params%cpp%Bo
                   call rsave_3d_array_to_hdf5(h5file_id, dset, &
                        units*F%B_3D%Z)
                   
                   if (ALLOCATED(F%FLAG3D)) then
                      dset = TRIM(gname) // "/Flag"
                      call rsave_3d_array_to_hdf5(h5file_id, dset, &
                           F%FLAG3D)
                   end if
                
                   if  (params%orbit_model(3:5).EQ.'pre') then

                      dset = TRIM(gname) // "/gradBR"
                      units = params%cpp%Bo/params%cpp%length
                      call rsave_3d_array_to_hdf5(h5file_id, dset, &
                           units*F%gradB_3D%R)

                      dset = TRIM(gname) // "/gradBPHI"
                      units = params%cpp%Bo/params%cpp%length
                      call rsave_3d_array_to_hdf5(h5file_id, dset, &
                           units*F%gradB_3D%PHI)

                      dset = TRIM(gname) // "/gradBZ"
                      units = params%cpp%Bo/params%cpp%length
                      call rsave_3d_array_to_hdf5(h5file_id, dset, &
                           units*F%gradB_3D%Z)

                      dset = TRIM(gname) // "/curlbR"
                      units = 1./params%cpp%length
                      call rsave_3d_array_to_hdf5(h5file_id, dset, &
                           units*F%curlb_3D%R)

                      dset = TRIM(gname) // "/curlbPHI"
                      units = 1./params%cpp%length
                      call rsave_3d_array_to_hdf5(h5file_id, dset, &
                           units*F%curlb_3D%PHI)

                      dset = TRIM(gname) // "/curlbZ"
                      units =1./params%cpp%length
                      call rsave_3d_array_to_hdf5(h5file_id, dset, &
                           units*F%curlb_3D%Z)

                   end if
                end if
                   
                DEALLOCATE(attr_array)
             end if

          else if (params%field_model(1:8) .EQ. 'EXTERNAL') then
             ALLOCATE(attr_array(1))
             
             dset = TRIM(gname) // "/dims"
             attr_array(1) = "Mesh dimension of the magnetic  &
                  field (NR,NPHI,NZ)"
             call save_1d_array_to_hdf5(h5file_id,dset,F%dims,attr_array)

             dset = TRIM(gname) // "/R"
             attr_array(1) = "Radial position of the magnetic field grid nodes"
             call save_1d_array_to_hdf5(h5file_id,dset, &
                  F%X%R*params%cpp%length,attr_array)

             if (ALLOCATED(F%X%PHI)) then
                if (F%Dim2x1t) then
                   dset = TRIM(gname) // "/PHI"
                   attr_array(1) = "Azimuthal angle of the magnetic &
                        field grid nodes"
                   call save_1d_array_to_hdf5(h5file_id,dset, &
                        F%X%PHI*params%cpp%time,attr_array)
                else
                   dset = TRIM(gname) // "/PHI"
                   attr_array(1) = "Azimuthal angle of the magnetic &
                        field grid nodes"
                   call save_1d_array_to_hdf5(h5file_id,dset,F%X%PHI,attr_array)
                end if
             end if

             dset = TRIM(gname) // "/Z"
             attr_array(1) = "Z position of the magnetic field grid nodes"
             call save_1d_array_to_hdf5(h5file_id,dset,F%X%Z* &
                  params%cpp%length,attr_array)

             dset = TRIM(gname) // "/Bo"
             attr = "Toroidal field at the magnetic axis in T"
             call save_to_hdf5(h5file_id,dset,F%Bo*params%cpp%Bo,attr)

             dset = TRIM(gname) // "/Eo"
             attr = "Electric field at the magnetic axis in V/m"
             call save_to_hdf5(h5file_id,dset,F%Eo*params%cpp%Eo,attr)

             dset = TRIM(gname) // "/E_dyn"
             attr = "Magnitude of dynamic E"
             call save_to_hdf5(h5file_id,dset,F%E_dyn*params%cpp%Eo,attr)
             
             dset = TRIM(gname) // "/E_pulse"
             attr = "Magnitude of dynamic E"
             call save_to_hdf5(h5file_id,dset,F%E_pulse*params%cpp%time,attr)

             dset = TRIM(gname) // "/E_width"
             attr = "Magnitude of dynamic E"
             call save_to_hdf5(h5file_id,dset,F%E_width*params%cpp%time,attr)
             
             dset = TRIM(gname) // "/Ro"
             attr = "Radial position of magnetic axis"
             call save_to_hdf5(h5file_id,dset,F%Ro*params%cpp%length,attr)

             dset = TRIM(gname) // "/Zo"
             attr = "Radial position of magnetic axis"
             call save_to_hdf5(h5file_id,dset,F%Zo*params%cpp%length,attr)

             dset = TRIM(gname) // "/Axisymmetric"
             attr = "Radial position of magnetic axis"
             if(F%axisymmetric_fields) then
                call save_to_hdf5(h5file_id,dset,1_idef,attr)
             else
                call save_to_hdf5(h5file_id,dset,0_idef,attr)
             end if

             if (ALLOCATED(F%PSIp)) then
                dset = TRIM(gname) // "/psi_p"
                units = params%cpp%Bo*params%cpp%length**2
                call rsave_2d_array_to_hdf5(h5file_id, dset, &
                     units*F%PSIp)
             end if

             if (ALLOCATED(F%E_3D%R)) then
                dset = TRIM(gname) // "/ER3D"
                units = params%cpp%Eo
                call rsave_3d_array_to_hdf5(h5file_id, dset, &
                     units*F%E_3D%R)
             end if
             
             if (ALLOCATED(F%E_3D%PHI)) then
                dset = TRIM(gname) // "/EPHI3D"
                units = params%cpp%Eo
                call rsave_3d_array_to_hdf5(h5file_id, dset, &
                     units*F%E_3D%PHI)
             end if
             
             if (ALLOCATED(F%PSIp3D)) then
                dset = TRIM(gname) // "/psi_p3D"
                units = params%cpp%Bo*params%cpp%length**2
                call rsave_3d_array_to_hdf5(h5file_id, dset, &
                     units*F%PSIp3D)
             end if

             if (ALLOCATED(F%FLAG2D)) then
                dset = TRIM(gname) // "/Flag2D"
                call rsave_2d_array_to_hdf5(h5file_id, dset, &
                     F%FLAG2D)
             end if

             if (ALLOCATED(F%FLAG3D)) then
                dset = TRIM(gname) // "/Flag3D"
                call rsave_3d_array_to_hdf5(h5file_id, dset, &
                     F%FLAG3D)
             end if


             if (params%SC_E) then
                dset = TRIM(gname) // "/dt_E_SC"
                attr = "Time step for self-consistent E calculation"
                call save_to_hdf5(h5file_id,dset,F%dt_E_SC,attr)

                dset = TRIM(gname) // "/Ip_exp"
                attr = "Scaling for self-consistent current density"
                call save_to_hdf5(h5file_id,dset,F%Ip_exp,attr)

                dset = TRIM(gname) // "/PSIP_1D"
                attr_array(1) = "1D minor radial mesh for &
                     self-consistent fields"
                call save_1d_array_to_hdf5(h5file_id,dset,F%PSIP_1D,attr_array)
             end if
             
             if  (F%axisymmetric_fields.and. &
                  .not.(params%field_model(10:12).eq.'PSI')) then

                dset = TRIM(gname) // "/BR"
                units = params%cpp%Bo
                call rsave_2d_array_to_hdf5(h5file_id, dset, &
                     units*F%B_2D%R)

                dset = TRIM(gname) // "/BPHI"
                units = params%cpp%Bo
                call rsave_2d_array_to_hdf5(h5file_id, dset, &
                     units*F%B_2D%PHI)

                dset = TRIM(gname) // "/BZ"
                units = params%cpp%Bo
                call rsave_2d_array_to_hdf5(h5file_id, dset, &
                     units*F%B_2D%Z)

                if  (params%orbit_model(3:5).EQ.'pre') then

                   dset = TRIM(gname) // "/gradBR"
                   units = params%cpp%Bo/params%cpp%length
                   call rsave_2d_array_to_hdf5(h5file_id, dset, &
                        units*F%gradB_2D%R)

                   dset = TRIM(gname) // "/gradBPHI"
                   units = params%cpp%Bo/params%cpp%length
                   call rsave_2d_array_to_hdf5(h5file_id, dset, &
                        units*F%gradB_2D%PHI)

                   dset = TRIM(gname) // "/gradBZ"
                   units = params%cpp%Bo/params%cpp%length
                   call rsave_2d_array_to_hdf5(h5file_id, dset, &
                        units*F%gradB_2D%Z)

                   dset = TRIM(gname) // "/curlbR"
                   units = 1./params%cpp%length
                   call rsave_2d_array_to_hdf5(h5file_id, dset, &
                        units*F%curlb_2D%R)

                   dset = TRIM(gname) // "/curlbPHI"
                   units = 1./params%cpp%length
                   call rsave_2d_array_to_hdf5(h5file_id, dset, &
                        units*F%curlb_2D%PHI)

                   dset = TRIM(gname) // "/curlbZ"
                   units = 1./params%cpp%length
                   call rsave_2d_array_to_hdf5(h5file_id, dset, &
                        units*F%curlb_2D%Z)

                end if
                
             else if ((.not.F%axisymmetric_fields).and. &
                  .not.(params%field_model(10:12).eq.'PSI')) then

                dset = TRIM(gname) // "/BR"
                units = params%cpp%Bo
                call rsave_3d_array_to_hdf5(h5file_id, dset, &
                     units*F%B_3D%R)

                dset = TRIM(gname) // "/BPHI"
                units = params%cpp%Bo
                call rsave_3d_array_to_hdf5(h5file_id, dset, &
                     units*F%B_3D%PHI)

                dset = TRIM(gname) // "/BZ"
                units = params%cpp%Bo
                call rsave_3d_array_to_hdf5(h5file_id, dset, &
                     units*F%B_3D%Z)

                if  (params%orbit_model(3:5).EQ.'pre') then

                   dset = TRIM(gname) // "/gradBR"
                   units = params%cpp%Bo/params%cpp%length
                   call rsave_3d_array_to_hdf5(h5file_id, dset, &
                        units*F%gradB_3D%R)

                   dset = TRIM(gname) // "/gradBPHI"
                   units = params%cpp%Bo/params%cpp%length
                   call rsave_3d_array_to_hdf5(h5file_id, dset, &
                        units*F%gradB_3D%PHI)

                   dset = TRIM(gname) // "/gradBZ"
                   units = params%cpp%Bo/params%cpp%length
                   call rsave_3d_array_to_hdf5(h5file_id, dset, &
                        units*F%gradB_3D%Z)

                   dset = TRIM(gname) // "/curlbR"
                   units = 1./params%cpp%length
                   call rsave_3d_array_to_hdf5(h5file_id, dset, &
                        units*F%curlb_3D%R)

                   dset = TRIM(gname) // "/curlbPHI"
                   units = 1./params%cpp%length
                   call rsave_3d_array_to_hdf5(h5file_id, dset, &
                        units*F%curlb_3D%PHI)

                   dset = TRIM(gname) // "/curlbZ"
                   units = 1./params%cpp%length
                   call rsave_3d_array_to_hdf5(h5file_id, dset, &
                        units*F%curlb_3D%Z)

                end if
             end if
             
             DEALLOCATE(attr_array)
          else if (params%field_model .EQ. 'UNIFORM') then
             dset = TRIM(gname) // "/Bo"
             attr = "Magnetic field in T"
             call save_to_hdf5(h5file_id,dset,F%Bo*params%cpp%Bo,attr)

             dset = TRIM(gname) // "/Eo"
             attr = "Electric field in V/m"
             call save_to_hdf5(h5file_id,dset,F%Eo*params%cpp%Eo,attr)
          end if

          call h5gclose_f(group_id, h5error)


          ! Characteristic scales
          gname = "scales"
          call h5gcreate_f(h5file_id, TRIM(gname), group_id, h5error)

          dset = TRIM(gname) // "/t"
          attr = "Characteristic time in secs"
          call save_to_hdf5(h5file_id,dset,params%cpp%time,attr)

          dset = TRIM(gname) // "/m"
          attr = "Characteristic mass in kg"
          call save_to_hdf5(h5file_id,dset,params%cpp%mass,attr)

          dset = TRIM(gname) // "/q"
          attr = "Characteristic charge in Coulombs"
          call save_to_hdf5(h5file_id,dset,params%cpp%charge,attr)

          dset = TRIM(gname) // "/l"
          attr = "Characteristic length in m"
          call save_to_hdf5(h5file_id,dset,params%cpp%length,attr)

          dset = TRIM(gname) // "/v"
          attr = "Characteristic velocity in m"
          call save_to_hdf5(h5file_id,dset,params%cpp%velocity,attr)

          dset = TRIM(gname) // "/K"
          attr = "Characteristic kinetic energy in J"
          call save_to_hdf5(h5file_id,dset,params%cpp%energy,attr)

          dset = TRIM(gname) // "/n"
          attr = "Characteristic plasma density in m^-3"
          call save_to_hdf5(h5file_id,dset,params%cpp%density,attr)

          dset = TRIM(gname) // "/E"
          attr = "Characteristic electric field in V/m"
          call save_to_hdf5(h5file_id,dset,params%cpp%Eo,attr)

          dset = TRIM(gname) // "/B"
          attr = "Characteristic magnetic field in T"
          call save_to_hdf5(h5file_id,dset,params%cpp%Bo,attr)

          dset = TRIM(gname) // "/P"
          attr = "Characteristic pressure in Pa"
          call save_to_hdf5(h5file_id,dset,params%cpp%pressure,attr)

          dset = TRIM(gname) // "/T"
          attr = "Characteristic plasma temperature in J"
          call save_to_hdf5(h5file_id,dset,params%cpp%temperature,attr)

          call h5gclose_f(group_id, h5error)

          call h5fclose_f(h5file_id, h5error)
       end if

    end if
  end subroutine save_simulation_parameters