submodule(init_zhdanov_m) init_zhdanov_profile_s !! Routines for initialisation of variables with profile implicit none contains module subroutine init_zhdanov_state_profile(comm_handler, & equi, mesh, map, & penalisation, params_species, & paramfile_zhdanov_species, & dens_species) !! Builds profile state type(comm_handler_t), intent(in) :: comm_handler !! Communicators class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(inout) :: mesh !! Mesh within poloidal plane type(parallel_map_t), intent(in) :: map !! Parallel map type(penalisation_t), intent(in) :: penalisation !! Penalisation type(params_zhdanov_species_t), intent(in) :: params_species !! Species paramters for the Zhdanov model character(len=*), intent(in) :: paramfile_zhdanov_species !! Filename, where parameters are read from type(variable_t), intent(inout) :: dens_species !! Species density integer :: io_error character(len=256) :: io_errmsg integer :: l real(GP) :: phi, x, y, rho real(GP) :: rho_ped_dens_species, rho_sol_dens_species, & val_ped_dens_species, & val_sol_dens_species, flucamp_dens_species namelist / init_profile / rho_ped_dens_species, rho_sol_dens_species, & val_ped_dens_species, & val_sol_dens_species, flucamp_dens_species ! Read parameters of init_profile ------------------------------------- open(unit = 20, file = paramfile_zhdanov_species, status = 'old', & action = 'read',& iostat = io_error, iomsg = io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) endif ! default values rho_ped_dens_species = equi%rhomin rho_sol_dens_species = 1.0_GP val_ped_dens_species = 1.0_GP val_sol_dens_species = 1.0E-1_GP flucamp_dens_species = 1.0E-5_GP read(20, nml = init_profile, iostat = io_error, iomsg = io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) endif if (is_rank_info_writer) then write(get_stdout(),*)'' write(get_stdout(),*)'Creating profile as initial state --------------------------------------' write(get_stdout(),202)rho_ped_dens_species, rho_sol_dens_species, & val_ped_dens_species, val_sol_dens_species, flucamp_dens_species 202 FORMAT(1X,'rho_ped_dens_species = ',ES14.6E3 /, & 1X,'rho_sol_dens_species = ',ES14.6E3 /, & 1X,'val_ped_dens_species = ',ES14.6E3 /, & 1X,'val_sol_dens_species = ',ES14.6E3 /, & 1X,'flucamp_dens_species = ',ES14.6E3 ) endif close(20) phi = mesh%get_phi() ! Set profiles dens_species !$OMP PARALLEL PRIVATE(l, x, y, rho) !$OMP DO do l = 1, mesh%get_n_points() x = mesh%get_x(l) y = mesh%get_y(l) rho = equi%rho(x, y, phi) dens_species%vals(l) = profile_sin(rho_ped_dens_species, & rho_sol_dens_species, val_ped_dens_species, & val_sol_dens_species, rho) enddo !$OMP END DO !$OMP END PARALLEL end subroutine pure real(GP) function profile_sin(rho_ped, rho_sol, val_ped, val_sol, rho) !! Returns profile based on sin function, i.e. real(GP), intent(in) :: rho_ped !! Radial position of pedestal top real(GP), intent(in) :: rho_sol !! Radial position of pedestal bottom real(GP), intent(in) :: val_ped !! Value at pedestal top real(GP), intent(in) :: val_sol !! Value at pedestal bottom real(GP), intent(in) :: rho !! Radial position real(GP) :: amp, offset, omega, phase if (rho < rho_ped) then profile_sin = val_ped elseif (rho > rho_sol) then profile_sin = val_sol else amp = (val_ped - val_sol) / 2.0_GP offset = val_ped - amp omega = PI / ( rho_sol - rho_ped ) phase = - 0.5_GP*PI * ( rho_sol + rho_ped ) / ( rho_sol - rho_ped ) profile_sin = -amp * sin(omega*rho + phase) + offset endif end function end submodule