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), domain%z) , trim(filename)//":Z")
  +            call check( nf90_put_var(ncid, varid(20), reshape(domain%z, output_shape, order=zlast)) , trim(filename)//":Z")
           endif
           
           ! write the actual data
           call check( nf90_put_var(ncid, time_id,   domain%model_time/86400.0+50000, start_scalar ), trim(filename)//":Time" )
           
           if (.not.surface_io_only) then
  -            call check( nf90_put_var(ncid, varid(1),  domain%qv,    start_three_D),    trim(filename)//":qv" )
  -            call check( nf90_put_var(ncid, varid(2),  domain%cloud, start_three_D),    trim(filename)//":cloud" )
  -            call check( nf90_put_var(ncid, varid(3),  domain%ice,   start_three_D),    trim(filename)//":ice" )
  -            call check( nf90_put_var(ncid, varid(4),  domain%qrain, start_three_D),    trim(filename)//":qrain" )
  -            call check( nf90_put_var(ncid, varid(5),  domain%qsnow, start_three_D),    trim(filename)//":qsnow" )
  +            call check( nf90_put_var(ncid, varid(1),  reshape(domain%qv,    output_shape, order=zlast), start_three_D),    trim(filename)//":qv" )
  +            call check( nf90_put_var(ncid, varid(2),  reshape(domain%cloud, output_shape, order=zlast), start_three_D),    trim(filename)//":cloud" )
  +            call check( nf90_put_var(ncid, varid(3),  reshape(domain%ice,   output_shape, order=zlast), start_three_D),    trim(filename)//":ice" )
  +            call check( nf90_put_var(ncid, varid(4),  reshape(domain%qrain, output_shape, order=zlast), start_three_D),    trim(filename)//":qrain" )
  +            call check( nf90_put_var(ncid, varid(5),  reshape(domain%qsnow, output_shape, order=zlast), start_three_D),    trim(filename)//":qsnow" )
               ! these should only be output for thompson microphysics
               if (options%physics%microphysics==kMP_THOMPSON) then
  -                call check( nf90_put_var(ncid, varid(6),  domain%qgrau, start_three_D),trim(filename)//":qgraupel" )
  -                call check( nf90_put_var(ncid, varid(7),  domain%nrain, start_three_D),trim(filename)//":nrain" )
  -                call check( nf90_put_var(ncid, varid(8),  domain%nice,  start_three_D),trim(filename)//":nice" )
  -            endif
  -            call check( nf90_put_var(ncid, varid(9),  domain%u,     start_three_D),    trim(filename)//":u" )
  -            call check( nf90_put_var(ncid, varid(10), domain%v,     start_three_D),    trim(filename)//":v" )
  -            call check( nf90_put_var(ncid, varid(11), domain%w_real,     start_three_D),    trim(filename)//":w" )
  -            call check( nf90_put_var(ncid, varid(12), domain%p,     start_three_D),    trim(filename)//":p" )
  -            call check( nf90_put_var(ncid, varid(13), domain%th,    start_three_D),    trim(filename)//":th" )
  -            call check( nf90_put_var(ncid, varid(21), domain%rho,   start_three_D),    trim(filename)//":rho" )
  +                call check( nf90_put_var(ncid, varid(6),  reshape(domain%qgrau, output_shape, order=zlast), start_three_D),trim(filename)//":qgraupel" )
  +                call check( nf90_put_var(ncid, varid(7),  reshape(domain%nrain, output_shape, order=zlast), start_three_D),trim(filename)//":nrain" )
  +                call check( nf90_put_var(ncid, varid(8),  reshape(domain%nice,  output_shape, order=zlast), start_three_D),trim(filename)//":nice" )
  +            endif
  +            output_shape(1)=output_shape(1)+1
  +            call check( nf90_put_var(ncid, varid(9),  reshape(domain%u,     output_shape, order=zlast), start_three_D),    trim(filename)//":u" )
  +            output_shape(1)=output_shape(1)-1
  +            output_shape(2)=output_shape(2)+1
  +            call check( nf90_put_var(ncid, varid(10), reshape(domain%v,     output_shape, order=zlast), start_three_D),    trim(filename)//":v" )
  +            output_shape(2)=output_shape(2)-1
  +            call check( nf90_put_var(ncid, varid(11), reshape(domain%w_real,output_shape, order=zlast), start_three_D),    trim(filename)//":w" )
  +            call check( nf90_put_var(ncid, varid(12), reshape(domain%p,     output_shape, order=zlast), start_three_D),    trim(filename)//":p" )
  +            call check( nf90_put_var(ncid, varid(13), reshape(domain%th,    output_shape, order=zlast), start_three_D),    trim(filename)//":th" )
  +            call check( nf90_put_var(ncid, varid(21), reshape(domain%rho,   output_shape, order=zlast), start_three_D),    trim(filename)//":rho" )
               if (options%physics%windtype==kWIND_LINEAR) then
  -                call check( nf90_put_var(ncid, varid(33), domain%nsquared, start_three_D), trim(filename)//":nsquared" )
  +                call check( nf90_put_var(ncid, varid(33), reshape(domain%nsquared, output_shape, order=zlast), start_three_D), trim(filename)//":nsquared" )
               endif
           else
               call check( nf90_put_var(ncid, varid(1),  domain%qv(:,1,:), start_two_D),  trim(filename)//":qv" )
  
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:             varid(9)=temp_id
              dimids(1)=x_id
              
              dimids(2)=yv_id
              call check( nf90_def_var(ncid, "v", NF90_REAL, dimids, temp_id), trim(err)//"v" )
              call check( nf90_put_att(ncid,temp_id,"standard_name","grid_northward_wind"))
              call check( nf90_put_att(ncid,temp_id,"long_name","Grid relative northward wind"))
              call check( nf90_put_att(ncid,temp_id,"units","m s-1"))
              varid(10)=temp_id
              dimids(2)=y_id
              
              call check( nf90_def_var(ncid, "w", NF90_REAL, dimids, temp_id), trim(err)//"w" )
              call check( nf90_put_att(ncid,temp_id,"standard_name","upward_air_velocity"))
  