readcheckpoint.F90

Go to the documentation of this file.

 
module checkpointer_read_checkpoint_mod
#ifndef TEST_MODE
  use netcdf, only : nf90_global, nf90_nowrite, nf90_inquire_attribute, nf90_open, nf90_inq_dimid, nf90_inquire_dimension, &
       nf90_inq_varid, nf90_get_var, nf90_get_att, nf90_close
#else
  use dummy_netcdf_mod, only : nf90_global, nf90_nowrite, nf90_inquire_attribute, nf90_open, nf90_inq_dimid, &
       nf90_inquire_dimension, nf90_inq_varid, nf90_get_var, nf90_get_att, nf90_close
#endif
  use datadefn_mod, only : string_length
  use state_mod, only : model_state_type
  use grids_mod, only : local_grid_type, global_grid_type, x_index, y_index, z_index, primal_grid, dual_grid
  use prognostics_mod, only : prognostic_field_type
  use logging_mod, only : log_info, log_error, log_log, log_master_log
  use conversions_mod, only : conv_is_integer, conv_to_integer, conv_is_real, conv_to_real, conv_is_logical, conv_to_logical, &
       conv_to_string
  use optionsdatabase_mod, only : options_add
  use checkpointer_common_mod, only : empty_dim_key, string_dim_key, x_dim_key, y_dim_key, &
       z_dim_key, q_dim_key, q_key, zq_key, th_key, zth_key, p_key, u_key, v_key, w_key, zu_key, zv_key, zw_key, x_key, y_key, &
       z_key, zn_key, nqfields, ugal, vgal, time_key, timestep, created_attribute_key, title_attribute_key, absolute_new_dtm_key, &
       dtm_key, dtm_new_key, q_indices_dim_key, q_indices_key, q_field_anonymous_name, zq_field_anonymous_name, &
       max_string_length, thref, olubar, olzubar, olvbar, olzvbar, olthbar, olzthbar, olqbar, olzqbar, olqbar_anonymous_name, &
       olzqbar_anonymous_name, rad_last_time_key, sth_lw_key, sth_sw_key, check_status, remove_null_terminator_from_string, &
       wup, wdwn
  use datadefn_mod, only : default_precision
  use q_indices_mod, only : q_metadata_type, set_q_index, get_q_index, get_indices_descriptor, standard_q_names
  implicit none
 
#ifndef TEST_MODE
  private
#endif
 
  public read_checkpoint_file
 
