variable_netcdf_s.f90 Source File

submodule (variable_m) variable_netcdf_s
    !! Routines related with NETCDF I/O of variables
    implicit none
contains

    module subroutine write_netcdf_variable(self, comm_handler, fgid, varname)
        class(variable_t), intent(in) :: self
        type(comm_handler_t), intent(in) :: comm_handler
        integer, intent(in) :: fgid
        character(len=*), intent(in) :: varname
        
        integer :: stag, npoints_max, l, ierr
        integer :: nf90_stat, id_npoints_max, id_nplanes, id_species, id_var
        integer, dimension(3) :: starts, counts
        real(GP), allocatable, dimension(:) :: var_with_fillers
    
        if (self%staggered) then
            stag = 1
        else
            stag = 0
        endif
        
        ! Define dimensions (Check if previously already defined)                
        nf90_stat =  nf90_inq_dimid(fgid, 'npoints_max', id_npoints_max)
        if (nf90_stat == 0) then
            nf90_stat = nf90_inquire_dimension(fgid, id_npoints_max, len=npoints_max)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        else
            call MPI_allreduce(self%ndim, npoints_max, 1, MPI_INTEGER, MPI_MAX, &
                           comm_handler%get_comm_planes(), ierr)
            nf90_stat = nf90_def_dim(fgid, 'npoints_max', npoints_max, id_npoints_max)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) 
        endif
        
        nf90_stat =  nf90_inq_dimid(fgid, 'planes', id_nplanes)
        if (nf90_stat /= 0) then
            nf90_stat = nf90_def_dim(fgid, 'planes', &
                                     comm_handler%get_nplanes(), id_nplanes)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) 
        endif
        
        nf90_stat =  nf90_inq_dimid(fgid, 'species', id_species)
        if (nf90_stat /= 0) then
            nf90_stat = nf90_def_dim(fgid, 'species', &
                                     comm_handler%get_nspecies(), id_species)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) 
        endif
        
        ! Define variables
        nf90_stat = nf90_def_var(fgid, varname, NF90_DOUBLE, &
                                 [id_npoints_max, id_nplanes, id_species], id_var)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)

        ! Put attribute
        nf90_stat = nf90_put_att(fgid, id_var, 'staggered', stag )
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)

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

        ! Put variables
        starts = [1, comm_handler%get_plane()+1, comm_handler%get_species()+1]
        counts = [npoints_max, 1, 1]
        if (self%ndim == npoints_max) then     
            nf90_stat = nf90_put_var(fgid, id_var, self%vals, &
                                     start=starts, count=counts)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        else
            ! For 3D equilibria create filler points
            allocate(var_with_fillers(npoints_max), source=GP_NAN)
            !$omp parallel private(l)
            !$omp do
            do l = 1, self%ndim
                var_with_fillers(l) = self%vals(l)
            enddo
            !$omp end do
            !$omp end parallel
            nf90_stat = nf90_put_var(fgid, id_var, var_with_fillers, &
                                     start=starts, count=counts)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
            deallocate(var_with_fillers)
        endif
        
        ! Set file into definition state
        nf90_stat = nf90_redef(fgid)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
    
    end subroutine
    
    module subroutine read_netcdf_variable(self, comm_handler, fgid, varname)
        class(variable_t), intent(out) :: self 
        type(comm_handler_t), intent(in) :: comm_handler
        integer, intent(in) :: fgid
        character(len=*), intent(in) :: varname

        integer :: nf90_stat, stag
        integer :: id_var, id_npoints
        integer, dimension(3) :: starts, counts
        integer, dimension(comm_handler%get_nplanes()) :: npoints
        
        nf90_stat = nf90_inq_varid(fgid, varname, id_var)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        
        ! Get staggered
        nf90_stat = nf90_get_att(fgid, id_var, 'staggered', stag )
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        if (stag == 1) then
            self%staggered = .true.
        else
            self%staggered = .false.
        endif
        
        ! Get dimension of variable via npoints_cano or npoints_stag variable
        if (self%staggered) then
            nf90_stat = nf90_inq_varid(fgid, 'npoints_stag', id_npoints)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        else
            nf90_stat = nf90_inq_varid(fgid, 'npoints_cano', id_npoints)
            call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        endif
        nf90_stat = nf90_get_var(fgid, id_npoints, npoints)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
        self%ndim = npoints(comm_handler%get_plane()+1)
        
        ! Get vals
        if (allocated(self%vals)) deallocate(self%vals)
        allocate(self%vals(self%ndim))
        
        starts = [1, comm_handler%get_plane()+1, comm_handler%get_species()+1]
        counts = [self%ndim, 1, 1]
        nf90_stat = nf90_get_var(fgid, id_var, self%vals, &
                                 start=starts, count=counts)
        call handle_error_netcdf(nf90_stat, __LINE__, __FILE__)
     
    end subroutine

end submodule