anchor.diff:          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
  
positive.code:         write(*,*) "Reading atmospheric restart data"
          write(*,*) "   timestep:",trim(str(timeslice))," from file:",trim(restart_file)
          call io_read3d(restart_file,"qv",inputdata,timeslice)
          call check_shapes_3d(inputdata,domain%qv)
          domain%qv=inputdata
          deallocate(inputdata)
          call io_read3d(restart_file,"qc",inputdata,timeslice)
          domain%cloud=inputdata
          deallocate(inputdata)
          call io_read3d(restart_file,"qr",inputdata,timeslice)
          domain%qrain=inputdata
          deallocate(inputdata)
          call io_read3d(restart_file,"qi",inputdata,timeslice)
          domain%ice=inputdata
          deallocate(inputdata)
          call io_read3d(restart_file,"qs",inputdata,timeslice)
          domain%qsnow=inputdata
          deallocate(inputdata)
          if (io_variable_is_present(restart_file,"qg")) then
              call io_read3d(restart_file,"qg",inputdata,timeslice)
              domain%qgrau=inputdata
              deallocate(inputdata)
          endif
          if (io_variable_is_present(restart_file,"nr")) then
              call io_read3d(restart_file,"nr",inputdata,timeslice)
              domain%nrain=inputdata
              deallocate(inputdata)
          endif
          if (io_variable_is_present(restart_file,"ni")) then
              call io_read3d(restart_file,"ni",inputdata,timeslice)
              domain%nice=inputdata
              deallocate(inputdata)
          endif
          call io_read3d(restart_file,"p",inputdata,timeslice)
          domain%p=inputdata
          deallocate(inputdata)
          call io_read3d(restart_file,"th",inputdata,timeslice)
          domain%th=inputdata
          deallocate(inputdata)
          call io_read3d(restart_file,"rho",inputdata,timeslice)
          domain%rho=inputdata
          deallocate(inputdata)
          
  
negative.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
          if (io_variable_is_present(restart_file,"crain")) then
              call io_read2d(restart_file,"crain",inputdata_2d,timeslice)
              domain%crain=inputdata_2d
              deallocate(inputdata_2d)
          endif
          
  