iol_source_io_s.f90 Source File


Contents

Source Code


Source Code

submodule(iol_source_m) iol_source_io_s
    !! Routines related with I/O of IOL sources
    implicit none
    
contains
    
    module subroutine display_iol_source(self)
        class(iol_source_t), intent(inout) :: self

        if (is_rank_info_writer) then
            write(get_stdout(),*)''            
            write(get_stdout(),*)'Parameters for ion-orbit-loss ---------------------'
                                           
            write(get_stdout(),202)self%signBt, self%signIp, &
                             self%R_magaxis, self%RX, self%zX, &
                             self%BphiX, self%psix, self%psirange_conf_red, &
                             self%invs_aspc, &
                             self%mass, self%charge, &
                             self%kmax, self%mmax, self%lmax, &
                             self%nsep_corr_cw, self%nsep_corr_ccw, &
                             self%er_mult, self%ti_mult, &
                             self%imax, self%jmax, &
                             self%idx_X, self%idx_imp, &
                             self%nrho_sp, self%xplane_n, &
                             self%bnrm, self%Rmin_m, self%Rmin_b, &
                             self%tau_av_frame, &
                             self%init_continue
                             
202         FORMAT(3X,'signBt              = ',I8       /, &
                   3X,'signIp              = ',I8       /, &
                   3X,'R_magaxis           = ',ES14.6E3 /, &
                   3X,'RX                  = ',ES14.6E3 /, &
                   3X,'zX                  = ',ES14.6E3 /, &
                   3X,'BphiX               = ',ES14.6E3 /, &
                   3X,'psix                = ',ES14.6E3 /, &
                   3X,'psirange_conf_red   = ',ES14.6E3 /, &
                   3X,'invs_aspc           = ',ES14.6E3 /, &
                   3X,'mass                = ',ES14.6E3 /, &
                   3X,'charge              = ',ES14.6E3 /, &
                   3X,'kmax                = ',I8       /, &
                   3X,'mmax                = ',I8       /, &
                   3X,'lmax                = ',I8       /, &
                   3X,'nsep_corr_cw        = ',I8       /, &
                   3X,'nsep_corr_ccw       = ',I8       /, &
                   3X,'er_mult             = ',ES14.6E3 /, &
                   3X,'ti_mult             = ',ES14.6E3 /, &
                   3X,'imax                = ',I8       /, &
                   3X,'jmax                = ',I8       /, &
                   3X,'idx_X               = ',I8       /, &
                   3X,'idx_imp             = ',I8       /, &
                   3X,'nrho_sp             = ',I8       /, &
                   3X,'xplane_n            = ',I8       /, &
                   3X,'bnrm                = ',ES14.6E3 /, &
                   3X,'Rmin_m              = ',ES14.6E3 /, &
                   3X,'Rmin_b              = ',ES14.6E3 /, &
                   3X,'tau_av_frame        = ',ES14.6E3 /, &
                   3X,'init_continue       = ',L1          )
            write(get_stdout(),*)'-------------------------------------------------------'
            write(get_stdout(),*)''            

        endif
        
    end subroutine  
    
    module subroutine write_netcdf_iol_source(self, tau, fexist, fgid)
        class(iol_source_t), intent(in) :: self
        real(GP), intent(in) :: tau
        logical, intent(in) :: fexist
        integer, intent(in) :: fgid
        
        integer :: npoints_inner, nsnaps
        integer :: nf90_stat
        integer :: iddim_snaps, iddim_npoints_inner, id_tau, id_src_iol
        integer :: offset_1d(1), offset_2d(2)
    
        npoints_inner = size(self%src_iol)
        
        if (fexist) then
            ! Inquire dimensions and variables
            nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'nsnaps', nsnaps )
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            nsnaps = nsnaps+1

            nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'nsnaps', nsnaps )
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
                        
            nf90_stat = nf90_inq_varid(fgid, 'tau', id_tau)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)

            nf90_stat = nf90_inq_varid(fgid, 'src_iol', id_src_iol)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)           
          
        else
            nsnaps = 1
            nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'nsnaps', nsnaps )
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        
            nf90_stat = nf90_def_dim(fgid, 'snaps', NF90_UNLIMITED, iddim_snaps)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) 
            
            nf90_stat = nf90_def_dim(fgid, 'npoints_inner', &
                npoints_inner, iddim_npoints_inner)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) 

            ! Define variables
            nf90_stat = nf90_def_var(fgid, 'tau', NF90_DOUBLE, &
                iddim_snaps, id_tau)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)

            nf90_stat = nf90_def_var(fgid, 'src_iol', NF90_DOUBLE, &
                (/iddim_npoints_inner, iddim_snaps/), id_src_iol)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)

            nf90_stat = nf90_enddef(fgid)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)

        endif
        
        offset_1d = (/nsnaps/)
        nf90_stat = nf90_put_var(fgid, id_tau, tau, offset_1d )
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)

        offset_2d = (/1, nsnaps/)        
        nf90_stat = nf90_put_var(fgid, id_src_iol, self%src_iol, offset_2d )
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
                
    end subroutine
    
    module subroutine read_netcdf_iol_source(self, fgid, nsnaps_req, nsnaps, tau)
        class(iol_source_t), intent(inout) :: self
        integer, intent(in) :: fgid
        integer, intent(in), optional :: nsnaps_req
        integer, intent(out), optional :: nsnaps
        real(GP), intent(out),  optional :: tau
        
        integer :: nf90_stat
        integer :: nsnaps_ret, npoints_inner
        real(GP) :: tau_ret
        integer :: iddim_npoints_inner, id_tau, id_src_iol, offset_1d(1), offset_2d(2)
        
        
        ! Determine snapshot to be returned
        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'nsnaps', nsnaps_ret)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        if (present(nsnaps_req)) then
            if (nsnaps_req > nsnaps_ret) then
                call handle_error('Nsnaps_req not available', &
                                  GRILLIX_ERR_OTHER, __LINE__, __FILE__, &
                                  error_info_t('nsnasp_req, nsnaps_ret: ', &
                                               [nsnaps_req, nsnaps_ret]))
            else
                nsnaps_ret = nsnaps_req
            endif
        endif
        
        if (present(nsnaps)) then
            nsnaps = nsnaps_ret
        endif
        
        ! Check consistency of dimension
        nf90_stat =  nf90_inq_dimid(fgid, 'npoints_inner', iddim_npoints_inner)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        nf90_stat =  nf90_inquire_dimension(fgid, iddim_npoints_inner, len = npoints_inner)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
    
        if (npoints_inner /= size(self%src_iol)) then
            call handle_error('Dimension mismatch', &
                              GRILLIX_ERR_OTHER, __LINE__, __FILE__, &
                              error_info_t('npoints_inner, size(src_iol): ', &
                                           [npoints_inner, size(self%src_iol)]))
        endif
        
        ! Get tau
        nf90_stat = nf90_inq_varid(fgid, 'tau', id_tau)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        offset_1d = (/nsnaps_ret/)
        nf90_stat = nf90_get_var(fgid, id_tau, tau_ret, offset_1d )
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        if (present(tau)) then
            tau = tau_ret
        endif
        
        ! Get src_iol
        nf90_stat = nf90_inq_varid(fgid, 'src_iol', id_src_iol)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        offset_2d = (/1, nsnaps_ret/)
        nf90_stat = nf90_get_var(fgid, id_src_iol, self%src_iol, offset_2d )
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
    end subroutine    

    module subroutine write_ascii_iol_source(self, dirpath)
        class(iol_source_t), intent(inout) :: self
        character(len=*), intent(in) :: dirpath
        
        integer :: io_error
        character(len=256) :: io_errmsg 
        
        character(len=32) :: fmt_nrho_sp, fmt_jmax, fmt_imax, fmt_xplane_n
        integer :: i
        
        ! Format descriptors
        write(fmt_nrho_sp,'("(",I0,"(ES22.14E3,X))")')self%nrho_sp
        write(fmt_jmax,'("(",I0,"(ES22.14E3,X))")')self%jmax
        write(fmt_imax,'("(",I0,"(ES22.14E3,X))")')self%imax
        write(fmt_xplane_n,'("(",I0,"(ES22.14E3,X))")')self%xplane_n

        if (is_rank_info_writer) then
        
            open (52, file = dirpath//'/rhopts_sp.txt', status = 'new', iostat = io_error, iomsg = io_errmsg)
            if (io_error /= 0) then
                write(*,*)io_errmsg
                error stop
            endif        
            write(52,fmt_nrho_sp) self%rhopts_sp
            close(52)
                 
            open (52, file = dirpath//'/Rmincontour_sp.txt', status = 'new', iostat = io_error, iomsg = io_errmsg)
            if (io_error /= 0) then
                 write(*,*)io_errmsg
                error stop
            endif      
            write(52,fmt_nrho_sp) self%Rmincontour_sp(:)
            close(52)
          
            open (52, file = dirpath//'/phi_sp.txt', status = 'new', iostat = io_error, iomsg = io_errmsg)
            if (io_error /= 0) then
                 write(*,*)io_errmsg
                error stop
            endif      
            write(52,fmt_nrho_sp) self%pot_sp
            close(52)
          
            open (52, file = dirpath//'/sep_Rcoord.txt', status = 'new', iostat = io_error, iomsg = io_errmsg)
            if (io_error /= 0) then
                 write(*,*)io_errmsg
                error stop
            endif     
            write(52,fmt_jmax) self%sep_Rcoord
            close(52)
          
            open (52, file = dirpath//'/sep_cw.txt', status = 'new', iostat = io_error, iomsg = io_errmsg)
            if (io_error /= 0) then
                 write(*,*)io_errmsg
                error stop
            endif     
            write(52,fmt_jmax) self%sep_cw
            close(52)
          
            open (52, file = dirpath//'/sep_cw_mod.txt', status = 'new', iostat = io_error, iomsg = io_errmsg)
            if (io_error /= 0) then
                 write(*,*)io_errmsg
                error stop
            endif     
            write(52,fmt_jmax) self%sep_cw_mod
            close(52)
          
            open (52, file = dirpath//'/sep_ccw.txt', status = 'new', iostat = io_error, iomsg = io_errmsg)
            if (io_error /= 0) then
                 write(*,*)io_errmsg
                error stop
            endif     
            write(52,fmt_jmax) self%sep_ccw
            close(52)
          
            open (52, file = dirpath//'/sep_ccw_mod.txt', status = 'new', iostat = io_error, iomsg = io_errmsg)
            if (io_error /= 0) then
                 write(*,*)io_errmsg
                error stop
            endif     
            write(52,fmt_jmax) self%sep_ccw_mod
            close(52)
          
            open (52, file = dirpath//'/Rcoord.txt', status = 'new', iostat = io_error, iomsg = io_errmsg)
            if (io_error /= 0) then
                 write(*,*)io_errmsg
                error stop
            endif     
            do i = 1,self%imax
                write(52,fmt_jmax) self%Rcoord(:,i)
            end do
            close(52)
          
            open (52, file = dirpath//'/zcoord.txt', status = 'new', iostat = io_error, iomsg = io_errmsg)
            if (io_error /= 0) then
                 write(*,*)io_errmsg
                error stop
            endif     
            do i = 1,self%imax
                write(52,fmt_jmax) self%zcoord(:,i)
            end do
            close(52)
          
            open (52, file = dirpath//'/Babs.txt', status = 'new', iostat = io_error, iomsg = io_errmsg)
            if (io_error /= 0) then
                 write(*,*)io_errmsg
                error stop
            endif     
            do i = 1,self%imax
                write(52,fmt_jmax) self%Babs(:,i)
            end do
            close(52)
          
            open (52, file = dirpath//'/Bphi.txt', status = 'new', iostat = io_error, iomsg = io_errmsg)
            if (io_error /= 0) then
                 write(*,*)io_errmsg
                error stop
            endif     
            do i = 1,self%imax
                write(52,fmt_jmax) self%Bphi(:,i)
            end do
            close(52)
          
            open (52, file = dirpath//'/psi.txt', status = 'new', iostat = io_error, iomsg = io_errmsg)
            if (io_error /= 0) then
                 write(*,*)io_errmsg
                error stop
            endif     
            do i = 1,self%imax
                write(52,fmt_jmax) self%psi(:,i)
            end do
            close(52)
          
            open (52, file = dirpath//'/temp_ion.txt', status = 'new', iostat = io_error, iomsg = io_errmsg)
            if (io_error /= 0) then
                 write(*,*)io_errmsg
                error stop
            endif     
            do i = 1,self%imax
                write(52,fmt_jmax) self%temp_ion(:,i)
            end do
            close(2)
          
            open (52, file = dirpath//'/phi.txt', status = 'new', iostat = io_error, iomsg = io_errmsg)
            if (io_error /= 0) then
                 write(*,*)io_errmsg
                error stop
            endif     
            do i = 1,self%imax
                write(52,fmt_jmax) self%pot(:,i)
            end do
            close(52)

            open (52, file = dirpath//'/safetyf.txt', status = 'new', iostat = io_error, iomsg = io_errmsg)
            if (io_error /= 0) then
                 write(*,*)io_errmsg
                error stop
            endif     
            write(52,fmt_imax) self%safetyf(:)
            close(52)

            open (52, file = dirpath//'/dens.txt', status = 'new', iostat = io_error, iomsg = io_errmsg)
            if (io_error /= 0) then
                 write(*,*)io_errmsg
                error stop
            endif     
            do i = 1,self%imax
                write(52,fmt_jmax) self%dens(:,i)
            end do
            close(52)

            open (52, file = dirpath//'/Xplane_R.txt', status = 'new', iostat = io_error, iomsg = io_errmsg)
            if (io_error /= 0) then
                 write(*,*)io_errmsg
                error stop
            endif     
            write(52,fmt_xplane_n) self%Xplane_R(:)
            close(52)

            open (52, file = dirpath//'/Xplane_psi.txt', status = 'new', iostat = io_error, iomsg = io_errmsg)
            if (io_error /= 0) then
                 write(*,*)io_errmsg
                error stop
            endif     
            write(52,fmt_xplane_n) self%Xplane_psi(:)
            close(52)

            open (52, file = dirpath//'/Xplane_Babs.txt', status = 'new', iostat = io_error, iomsg = io_errmsg)
            if (io_error /= 0) then
                 write(*,*)io_errmsg
                error stop
            endif     
            write(52,fmt_xplane_n) self%Xplane_Babs(:)
            close(52)

            open (52, file = dirpath//'/Xplane_Bphi.txt', status = 'new', iostat = io_error, iomsg = io_errmsg)
            if (io_error /= 0) then
                 write(*,*)io_errmsg
                error stop
            endif     
            write(52,fmt_xplane_n) self%Xplane_Bphi(:)
            close(52)

        endif
        
    end subroutine
         
end submodule