anchor.diff:          write(*,*) "Reading atmospheric restart data"
           write(*,*) "   timestep:",trim(str(timeslice))," from file:",trim(restart_file)
           call io_read3d(restart_file,"qv",inputdata,timeslice)
  +        call swap_y_z_dimensions(inputdata)
           call check_shapes_3d(inputdata,domain%qv)
           domain%qv=inputdata
           deallocate(inputdata)
  +        
           call io_read3d(restart_file,"qc",inputdata,timeslice)
  +        call swap_y_z_dimensions(inputdata)
           domain%cloud=inputdata
           deallocate(inputdata)
  +
           call io_read3d(restart_file,"qr",inputdata,timeslice)
  +        call swap_y_z_dimensions(inputdata)
           domain%qrain=inputdata
           deallocate(inputdata)
           call io_read3d(restart_file,"qi",inputdata,timeslice)
  +        call swap_y_z_dimensions(inputdata)
           domain%ice=inputdata
           deallocate(inputdata)
           call io_read3d(restart_file,"qs",inputdata,timeslice)
  +        call swap_y_z_dimensions(inputdata)
           domain%qsnow=inputdata
           deallocate(inputdata)
           if (io_variable_is_present(restart_file,"qg")) then
               call io_read3d(restart_file,"qg",inputdata,timeslice)
  +            call swap_y_z_dimensions(inputdata)
               domain%qgrau=inputdata
               deallocate(inputdata)
           endif
           if (io_variable_is_present(restart_file,"nr")) then
               call io_read3d(restart_file,"nr",inputdata,timeslice)
  +            call swap_y_z_dimensions(inputdata)
               domain%nrain=inputdata
               deallocate(inputdata)
           endif
           if (io_variable_is_present(restart_file,"ni")) then
               call io_read3d(restart_file,"ni",inputdata,timeslice)
  +            call swap_y_z_dimensions(inputdata)
               domain%nice=inputdata
               deallocate(inputdata)
           endif
           call io_read3d(restart_file,"p",inputdata,timeslice)
  +        call swap_y_z_dimensions(inputdata)
           domain%p=inputdata
           deallocate(inputdata)
           call io_read3d(restart_file,"th",inputdata,timeslice)
  +        call swap_y_z_dimensions(inputdata)
           domain%th=inputdata
           deallocate(inputdata)
           call io_read3d(restart_file,"rho",inputdata,timeslice)
  +        call swap_y_z_dimensions(inputdata)
           domain%rho=inputdata
           deallocate(inputdata)
           
  
positive.code:         if (io_variable_is_present(restart_file,"graupel")) then
              call io_read2d(restart_file,"graupel",inputdata_2d,timeslice)
              domain%graupel=inputdata_2d
              deallocate(inputdata_2d)
          endif
          
  
negative.code:         do i=1,3
              if (size(data1,i).ne.size(data2,i)) then
                  write(*,*) "Restart file 3D dimensions don't match domain"
                  write(*,*) shape(data1)
                  write(*,*) shape(data2)
                  stop
              endif
          enddo
      end subroutine check_shapes_3d
      
      !>------------------------------------------------------------
      !! Swap the last two dimensions of an array
      !!
      !! Call reshape after finding the nx,ny,nz values
      !!
      !! @param data     3D array to be reshaped
      !!
      !!------------------------------------------------------------
      subroutine swap_y_z_dimensions(data)
          implicit none
          real,dimension(:,:,:),intent(inout),allocatable :: data
          real,dimension(:,:,:), allocatable :: temporary_data
          integer :: nx,ny,nz
          
          nx=size(data,1)
          ny=size(data,2)
          nz=size(data,3)
          allocate(temporary_data(nx,nz,ny))
          temporary_data = reshape(data, [nx,nz,ny], order=[1,3,2])
  
          deallocate(data)
          allocate(data(nx,nz,ny))
          data=temporary_data
          
      end subroutine swap_y_z_dimensions
      
      !>------------------------------------------------------------
      !!  Load restart file
  