anchor.diff:          call swap_y_z_dimensions(inputdata)
           domain%qrain=inputdata
           deallocate(inputdata)
  -        call io_read3d(restart_file,"qi",inputdata,timeslice)
  +        
  +        call io_read3d(restart_file,"p",inputdata,timeslice)
           call swap_y_z_dimensions(inputdata)
  -        domain%ice=inputdata
  +        domain%p=inputdata
           deallocate(inputdata)
  -        call io_read3d(restart_file,"qs",inputdata,timeslice)
  +        
  +        call io_read3d(restart_file,"th",inputdata,timeslice)
           call swap_y_z_dimensions(inputdata)
  -        domain%qsnow=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)
  +        
  +
  +        ! The variables below are optional depending on the physics options used so test for their presence
  +        if (io_variable_is_present(restart_file,"qi")) then
  +            call io_read3d(restart_file,"qi",inputdata,timeslice)
  +            call swap_y_z_dimensions(inputdata)
  +            domain%ice=inputdata
  +            deallocate(inputdata)
  +        endif
  +        if (io_variable_is_present(restart_file,"qs")) then
  +            call io_read3d(restart_file,"qs",inputdata,timeslice)
  +            call swap_y_z_dimensions(inputdata)
  +            domain%qsnow=inputdata
  +            deallocate(inputdata)
  +        endif
           if (io_variable_is_present(restart_file,"qg")) then
               call io_read3d(restart_file,"qg",inputdata,timeslice)
               call swap_y_z_dimensions(inputdata)
  
positive.code:             call io_read2d(restart_file,"rsds",inputdata_2d,timeslice)
              domain%swdown=inputdata_2d
              deallocate(inputdata_2d)
              
              call io_read2d(restart_file,"rlds",inputdata_2d,timeslice)
              domain%lwdown=inputdata_2d
              deallocate(inputdata_2d)
          endif
          
          
      end subroutine load_restart_file
      
  
negative.code:         endif
          
          if (io_variable_is_present(restart_file,"soil_t")) then
              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)
          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
          
  