anchor.diff:              ! and write constant (in time) variables
               call check( nf90_put_var(ncid, lat_id,    domain%lat), trim(filename)//":Latitude" )
               call check( nf90_put_var(ncid, lon_id,    domain%lon), trim(filename)//":Longitude" )
  -            call check( nf90_put_var(ncid, varid(20), reshape(domain%z, output_shape, order=zlast)) , trim(filename)//":Z")
  +            call check( nf90_put_var(ncid, lev_id,    domain%z(1,:,1)-domain%terrain(1,1)), trim(filename)//":Level" )
  +
  +            if (.not.surface_io_only) then
  +                call check( nf90_put_var(ncid, varid(20), reshape(domain%z, output_shape, order=zlast)) , trim(filename)//":Z")
  +            else
  +                call check( nf90_put_var(ncid, varid(20), domain%z(:,1,:)) , trim(filename)//":Z")
  +            endif
           endif
   
           ! write the actual data
  
positive.code:         call check( nf90_put_att(ncid,NF90_GLOBAL,"contact","Ethan Gutmann : gutmann@ucar.edu"), trim(err))
          call check( nf90_put_att(ncid,NF90_GLOBAL,"git",VERSION), trim(err))
  
          ! call check( nf90_put_att(ncid,NF90_GLOBAL,"bucket_size",kPRECIP_BUCKET_SIZE), trim(err))
          call check( nf90_put_att(ncid,NF90_GLOBAL,"dx",options%dx), trim(err))
          call check( nf90_put_att(ncid,NF90_GLOBAL,"wind_smoothing",options%smooth_wind_distance), trim(err))
  
  
negative.code:     integer :: dimtwo(2)
      integer :: dimtwo_time(3)
      !> variable IDs
      integer :: lat_id,lon_id,time_id, lev_id
      !> array to store ALL var ids
      integer :: varid(nvars)
  
  