anchor.diff:          endif
           
           if (io_variable_is_present(restart_file,"soil_t")) then
  -            write(*,*) "Reading land surface restart data"
               call io_read3d(restart_file,"soil_t",inputdata,timeslice)
               call swap_y_z_dimensions(inputdata)
               call check_shapes_3d(inputdata,domain%soil_t)
               domain%soil_t=inputdata
               deallocate(inputdata)
  +        endif
  +        if (io_variable_is_present(restart_file, "soil_w")) then
               call io_read3d(restart_file,"soil_w",inputdata,timeslice)
               call swap_y_z_dimensions(inputdata)
               domain%soil_vwc=inputdata
               deallocate(inputdata)
  -        
  -            call io_read2d(restart_file,"ts",inputdata_2d,timeslice)
  -            domain%skin_t=inputdata_2d
  -            deallocate(inputdata_2d)
  +        endif
  +        
  +        if (io_variable_is_present(restart_file, "hfgs")) then
               call io_read2d(restart_file,"hfgs",inputdata_2d,timeslice)
               domain%ground_heat=inputdata_2d
               deallocate(inputdata_2d)
  +        endif
  +
  +        if (io_variable_is_present(restart_file, "snw")) then
               call io_read2d(restart_file,"snw",inputdata_2d,timeslice)
               domain%snow_swe=inputdata_2d
               deallocate(inputdata_2d)
  +        endif
  +        
  +        if (io_variable_is_present(restart_file, "canwat")) then
               call io_read2d(restart_file,"canwat",inputdata_2d,timeslice)
               domain%canopy_water=inputdata_2d
  +            deallocate(inputdata_2d)
  +        endif
  +
  +        if (io_variable_is_present(restart_file, "ts")) then
  +            call io_read2d(restart_file,"ts",inputdata_2d,timeslice)
  +            domain%skin_t=inputdata_2d
  +            domain%sst = inputdata_2d
               deallocate(inputdata_2d)
           endif
           
  
positive.code:             domain%nsnow=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)
          
          call io_read2d(restart_file,"rain",inputdata_2d,timeslice)
          domain%rain=inputdata_2d
          deallocate(inputdata_2d)
          call io_read2d(restart_file,"snow",inputdata_2d,timeslice)
          domain%snow=inputdata_2d
          deallocate(inputdata_2d)
          if (io_variable_is_present(restart_file,"graupel")) then
              call io_read2d(restart_file,"graupel",inputdata_2d,timeslice)
              domain%graupel=inputdata_2d
  
negative.code:         
          write(*,*) "Reading atmospheric restart data"
          write(*,*) "   timestep:",trim(str(timeslice))," from file:",trim(restart_file)
          
          ! The first variables here are required. If they do not exist, ICAR *should* exit with a netCDF error
          call io_read3d(restart_file,"qv",inputdata,timeslice)
          call swap_y_z_dimensions(inputdata)
          call check_shapes_3d(inputdata,domain%qv)
  