contains
 
  subroutine read_checkpoint_file(current_state, filename)
    type(model_state_type), intent(inout) :: current_state
    character(len=*), intent(in) :: filename
 
    integer :: ncid, z_dim, y_dim, x_dim
    logical :: z_found, y_found, x_found
    character(len=:), allocatable :: attribute_value
 
    call check_status(nf90_open(path = filename, mode = nf90_nowrite, ncid = ncid))
    attribute_value=read_specific_global_attribute(ncid, "created")
    call read_dimensions(ncid, z_dim, y_dim, x_dim, z_found, y_found, x_found)
    call load_global_grid(current_state, ncid, z_dim, y_dim, x_dim, z_found, y_found, x_found)    
 
    call decompose_grid(current_state)
 
    call load_q_indices(ncid)    
    call load_misc(current_state, ncid)
    call load_all_fields(current_state, ncid)
    call initalise_source_and_sav_fields(current_state)  
    call load_mean_profiles(current_state, ncid, z_dim)
    call load_pdf_profiles(current_state, ncid, z_dim)
    call check_status(nf90_close(ncid))
 
    if (current_state%number_q_fields .gt. 0) then
      ! Retrieve standard q indices
      current_state%water_vapour_mixing_ratio_index=get_q_index(standard_q_names%VAPOUR, 'checkpoint')
      current_state%liquid_water_mixing_ratio_index=get_q_index(standard_q_names%CLOUD_LIQUID_MASS, 'checkpoint')
    end if
 
    call log_master_log(log_info, "Restarted configuration from checkpoint file `"//trim(filename)//"` created at "&
         //attribute_value)
    deallocate(attribute_value)
    current_state%initialised=.true.
  end subroutine read_checkpoint_file
 
  subroutine decompose_grid(current_state)
    type(model_state_type), intent(inout) :: current_state
 
    if (associated(current_state%parallel%decomposition_procedure)) then
      call current_state%parallel%decomposition_procedure(current_state)
    else
      call log_log(log_error, "No decomposition specified")
    end if
  end subroutine decompose_grid
 
  subroutine initalise_source_and_sav_fields(current_state)
    type(model_state_type), intent(inout) :: current_state
 
    integer :: x_size, y_size, z_size, i
 
    z_size = current_state%local_grid%size(z_index) + current_state%local_grid%halo_size(z_index) * 2
    y_size = current_state%local_grid%size(y_index) + current_state%local_grid%halo_size(y_index) * 2
    x_size = current_state%local_grid%size(x_index) + current_state%local_grid%halo_size(x_index) * 2
#ifdef U_ACTIVE
    allocate(current_state%su%data(z_size, y_size, x_size))
    current_state%su%data(:,:,:) = 0.
    current_state%su%active=.true.
    allocate(current_state%savu%data(z_size, y_size, x_size))
    current_state%savu%data(:,:,:) = 0.
    current_state%savu%active=.true.
#endif
#ifdef V_ACTIVE
    allocate(current_state%sv%data(z_size, y_size, x_size))
    current_state%sv%data(:,:,:) = 0.
    current_state%sv%active=.true.
    allocate(current_state%savv%data(z_size, y_size, x_size))
    current_state%savv%data(:,:,:) = 0.
    current_state%savv%active=.true.
#endif
#ifdef W_ACTIVE
    allocate(current_state%sw%data(z_size, y_size, x_size))
    current_state%sw%data(:,:,:) = 0.
    current_state%sw%active=.true.
    allocate(current_state%savw%data(z_size, y_size, x_size))
    current_state%savw%data(:,:,:) = 0.
    current_state%savw%active=.true.
#endif
    if (current_state%th%active) then
      allocate(current_state%sth%data(z_size, y_size, x_size))
      current_state%sth%data(:,:,:) = 0.
      current_state%sth%active=.true.
    end if
    do i=1,current_state%number_q_fields
      current_state%sq(i)%active=.true.
      allocate(current_state%sq(i)%data(z_size, y_size, x_size))
      current_state%sq(i)%data(:,:,:) = 0.
    end do
  end subroutine initalise_source_and_sav_fields
 
  subroutine load_misc(current_state, ncid)
    type(model_state_type), intent(inout) :: current_state
    integer, intent(in) :: ncid
 
    integer :: i_data(1)  ! Procedure requires a vector rather than scalar
    real(kind=default_precision) :: r_data(1)
 
    call read_single_variable(ncid, timestep, integer_data_1d=i_data)
    current_state%timestep = i_data(1)+1 ! plus one to increment for next timestep
    !current_state%start_timestep = current_state%timestep 
    call read_single_variable(ncid, ugal, real_data_1d_double=r_data)
    current_state%ugal = r_data(1)
    call read_single_variable(ncid, vgal, real_data_1d_double=r_data)
    current_state%vgal = r_data(1)
    call read_single_variable(ncid, nqfields, integer_data_1d=i_data)
    current_state%number_q_fields = i_data(1)
    call read_single_variable(ncid, dtm_key, real_data_1d_double=r_data)
    current_state%dtm = r_data(1)
    call read_single_variable(ncid, dtm_new_key, real_data_1d_double=r_data)
    current_state%dtm_new = r_data(1)
    current_state%update_dtm = current_state%dtm .ne. current_state%dtm_new
    call read_single_variable(ncid, absolute_new_dtm_key, real_data_1d_double=r_data)
    current_state%absolute_new_dtm = r_data(1)
    call read_single_variable(ncid, time_key, real_data_1d_double=r_data)
    ! The time is written into checkpoint as time+dtm, therefore the time as read in has been correctly advanced
    current_state%time = r_data(1)
    call read_single_variable(ncid, rad_last_time_key, real_data_1d_double=r_data)
    current_state%rad_last_time = r_data(1)
  end subroutine load_misc
 
  function read_specific_global_attribute(ncid, key)
    integer, intent(in) :: ncid
    character(len=*), intent(in) :: key
 
    integer :: length
    character(len=:),allocatable,target :: read_specific_global_attribute
 
    call check_status(nf90_inquire_attribute(ncid, nf90_global, key, len = length))
    allocate(character(length) :: read_specific_global_attribute)
    call check_status(nf90_get_att(ncid, nf90_global, key, read_specific_global_attribute))
  end function read_specific_global_attribute
 
  subroutine load_all_fields(current_state, ncid)
    type(model_state_type), intent(inout) :: current_state
    integer, intent(in) :: ncid
 
    integer :: i
    logical :: multi_process
    type(q_metadata_type) :: q_metadata
    character(len=STRING_LENGTH) :: q_field_name, zq_field_name
 
    multi_process = current_state%parallel%processes .gt. 1
 
    if (does_field_exist(ncid, u_key)) then
      call load_single_3d_field(ncid, current_state%local_grid, current_state%u, dual_grid, &
           dual_grid, primal_grid, u_key, multi_process)
      call load_single_3d_field(ncid, current_state%local_grid, current_state%zu, dual_grid, &
           dual_grid, primal_grid, zu_key, multi_process)
    end if
    if (does_field_exist(ncid, v_key)) then
      call load_single_3d_field(ncid, current_state%local_grid, current_state%v, dual_grid, &
           primal_grid, dual_grid, v_key, multi_process)
      call load_single_3d_field(ncid, current_state%local_grid, current_state%zv, dual_grid, &
           primal_grid, dual_grid, zv_key, multi_process)
    end if
    if (does_field_exist(ncid, w_key)) then
      call load_single_3d_field(ncid, current_state%local_grid, current_state%w, primal_grid, &
           dual_grid, dual_grid, w_key, multi_process)
      call load_single_3d_field(ncid, current_state%local_grid, current_state%zw, primal_grid, &
           dual_grid, dual_grid, zw_key, multi_process)
    end if
    if (does_field_exist(ncid, th_key)) then
      call load_single_3d_field(ncid, current_state%local_grid, current_state%th, dual_grid, &
           dual_grid, dual_grid, th_key, multi_process)
      call load_single_3d_field(ncid, current_state%local_grid, current_state%zth, dual_grid, &
           dual_grid, dual_grid, zth_key, multi_process)
    end if
    if (does_field_exist(ncid, p_key)) then
      call load_single_3d_field(ncid, current_state%local_grid, current_state%p, dual_grid, &
           dual_grid, dual_grid, p_key, multi_process)
    end if    
    allocate(current_state%q(current_state%number_q_fields), current_state%zq(current_state%number_q_fields), &
         current_state%sq(current_state%number_q_fields))
    if (does_field_exist(ncid, q_key)) then
      do i=1,current_state%number_q_fields
        call load_single_3d_field(ncid, current_state%local_grid, current_state%q(i), dual_grid, &
             dual_grid, dual_grid, q_key, multi_process, i)
        call load_single_3d_field(ncid, current_state%local_grid, current_state%zq(i), dual_grid, &
             dual_grid, dual_grid, zq_key, multi_process, i)
      end do
    else
      do i=1,current_state%number_q_fields
        q_metadata=get_indices_descriptor(i)
        q_field_name=trim(q_key)//"_"//trim(q_metadata%name)
        zq_field_name=trim(zq_key)//"_"//trim(q_metadata%name)
        if (.not. does_field_exist(ncid, q_field_name)) then
          q_field_name=trim(q_field_anonymous_name)//"_"//trim(conv_to_string(i))
          zq_field_name=trim(zq_field_anonymous_name)//"_"//trim(conv_to_string(i))
          if (.not. does_field_exist(ncid, q_field_name)) then
            call log_log(log_error, "No entry in checkpoint file for Q field "//trim(conv_to_string(i)))
          end if          
        end if
        if (.not. does_field_exist(ncid, zq_field_name)) then
          call log_log(log_error, "Missmatch between q and zq field name in the checkpoint file")
        end if
        call load_single_3d_field(ncid, current_state%local_grid, current_state%q(i), dual_grid, &
             dual_grid, dual_grid, q_field_name, multi_process)
        call load_single_3d_field(ncid, current_state%local_grid, current_state%zq(i), dual_grid, &
             dual_grid, dual_grid, zq_field_name, multi_process)
      end do      
    end if
    if (does_field_exist(ncid, sth_lw_key)) then
       call load_single_3d_field(ncid, current_state%local_grid, current_state%sth_lw, dual_grid, &
            dual_grid, dual_grid, sth_lw_key, multi_process)
    end if
    if (does_field_exist(ncid, sth_lw_key)) then
       call load_single_3d_field(ncid, current_state%local_grid, current_state%sth_sw, dual_grid, &
            dual_grid, dual_grid, sth_sw_key, multi_process)
    endif
  end subroutine load_all_fields
 
  logical function does_field_exist(ncid, variable_key)
    integer, intent(in) :: ncid
    character(len=*), intent(in) :: variable_key
 
    integer :: variable_id
 
    call check_status(nf90_inq_varid(ncid, variable_key, variable_id), does_field_exist)
  end function does_field_exist
 
  subroutine load_q_indices(ncid)
    integer, intent(in) :: ncid
 
    integer :: number_q_indices, i, q_indices_id
    character(len=MAX_STRING_LENGTH) :: key, value
 
    number_q_indices=get_number_q_indices(ncid)
    if (number_q_indices .gt. 0) then
      call check_status(nf90_inq_varid(ncid, q_indices_key, q_indices_id))
      do i=1, number_q_indices
        call check_status(nf90_get_var(ncid, q_indices_id, key, (/ 1, 1, i /)))
        call check_status(nf90_get_var(ncid, q_indices_id, value, (/ 1, 2, i /)))
        call remove_null_terminator_from_string(key)
        call remove_null_terminator_from_string(value)
        call set_q_index(conv_to_integer(trim(value)), key)
      end do
    end if
  end subroutine load_q_indices
 
  integer function get_number_q_indices(ncid)
    integer, intent(in) :: ncid
 
    integer :: q_indices_dimid, q_indices_dim
    logical :: found_flag
 
    call check_status(nf90_inq_dimid(ncid, q_indices_dim_key, q_indices_dimid), found_flag)
    if (found_flag) then
      call check_status(nf90_inquire_dimension(ncid, q_indices_dimid, len=q_indices_dim))
      get_number_q_indices=q_indices_dim
    else
      get_number_q_indices=0
    end if
  end function get_number_q_indices  
 
  subroutine load_single_3d_field(ncid, local_grid, field, z_grid, y_grid, x_grid, variable_key, multi_process, fourth_dim_loc)
    character(len=*), intent(in) :: variable_key
    integer, intent(in) :: z_grid, y_grid, x_grid, ncid
    type(prognostic_field_type), intent(inout) :: field
    type(local_grid_type), intent(inout) :: local_grid
    logical, intent(in) :: multi_process
    integer, optional, intent(in) :: fourth_dim_loc
 
    integer :: start(5), count(5), i, map(5)
 integer :: variable_id, nd
 
    if (allocated(field%data)) deallocate(field%data)
    allocate(field%data(local_grid%size(z_index) + local_grid%halo_size(z_index) * 2, local_grid%size(y_index) + &
         local_grid%halo_size(y_index) * 2, local_grid%size(x_index) + local_grid%halo_size(x_index) * 2))
    field%data(:,:,:)=0.
 
    if (multi_process .or. present(fourth_dim_loc)) then
      do i=1,3
        if (i==1) then
          map(i)=1
        else
          map(i)=map(i-1)*local_grid%size(i-1)
        end if
        start(i) = local_grid%start(i)
        count(i) = local_grid%size(i)
      end do
      if (present(fourth_dim_loc)) then
        start(4) = fourth_dim_loc
        count(4) = 1
        map(4)=map(3)*local_grid%size(3)
 
        start(5)=1
        count(5)=1
        map(5)=map(4)
      else
        start(4)=1
        map(4)=map(3)*local_grid%size(3)
        count(4)=1
      end if
 
      call read_single_variable(ncid, variable_key, real_data_3d=field%data(local_grid%local_domain_start_index(z_index):&
           local_grid%local_domain_end_index(z_index),local_grid%local_domain_start_index(y_index):&
           local_grid%local_domain_end_index(y_index), local_grid%local_domain_start_index(x_index):&
           local_grid%local_domain_end_index(x_index)), start=start, count=count, map=map)
    else
      call read_single_variable(ncid, variable_key, real_data_3d=field%data(local_grid%local_domain_start_index(z_index):&
           local_grid%local_domain_end_index(z_index),local_grid%local_domain_start_index(y_index):&
           local_grid%local_domain_end_index(y_index), local_grid%local_domain_start_index(x_index):&
           local_grid%local_domain_end_index(x_index)))
    end if
 
    field%grid(z_index) = z_grid
    field%grid(y_index) = y_grid
    field%grid(x_index) = x_grid
    field%active = .true.
  end subroutine load_single_3d_field
 
  subroutine load_global_grid(current_state, ncid, z_dim, y_dim, x_dim, z_found, y_found, x_found)
    type(model_state_type), intent(inout) :: current_state
    integer, intent(in) :: ncid, z_dim, y_dim, x_dim
    logical, intent(in) :: z_found, y_found, x_found
 
    current_state%global_grid%dimensions = 0
 
    if (z_found) then
      call read_dimension_of_grid(ncid, current_state%global_grid, z_key, z_index, z_dim)
      call define_vertical_levels(ncid, current_state, z_key, z_dim)
      
    end if
    if (y_found) call read_dimension_of_grid(ncid, current_state%global_grid, y_key, y_index, y_dim)
    if (x_found) call read_dimension_of_grid(ncid, current_state%global_grid, x_key, x_index, x_dim)
  end subroutine load_global_grid
 
  subroutine load_mean_profiles(current_state, ncid, z_dim_id)
    type(model_state_type), intent(inout) :: current_state
    integer, intent(in) :: ncid, z_dim_id
 
    integer :: z_size, i
    type(q_metadata_type) :: q_metadata
    character(len=STRING_LENGTH) :: q_field_name, zq_field_name
 
    call check_status(nf90_inquire_dimension(ncid, z_dim_id, len=z_size))
    if (does_field_exist(ncid, olubar)) then
      allocate(current_state%global_grid%configuration%vertical%olubar(z_size))
      call read_single_variable(ncid, olubar, real_data_1d_double=current_state%global_grid%configuration%vertical%olubar)
    end if
    if (does_field_exist(ncid, olzubar)) then
      allocate(current_state%global_grid%configuration%vertical%olzubar(z_size))
      call read_single_variable(ncid, olzubar, real_data_1d_double=current_state%global_grid%configuration%vertical%olzubar)
    end if
    if (does_field_exist(ncid, olvbar)) then
      allocate(current_state%global_grid%configuration%vertical%olvbar(z_size))
      call read_single_variable(ncid, olvbar, real_data_1d_double=current_state%global_grid%configuration%vertical%olvbar)
    end if
    if (does_field_exist(ncid, olzvbar)) then
      allocate(current_state%global_grid%configuration%vertical%olzvbar(z_size))
      call read_single_variable(ncid, olzvbar, real_data_1d_double=current_state%global_grid%configuration%vertical%olzvbar)
    end if
    if (does_field_exist(ncid, olthbar)) then
      allocate(current_state%global_grid%configuration%vertical%olthbar(z_size))
      call read_single_variable(ncid, olthbar, real_data_1d_double=current_state%global_grid%configuration%vertical%olthbar)
    end if
    if (does_field_exist(ncid, olzthbar)) then
      allocate(current_state%global_grid%configuration%vertical%olzthbar(z_size))
      call read_single_variable(ncid, olzthbar, real_data_1d_double=current_state%global_grid%configuration%vertical%olzthbar)
    end if
   if (does_field_exist(ncid, olqbar)) then
     allocate(current_state%global_grid%configuration%vertical%olqbar(z_size, current_state%number_q_fields))
     call read_single_variable(ncid, olqbar, real_data_2d_double=current_state%global_grid%configuration%vertical%olqbar)
   else if (current_state%number_q_fields .gt. 0) then
     do i=1,current_state%number_q_fields
        q_metadata=get_indices_descriptor(i)
        q_field_name=trim(olqbar)//"_"//trim(q_metadata%name)
        if (.not. does_field_exist(ncid, q_field_name)) then
          q_field_name=trim(olqbar_anonymous_name)//"_"//trim(conv_to_string(i))          
          if (.not. does_field_exist(ncid, q_field_name)) then
            cycle
          end if        
        end if
        if (.not. allocated(current_state%global_grid%configuration%vertical%olqbar)) then
          allocate(current_state%global_grid%configuration%vertical%olqbar(z_size, current_state%number_q_fields))
        end if
        call read_single_variable(ncid, q_field_name, &
             real_data_1d_double=current_state%global_grid%configuration%vertical%olqbar(:, i))        
      end do
    end if
   if (does_field_exist(ncid, olzqbar)) then
     allocate(current_state%global_grid%configuration%vertical%olzqbar(z_size, current_state%number_q_fields))
     call read_single_variable(ncid, olzqbar, real_data_2d_double=current_state%global_grid%configuration%vertical%olzqbar)
   else if (current_state%number_q_fields .gt. 0) then
     do i=1,current_state%number_q_fields
        q_metadata=get_indices_descriptor(i)
        q_field_name=trim(olzqbar)//"_"//trim(q_metadata%name)
        if (.not. does_field_exist(ncid, q_field_name)) then
          q_field_name=trim(olzqbar_anonymous_name)//"_"//trim(conv_to_string(i))          
          if (.not. does_field_exist(ncid, q_field_name)) then
            cycle
          end if        
        end if
        if (.not. allocated(current_state%global_grid%configuration%vertical%olzqbar)) then
          allocate(current_state%global_grid%configuration%vertical%olzqbar(z_size, current_state%number_q_fields))
        end if
        call read_single_variable(ncid, q_field_name, &
             real_data_1d_double=current_state%global_grid%configuration%vertical%olzqbar(:, i))        
      end do
    end if
  end subroutine load_mean_profiles  
 
 
  subroutine load_pdf_profiles(current_state, ncid, z_dim_id)
    type(model_state_type), intent(inout) :: current_state
    integer, intent(in) :: ncid, z_dim_id
 
    integer :: z_size, i
 
    call check_status(nf90_inquire_dimension(ncid, z_dim_id, len=z_size))
    if (does_field_exist(ncid, wup)) then
      allocate(current_state%global_grid%configuration%vertical%w_up(z_size))
      call read_single_variable(ncid, wup, real_data_1d_double=current_state%global_grid%configuration%vertical%w_up)
    end if
    if (does_field_exist(ncid, wdwn)) then
      allocate(current_state%global_grid%configuration%vertical%w_dwn(z_size))
      call read_single_variable(ncid, wdwn, real_data_1d_double=current_state%global_grid%configuration%vertical%w_dwn)
    end if
  end subroutine load_pdf_profiles
 
 
  subroutine define_vertical_levels(ncid, current_state, z_key, z_dim_id)
    type(model_state_type), intent(inout) :: current_state
    character(len=*), intent(in) :: z_key
    integer, intent(in) :: z_dim_id, ncid
 
    integer :: i, z_size
    real, dimension(:), allocatable :: data
 
    call check_status(nf90_inquire_dimension(ncid, z_dim_id, len=z_size))
    allocate(data(z_size))
    call read_single_variable(ncid, z_key, real_data_1d=data)
 
    allocate(current_state%global_grid%configuration%vertical%kgd(z_size), &
         current_state%global_grid%configuration%vertical%hgd(z_size), &
         current_state%global_grid%configuration%vertical%thref(z_size))
 
    call read_single_variable(ncid, thref, real_data_1d_double=current_state%global_grid%configuration%vertical%thref)
 
    do i=1,z_size
      current_state%global_grid%configuration%vertical%kgd(i) = i
      current_state%global_grid%configuration%vertical%hgd(i) = real(data(i))
    end do
    deallocate(data)    
  end subroutine define_vertical_levels
 
  subroutine read_dimension_of_grid(ncid, grid, variable_key, dimension, dimension_id)
    integer, intent(in) :: ncid, dimension_id, dimension
    character(len=*), intent(in) :: variable_key
    type(global_grid_type), intent(inout) :: grid
 
    integer :: dim_size
 
    call check_status(nf90_inquire_dimension(ncid, dimension_id, len=dim_size))
 
    if (variable_key .eq. "x" .or. variable_key .eq. "y") then
      call read_single_variable(ncid, trim(variable_key)//"_resolution", real_data=grid%resolution(dimension))
      call read_single_variable(ncid, trim(variable_key)//"_top", real_data=grid%top(dimension))
      call read_single_variable(ncid, trim(variable_key)//"_bottom", real_data=grid%bottom(dimension))
    else if (variable_key .eq. "z") then
      allocate(grid%configuration%vertical%z(dim_size), grid%configuration%vertical%zn(dim_size))
      call read_single_variable(ncid, z_key, real_data_1d_double=grid%configuration%vertical%z)
      call read_single_variable(ncid, zn_key, real_data_1d_double=grid%configuration%vertical%zn)
      grid%top(dimension) = int(grid%configuration%vertical%z(dim_size))
      grid%resolution(dimension) = int(grid%configuration%vertical%z(2) - grid%configuration%vertical%z(1))
      ! For now hard code the bottom of the grid in each dimension as 0, first element is the 0th + size
      ! I.e. if dxx=100, start=0 then first point is 100 rather than 0 (in x and y), is correct in z in CP file
      grid%bottom(dimension) = 0
    end if            
    grid%size(dimension) = dim_size
    grid%dimensions = grid%dimensions + 1
    grid%active(dimension) = .true.    
  end subroutine read_dimension_of_grid
 
  subroutine read_single_variable(ncid, key, int_data, real_data, real_data_1d, real_data_1d_double, real_data_2d_double, &
       real_data_3d, integer_data_1d, start, count, map)
    integer, intent(in) :: ncid
    character(len=*), intent(in) :: key
    integer, intent(inout), optional :: int_data
    real(kind=default_precision) , intent(inout), optional :: real_data
    real, dimension(:), intent(inout), optional :: real_data_1d
    real(kind=default_precision), dimension(:), intent(inout), optional :: real_data_1d_double
    real(kind=default_precision), dimension(:,:), intent(inout), optional :: real_data_2d_double
    real(kind=default_precision), dimension(:,:,:), intent(inout), optional :: real_data_3d
    integer, dimension(:), intent(inout), optional :: integer_data_1d
    integer, dimension(:), intent(in), optional :: start, count, map
 
    integer :: variable_id
 
    call check_status(nf90_inq_varid(ncid, key, variable_id))
 
    if (.not. present(int_data) .and. .not. present(real_data) .and. .not. present(real_data_1d) .and. &
         .not. present(real_data_1d_double) .and. .not. present(real_data_2d_double) .and. &
         .not. present(real_data_3d) .and. .not. present(integer_data_1d)) return
 
    if (present(int_data)) then
      call check_status(nf90_get_var(ncid, variable_id, int_data))
      return
    end if  
 
    if (present(real_data)) then
      call check_status(nf90_get_var(ncid, variable_id, real_data))
      return
    end if
 
    if (present(real_data_1d)) then
      if (present(start) .and. present(count) .and. present(map)) then
        call check_status(nf90_get_var(ncid, variable_id, real_data_1d, start=start, count=count, map=map))
      else if (present(start) .and. present(count)) then
        call check_status(nf90_get_var(ncid, variable_id, real_data_1d, start=start, count=count))
      else
        call check_status(nf90_get_var(ncid, variable_id, real_data_1d))
      end if
    else if (present(real_data_1d_double)) then
      if (present(start) .and. present(count) .and. present(map)) then
        call check_status(nf90_get_var(ncid, variable_id, real_data_1d_double, start=start, count=count, map=map))
      else if (present(start) .and. present(count)) then
        call check_status(nf90_get_var(ncid, variable_id, real_data_1d_double, start=start, count=count))
      else
        call check_status(nf90_get_var(ncid, variable_id, real_data_1d_double))
      end if
    else if (present(real_data_2d_double)) then
      if (present(start) .and. present(count) .and. present(map)) then
        call check_status(nf90_get_var(ncid, variable_id, real_data_2d_double, start=start, count=count, map=map))
      else if (present(start) .and. present(count)) then
        call check_status(nf90_get_var(ncid, variable_id, real_data_2d_double, start=start, count=count))
      else
        call check_status(nf90_get_var(ncid, variable_id, real_data_2d_double))
      end if
    else if (present(real_data_3d)) then
      if (present(start) .and. present(count) .and. present(map)) then
        call check_status(nf90_get_var(ncid, variable_id, real_data_3d, start=start, count=count, map=map))
      else if (present(start) .and. present(count)) then
        call check_status(nf90_get_var(ncid, variable_id, real_data_3d, start=start, count=count))
      else
        call check_status(nf90_get_var(ncid, variable_id, real_data_3d))
      end if
    else if (present(integer_data_1d)) then
      if (present(start) .and. present(count) .and. present(map)) then
        call check_status(nf90_get_var(ncid, variable_id, integer_data_1d, start, count, map))
      else if (present(start) .and. present(count)) then
        call check_status(nf90_get_var(ncid, variable_id, integer_data_1d, start, count))
      else
        call check_status(nf90_get_var(ncid, variable_id, integer_data_1d))
      end if
    end if
  end subroutine read_single_variable
 
  subroutine read_dimensions(ncid, z_dim, y_dim, x_dim, z_found, y_found, x_found)
    integer, intent(in) :: ncid
    integer, intent(out) :: z_dim, y_dim, x_dim
    logical, intent(out) :: z_found, y_found, x_found
 
    call check_status(nf90_inq_dimid(ncid, z_dim_key, z_dim), z_found)
    call check_status(nf90_inq_dimid(ncid, y_dim_key, y_dim), y_found)
    call check_status(nf90_inq_dimid(ncid, x_dim_key, x_dim), x_found)
  end subroutine read_dimensions
end module checkpointer_read_checkpoint_mod