mesh_cart_netcdfio_s.f90 Source File


Source Code

submodule(mesh_cart_m) mesh_2d_netcdfio_s
    use screen_io_m, only : nf90_handle_err

    implicit none
    !! Submodule related with I/O to netcdf of mesh

contains

    module subroutine write_netcdf_mesh_cart(self, fgid)
        class(mesh_cart_t), intent(inout) :: self
        integer, intent(in) :: fgid

        integer :: nf90_stat, yperiodic_int, extend_beyond_wall_int
        integer :: iddim_n_points, iddim_n_points_inner, iddim_n_points_boundary, iddim_n_points_ghost, iddim_size_neighbor
        integer :: id_cart_i, id_cart_j, id_pinfo, id_district
        integer :: id_inner_indices, id_boundary_indices, id_ghost_indices, id_index_neighbor, id_redblack_indices

        character(len=*), parameter :: cname = "write_netcdf_mesh_cart"
        ! Name used for logging errors

        ! Write global attributes (metadata) ------------------------------------------------------
        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'phi', self%phi)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'lvl', self%lvl)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'lvst', self%lvst)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'spacing_f', self%spacing_f)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'spacing_c', self%spacing_c)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'xmin', self%xmin)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'ymin', self%ymin)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'nx_f', self%nx_f)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'ny_f', self%ny_f)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'size_ghost_layer', self%size_ghost_layer)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        if (self%yperiodic) then
            yperiodic_int = 1
        else
            yperiodic_int = 0
        endif
        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'yperiodic', yperiodic_int)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        if (self%extend_beyond_wall) then
            extend_beyond_wall_int = 1
        else
            extend_beyond_wall_int = 0
        endif
        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'extend_beyond_wall', &
                                 extend_beyond_wall_int)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'n_points_red', self%n_points_red)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'n_points_black', self%n_points_black)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        ! Define dimensions -----------------------------------------------------------------------
        nf90_stat = nf90_def_dim(fgid, 'n_points', self%n_points, iddim_n_points)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_def_dim(fgid, 'n_points_inner', self%n_points_inner, iddim_n_points_inner)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_def_dim(fgid, 'n_points_boundary', self%n_points_boundary, iddim_n_points_boundary)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_def_dim(fgid, 'n_points_ghost', self%n_points_ghost, iddim_n_points_ghost)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_def_dim(fgid, 'size_neighbor', 2*self%size_neighbor+1, iddim_size_neighbor)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        ! Define variables ------------------------------------------------------------------------
        nf90_stat = nf90_def_var(fgid, 'cart_i', NF90_INT, iddim_n_points, id_cart_i )
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_def_var(fgid, 'cart_j', NF90_INT, iddim_n_points, id_cart_j )
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_def_var(fgid, 'pinfo', NF90_INT, iddim_n_points, id_pinfo )
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_def_var(fgid, 'district', NF90_INT, iddim_n_points, id_district )
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_def_var(fgid, 'inner_indices', NF90_INT, iddim_n_points_inner, id_inner_indices)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_def_var(fgid, 'boundary_indices', NF90_INT, iddim_n_points_boundary, id_boundary_indices)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_def_var(fgid, 'ghost_indices', NF90_INT, iddim_n_points_ghost, id_ghost_indices)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_def_var(fgid, 'index_neighbor', NF90_INT, &
                                   (/ iddim_size_neighbor, iddim_size_neighbor, iddim_n_points /), &
                                    id_index_neighbor )
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_def_var(fgid, 'redblack_indices', NF90_INT, iddim_n_points, id_redblack_indices)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        ! end of definition
        nf90_stat = nf90_enddef(fgid)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        ! Put variables ---------------------------------------------------------------------------
        nf90_stat = nf90_put_var(fgid, id_cart_i, self%cart_i )
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_var(fgid, id_cart_j, self%cart_j )
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_var(fgid, id_pinfo, self%pinfo )
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_var(fgid, id_district, self%district )
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_var(fgid, id_inner_indices, self%inner_indices )
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_var(fgid, id_boundary_indices, self%boundary_indices )
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_var(fgid, id_ghost_indices, self%ghost_indices )
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_var(fgid, id_index_neighbor, self%index_neighbor )
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_put_var(fgid, id_redblack_indices, self%redblack_indices )
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

    end subroutine


    module subroutine read_netcdf_mesh_cart(self, fgid)
        class(mesh_cart_t), intent(inout) :: self
        integer, intent(in) :: fgid

        integer :: nf90_stat, yperiodic_int, extend_beyond_wall_int, ngb
        integer :: iddim_n_points, iddim_n_points_inner, iddim_n_points_boundary, iddim_n_points_ghost, iddim_size_neighbor
        integer :: id_cart_i, id_cart_j, id_pinfo, id_district
        integer :: id_inner_indices, id_boundary_indices, id_ghost_indices, id_index_neighbor, id_redblack_indices

        character(len=*), parameter :: cname = "read_netcdf_mesh_cart"
        ! Name used for logging errors

        ! Read global attributes (metadata) -------------------------------------------------------
        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'phi', self%phi)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'lvl', self%lvl)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'lvst', self%lvst)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'spacing_f', self%spacing_f)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'spacing_c', self%spacing_c)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'xmin', self%xmin)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'ymin', self%ymin)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'nx_f', self%nx_f)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'ny_f', self%ny_f)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'size_ghost_layer', self%size_ghost_layer)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'yperiodic', yperiodic_int)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)
        if (yperiodic_int == 0) then
            self%yperiodic = .false.
        else
            self%yperiodic = .true.
        endif
        
        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'extend_beyond_wall', &
                                 extend_beyond_wall_int)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)
        if (extend_beyond_wall_int == 0) then
            self%extend_beyond_wall = .false.
        else
            self%extend_beyond_wall = .true.
        endif

        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'n_points_red', self%n_points_red)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'n_points_black', self%n_points_black)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        ! Get dimensions --------------------------------------------------------------------------
        nf90_stat = nf90_inq_dimid(fgid, 'n_points', iddim_n_points)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)
        nf90_stat =  nf90_inquire_dimension(fgid, iddim_n_points, len = self%n_points)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_inq_dimid(fgid, 'n_points_inner', iddim_n_points_inner)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)
        nf90_stat =  nf90_inquire_dimension(fgid, iddim_n_points_inner, len = self%n_points_inner)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_inq_dimid(fgid, 'n_points_boundary', iddim_n_points_boundary)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)
        nf90_stat =  nf90_inquire_dimension(fgid, iddim_n_points_boundary, len = self%n_points_boundary)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_inq_dimid(fgid, 'n_points_ghost', iddim_n_points_ghost)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)
        nf90_stat =  nf90_inquire_dimension(fgid, iddim_n_points_ghost, len = self%n_points_ghost)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_inq_dimid(fgid, 'size_neighbor', iddim_size_neighbor)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)
        nf90_stat =  nf90_inquire_dimension(fgid, iddim_size_neighbor, len = ngb)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)
        if (ngb > 0) then
            ngb = (ngb-1)/2
        else
            ngb = 0
        endif
        self%size_neighbor = ngb

       ! Allocate fields -------------------------------------------------------------------------
        if (allocated(self%cart_i)) deallocate(self%cart_i)
        allocate( self%cart_i( self%n_points ) )

        if (allocated(self%cart_j)) deallocate(self%cart_j)
        allocate( self%cart_j( self%n_points ) )

        if (allocated(self%pinfo)) deallocate(self%pinfo)
        allocate( self%pinfo( self%n_points ) )

        if (allocated(self%district)) deallocate(self%district)
        allocate( self%district( self%n_points ) )

        if (allocated(self%inner_indices)) deallocate(self%inner_indices)
        allocate( self%inner_indices( self%n_points_inner ) )

        if (allocated(self%boundary_indices)) deallocate(self%boundary_indices)
        allocate( self%boundary_indices( self%n_points_boundary ) )

        if (allocated(self%ghost_indices)) deallocate(self%ghost_indices)
        allocate( self%ghost_indices( self%n_points_ghost ) )

        if (allocated(self%index_neighbor)) deallocate(self%index_neighbor)
        allocate( self%index_neighbor( -ngb:ngb, -ngb:ngb, self%n_points ) )

        if (allocated(self%redblack_indices)) deallocate(self%redblack_indices)
        allocate( self%redblack_indices( self%n_points ) )

        ! Get variables ---------------------------------------------------------------------------
        nf90_stat = nf90_inq_varid(fgid, 'cart_i', id_cart_i)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)
        nf90_stat = nf90_get_var(fgid, id_cart_i, self%cart_i)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_inq_varid(fgid, 'cart_j', id_cart_j)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)
        nf90_stat = nf90_get_var(fgid, id_cart_j, self%cart_j)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_inq_varid(fgid, 'pinfo', id_pinfo)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)
        nf90_stat = nf90_get_var(fgid, id_pinfo, self%pinfo)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_inq_varid(fgid, 'district', id_district)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)
        nf90_stat = nf90_get_var(fgid, id_district, self%district)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_inq_varid(fgid, 'inner_indices', id_inner_indices)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)
        nf90_stat = nf90_get_var(fgid, id_inner_indices, self%inner_indices)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_inq_varid(fgid, 'boundary_indices', id_boundary_indices)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)
        nf90_stat = nf90_get_var(fgid, id_boundary_indices, self%boundary_indices)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_inq_varid(fgid, 'ghost_indices', id_ghost_indices)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)
        nf90_stat = nf90_get_var(fgid, id_ghost_indices, self%ghost_indices)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_inq_varid(fgid, 'index_neighbor', id_index_neighbor)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)
        nf90_stat = nf90_get_var(fgid, id_index_neighbor, self%index_neighbor)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

        nf90_stat = nf90_inq_varid(fgid, 'redblack_indices', id_redblack_indices)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)
        nf90_stat = nf90_get_var(fgid, id_redblack_indices, self%redblack_indices)
        if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname)

    end subroutine

end submodule