parallel_map_netcdf_s.f90 Source File

submodule (parallel_map_m) parallel_map_netcdf_s
    !! Rotines related with netcdf I/O
    implicit none
    
contains

    module subroutine write_netcdf_parallel_map(self, axisymmetric, fgid)
        class(parallel_map_t), intent(in) :: self
        logical, intent(in) :: axisymmetric
        integer, intent(in) :: fgid

        integer :: itemp, n_points_cano, n_points_stag
        integer :: nf90_stat
        integer :: iddim_n_points_cano, iddim_n_points_stag
        integer :: id_dpar_cano_stag_fwd, id_dpar_cano_stag_bwd, &
                   id_dpar_cano_cano_fwd, id_dpar_cano_cano_bwd, &
                   id_fluxbox_vol_cano
        integer :: id_dpar_stag_cano_fwd, id_dpar_stag_cano_bwd, &
                   id_dpar_stag_stag_fwd, id_dpar_stag_stag_bwd, &
                   id_fluxbox_vol_stag 
        integer :: id_map_cano_stag_fwd, id_map_cano_stag_bwd, &
                   id_map_cano_cano_fwd, id_map_cano_cano_bwd, &
                   id_map_cano_cano_fwd2, id_map_cano_cano_bwd2, &
                   id_map_stag_cano_fwd, id_map_stag_cano_bwd, &
                   id_map_stag_stag_fwd, id_map_stag_stag_bwd
        integer :: id_grad_stag_cano_fwd, id_grad_stag_cano_bwd, &
                   id_pdiv_cano_stag_fwd, id_pdiv_cano_stag_bwd
                   
        ! Define dimensions 
        n_points_cano = self%map_cano_stag_fwd%ndim
        n_points_stag = self%map_stag_cano_fwd%ndim
        
        nf90_stat = nf90_def_dim(fgid, 'n_points_cano', &
                                 n_points_cano, iddim_n_points_cano)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        nf90_stat = nf90_def_dim(fgid, 'n_points_stag', &
                                 n_points_stag, iddim_n_points_stag)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        ! Write global attributes (metadata)
        itemp = logical_to_integer(axisymmetric)
        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'axisymmetric', itemp)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'nplanes', self%nplanes)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'iplane', self%iplane)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'phi_cano', self%phi_cano)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'phi_stag', self%phi_stag)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'dphi', self%dphi)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'intorder', self%intorder)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'xorder', self%xorder)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        itemp = logical_to_integer(self%use_fixed_stencil)
        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'use_fixed_stencil', itemp)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
  
        itemp = logical_to_integer(self%use_gauss_quadrature)
        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'use_gauss_quadrature', itemp)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        ! Define cano-variables
        nf90_stat = nf90_def_var(fgid, 'dpar_cano_stag_fwd', NF90_DOUBLE, &
                                 iddim_n_points_cano, id_dpar_cano_stag_fwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
     
        nf90_stat = nf90_def_var(fgid, 'dpar_cano_stag_bwd', NF90_DOUBLE, &
                                 iddim_n_points_cano, id_dpar_cano_stag_bwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
     
        nf90_stat = nf90_def_var(fgid, 'dpar_cano_cano_fwd', NF90_DOUBLE, &
                                 iddim_n_points_cano, id_dpar_cano_cano_fwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
     
        nf90_stat = nf90_def_var(fgid, 'dpar_cano_cano_bwd', NF90_DOUBLE, &
                                 iddim_n_points_cano, id_dpar_cano_cano_bwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
     
        nf90_stat = nf90_def_var(fgid, 'fluxbox_vol_cano', NF90_DOUBLE, &
                                 iddim_n_points_cano, id_fluxbox_vol_cano)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
     
        ! Define stag-variables (only for non-axisymmetric case)
        if (.not.axisymmetric) then
            nf90_stat = nf90_def_var(fgid, 'dpar_stag_cano_fwd', NF90_DOUBLE, &
                                     iddim_n_points_stag, id_dpar_stag_cano_fwd)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            
            nf90_stat = nf90_def_var(fgid, 'dpar_stag_cano_bwd', NF90_DOUBLE, &
                                     iddim_n_points_stag, id_dpar_stag_cano_bwd)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            
            nf90_stat = nf90_def_var(fgid, 'dpar_stag_stag_fwd', NF90_DOUBLE, &
                                     iddim_n_points_stag, id_dpar_stag_stag_fwd)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            
            nf90_stat = nf90_def_var(fgid, 'dpar_stag_stag_bwd', NF90_DOUBLE, &
                                     iddim_n_points_stag, id_dpar_stag_stag_bwd)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            
            nf90_stat = nf90_def_var(fgid, 'fluxbox_vol_stag', NF90_DOUBLE, &
                                     iddim_n_points_stag, id_fluxbox_vol_stag)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        endif       
        
        ! Define netcdf groups for map matrices
        nf90_stat = nf90_def_grp(fgid, 'map_stag_cano_fwd', &
                                       id_map_stag_cano_fwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        nf90_stat = nf90_def_grp(fgid, 'map_stag_cano_bwd', &
                                       id_map_stag_cano_bwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        nf90_stat = nf90_def_grp(fgid, 'map_cano_cano_fwd', &
                                       id_map_cano_cano_fwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)

        nf90_stat = nf90_def_grp(fgid, 'map_cano_cano_fwd2', &
                                        id_map_cano_cano_fwd2)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        nf90_stat = nf90_def_grp(fgid, 'map_cano_cano_bwd', &
                                       id_map_cano_cano_bwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)

        nf90_stat = nf90_def_grp(fgid, 'map_cano_cano_bwd2', &
                                       id_map_cano_cano_bwd2)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        if (.not.axisymmetric) then
            nf90_stat = nf90_def_grp(fgid, 'map_cano_stag_fwd', &
                                           id_map_cano_stag_fwd)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            
            nf90_stat = nf90_def_grp(fgid, 'map_cano_stag_bwd', &
                                           id_map_cano_stag_bwd)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            
            nf90_stat = nf90_def_grp(fgid, 'map_stag_stag_fwd', &
                                           id_map_stag_stag_fwd)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            
            nf90_stat = nf90_def_grp(fgid, 'map_stag_stag_bwd', &
                                           id_map_stag_stag_bwd)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        endif
        
        ! Define groups for operator matrices (gradpar and divpar)
        nf90_stat = nf90_def_grp(fgid, 'grad_stag_cano_fwd', &
                                       id_grad_stag_cano_fwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            
        nf90_stat = nf90_def_grp(fgid, 'grad_stag_cano_bwd', &
                                       id_grad_stag_cano_bwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            
        nf90_stat = nf90_def_grp(fgid, 'pdiv_cano_stag_fwd', &
                                       id_pdiv_cano_stag_fwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        nf90_stat = nf90_def_grp(fgid, 'pdiv_cano_stag_bwd', &
                                       id_pdiv_cano_stag_bwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
              
        ! End of definition
        nf90_stat = nf90_enddef(fgid)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        ! Put variables
        nf90_stat = nf90_put_var(fgid, id_dpar_cano_stag_fwd, &
                                 self%dpar_cano_stag_fwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        nf90_stat = nf90_put_var(fgid, id_dpar_cano_stag_bwd, &
                                 self%dpar_cano_stag_bwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        nf90_stat = nf90_put_var(fgid, id_dpar_cano_cano_fwd, &
                                 self%dpar_cano_cano_fwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        nf90_stat = nf90_put_var(fgid, id_dpar_cano_cano_bwd, &
                                 self%dpar_cano_cano_bwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        nf90_stat = nf90_put_var(fgid, id_fluxbox_vol_cano, &
                                 self%fluxbox_vol_cano)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        if (.not.axisymmetric) then
            nf90_stat = nf90_put_var(fgid, id_dpar_stag_cano_fwd, &
                                 self%dpar_stag_cano_fwd)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            
            nf90_stat = nf90_put_var(fgid, id_dpar_stag_cano_bwd, &
                                     self%dpar_stag_cano_bwd)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            
            nf90_stat = nf90_put_var(fgid, id_dpar_stag_stag_fwd, &
                                     self%dpar_stag_stag_fwd)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            
            nf90_stat = nf90_put_var(fgid, id_dpar_stag_stag_bwd, &
                                     self%dpar_stag_stag_bwd)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            
            nf90_stat = nf90_put_var(fgid, id_fluxbox_vol_stag, &
                                     self%fluxbox_vol_stag)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        endif

        ! Write matrices
        call self%map_stag_cano_fwd%write_netcdf(id_map_stag_cano_fwd)
        call self%map_stag_cano_bwd%write_netcdf(id_map_stag_cano_bwd)

        call self%map_cano_cano_fwd%write_netcdf(id_map_cano_cano_fwd)
        call self%map_cano_cano_fwd2%write_netcdf(id_map_cano_cano_fwd2)
        call self%map_cano_cano_bwd%write_netcdf(id_map_cano_cano_bwd)
        call self%map_cano_cano_bwd2%write_netcdf(id_map_cano_cano_bwd2)

        if (.not.axisymmetric) then
            call self%map_cano_stag_fwd%write_netcdf(id_map_cano_stag_fwd)
            call self%map_cano_stag_bwd%write_netcdf(id_map_cano_stag_bwd)

            call self%map_stag_stag_fwd%write_netcdf(id_map_stag_stag_fwd)
            call self%map_stag_stag_bwd%write_netcdf(id_map_stag_stag_bwd)
        endif

        call self%grad_stag_cano_fwd%write_netcdf(id_grad_stag_cano_fwd)
        call self%grad_stag_cano_bwd%write_netcdf(id_grad_stag_cano_bwd)

        call self%pdiv_cano_stag_fwd%write_netcdf(id_pdiv_cano_stag_fwd)
        call self%pdiv_cano_stag_bwd%write_netcdf(id_pdiv_cano_stag_bwd)

    end subroutine

    module subroutine read_netcdf_parallel_map(self, axisymmetric, fgid)
        class(parallel_map_t), intent(inout), target :: self
        logical, intent(in) :: axisymmetric
        integer, intent(in) :: fgid
        
        integer :: itemp, n_points_cano, n_points_stag
        logical :: axsym_from_file
        integer :: nf90_stat
        integer :: iddim_n_points_cano, iddim_n_points_stag
        integer :: id_dpar_cano_stag_fwd, id_dpar_cano_stag_bwd, &
                   id_dpar_cano_cano_fwd, id_dpar_cano_cano_bwd, &
                   id_fluxbox_vol_cano
        integer :: id_dpar_stag_cano_fwd, id_dpar_stag_cano_bwd, &
                   id_dpar_stag_stag_fwd, id_dpar_stag_stag_bwd, &
                   id_fluxbox_vol_stag 
        integer :: id_map_cano_stag_fwd, id_map_cano_stag_bwd, &
                   id_map_cano_cano_fwd, id_map_cano_cano_bwd, &
                   id_map_cano_cano_fwd2, id_map_cano_cano_bwd2, &
                   id_map_stag_cano_fwd, id_map_stag_cano_bwd, &
                   id_map_stag_stag_fwd, id_map_stag_stag_bwd
        integer :: id_grad_stag_cano_fwd, id_grad_stag_cano_bwd, &
                   id_pdiv_cano_stag_fwd, id_pdiv_cano_stag_bwd

        ! Check if data matches axiymmetry parameter
        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'axisymmetric', itemp)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        axsym_from_file = integer_to_logical(itemp)
        if (.not.(axisymmetric.eqv.axsym_from_file)) then
            call handle_error('Axisymmetric does not match file', GRILLIX_ERR_OTHER, __LINE__, __FILE__)
        endif
                
        ! Read attributes
        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'nplanes', self%nplanes)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'iplane', self%iplane)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'phi_cano', self%phi_cano)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'phi_stag', self%phi_stag)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'dphi', self%dphi)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'intorder', self%intorder)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'xorder', self%xorder)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'use_fixed_stencil', itemp)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        self%use_fixed_stencil = integer_to_logical(itemp)
  
        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'use_gauss_quadrature', itemp)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        self%use_gauss_quadrature = integer_to_logical(itemp)
        
        ! Read dimensions
        nf90_stat = nf90_inq_dimid(fgid, 'n_points_cano', iddim_n_points_cano)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        nf90_stat =  nf90_inquire_dimension(fgid, iddim_n_points_cano, len=n_points_cano)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        nf90_stat = nf90_inq_dimid(fgid, 'n_points_stag', iddim_n_points_stag)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        nf90_stat =  nf90_inquire_dimension(fgid, iddim_n_points_stag, len=n_points_stag)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        ! Read cano-variables
        if (.not. allocated(self%dpar_cano_stag_fwd)) then
            allocate(self%dpar_cano_stag_fwd(n_points_cano)) 
        else
            call handle_error('Already allocated', &
                              GRILLIX_ERR_ALLOC, __LINE__, __FILE__)
        endif        
        nf90_stat = nf90_inq_varid(fgid, 'dpar_cano_stag_fwd', id_dpar_cano_stag_fwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        nf90_stat = nf90_get_var(fgid, id_dpar_cano_stag_fwd, self%dpar_cano_stag_fwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        if (.not. allocated(self%dpar_cano_stag_bwd)) then
            allocate(self%dpar_cano_stag_bwd(n_points_cano)) 
        else
            call handle_error('Already allocated', &
                              GRILLIX_ERR_ALLOC, __LINE__, __FILE__)
        endif        
        nf90_stat = nf90_inq_varid(fgid, 'dpar_cano_stag_bwd', id_dpar_cano_stag_bwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        nf90_stat = nf90_get_var(fgid, id_dpar_cano_stag_bwd, self%dpar_cano_stag_bwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)

        if (.not. allocated(self%dpar_cano_cano_fwd)) then
            allocate(self%dpar_cano_cano_fwd(n_points_cano)) 
        else
            call handle_error('Already allocated', &
                              GRILLIX_ERR_ALLOC, __LINE__, __FILE__)
        endif        
        nf90_stat = nf90_inq_varid(fgid, 'dpar_cano_cano_fwd', id_dpar_cano_cano_fwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        nf90_stat = nf90_get_var(fgid, id_dpar_cano_cano_fwd, self%dpar_cano_cano_fwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        if (.not. allocated(self%dpar_cano_cano_bwd)) then
            allocate(self%dpar_cano_cano_bwd(n_points_cano)) 
        else
            call handle_error('Already allocated', &
                              GRILLIX_ERR_ALLOC, __LINE__, __FILE__)
        endif        
        nf90_stat = nf90_inq_varid(fgid, 'dpar_cano_cano_bwd', id_dpar_cano_cano_bwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        nf90_stat = nf90_get_var(fgid, id_dpar_cano_cano_bwd, self%dpar_cano_cano_bwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        if (.not. allocated(self%fluxbox_vol_cano)) then
            allocate(self%fluxbox_vol_cano(n_points_cano)) 
        else
            call handle_error('Already allocated', &
                              GRILLIX_ERR_ALLOC, __LINE__, __FILE__)
        endif        
        nf90_stat = nf90_inq_varid(fgid, 'fluxbox_vol_cano', id_fluxbox_vol_cano)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        nf90_stat = nf90_get_var(fgid, id_fluxbox_vol_cano, self%fluxbox_vol_cano)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)

        if (axisymmetric) then
            ! Associate stag-variables
            ! For axisymmetric case quantities for staggered mesh 
            ! are the same asfor canonical mesh
            self%dpar_stag_cano_fwd => self%dpar_cano_stag_fwd
            self%dpar_stag_cano_bwd => self%dpar_cano_stag_bwd
            self%dpar_stag_stag_fwd => self%dpar_cano_cano_fwd
            self%dpar_stag_stag_bwd => self%dpar_cano_cano_bwd
            self%fluxbox_vol_stag => self%fluxbox_vol_cano
        else 
            ! Read stag-variables   
            allocate(self%dpar_stag_cano_fwd_mem(n_points_stag))
            self%dpar_stag_cano_fwd => self%dpar_stag_cano_fwd_mem            
            nf90_stat = nf90_inq_varid(fgid, 'dpar_stag_cano_fwd', id_dpar_stag_cano_fwd)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            nf90_stat = nf90_get_var(fgid, id_dpar_stag_cano_fwd, self%dpar_stag_cano_fwd)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            
            allocate(self%dpar_stag_cano_bwd_mem(n_points_stag))
            self%dpar_stag_cano_bwd => self%dpar_stag_cano_bwd_mem            
            nf90_stat = nf90_inq_varid(fgid, 'dpar_stag_cano_bwd', id_dpar_stag_cano_bwd)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            nf90_stat = nf90_get_var(fgid, id_dpar_stag_cano_bwd, self%dpar_stag_cano_bwd)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            
            allocate(self%dpar_stag_stag_fwd_mem(n_points_stag))
            self%dpar_stag_stag_fwd => self%dpar_stag_stag_fwd_mem
            nf90_stat = nf90_inq_varid(fgid, 'dpar_stag_stag_fwd', id_dpar_stag_stag_fwd)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            nf90_stat = nf90_get_var(fgid, id_dpar_stag_stag_fwd, self%dpar_stag_stag_fwd)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            
            allocate(self%dpar_stag_stag_bwd_mem(n_points_stag))
            self%dpar_stag_stag_bwd => self%dpar_stag_stag_bwd_mem
            nf90_stat = nf90_inq_varid(fgid, 'dpar_stag_stag_bwd', id_dpar_stag_stag_bwd)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            nf90_stat = nf90_get_var(fgid, id_dpar_stag_stag_bwd, self%dpar_stag_stag_bwd)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            
            allocate(self%fluxbox_vol_stag_mem(n_points_stag))
            self%fluxbox_vol_stag => self%fluxbox_vol_stag_mem
            nf90_stat = nf90_inq_varid(fgid, 'fluxbox_vol_stag', id_fluxbox_vol_stag)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            nf90_stat = nf90_get_var(fgid, id_fluxbox_vol_stag, self%fluxbox_vol_stag)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            
        endif

        ! Read map matrices
        if (.not. allocated(self%map_stag_cano_fwd)) then
            allocate(self%map_stag_cano_fwd) 
        else
            call handle_error('Already allocated', &
                              GRILLIX_ERR_ALLOC, __LINE__, __FILE__)
        endif
        nf90_stat = nf90_inq_grp_ncid(fgid, 'map_stag_cano_fwd', id_map_stag_cano_fwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        call self%map_stag_cano_fwd%read_netcdf(id_map_stag_cano_fwd)
        
        if (.not. allocated(self%map_stag_cano_bwd)) then
            allocate(self%map_stag_cano_bwd) 
        else
            call handle_error('Already allocated', &
                              GRILLIX_ERR_ALLOC, __LINE__, __FILE__)
        endif
        nf90_stat = nf90_inq_grp_ncid(fgid, 'map_stag_cano_bwd', id_map_stag_cano_bwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        call self%map_stag_cano_bwd%read_netcdf(id_map_stag_cano_bwd)
        
        if (.not. allocated(self%map_cano_cano_fwd)) then
            allocate(self%map_cano_cano_fwd) 
        else
            call handle_error('Already allocated', &
                              GRILLIX_ERR_ALLOC, __LINE__, __FILE__)
        endif
        nf90_stat = nf90_inq_grp_ncid(fgid, 'map_cano_cano_fwd', id_map_cano_cano_fwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        call self%map_cano_cano_fwd%read_netcdf(id_map_cano_cano_fwd)

        if (.not. allocated(self%map_cano_cano_fwd2)) then
            allocate(self%map_cano_cano_fwd2) 
        else
            call handle_error('Already allocated', &
                              GRILLIX_ERR_ALLOC, __LINE__, __FILE__)
        endif
        nf90_stat = nf90_inq_grp_ncid(fgid, 'map_cano_cano_fwd2', id_map_cano_cano_fwd2)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        call self%map_cano_cano_fwd2%read_netcdf(id_map_cano_cano_fwd2)
        
        if (.not. allocated(self%map_cano_cano_bwd)) then
            allocate(self%map_cano_cano_bwd) 
        else
            call handle_error('Already allocated', &
                              GRILLIX_ERR_ALLOC, __LINE__, __FILE__)
        endif
        nf90_stat = nf90_inq_grp_ncid(fgid, 'map_cano_cano_bwd', id_map_cano_cano_bwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        call self%map_cano_cano_bwd%read_netcdf(id_map_cano_cano_bwd)

        if (.not. allocated(self%map_cano_cano_bwd2)) then
            allocate(self%map_cano_cano_bwd2) 
        else
            call handle_error('Already allocated', &
                              GRILLIX_ERR_ALLOC, __LINE__, __FILE__)
        endif
        nf90_stat = nf90_inq_grp_ncid(fgid, 'map_cano_cano_bwd2', id_map_cano_cano_bwd2)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        call self%map_cano_cano_bwd2%read_netcdf(id_map_cano_cano_bwd2)

        if (axisymmetric) then
            ! Associate further map matrices
            self%map_cano_stag_fwd => self%map_stag_cano_fwd
            self%map_cano_stag_bwd => self%map_stag_cano_bwd
            self%map_stag_stag_fwd => self%map_cano_cano_fwd
            self%map_stag_stag_bwd => self%map_cano_cano_bwd
        else 
            
            allocate(self%map_cano_stag_fwd_mem)
            self%map_cano_stag_fwd => self%map_cano_stag_fwd_mem
            nf90_stat = nf90_inq_grp_ncid(fgid, 'map_cano_stag_fwd', id_map_cano_stag_fwd)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            call self%map_cano_stag_fwd%read_netcdf(id_map_cano_stag_fwd)
            
            allocate(self%map_cano_stag_bwd_mem)
            self%map_cano_stag_bwd => self%map_cano_stag_bwd_mem
            nf90_stat = nf90_inq_grp_ncid(fgid, 'map_cano_stag_bwd', id_map_cano_stag_bwd)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            call self%map_cano_stag_bwd%read_netcdf(id_map_cano_stag_bwd)
            
            allocate(self%map_stag_stag_fwd_mem)
            self%map_stag_stag_fwd => self%map_stag_stag_fwd_mem
            nf90_stat = nf90_inq_grp_ncid(fgid, 'map_stag_stag_fwd', id_map_stag_stag_fwd)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            call self%map_stag_stag_fwd%read_netcdf(id_map_stag_stag_fwd)
            
            allocate(self%map_stag_stag_bwd_mem)
            self%map_stag_stag_bwd => self%map_stag_stag_bwd_mem
            nf90_stat = nf90_inq_grp_ncid(fgid, 'map_stag_stag_bwd', id_map_stag_stag_bwd)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            call self%map_stag_stag_bwd%read_netcdf(id_map_stag_stag_bwd)
            
        endif
        
        ! Read operator matrices
        if (.not. allocated(self%grad_stag_cano_fwd)) then
            allocate(self%grad_stag_cano_fwd) 
        else
            call handle_error('Already allocated', &
                              GRILLIX_ERR_ALLOC, __LINE__, __FILE__)
        endif
        nf90_stat = nf90_inq_grp_ncid(fgid, 'grad_stag_cano_fwd', id_grad_stag_cano_fwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        call self%grad_stag_cano_fwd%read_netcdf(id_grad_stag_cano_fwd)
        
        if (.not. allocated(self%grad_stag_cano_bwd)) then
            allocate(self%grad_stag_cano_bwd) 
        else
            call handle_error('Already allocated', &
                              GRILLIX_ERR_ALLOC, __LINE__, __FILE__)
        endif
        nf90_stat = nf90_inq_grp_ncid(fgid, 'grad_stag_cano_bwd', id_grad_stag_cano_bwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        call self%grad_stag_cano_bwd%read_netcdf(id_grad_stag_cano_bwd)
        
        if (.not. allocated(self%pdiv_cano_stag_fwd)) then
            allocate(self%pdiv_cano_stag_fwd) 
        else
            call handle_error('Already allocated', &
                              GRILLIX_ERR_ALLOC, __LINE__, __FILE__)
        endif
        nf90_stat = nf90_inq_grp_ncid(fgid, 'pdiv_cano_stag_fwd', id_pdiv_cano_stag_fwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        call self%pdiv_cano_stag_fwd%read_netcdf(id_pdiv_cano_stag_fwd)
        
        if (.not. allocated(self%pdiv_cano_stag_bwd)) then
            allocate(self%pdiv_cano_stag_bwd) 
        else
            call handle_error('Already allocated', &
                              GRILLIX_ERR_ALLOC, __LINE__, __FILE__)
        endif
        nf90_stat = nf90_inq_grp_ncid(fgid, 'pdiv_cano_stag_bwd', id_pdiv_cano_stag_bwd)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        call self%pdiv_cano_stag_bwd%read_netcdf(id_pdiv_cano_stag_bwd)
 
    end subroutine
    
end submodule