anchor.diff:              call check( nf90_put_att(ncid,temp_id,"units","W m-2"))
               varid(25)=temp_id
       
  -            dimids(2)=soil_id
  +            dimids(3)=soil_id
               call check( nf90_def_var(ncid, "soil_w", NF90_REAL, dimids, temp_id), trim(err)//"soil_w" )
               call check( nf90_put_att(ncid,temp_id,"standard_name","soil_moisture_content"))
               call check( nf90_put_att(ncid,temp_id,"long_name","Column Soil Moisture"))
  
positive.code:             varid(29)=temp_id
              call check( nf90_inq_varid(ncid, "hfgs",  temp_id), trim(err)//"hfgs" )
              varid(25)=temp_id
              dimids(2)=soil_id
              call check( nf90_inq_varid(ncid, "soil_w", temp_id), trim(err)//"soil_w" )
              varid(26)=temp_id
              call check( nf90_inq_varid(ncid, "soil_t", temp_id), trim(err)//"soil_t" )
  
negative.code:         surface_io_only = options%surface_io_only
          
          ! these are module level variables
          nx = size(domain%qv,1)
          nz = size(domain%qv,2)
          ny = size(domain%qv,3)
          nsoil = size(domain%soil_t,2)
          output_shape = [nx,ny,nz]
          zlast = [1,3,2]
          
          if (.not.allocated(last_rain)) then
              output_rain_rate=.False.
              allocate(last_rain(nx,ny))
  