Write dedicated enthalpy/energy budget diagnostics.
This diagnostics-only routine separates the terms that are currently available from the missing conductive-boundary-flux closure term.
Sign convention: boundary fluxes are positive outward from the domain.
Reported balance without conduction: d/dt integral(rho h dV) + net_boundary_advective_rho_h_flux - qrad_integral - rho_species_enthalpy_diffusion_integral
If the full enthalpy equation is written with an outward conductive flux on the left-hand side, the conductive boundary flux required for closure is: required_conductive_out = -balance_without_conduction
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(mesh_t), | intent(in) | :: | mesh | |||
| type(flow_mpi_t), | intent(in) | :: | flow | |||
| type(case_params_t), | intent(in) | :: | params | |||
| type(flow_fields_t), | intent(in) | :: | fields | |||
| type(energy_fields_t), | intent(in) | :: | energy | |||
| type(transport_properties_t), | intent(in) | :: | transport | |||
| integer, | intent(in) | :: | step | |||
| real(kind=rk), | intent(in) | :: | time |
subroutine write_enthalpy_energy_budget_diagnostics(mesh, flow, params, fields, energy, transport, step, time) type(mesh_t), intent(in) :: mesh type(flow_mpi_t), intent(in) :: flow type(case_params_t), intent(in) :: params type(flow_fields_t), intent(in) :: fields type(energy_fields_t), intent(in) :: energy type(transport_properties_t), intent(in) :: transport integer, intent(in) :: step real(rk), intent(in) :: time integer :: c, f, owner, nb, ierr, unit_id logical, save :: enthalpy_budget_initialized = .false. logical, save :: previous_initialized = .false. character(len=1024) :: filename real(rk), save :: previous_time = 0.0_rk real(rk), save :: previous_rho_h_integral = 0.0_rk real(rk) :: vol, rho_c, h_c, T_c, cp_c, lambda_c, mflux real(rk) :: dt_since_previous, delta_rho_h, rho_h_rate real(rk) :: balance_without_conduction, balance_with_reported_conduction real(rk) :: required_conductive_boundary_flux_out real(rk) :: relative_balance_without_conduction real(rk) :: denom real(rk) :: local_volume, volume_total real(rk) :: local_rho_h_integral, rho_h_integral real(rk) :: local_boundary_rho_h_advective_flux, net_boundary_rho_h_advective_flux real(rk) :: local_boundary_rho_h_advective_outflow, boundary_rho_h_advective_outflow real(rk) :: local_boundary_rho_h_advective_inflow, boundary_rho_h_advective_inflow real(rk) :: local_qrad_integral, qrad_integral real(rk) :: local_rho_species_hdiff_integral, rho_species_hdiff_integral real(rk) :: local_h_min, h_min real(rk) :: local_h_max, h_max real(rk) :: local_T_min, T_min real(rk) :: local_T_max, T_max real(rk) :: local_rho_min, rho_min real(rk) :: local_rho_max, rho_max real(rk) :: local_cp_min, cp_min real(rk) :: local_cp_max, cp_max real(rk) :: local_lambda_min, lambda_min real(rk) :: local_lambda_max, lambda_max real(rk) :: local_h_int, h_mean real(rk) :: local_T_int, T_mean real(rk) :: local_rho_int, rho_mean real(rk) :: local_cp_int, cp_mean real(rk) :: local_lambda_int, lambda_mean real(rk) :: local_rho_h_l2_num, rho_h_l2 real(rk) :: sum_send(24), sum_recv(24) real(rk) :: min_send(5), min_recv(5) real(rk) :: max_send(5), max_recv(5) real(rk) :: reported_conductive_boundary_flux_out integer :: conductive_boundary_flux_available real(rk), save :: previous_cumulative_advective_boundary_rho_h_flux = 0.0_rk real(rk), save :: previous_cumulative_conductive_boundary_flux = 0.0_rk real(rk), save :: previous_cumulative_qrad_integral = 0.0_rk real(rk), save :: previous_cumulative_rho_species_hdiff_integral = 0.0_rk real(rk) :: cumulative_advective_boundary_rho_h_flux real(rk) :: cumulative_conductive_boundary_flux real(rk) :: cumulative_qrad_integral_total real(rk) :: cumulative_rho_species_hdiff_integral_total integer :: cumulative_energy_budget_available real(rk) :: cumulative_advective_boundary_rho_h_flux_average real(rk) :: cumulative_conductive_boundary_flux_average real(rk) :: cumulative_qrad_integral_average real(rk) :: cumulative_rho_species_hdiff_integral_average real(rk) :: cumulative_operator_balance_defect real(rk) :: relative_cumulative_operator_balance_defect real(rk), save :: previous_operator_consistent_rho_h_integral = 0.0_rk real(rk) :: operator_consistent_rho_h_integral real(rk) :: operator_consistent_delta_rho_h real(rk) :: operator_consistent_rho_h_rate real(rk) :: operator_consistent_cumulative_balance_defect real(rk) :: relative_operator_consistent_cumulative_balance_defect real(rk), save :: previous_cumulative_energy_update_delta_integral = 0.0_rk real(rk), save :: previous_cumulative_energy_update_rhs_integral = 0.0_rk real(rk) :: cumulative_energy_update_delta_rate_average real(rk) :: cumulative_energy_update_rhs_rate_average real(rk) :: output_state_density_reconciliation_rate real(rk) :: operator_consistent_density_reconciliation_rate real(rk) :: output_state_budget_defect_after_density_reconciliation real(rk) :: operator_consistent_budget_defect_after_density_reconciliation real(rk) :: rel_output_recon_defect real(rk) :: rel_operator_recon_defect if (.not. params%write_diagnostics) return if (.not. params%enable_energy) return if (.not. allocated(energy%h)) return if (.not. allocated(transport%rho)) return local_volume = 0.0_rk local_rho_h_integral = 0.0_rk local_boundary_rho_h_advective_flux = 0.0_rk local_boundary_rho_h_advective_outflow = 0.0_rk local_boundary_rho_h_advective_inflow = 0.0_rk local_qrad_integral = 0.0_rk local_rho_species_hdiff_integral = 0.0_rk local_h_min = huge(1.0_rk) local_h_max = -huge(1.0_rk) local_T_min = huge(1.0_rk) local_T_max = -huge(1.0_rk) local_rho_min = huge(1.0_rk) local_rho_max = -huge(1.0_rk) local_cp_min = huge(1.0_rk) local_cp_max = -huge(1.0_rk) local_lambda_min = huge(1.0_rk) local_lambda_max = -huge(1.0_rk) local_h_int = 0.0_rk local_T_int = 0.0_rk local_rho_int = 0.0_rk local_cp_int = 0.0_rk local_lambda_int = 0.0_rk local_rho_h_l2_num = 0.0_rk do c = 1, mesh%ncells if (.not. flow%owned(c)) cycle vol = mesh%cells(c)%volume if (vol <= 0.0_rk) cycle rho_c = transport%rho(c) h_c = energy%h(c) T_c = 0.0_rk if (allocated(energy%T)) then if (size(energy%T) >= c) T_c = energy%T(c) end if cp_c = 0.0_rk if (allocated(energy%cp)) then if (size(energy%cp) >= c) cp_c = energy%cp(c) end if lambda_c = 0.0_rk if (allocated(energy%lambda)) then if (size(energy%lambda) >= c) lambda_c = energy%lambda(c) end if local_volume = local_volume + vol local_rho_h_integral = local_rho_h_integral + rho_c * h_c * vol local_rho_h_l2_num = local_rho_h_l2_num + (rho_c * h_c) * (rho_c * h_c) * vol local_h_min = min(local_h_min, h_c) local_h_max = max(local_h_max, h_c) local_T_min = min(local_T_min, T_c) local_T_max = max(local_T_max, T_c) local_rho_min = min(local_rho_min, rho_c) local_rho_max = max(local_rho_max, rho_c) local_cp_min = min(local_cp_min, cp_c) local_cp_max = max(local_cp_max, cp_c) local_lambda_min = min(local_lambda_min, lambda_c) local_lambda_max = max(local_lambda_max, lambda_c) local_h_int = local_h_int + h_c * vol local_T_int = local_T_int + T_c * vol local_rho_int = local_rho_int + rho_c * vol local_cp_int = local_cp_int + cp_c * vol local_lambda_int = local_lambda_int + lambda_c * vol if (allocated(energy%qrad)) then if (size(energy%qrad) >= c) local_qrad_integral = local_qrad_integral + energy%qrad(c) * vol end if if (allocated(energy%species_enthalpy_diffusion)) then if (size(energy%species_enthalpy_diffusion) >= c) then local_rho_species_hdiff_integral = local_rho_species_hdiff_integral + & rho_c * energy%species_enthalpy_diffusion(c) * vol end if end if end do if (allocated(fields%mass_flux)) then do f = 1, mesh%nfaces owner = mesh%faces(f)%owner nb = mesh%faces(f)%neighbor if (nb <= 0) nb = mesh%faces(f)%periodic_neighbor if (nb <= 0 .and. owner >= 1 .and. owner <= mesh%ncells) then if (.not. flow%owned(owner)) cycle if (size(energy%h) < owner) cycle mflux = fields%mass_flux(f) h_c = energy%h(owner) local_boundary_rho_h_advective_flux = local_boundary_rho_h_advective_flux + mflux * h_c if (mflux >= 0.0_rk) then local_boundary_rho_h_advective_outflow = local_boundary_rho_h_advective_outflow + mflux * h_c else local_boundary_rho_h_advective_inflow = local_boundary_rho_h_advective_inflow + mflux * h_c end if end if end do end if sum_send = 0.0_rk sum_send(1) = local_volume sum_send(2) = local_rho_h_integral sum_send(3) = local_boundary_rho_h_advective_flux sum_send(4) = local_boundary_rho_h_advective_outflow sum_send(5) = local_boundary_rho_h_advective_inflow sum_send(6) = local_qrad_integral sum_send(7) = local_rho_species_hdiff_integral sum_send(8) = local_h_int sum_send(9) = local_T_int sum_send(10) = local_rho_int sum_send(11) = local_cp_int sum_send(12) = local_lambda_int sum_send(13) = local_rho_h_l2_num sum_send(14) = energy%last_conductive_boundary_flux_out sum_send(15) = real(energy%last_conductive_boundary_flux_available, rk) sum_send(16) = energy%cumulative_boundary_rho_h_advective_flux_out sum_send(17) = energy%cumulative_boundary_rho_h_conductive_flux_out sum_send(18) = energy%cumulative_qrad_integral sum_send(19) = energy%cumulative_rho_species_hdiff_integral sum_send(20) = real(energy%cumulative_energy_budget_available, rk) call MPI_Allreduce(sum_send, sum_recv, 24, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing enthalpy budget sums') min_send(1) = local_h_min min_send(2) = local_T_min min_send(3) = local_rho_min min_send(4) = local_cp_min min_send(5) = local_lambda_min call MPI_Allreduce(min_send, min_recv, 5, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing enthalpy budget minima') max_send(1) = local_h_max max_send(2) = local_T_max max_send(3) = local_rho_max max_send(4) = local_cp_max max_send(5) = local_lambda_max call MPI_Allreduce(max_send, max_recv, 5, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing enthalpy budget maxima') volume_total = sum_recv(1) rho_h_integral = sum_recv(2) net_boundary_rho_h_advective_flux = sum_recv(3) boundary_rho_h_advective_outflow = sum_recv(4) boundary_rho_h_advective_inflow = sum_recv(5) qrad_integral = sum_recv(6) rho_species_hdiff_integral = sum_recv(7) if (volume_total > 0.0_rk) then h_mean = sum_recv(8) / volume_total T_mean = sum_recv(9) / volume_total rho_mean = sum_recv(10) / volume_total cp_mean = sum_recv(11) / volume_total lambda_mean = sum_recv(12) / volume_total rho_h_l2 = sqrt(max(0.0_rk, sum_recv(13) / volume_total)) else h_mean = 0.0_rk T_mean = 0.0_rk rho_mean = 0.0_rk cp_mean = 0.0_rk lambda_mean = 0.0_rk rho_h_l2 = 0.0_rk end if h_min = min_recv(1) T_min = min_recv(2) rho_min = min_recv(3) cp_min = min_recv(4) lambda_min = min_recv(5) h_max = max_recv(1) T_max = max_recv(2) rho_max = max_recv(3) cp_max = max_recv(4) lambda_max = max_recv(5) if (previous_initialized) then dt_since_previous = time - previous_time else dt_since_previous = 0.0_rk end if if (previous_initialized .and. dt_since_previous > tiny(1.0_rk)) then delta_rho_h = rho_h_integral - previous_rho_h_integral rho_h_rate = delta_rho_h / dt_since_previous else delta_rho_h = 0.0_rk rho_h_rate = 0.0_rk end if reported_conductive_boundary_flux_out = sum_recv(14) if (sum_recv(15) > 0.5_rk) then conductive_boundary_flux_available = 1 else conductive_boundary_flux_available = 0 end if cumulative_advective_boundary_rho_h_flux = sum_recv(16) cumulative_conductive_boundary_flux = sum_recv(17) cumulative_qrad_integral_total = sum_recv(18) cumulative_rho_species_hdiff_integral_total = sum_recv(19) if (sum_recv(20) > 0.5_rk) then cumulative_energy_budget_available = 1 else cumulative_energy_budget_available = 0 end if if (previous_initialized .and. dt_since_previous > tiny(1.0_rk) .and. & cumulative_energy_budget_available == 1) then cumulative_advective_boundary_rho_h_flux_average = & (cumulative_advective_boundary_rho_h_flux - previous_cumulative_advective_boundary_rho_h_flux) / & dt_since_previous cumulative_conductive_boundary_flux_average = & (cumulative_conductive_boundary_flux - previous_cumulative_conductive_boundary_flux) / & dt_since_previous cumulative_qrad_integral_average = & (cumulative_qrad_integral_total - previous_cumulative_qrad_integral) / dt_since_previous cumulative_rho_species_hdiff_integral_average = & (cumulative_rho_species_hdiff_integral_total - previous_cumulative_rho_species_hdiff_integral) / & dt_since_previous cumulative_operator_balance_defect = rho_h_rate + & cumulative_advective_boundary_rho_h_flux_average + & cumulative_conductive_boundary_flux_average - & cumulative_qrad_integral_average - & cumulative_rho_species_hdiff_integral_average denom = max(abs(rho_h_rate), abs(cumulative_advective_boundary_rho_h_flux_average), & abs(cumulative_conductive_boundary_flux_average), abs(cumulative_qrad_integral_average), & abs(cumulative_rho_species_hdiff_integral_average), tiny(1.0_rk)) relative_cumulative_operator_balance_defect = abs(cumulative_operator_balance_defect) / denom else cumulative_advective_boundary_rho_h_flux_average = 0.0_rk cumulative_conductive_boundary_flux_average = 0.0_rk cumulative_qrad_integral_average = 0.0_rk cumulative_rho_species_hdiff_integral_average = 0.0_rk cumulative_operator_balance_defect = 0.0_rk relative_cumulative_operator_balance_defect = 0.0_rk end if operator_consistent_rho_h_integral = energy%last_operator_consistent_rho_h_integral if (previous_initialized .and. dt_since_previous > tiny(1.0_rk) .and. & cumulative_energy_budget_available == 1) then operator_consistent_delta_rho_h = operator_consistent_rho_h_integral - & previous_operator_consistent_rho_h_integral operator_consistent_rho_h_rate = operator_consistent_delta_rho_h / dt_since_previous operator_consistent_cumulative_balance_defect = operator_consistent_rho_h_rate + & cumulative_advective_boundary_rho_h_flux_average + & cumulative_conductive_boundary_flux_average - & cumulative_qrad_integral_average - & cumulative_rho_species_hdiff_integral_average denom = max(abs(operator_consistent_rho_h_rate), & abs(cumulative_advective_boundary_rho_h_flux_average), & abs(cumulative_conductive_boundary_flux_average), & abs(cumulative_qrad_integral_average), & abs(cumulative_rho_species_hdiff_integral_average), tiny(1.0_rk)) relative_operator_consistent_cumulative_balance_defect = & abs(operator_consistent_cumulative_balance_defect) / denom else operator_consistent_delta_rho_h = 0.0_rk operator_consistent_rho_h_rate = 0.0_rk operator_consistent_cumulative_balance_defect = 0.0_rk relative_operator_consistent_cumulative_balance_defect = 0.0_rk end if if (previous_initialized .and. dt_since_previous > tiny(1.0_rk) .and. & cumulative_energy_budget_available == 1) then cumulative_energy_update_delta_rate_average = & (energy%cumulative_energy_update_delta_integral - previous_cumulative_energy_update_delta_integral) / & dt_since_previous cumulative_energy_update_rhs_rate_average = & (energy%cumulative_energy_update_rhs_integral - previous_cumulative_energy_update_rhs_integral) / & dt_since_previous output_state_density_reconciliation_rate = rho_h_rate - cumulative_energy_update_delta_rate_average operator_consistent_density_reconciliation_rate = operator_consistent_rho_h_rate - & cumulative_energy_update_delta_rate_average output_state_budget_defect_after_density_reconciliation = & cumulative_operator_balance_defect - output_state_density_reconciliation_rate operator_consistent_budget_defect_after_density_reconciliation = & operator_consistent_cumulative_balance_defect - operator_consistent_density_reconciliation_rate denom = max(abs(cumulative_energy_update_delta_rate_average), & abs(cumulative_energy_update_rhs_rate_average), & abs(cumulative_advective_boundary_rho_h_flux_average), & abs(cumulative_conductive_boundary_flux_average), & abs(cumulative_qrad_integral_average), & abs(cumulative_rho_species_hdiff_integral_average), tiny(1.0_rk)) rel_output_recon_defect = & abs(output_state_budget_defect_after_density_reconciliation) / denom rel_operator_recon_defect = & abs(operator_consistent_budget_defect_after_density_reconciliation) / denom else cumulative_energy_update_delta_rate_average = 0.0_rk cumulative_energy_update_rhs_rate_average = 0.0_rk output_state_density_reconciliation_rate = 0.0_rk operator_consistent_density_reconciliation_rate = 0.0_rk output_state_budget_defect_after_density_reconciliation = 0.0_rk operator_consistent_budget_defect_after_density_reconciliation = 0.0_rk rel_output_recon_defect = 0.0_rk rel_operator_recon_defect = 0.0_rk end if balance_without_conduction = rho_h_rate + net_boundary_rho_h_advective_flux - & qrad_integral - rho_species_hdiff_integral required_conductive_boundary_flux_out = -balance_without_conduction balance_with_reported_conduction = balance_without_conduction + reported_conductive_boundary_flux_out denom = max(abs(rho_h_rate), abs(net_boundary_rho_h_advective_flux), abs(qrad_integral), & abs(rho_species_hdiff_integral), tiny(1.0_rk)) relative_balance_without_conduction = abs(balance_without_conduction) / denom if (flow%rank == 0) then call execute_command_line('mkdir -p ' // trim(params%output_dir) // '/diagnostics') filename = trim(params%output_dir) // '/diagnostics/enthalpy_energy_budget.csv' if (.not. enthalpy_budget_initialized) then open(newunit=unit_id, file=trim(filename), status='replace', action='write') enthalpy_budget_initialized = .true. write(unit_id,'(a)') 'step,time,volume_total,rho_h_integral,delta_time_since_previous,' // & 'delta_rho_h_since_previous,rho_h_rate_since_previous,' // & 'net_boundary_rho_h_advective_flux,boundary_rho_h_advective_outflow,' // & 'boundary_rho_h_advective_inflow,qrad_integral,' // & 'rho_species_enthalpy_diffusion_integral,reported_conductive_boundary_flux_out,' // & 'conductive_boundary_flux_available,required_conductive_boundary_flux_out,' // & 'balance_defect_without_conduction,balance_defect_with_reported_conduction,' // & 'relative_balance_defect_without_conduction,' // & 'cumulative_advective_boundary_rho_h_flux_average,' // & 'cumulative_conductive_boundary_flux_average,cumulative_qrad_integral_average,' // & 'cumulative_rho_species_enthalpy_diffusion_integral_average,' // & 'cumulative_energy_budget_available,cumulative_operator_balance_defect,' // & 'relative_cumulative_operator_balance_defect,' // & 'last_energy_update_delta_rate_integral,last_energy_update_rhs_integral,' // & 'last_energy_update_balance_defect,relative_last_energy_update_balance_defect,' // & 'operator_consistent_rho_h_integral,' // & 'operator_consistent_delta_rho_h_since_previous,' // & 'operator_consistent_rho_h_rate_since_previous,' // & 'operator_consistent_cumulative_balance_defect,' // & 'relative_operator_consistent_cumulative_balance_defect,' // & 'cumulative_energy_update_delta_rate_average,' // & 'cumulative_energy_update_rhs_rate_average,' // & 'output_state_density_reconciliation_rate,' // & 'operator_consistent_density_reconciliation_rate,' // & 'output_state_budget_defect_after_density_reconciliation,' // & 'operator_consistent_budget_defect_after_density_reconciliation,' // & 'rel_output_recon_defect,' // & 'rel_operator_recon_defect,' // & 'h_min,h_max,h_mean,T_min,T_max,T_mean,' // & 'rho_min,rho_max,rho_mean,cp_min,cp_max,cp_mean,lambda_min,lambda_max,lambda_mean,rho_h_l2' else open(newunit=unit_id, file=trim(filename), status='unknown', position='append', action='write') end if call write_csv_value_i(unit_id, step) call write_csv_value_r(unit_id, time) call write_csv_value_r(unit_id, volume_total) call write_csv_value_r(unit_id, rho_h_integral) call write_csv_value_r(unit_id, dt_since_previous) call write_csv_value_r(unit_id, delta_rho_h) call write_csv_value_r(unit_id, rho_h_rate) call write_csv_value_r(unit_id, net_boundary_rho_h_advective_flux) call write_csv_value_r(unit_id, boundary_rho_h_advective_outflow) call write_csv_value_r(unit_id, boundary_rho_h_advective_inflow) call write_csv_value_r(unit_id, qrad_integral) call write_csv_value_r(unit_id, rho_species_hdiff_integral) call write_csv_value_r(unit_id, reported_conductive_boundary_flux_out) write(unit_id,'(",",i0)', advance='no') conductive_boundary_flux_available call write_csv_value_r(unit_id, required_conductive_boundary_flux_out) call write_csv_value_r(unit_id, balance_without_conduction) call write_csv_value_r(unit_id, balance_with_reported_conduction) call write_csv_value_r(unit_id, relative_balance_without_conduction) call write_csv_value_r(unit_id, cumulative_advective_boundary_rho_h_flux_average) call write_csv_value_r(unit_id, cumulative_conductive_boundary_flux_average) call write_csv_value_r(unit_id, cumulative_qrad_integral_average) call write_csv_value_r(unit_id, cumulative_rho_species_hdiff_integral_average) write(unit_id,'(",",i0)', advance='no') cumulative_energy_budget_available call write_csv_value_r(unit_id, cumulative_operator_balance_defect) call write_csv_value_r(unit_id, relative_cumulative_operator_balance_defect) call write_csv_value_r(unit_id, energy%last_energy_update_delta_rate_integral) call write_csv_value_r(unit_id, energy%last_energy_update_rhs_integral) call write_csv_value_r(unit_id, energy%last_energy_update_balance_defect) call write_csv_value_r(unit_id, energy%relative_last_energy_update_balance_defect) call write_csv_value_r(unit_id, operator_consistent_rho_h_integral) call write_csv_value_r(unit_id, operator_consistent_delta_rho_h) call write_csv_value_r(unit_id, operator_consistent_rho_h_rate) call write_csv_value_r(unit_id, operator_consistent_cumulative_balance_defect) call write_csv_value_r(unit_id, relative_operator_consistent_cumulative_balance_defect) call write_csv_value_r(unit_id, cumulative_energy_update_delta_rate_average) call write_csv_value_r(unit_id, cumulative_energy_update_rhs_rate_average) call write_csv_value_r(unit_id, output_state_density_reconciliation_rate) call write_csv_value_r(unit_id, operator_consistent_density_reconciliation_rate) call write_csv_value_r(unit_id, output_state_budget_defect_after_density_reconciliation) call write_csv_value_r(unit_id, operator_consistent_budget_defect_after_density_reconciliation) call write_csv_value_r(unit_id, rel_output_recon_defect) call write_csv_value_r(unit_id, rel_operator_recon_defect) call write_csv_value_r(unit_id, h_min) call write_csv_value_r(unit_id, h_max) call write_csv_value_r(unit_id, h_mean) call write_csv_value_r(unit_id, T_min) call write_csv_value_r(unit_id, T_max) call write_csv_value_r(unit_id, T_mean) call write_csv_value_r(unit_id, rho_min) call write_csv_value_r(unit_id, rho_max) call write_csv_value_r(unit_id, rho_mean) call write_csv_value_r(unit_id, cp_min) call write_csv_value_r(unit_id, cp_max) call write_csv_value_r(unit_id, cp_mean) call write_csv_value_r(unit_id, lambda_min) call write_csv_value_r(unit_id, lambda_max) call write_csv_value_r(unit_id, lambda_mean) call write_csv_value_r(unit_id, rho_h_l2) write(unit_id,*) close(unit_id) end if previous_time = time previous_rho_h_integral = rho_h_integral previous_cumulative_advective_boundary_rho_h_flux = cumulative_advective_boundary_rho_h_flux previous_cumulative_conductive_boundary_flux = cumulative_conductive_boundary_flux previous_cumulative_qrad_integral = cumulative_qrad_integral_total previous_cumulative_rho_species_hdiff_integral = cumulative_rho_species_hdiff_integral_total previous_operator_consistent_rho_h_integral = operator_consistent_rho_h_integral previous_cumulative_energy_update_delta_integral = energy%cumulative_energy_update_delta_integral previous_cumulative_energy_update_rhs_integral = energy%cumulative_energy_update_rhs_integral previous_initialized = .true. end subroutine write_enthalpy_energy_budget_diagnostics