anchor.diff:          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
           
  
positive.code:         do i=1,3
              if (size(data1,i).ne.size(data2,i)) then
                  write(*,*) "Restart file 3D dimensions don't match domain"
                  stop
              endif
          enddo
      end subroutine check_shapes_3d
      
      !>------------------------------------------------------------
      !!  Load restart file
  
negative.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 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)
          
  