arg-reduced solve_rate_cool_g-cpp.C

#include <cstdio>

#include <vector>

#include "grackle.h"

#include "fortran_func_decls.h"

#include "utils-cpp.hpp"

#include "solve_rate_cool_g-cpp.h"

//_// TODO: ADD ANY OTHER INCLUDE DIRECTIVES

extern "C" {

void solve_rate_cool_g(

int* imetal, double* dt, double* utem, double* uxyz, double* urho, int* ierr,

chemistry_data* my_chemistry, chemistry_data_storage* my_rates,

code_units* my_units, grackle_field_data* my_fields,

photo_rate_storage* my_uvb_rates

)

{

// SOLVE MULTI-SPECIES RATE EQUATIONS AND RADIATIVE COOLING

// written by: Yu Zhang, Peter Anninos and Tom Abel

// date:

// modified1: January, 1996 by Greg Bryan; converted to KRONOS

// modified2: October, 1996 by GB; adapted to AMR

// modified3: May, 1999 by GB and Tom Abel, 3bodyH2, solver, HD

// modified4: June, 2005 by GB to solve rate & cool at same time

// modified5: April, 2009 by JHW to include radiative transfer

// modified6: September, 2009 by BDS to include cloudy cooling

// PURPOSE:

// Solve the multi-species rate and cool equations.

// INPUTS:

// icool - flag to update energy from radiative cooling

// in,jn,kn - dimensions of 3D fields

// d - total density field

// de - electron density field

// HI,HII - H density fields (neutral & ionized)

// HeI/II/III - He density fields

// DI/II - D density fields (neutral & ionized)

// HDI - neutral HD molecule density field

// HM - H- density field

// H2I - H_2 (molecular H) density field

// H2II - H_2+ density field

// metal - metal density field

// dust - dust density field

// kph* - photoionization fields

// gamma* - photoheating fields

// f_shield_custom - custom H2 shielding factor

// is,ie - start and end indices of active region (zero based)

// iexpand - comoving coordinates flag (0 = off, 1 = on)

// idim - dimensionality (rank) of problem

// ispecies - chemistry module (1 - H/He only, 2 - molecular H, 3 - D)

// imetal - flag if metal field is active (0 = no, 1 = yes)

// imcool - flag if there is metal cooling

// idust - flag for H2 formation on dust grains

// idustall - flag for dust (0 - none, 1 - heating/cooling + H2 form.)

// idustfield - flag if a dust density field is present

// iisrffield - flag if a field for the interstellar radiation field is present

// ih2co - flag to include H2 cooling (1 = on, 0 = off)

// ipiht - flag to include photoionization heating (1 = on, 0 = off)

// idustrec - flag to include dust recombination cooling (1 = on, -1 = off)

// iH2shield - flag for approximate self-shielding of H2 (Wolcott-Green+ 2011)

// iradshield - flag for approximate self-shielding of UV background

// avgsighi - spectrum averaged ionization crs for HI for use with shielding

// avgsighei - spectrum averaged ionization crs for HeI for use with shielding

// avgsigheii - spectrum averaged ionization crs for HeII for use with shielding

// iradtrans - flag to include radiative transfer (1 = on, 0 = off)

// iradcoupled - flag to indicate coupled radiative transfer

// iradstep - flag to indicate intermediate coupled radiative transfer timestep

// irt_honly - flag to indicate applying RT ionization and heating to HI only

// iH2shieldcustom - flag to indicate a custom H2 shielding factor is provided

// fh - Hydrogen mass fraction (typically 0.76)

// dtoh - Deuterium to H mass ratio

// z_solar - Solar metal mass fraction

// fgr - the local dust to gas ratio (by mass)

// dt - timestep to integrate over

// aye - expansion factor (in code units)

// utim - time units (i.e. code units to CGS conversion factor)

// uaye - expansion factor conversion factor (uaye = 1/(1+zinit))

// urho - density units

// uxyz - length units

// utem - temperature(-like) units

// temstart, temend - start and end of temperature range for rate table

// nratec - dimensions of chemical rate arrays (functions of temperature)

// dtemstart, dtemend - start and end of dust temperature range

// ndratec - extra dimension for H2 formation on dust rate (dust temperature)

// icmbTfloor - flag to include temperature floor from cmb

// iClHeat - flag to include cloudy heating

// priGridRank - rank of cloudy primordial cooling data grid

// priGridDim - array containing dimensions of cloudy primordial data

// priPar1, priPar2, priPar3 - arrays containing primordial grid parameter values

// priDataSize - total size of flattened 1D primordial cooling data array

// priCooling - primordial cooling data

// priHeating - primordial heating data

// priMMW - primordial mmw data

// metGridRank - rank of cloudy metal cooling data grid

// metGridDim - array containing dimensions of cloudy metal data

// metPar1, metPar2, metPar3 - arrays containing metal grid parameter values

// metDataSize - total size of flattened 1D metal cooling data array

// metCooling - metal cooling data

// metHeating - metal heating data

// iVheat - flag for using volumetric heating rate

// iMheat - flag for using specific heating rate

// Vheat - array of volumetric heating rates

// Mheat - array of specific heating rates

// iTfloor - flag for using temperature floor field

// Tfloor_scalar - scalar temperature floor value

// Tfloor - array of temperature floor values

// itmax - maximum allowed sub-cycle iterations

// exititmax - flag to exit if max iterations exceeded

// OUTPUTS:

// update chemical rate densities (HI, HII, etc)

// PARAMETERS:

// mh - H mass in cgs units

// -----------------------------------------------------------------------

// General Arguments

int ierror;

// -- removed line (previously just declared arg types) --

// Density, energy and velocity fields fields

grackle::impl::View<gr_float***> de(my_fields->e_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> HI(my_fields->HI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> HII(my_fields->HII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> HeI(my_fields->HeI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> HeII(my_fields->HeII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> HeIII(my_fields->HeIII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> HM(my_fields->HM_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> H2I(my_fields->H2I_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> H2II(my_fields->H2II_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> DI(my_fields->DI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> DII(my_fields->DII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> HDI(my_fields->HDI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> d(my_fields->density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> e(my_fields->internal_energy, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> u(my_fields->x_velocity, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> v(my_fields->y_velocity, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> w(my_fields->z_velocity, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> metal(my_fields->metal_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> dust(my_fields->dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> Vheat(my_fields->volumetric_heating_rate, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> Mheat(my_fields->specific_heating_rate, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> Tfloor(my_fields->temperature_floor, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> DM(my_fields->DM_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> HDII(my_fields->HDII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> HeHII(my_fields->HeHII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> CI(my_fields->CI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> CII(my_fields->CII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> CO(my_fields->CO_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> CO2(my_fields->CO2_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> OI(my_fields->OI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> OH(my_fields->OH_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> H2O(my_fields->H2O_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> O2(my_fields->O2_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> SiI(my_fields->SiI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> SiOI(my_fields->SiOI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> SiO2I(my_fields->SiO2I_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> CH(my_fields->CH_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> CH2(my_fields->CH2_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> COII(my_fields->COII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> OII(my_fields->OII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> OHII(my_fields->OHII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> H2OII(my_fields->H2OII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> H3OII(my_fields->H3OII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> O2II(my_fields->O2II_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> Mg(my_fields->Mg_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> Al(my_fields->Al_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> S(my_fields->S_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> Fe(my_fields->Fe_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> SiM(my_fields->SiM_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> FeM(my_fields->FeM_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> Mg2SiO4(my_fields->Mg2SiO4_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> MgSiO3(my_fields->MgSiO3_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> Fe3O4(my_fields->Fe3O4_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> AC(my_fields->AC_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> SiO2D(my_fields->SiO2_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> MgO(my_fields->MgO_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> FeS(my_fields->FeS_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> Al2O3(my_fields->Al2O3_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> reforg(my_fields->ref_org_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> volorg(my_fields->vol_org_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> H2Oice(my_fields->H2O_ice_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> metal_loc(my_fields->local_ISM_metal_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> metal_C13(my_fields->ccsn13_metal_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> metal_C20(my_fields->ccsn20_metal_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> metal_C25(my_fields->ccsn25_metal_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> metal_C30(my_fields->ccsn30_metal_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> metal_F13(my_fields->fsn13_metal_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> metal_F15(my_fields->fsn15_metal_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> metal_F50(my_fields->fsn50_metal_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> metal_F80(my_fields->fsn80_metal_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> metal_P170(my_fields->pisn170_metal_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> metal_P200(my_fields->pisn200_metal_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> metal_Y19(my_fields->y19_metal_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

// Radiative transfer fields

grackle::impl::View<gr_float***> kphHI(my_fields->RT_HI_ionization_rate, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> kphHeI(my_fields->RT_HeI_ionization_rate, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> kphHeII(my_fields->RT_HeII_ionization_rate, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> kdissH2I(my_fields->RT_H2_dissociation_rate, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> photogamma(my_fields->RT_heating_rate, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

// -- removed line (previously just declared arg types) --

grackle::impl::View<gr_float***> kdissHDI(my_fields->RT_HDI_dissociation_rate, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> kphCI(my_fields->RT_CI_ionization_rate, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> kphOI(my_fields->RT_OI_ionization_rate, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> kdissCO(my_fields->RT_CO_dissociation_rate, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> kdissOH(my_fields->RT_OH_dissociation_rate, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

grackle::impl::View<gr_float***> kdissH2O(my_fields->RT_H2O_dissociation_rate, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

// -- removed line (previously just declared arg types) --

// H2 self-shielding length-scale field

grackle::impl::View<gr_float***> xH2shield(my_fields->H2_self_shielding_length, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]);

// Interstellar radiation field for dust heating

// Custom H2 shielding factor

// Cooling tables (coolings rates as a function of temperature)

// -- removed line (previously just declared arg types) --

// opacity table

// -- removed line (previously just declared arg types) --

// Chemistry tables (rates as a function of temperature)

// -- removed line (previously just declared arg types) --

// Cloudy cooling data

// -- removed line (previously just declared arg types) --

// Parameters

#ifdef GRACKLE_FLOAT_4

const gr_float tolerance = (gr_float)(1.0e-05);

#else

const gr_float tolerance = (gr_float)(1.0e-10);

#endif

const double mh = mass_h;

const double pi = pi_val;

// Locals

int i, j, k, iter;

int t, dj, dk;

double ttmin, dom, energy, comp1, comp2;

double coolunit, dbase1, tbase1, xbase1, chunit, uvel;

double heq1, heq2, eqk221, eqk222, eqk131, eqk132, eqt1, eqt2, eqtdef, dheq, heq, dlogtem, dx_cgs, c_ljeans, min_metallicity;

gr_float factor;

// row temporaries

std::vector<long long> indixe(my_fields->grid_dimension[0]);

double olddtit;

std::vector<double> t1(my_fields->grid_dimension[0]);

std::vector<double> t2(my_fields->grid_dimension[0]);

std::vector<double> logtem(my_fields->grid_dimension[0]);

std::vector<double> tdef(my_fields->grid_dimension[0]);

std::vector<double> dtit(my_fields->grid_dimension[0]);

std::vector<double> ttot(my_fields->grid_dimension[0]);

std::vector<double> p2d(my_fields->grid_dimension[0]);

std::vector<double> tgas(my_fields->grid_dimension[0]);

std::vector<double> tgasold(my_fields->grid_dimension[0]);

std::vector<double> tdust(my_fields->grid_dimension[0]);

std::vector<double> metallicity(my_fields->grid_dimension[0]);

std::vector<double> dust2gas(my_fields->grid_dimension[0]);

std::vector<double> rhoH(my_fields->grid_dimension[0]);

std::vector<double> mmw(my_fields->grid_dimension[0]);

std::vector<double> mynh(my_fields->grid_dimension[0]);

std::vector<double> myde(my_fields->grid_dimension[0]);

std::vector<double> gammaha_eff(my_fields->grid_dimension[0]);

std::vector<double> gasgr_tdust(my_fields->grid_dimension[0]);

std::vector<double> regr(my_fields->grid_dimension[0]);

std::vector<double> ddom(my_fields->grid_dimension[0]);

// Rate equation row temporaries

std::vector<double> HIp(my_fields->grid_dimension[0]);

std::vector<double> HIIp(my_fields->grid_dimension[0]);

std::vector<double> HeIp(my_fields->grid_dimension[0]);

std::vector<double> HeIIp(my_fields->grid_dimension[0]);

std::vector<double> HeIIIp(my_fields->grid_dimension[0]);

std::vector<double> HMp(my_fields->grid_dimension[0]);

std::vector<double> H2Ip(my_fields->grid_dimension[0]);

std::vector<double> H2IIp(my_fields->grid_dimension[0]);

std::vector<double> dep(my_fields->grid_dimension[0]);

std::vector<double> dedot(my_fields->grid_dimension[0]);

std::vector<double> HIdot(my_fields->grid_dimension[0]);

std::vector<double> dedot_prev(my_fields->grid_dimension[0]);

std::vector<double> DIp(my_fields->grid_dimension[0]);

std::vector<double> DIIp(my_fields->grid_dimension[0]);

std::vector<double> HDIp(my_fields->grid_dimension[0]);

std::vector<double> HIdot_prev(my_fields->grid_dimension[0]);

std::vector<double> k24shield(my_fields->grid_dimension[0]);

std::vector<double> k25shield(my_fields->grid_dimension[0]);

std::vector<double> k26shield(my_fields->grid_dimension[0]);

std::vector<double> k28shield(my_fields->grid_dimension[0]);

std::vector<double> k29shield(my_fields->grid_dimension[0]);

std::vector<double> k30shield(my_fields->grid_dimension[0]);

std::vector<double> k31shield(my_fields->grid_dimension[0]);

std::vector<double> k1(my_fields->grid_dimension[0]);

std::vector<double> k2(my_fields->grid_dimension[0]);

std::vector<double> k3(my_fields->grid_dimension[0]);

std::vector<double> k4(my_fields->grid_dimension[0]);

std::vector<double> k5(my_fields->grid_dimension[0]);

std::vector<double> k6(my_fields->grid_dimension[0]);

std::vector<double> k7(my_fields->grid_dimension[0]);

std::vector<double> k8(my_fields->grid_dimension[0]);

std::vector<double> k9(my_fields->grid_dimension[0]);

std::vector<double> k10(my_fields->grid_dimension[0]);

std::vector<double> k11(my_fields->grid_dimension[0]);

std::vector<double> k12(my_fields->grid_dimension[0]);

std::vector<double> k13(my_fields->grid_dimension[0]);

std::vector<double> k14(my_fields->grid_dimension[0]);

std::vector<double> k15(my_fields->grid_dimension[0]);

std::vector<double> k16(my_fields->grid_dimension[0]);

std::vector<double> k17(my_fields->grid_dimension[0]);

std::vector<double> k18(my_fields->grid_dimension[0]);

std::vector<double> k19(my_fields->grid_dimension[0]);

std::vector<double> k22(my_fields->grid_dimension[0]);

std::vector<double> k50(my_fields->grid_dimension[0]);

std::vector<double> k51(my_fields->grid_dimension[0]);

std::vector<double> k52(my_fields->grid_dimension[0]);

std::vector<double> k53(my_fields->grid_dimension[0]);

std::vector<double> k54(my_fields->grid_dimension[0]);

std::vector<double> k55(my_fields->grid_dimension[0]);

std::vector<double> k56(my_fields->grid_dimension[0]);

std::vector<double> k57(my_fields->grid_dimension[0]);

std::vector<double> k58(my_fields->grid_dimension[0]);

std::vector<double> k13dd(my_fields->grid_dimension[0] * 14);

std::vector<double> h2dust(my_fields->grid_dimension[0]);

std::vector<double> ncrn(my_fields->grid_dimension[0]);

std::vector<double> ncrd1(my_fields->grid_dimension[0]);

std::vector<double> ncrd2(my_fields->grid_dimension[0]);

std::vector<double> DMp(my_fields->grid_dimension[0]);

std::vector<double> HDIIp(my_fields->grid_dimension[0]);

std::vector<double> HeHIIp(my_fields->grid_dimension[0]);

std::vector<double> CIp(my_fields->grid_dimension[0]);

std::vector<double> CIIp(my_fields->grid_dimension[0]);

std::vector<double> COp(my_fields->grid_dimension[0]);

std::vector<double> CO2p(my_fields->grid_dimension[0]);

std::vector<double> OIp(my_fields->grid_dimension[0]);

std::vector<double> OHp(my_fields->grid_dimension[0]);

std::vector<double> H2Op(my_fields->grid_dimension[0]);

std::vector<double> O2p(my_fields->grid_dimension[0]);

std::vector<double> SiIp(my_fields->grid_dimension[0]);

std::vector<double> SiOIp(my_fields->grid_dimension[0]);

std::vector<double> SiO2Ip(my_fields->grid_dimension[0]);

std::vector<double> CHp(my_fields->grid_dimension[0]);

std::vector<double> CH2p(my_fields->grid_dimension[0]);

std::vector<double> COIIp(my_fields->grid_dimension[0]);

std::vector<double> OIIp(my_fields->grid_dimension[0]);

std::vector<double> OHIIp(my_fields->grid_dimension[0]);

std::vector<double> H2OIIp(my_fields->grid_dimension[0]);

std::vector<double> H3OIIp(my_fields->grid_dimension[0]);

std::vector<double> O2IIp(my_fields->grid_dimension[0]);

std::vector<double> Mgp(my_fields->grid_dimension[0]);

std::vector<double> Alp(my_fields->grid_dimension[0]);

std::vector<double> Sp(my_fields->grid_dimension[0]);

std::vector<double> Fep(my_fields->grid_dimension[0]);

std::vector<gr_float> SiMp(my_fields->grid_dimension[0]);

std::vector<gr_float> FeMp(my_fields->grid_dimension[0]);

std::vector<gr_float> Mg2SiO4p(my_fields->grid_dimension[0]);

std::vector<gr_float> MgSiO3p(my_fields->grid_dimension[0]);

std::vector<gr_float> Fe3O4p(my_fields->grid_dimension[0]);

std::vector<gr_float> ACp(my_fields->grid_dimension[0]);

std::vector<gr_float> SiO2Dp(my_fields->grid_dimension[0]);

std::vector<gr_float> MgOp(my_fields->grid_dimension[0]);

std::vector<gr_float> FeSp(my_fields->grid_dimension[0]);

std::vector<gr_float> Al2O3p(my_fields->grid_dimension[0]);

std::vector<gr_float> reforgp(my_fields->grid_dimension[0]);

std::vector<gr_float> volorgp(my_fields->grid_dimension[0]);

std::vector<gr_float> H2Oicep(my_fields->grid_dimension[0]);

std::vector<double> k125(my_fields->grid_dimension[0]);

std::vector<double> k129(my_fields->grid_dimension[0]);

std::vector<double> k130(my_fields->grid_dimension[0]);

std::vector<double> k131(my_fields->grid_dimension[0]);

std::vector<double> k132(my_fields->grid_dimension[0]);

std::vector<double> k133(my_fields->grid_dimension[0]);

std::vector<double> k134(my_fields->grid_dimension[0]);

std::vector<double> k135(my_fields->grid_dimension[0]);

std::vector<double> k136(my_fields->grid_dimension[0]);

std::vector<double> k137(my_fields->grid_dimension[0]);

std::vector<double> k148(my_fields->grid_dimension[0]);

std::vector<double> k149(my_fields->grid_dimension[0]);

std::vector<double> k150(my_fields->grid_dimension[0]);

std::vector<double> k151(my_fields->grid_dimension[0]);

std::vector<double> k152(my_fields->grid_dimension[0]);

std::vector<double> k153(my_fields->grid_dimension[0]);

std::vector<double> kz15(my_fields->grid_dimension[0]);

std::vector<double> kz16(my_fields->grid_dimension[0]);

std::vector<double> kz17(my_fields->grid_dimension[0]);

std::vector<double> kz18(my_fields->grid_dimension[0]);

std::vector<double> kz19(my_fields->grid_dimension[0]);

std::vector<double> kz20(my_fields->grid_dimension[0]);

std::vector<double> kz21(my_fields->grid_dimension[0]);

std::vector<double> kz22(my_fields->grid_dimension[0]);

std::vector<double> kz23(my_fields->grid_dimension[0]);

std::vector<double> kz24(my_fields->grid_dimension[0]);

std::vector<double> kz25(my_fields->grid_dimension[0]);

std::vector<double> kz26(my_fields->grid_dimension[0]);

std::vector<double> kz27(my_fields->grid_dimension[0]);

std::vector<double> kz28(my_fields->grid_dimension[0]);

std::vector<double> kz29(my_fields->grid_dimension[0]);

std::vector<double> kz30(my_fields->grid_dimension[0]);

std::vector<double> kz31(my_fields->grid_dimension[0]);

std::vector<double> kz32(my_fields->grid_dimension[0]);

std::vector<double> kz33(my_fields->grid_dimension[0]);

std::vector<double> kz34(my_fields->grid_dimension[0]);

std::vector<double> kz35(my_fields->grid_dimension[0]);

std::vector<double> kz36(my_fields->grid_dimension[0]);

std::vector<double> kz37(my_fields->grid_dimension[0]);

std::vector<double> kz38(my_fields->grid_dimension[0]);

std::vector<double> kz39(my_fields->grid_dimension[0]);

std::vector<double> kz40(my_fields->grid_dimension[0]);

std::vector<double> kz41(my_fields->grid_dimension[0]);

std::vector<double> kz42(my_fields->grid_dimension[0]);

std::vector<double> kz43(my_fields->grid_dimension[0]);

std::vector<double> kz44(my_fields->grid_dimension[0]);

std::vector<double> kz45(my_fields->grid_dimension[0]);

std::vector<double> kz46(my_fields->grid_dimension[0]);

std::vector<double> kz47(my_fields->grid_dimension[0]);

std::vector<double> kz48(my_fields->grid_dimension[0]);

std::vector<double> kz49(my_fields->grid_dimension[0]);

std::vector<double> kz50(my_fields->grid_dimension[0]);

std::vector<double> kz51(my_fields->grid_dimension[0]);

std::vector<double> kz52(my_fields->grid_dimension[0]);

std::vector<double> kz53(my_fields->grid_dimension[0]);

std::vector<double> kz54(my_fields->grid_dimension[0]);

// -- removed line (previously just declared arg types) --

std::vector<double> kdSiM(my_fields->grid_dimension[0]);

std::vector<double> kdFeM(my_fields->grid_dimension[0]);

std::vector<double> kdMg2SiO4(my_fields->grid_dimension[0]);

std::vector<double> kdMgSiO3(my_fields->grid_dimension[0]);

std::vector<double> kdFe3O4(my_fields->grid_dimension[0]);

std::vector<double> kdAC(my_fields->grid_dimension[0]);

std::vector<double> kdSiO2D(my_fields->grid_dimension[0]);

std::vector<double> kdMgO(my_fields->grid_dimension[0]);

std::vector<double> kdFeS(my_fields->grid_dimension[0]);

std::vector<double> kdAl2O3(my_fields->grid_dimension[0]);

std::vector<double> kdreforg(my_fields->grid_dimension[0]);

std::vector<double> kdvolorg(my_fields->grid_dimension[0]);

std::vector<double> kdH2Oice(my_fields->grid_dimension[0]);

// grain temperature

std::vector<double> tSiM(my_fields->grid_dimension[0]);

std::vector<double> tFeM(my_fields->grid_dimension[0]);

std::vector<double> tMg2SiO4(my_fields->grid_dimension[0]);

std::vector<double> tMgSiO3(my_fields->grid_dimension[0]);

std::vector<double> tFe3O4(my_fields->grid_dimension[0]);

std::vector<double> tAC(my_fields->grid_dimension[0]);

std::vector<double> tSiO2D(my_fields->grid_dimension[0]);

std::vector<double> tMgO(my_fields->grid_dimension[0]);

std::vector<double> tFeS(my_fields->grid_dimension[0]);

std::vector<double> tAl2O3(my_fields->grid_dimension[0]);

std::vector<double> treforg(my_fields->grid_dimension[0]);

std::vector<double> tvolorg(my_fields->grid_dimension[0]);

std::vector<double> tH2Oice(my_fields->grid_dimension[0]);

// Cooling/heating row locals

std::vector<double> ceHI(my_fields->grid_dimension[0]);

std::vector<double> ceHeI(my_fields->grid_dimension[0]);

std::vector<double> ceHeII(my_fields->grid_dimension[0]);

std::vector<double> ciHI(my_fields->grid_dimension[0]);

std::vector<double> ciHeI(my_fields->grid_dimension[0]);

std::vector<double> ciHeIS(my_fields->grid_dimension[0]);

std::vector<double> ciHeII(my_fields->grid_dimension[0]);

std::vector<double> reHII(my_fields->grid_dimension[0]);

std::vector<double> reHeII1(my_fields->grid_dimension[0]);

std::vector<double> reHeII2(my_fields->grid_dimension[0]);

std::vector<double> reHeIII(my_fields->grid_dimension[0]);

std::vector<double> brem(my_fields->grid_dimension[0]);

std::vector<double> edot(my_fields->grid_dimension[0]);

std::vector<double> hyd01k(my_fields->grid_dimension[0]);

std::vector<double> h2k01(my_fields->grid_dimension[0]);

std::vector<double> vibh(my_fields->grid_dimension[0]);

std::vector<double> roth(my_fields->grid_dimension[0]);

std::vector<double> rotl(my_fields->grid_dimension[0]);

std::vector<double> gpldl(my_fields->grid_dimension[0]);

std::vector<double> gphdl(my_fields->grid_dimension[0]);

std::vector<double> hdlte(my_fields->grid_dimension[0]);

std::vector<double> hdlow(my_fields->grid_dimension[0]);

std::vector<double> cieco(my_fields->grid_dimension[0]);

// Iteration mask

gr_mask_type anydust;

std::vector<gr_mask_type> itmask(my_fields->grid_dimension[0]);

std::vector<gr_mask_type> itmask_tmp(my_fields->grid_dimension[0]);

std::vector<gr_mask_type> itmask_nr(my_fields->grid_dimension[0]);

std::vector<gr_mask_type> itmask_metal(my_fields->grid_dimension[0]);

int itr, itr_time;

std::vector<int> imp_eng(my_fields->grid_dimension[0]);

int nsp, isp, jsp, id;

double dspj, err, err_max;

const int i_eng = 52;

std::vector<double> dsp(i_eng);

std::vector<double> dsp0(i_eng);

std::vector<double> dsp1(i_eng);

std::vector<double> dspdot(i_eng);

std::vector<double> dspdot1(i_eng);

std::vector<double> ddsp(i_eng);

std::vector<double> der_data_(i_eng * i_eng);

grackle::impl::View<double**> der(der_data_.data(), i_eng, i_eng);

std::vector<int> idsp;

std::vector<double> mtrx_data_;

grackle::impl::View<double**> mtrx;

std::vector<double> vec;

const double eps = 1.e-4;

// \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\/////////////////////////////////

// =======================================================================

// Set error indicator

(*ierr) = 1;

// Set flag for dust-related options

if ((my_chemistry->h2_on_dust > 0) || (my_chemistry->dust_chemistry > 0)) {

anydust = MASK_TRUE;

} else {

anydust = MASK_FALSE;

}

// ignore metal chemistry/cooling below this metallicity

min_metallicity = 1.e-9 / my_chemistry->SolarMetalFractionByMass;

// Set units

dom = (*urho)*(std::pow(my_units->a_value,3))/mh;

tbase1 = my_units->time_units;

xbase1 = (*uxyz)/(my_units->a_value*my_units->a_units); // uxyz is [x]*a = [x]*[a]*a' '

dbase1 = (*urho)*std::pow((my_units->a_value*my_units->a_units),3); // urho is [dens]/a^3 = [dens]/([a]*a')^3 '

coolunit = (std::pow(my_units->a_units,5) * std::pow(xbase1,2) * std::pow(mh,2)) / (std::pow(tbase1,3) * dbase1);

uvel = ((*uxyz)/my_units->a_value) / my_units->time_units;

// chunit = (1.60218e-12_DKIND)/(2._DKIND*uvel*uvel*mh) ! 1 eV per H2 formed

chunit = (1.60218e-12)/(uvel*uvel*mh); // 1 eV per REACTION (Feb 2020, Gen Chiaki)

dx_cgs = my_fields->grid_dx * xbase1;

c_ljeans = std::sqrt((my_chemistry->Gamma * pi * kboltz) /

(GravConst * mh * dbase1));

dlogtem = (std::log(my_chemistry->TemperatureEnd) - std::log(my_chemistry->TemperatureStart))/(double)(my_chemistry->NumberOfTemperatureBins-1 );

// We better make consistent at first GC202002

if (my_chemistry->primordial_chemistry > 0) {

#define ABUNDANCE_CORRECTION

#ifdef ABUNDANCE_CORRECTION

FORTRAN_NAME(make_consistent_g)(de.data(), HI.data(), HII.data(), HeI.data(), HeII.data(), HeIII.data(),

HM.data(), H2I.data(), H2II.data(), DI.data(), DII.data(), HDI.data(), metal.data(), dust.data(),

d.data(), &my_fields->grid_start[0], &my_fields->grid_end[0], &my_fields->grid_start[1], &my_fields->grid_end[1], &my_fields->grid_start[2], &my_fields->grid_end[2],

&my_fields->grid_dimension[0], &my_fields->grid_dimension[1], &my_fields->grid_dimension[2], &my_chemistry->primordial_chemistry, imetal, &my_chemistry->HydrogenFractionByMass, &my_chemistry->DeuteriumToHydrogenRatio,

&my_chemistry->use_dust_density_field, &my_chemistry->metal_chemistry, &my_chemistry->grain_growth, &dom,

DM.data(), HDII.data(), HeHII.data(),

CI.data(), CII.data(), CO.data(), CO2.data(),

OI.data(), OH.data(), H2O.data(), O2.data(),

SiI.data(), SiOI.data(), SiO2I.data(),

CH.data(), CH2.data(), COII.data(), OII.data(),

OHII.data(), H2OII.data(), H3OII.data(), O2II.data(),

Mg.data(), Al.data(), S.data(), Fe.data(),

SiM.data(), FeM.data(), Mg2SiO4.data(), MgSiO3.data(), Fe3O4.data(),

AC.data(), SiO2D.data(), MgO.data(), FeS.data(), Al2O3.data(),

reforg.data(), volorg.data(), H2Oice.data(),

&my_chemistry->multi_metals, &my_chemistry->metal_abundances, &my_chemistry->dust_species, &my_chemistry->use_multiple_dust_temperatures, &my_chemistry->dust_sublimation,

metal_loc.data(),

metal_C13.data(), metal_C20.data(), metal_C25.data(), metal_C30.data(),

metal_F13.data(), metal_F15.data(), metal_F50.data(), metal_F80.data(),

metal_P170.data(), metal_P200.data(), metal_Y19.data(),

&my_rates->SN0_N,

my_rates->SN0_XC, my_rates->SN0_XO, my_rates->SN0_XMg, my_rates->SN0_XAl, my_rates->SN0_XSi,

my_rates->SN0_XS, my_rates->SN0_XFe,

my_rates->SN0_fC, my_rates->SN0_fO, my_rates->SN0_fMg, my_rates->SN0_fAl, my_rates->SN0_fSi,

my_rates->SN0_fS, my_rates->SN0_fFe

);

#endif

}

// Convert densities from comoving to proper

if (my_units->comoving_coordinates == 1) {

factor = (gr_float)(std::pow(my_units->a_value,(-3)) );

FORTRAN_NAME(scale_fields_g)(

&my_chemistry->primordial_chemistry, imetal, &my_chemistry->use_dust_density_field, &my_chemistry->metal_chemistry,

&my_chemistry->metal_abundances, &my_chemistry->dust_species, &my_chemistry->multi_metals, &my_chemistry->grain_growth, &my_chemistry->dust_sublimation, &factor,

&my_fields->grid_start[0], &my_fields->grid_end[0], &my_fields->grid_start[1], &my_fields->grid_end[1], &my_fields->grid_start[2], &my_fields->grid_end[2], &my_fields->grid_dimension[0], &my_fields->grid_dimension[1], &my_fields->grid_dimension[2],

d.data(), de.data(), HI.data(), HII.data(), HeI.data(), HeII.data(), HeIII.data(),

HM.data(), H2I.data(), H2II.data(), DI.data(), DII.data(), HDI.data(), DM.data(), HDII.data(), HeHII.data(),

metal.data(), dust.data(),

CI.data(), CII.data(), CO.data(), CO2.data(),

OI.data(), OH.data(), H2O.data(), O2.data(),

SiI.data(), SiOI.data(), SiO2I.data(),

CH.data(), CH2.data(), COII.data(), OII.data(),

OHII.data(), H2OII.data(), H3OII.data(), O2II.data(),

Mg.data(), Al.data(), S.data(), Fe.data(),

SiM.data(), FeM.data(), Mg2SiO4.data(), MgSiO3.data(), Fe3O4.data(),

AC.data(), SiO2D.data(), MgO.data(), FeS.data(), Al2O3.data(),

reforg.data(), volorg.data(), H2Oice.data(),

metal_loc.data(), metal_C13.data(), metal_C20.data(), metal_C25.data(), metal_C30.data(),

metal_F13.data(), metal_F15.data(), metal_F50.data(), metal_F80.data(),

metal_P170.data(), metal_P200.data(), metal_Y19.data());

}

#ifdef ABUNDANCE_CORRECTION

FORTRAN_NAME(ceiling_species_g)(d.data(), de.data(), HI.data(), HII.data(), HeI.data(), HeII.data(), HeIII.data(),

HM.data(), H2I.data(), H2II.data(), DI.data(), DII.data(), HDI.data(), metal.data(), dust.data(),

&my_fields->grid_start[0], &my_fields->grid_end[0], &my_fields->grid_start[1], &my_fields->grid_end[1], &my_fields->grid_start[2], &my_fields->grid_end[2],

&my_fields->grid_dimension[0], &my_fields->grid_dimension[1], &my_fields->grid_dimension[2], &my_chemistry->primordial_chemistry, imetal, &my_chemistry->use_dust_density_field,

DM.data(), HDII.data(), HeHII.data(),

&my_chemistry->metal_abundances, &my_chemistry->metal_chemistry, &my_chemistry->dust_species, &my_chemistry->multi_metals,

&my_chemistry->grain_growth, &my_chemistry->dust_sublimation,

CI.data(), CII.data(), CO.data(), CO2.data(),

OI.data(), OH.data(), H2O.data(), O2.data(),

SiI.data(), SiOI.data(), SiO2I.data(),

CH.data(), CH2.data(), COII.data(), OII.data(),

OHII.data(), H2OII.data(), H3OII.data(), O2II.data(),

Mg.data(), Al.data(), S.data(), Fe.data(),

SiM.data(), FeM.data(), Mg2SiO4.data(), MgSiO3.data(), Fe3O4.data(),

AC.data(), SiO2D.data(), MgO.data(), FeS.data(), Al2O3.data(),

reforg.data(), volorg.data(), H2Oice.data(),

metal_loc.data(),

metal_C13.data(), metal_C20.data(), metal_C25.data(), metal_C30.data(),

metal_F13.data(), metal_F15.data(), metal_F50.data(), metal_F80.data(),

metal_P170.data(), metal_P200.data(), metal_Y19.data());

#endif

// Loop over zones, and do an entire i-column in one go

dk = my_fields->grid_end[2] - my_fields->grid_start[2] + 1;

dj = my_fields->grid_end[1] - my_fields->grid_start[1] + 1;

// parallelize the k and j loops with OpenMP

// flat j and k loops for better parallelism

//_// PORT: #ifdef _OPENMP

//_// PORT: ! ierr is declared as shared and should be modified with atomic operation

//_// PORT: !$omp parallel do schedule(runtime) private(

//_// PORT: !$omp& i, j, k, iter,

//_// PORT: !$omp& ttmin, energy, comp1, comp2,

//_// PORT: !$omp& heq1, heq2, eqk221, eqk222, eqk131, eqk132,

//_// PORT: !$omp& eqt1, eqt2, eqtdef, dheq, heq,

//_// PORT: !$omp& indixe,

//_// PORT: !$omp& t1, t2, logtem, tdef,

//_// PORT: !$omp& dtit, ttot, p2d, tgas, tgasold,

//_// PORT: !$omp& tdust, metallicity, dust2gas, rhoH, mmw,

//_// PORT: !$omp& mynh, myde, gammaha_eff, gasgr_tdust, regr, ddom,

//_// PORT: !$omp& olddtit,

//_// PORT: !$omp& HIp, HIIp, HeIp, HeIIp, HeIIIp,

//_// PORT: !$omp& HMp, H2Ip, H2IIp,

//_// PORT: !$omp& dep, dedot,HIdot, dedot_prev,

//_// PORT: !$omp& DIp, DIIp, HDIp, HIdot_prev,

//_// PORT: !$omp& k24shield, k25shield, k26shield,

//_// PORT: !$omp& k28shield, k29shield, k30shield,

//_// PORT: !$omp& k31shield,

//_// PORT: !$omp& k1 , k2 , k3 , k4 , k5,

//_// PORT: !$omp& k6 , k7 , k8 , k9 , k10,

//_// PORT: !$omp& k11, k12, k13, k14, k15,

//_// PORT: !$omp& k16, k17, k18, k19, k22,

//_// PORT: !$omp& k50, k51, k52, k53, k54,

//_// PORT: !$omp& k55, k56, k57, k58, k13dd, h2dust,

//_// PORT: !$omp& ncrn, ncrd1, ncrd2,

//_// PORT: !$omp& ceHI, ceHeI, ceHeII,

//_// PORT: !$omp& ciHI, ciHeI, ciHeIS, ciHeII,

//_// PORT: !$omp& reHII, reHeII1, reHeII2, reHeIII,

//_// PORT: !$omp& brem, edot,

//_// PORT: !$omp& hyd01k, h2k01, vibh, roth, rotl,

//_// PORT: !$omp& gpldl, gphdl, hdlte, hdlow, cieco,

//_// PORT: !$omp& itmask, itmask_metal )

//_// PORT: #endif

for (t = 0; t<=(dk * dj - 1); t++) {

k = t/dj + my_fields->grid_start[2]+1;

j = grackle::impl::mod(t,dj) + my_fields->grid_start[1]+1;

// tolerance = 1.0e-06_DKIND * dt

// Initialize iteration mask to true for all cells.

for (i = my_fields->grid_start[0] + 1; i<=(my_fields->grid_end[0] + 1); i++) {

itmask[i-1] = MASK_TRUE;

}

// If we are using coupled radiation with intermediate stepping,

// set iteration mask to include only cells with radiation in the

// intermediate coupled chemistry / energy step

if (my_chemistry->use_radiative_transfer == 1) {

if (my_chemistry->radiative_transfer_coupled_rate_solver == 1 && my_chemistry->radiative_transfer_intermediate_step == 1) {

for (i = my_fields->grid_start[0] + 1; i<=(my_fields->grid_end[0] + 1); i++) {

if (kphHI(i-1,j-1,k-1) > 0) {

itmask[i-1] = MASK_TRUE;

} else {

itmask[i-1] = MASK_FALSE;

}

// Normal rate solver, but don't double count cells with radiation

if (my_chemistry->radiative_transfer_coupled_rate_solver == 1 && my_chemistry->radiative_transfer_intermediate_step == 0) {

for (i = my_fields->grid_start[0] + 1; i<=(my_fields->grid_end[0] + 1); i++) {

if (kphHI(i-1,j-1,k-1) > 0) {

itmask[i-1] = MASK_FALSE;

} else {

itmask[i-1] = MASK_TRUE;

}

// Set time elapsed to zero for each cell in 1D section

for (i = my_fields->grid_start[0] + 1; i<=(my_fields->grid_end[0] + 1); i++) {

ttot[i-1] = 0.;

}

// A useful slice variable since we do this a lot

for (i = my_fields->grid_start[0] + 1; i<=(my_fields->grid_end[0] + 1); i++) {

ddom[i-1] = d(i-1,j-1,k-1) * dom;

}

// ------------------ Loop over subcycles ----------------

for (iter = 1; iter<=(my_chemistry->max_iterations); iter++) {

for (i = my_fields->grid_start[0] + 1; i<=(my_fields->grid_end[0] + 1); i++) {

if (itmask[i-1] != MASK_FALSE) {

dtit[i-1] = huge8;

}

// Compute the cooling rate, tgas, tdust, and metallicity for this row// Look-up rates as a function of temperature for 1D set of zones

// (maybe should add itmask to this call)// Compute dedot and HIdot, the rates of change of de and HI

// (should add itmask to this call)

FORTRAN_NAME(rate_timestep_g)(

dedot.data(), HIdot.data(), &my_chemistry->primordial_chemistry, &anydust,

de.data(), HI.data(), HII.data(), HeI.data(), HeII.data(), HeIII.data(), d.data(),

HM.data(), H2I.data(), H2II.data(),

&my_fields->grid_dimension[0], &my_fields->grid_dimension[1], &my_fields->grid_dimension[2], &my_fields->grid_start[0], &my_fields->grid_end[0], &j, &k,

k1.data(), k2.data(), k3.data(), k4.data(), k5.data(), k6.data(), k7.data(), k8.data(), k9.data(), k10.data(), k11.data(),

k12.data(), k13.data(), k14.data(), k15.data(), k16.data(), k17.data(), k18.data(), k19.data(), k22.data(),

&my_uvb_rates->k24, &my_uvb_rates->k25, &my_uvb_rates->k26, &my_uvb_rates->k27, &my_uvb_rates->k28, &my_uvb_rates->k29, &my_uvb_rates->k30,

k50.data(), k51.data(), k52.data(), k53.data(), k54.data(), k55.data(), k56.data(), k57.data(), k58.data(),

h2dust.data(), ncrn.data(), ncrd1.data(), ncrd2.data(), rhoH.data(),

k24shield.data(), k25shield.data(), k26shield.data(),

k28shield.data(), k29shield.data(), k30shield.data(), k31shield.data(),

&my_chemistry->use_radiative_transfer, &my_chemistry->radiative_transfer_hydrogen_only,

kphHI.data(), kphHeI.data(), kphHeII.data(),

itmask.data(), edot.data(), &chunit, &dom, metal.data(),

HDI.data(), &my_chemistry->metal_chemistry, CI.data(), OI.data(), OH.data(), CO.data(), H2O.data(),

&my_chemistry->radiative_transfer_HDI_dissociation, kdissHDI.data(), &my_chemistry->radiative_transfer_metal_ionization, kphCI.data(), kphOI.data(),

&my_chemistry->radiative_transfer_metal_dissociation, kdissCO.data(), kdissOH.data(), kdissH2O.data()

);

// move itmask temporary array

// then split cells with low densities

// => Gauss-Seidel scheme

// and with high densities

// => Newton-Raphson scheme

std::memcpy(itmask_tmp.data(), itmask.data(), sizeof(gr_mask_type)*my_fields->grid_dimension[0]);

std::memcpy(itmask_nr.data(), itmask.data(), sizeof(gr_mask_type)*my_fields->grid_dimension[0]);

for (i = my_fields->grid_start[0] + 1; i<=(my_fields->grid_end[0] + 1); i++) {

if ( itmask_tmp[i-1] != MASK_FALSE ) {

if ( ( ((*imetal) == 0)

&& (ddom[i-1] < 1.e8) )

|| ( ((*imetal) == 1)

&& ( ( (metallicity[i-1] <= min_metallicity)

&& (ddom[i-1] < 1.e8) )

|| ( (metallicity[i-1] > min_metallicity)

&& (ddom[i-1] < 1.e6) ) ) ) ) {

itmask_nr[i-1] = MASK_FALSE;

} else {

itmask[i-1] = MASK_FALSE;

}

for (i = my_fields->grid_start[0] + 1; i<=(my_fields->grid_end[0] + 1); i++) {

if (itmask_nr[i-1] != MASK_FALSE) {

if ( (my_chemistry->with_radiative_cooling == 1) && (my_chemistry->primordial_chemistry > 1) &&

((ddom[i-1] > 1.e7)

&& (tgas[i-1] > 1650.e0)) ) {

imp_eng[i-1] = 1;

} else {

imp_eng[i-1] = 0;

}

// Find timestep that keeps relative chemical changes below 10%

for (i = my_fields->grid_start[0] + 1; i<=(my_fields->grid_end[0] + 1); i++) {

if (itmask[i-1] != MASK_FALSE) {

// Bound from below to prevent numerical errors

if (std::fabs(dedot[i-1]) < tiny8)

{ dedot[i-1] = std::fmin(tiny,de(i-1,j-1,k-1)); }

if (std::fabs(HIdot[i-1]) < tiny8)

{ HIdot[i-1] = std::fmin(tiny,HI(i-1,j-1,k-1)); }

// If the net rate is almost perfectly balanced then set

// it to zero (since it is zero to available precision)

if (std::fmin(std::fabs(k1[i-1]* de(i-1,j-1,k-1)*HI(i-1,j-1,k-1)),

std::fabs(k2[i-1]*HII(i-1,j-1,k-1)*de(i-1,j-1,k-1)))/

std::fmax(std::fabs(dedot[i-1]),std::fabs(HIdot[i-1])) >

1.0e6) {

dedot[i-1] = tiny8;

HIdot[i-1] = tiny8;

}

// If the iteration count is high then take the smaller of

// the calculated dedot and last time step's actual dedot.

// This is intended to get around the problem of a low

// electron or HI fraction which is in equilibrium with high

// individual terms (which all nearly cancel).

if (iter > 50) {

dedot[i-1] = std::fmin(std::fabs(dedot[i-1]), std::fabs(dedot_prev[i-1]));

HIdot[i-1] = std::fmin(std::fabs(HIdot[i-1]), std::fabs(HIdot_prev[i-1]));

}

// compute minimum rate timestep

olddtit = dtit[i-1];

dtit[i-1] = grackle::impl::fmin(std::fabs(0.1*de(i-1,j-1,k-1)/dedot[i-1]),

std::fabs(0.1*HI(i-1,j-1,k-1)/HIdot[i-1]),

(*dt)-ttot[i-1], 0.5*(*dt));

if (ddom[i-1] > 1.e8 &&

edot[i-1] > 0. &&

my_chemistry->primordial_chemistry > 1) {

// Equilibrium value for H is:

// H = (-1._DKIND / (4*k22)) * (k13 - sqrt(8 k13 k22 rho + k13^2))

// We now want this to change by 10% or less, but we're only

// differentiating by dT. We have de/dt. We need dT/de.

// T = (g-1)*p2d*utem/N; tgas == (g-1)(p2d*utem/N)

// dH_eq / dt = (dH_eq/dT) * (dT/de) * (de/dt)

// dH_eq / dT (see above; we can calculate the derivative here)

// dT / de = utem * (gamma - 1._DKIND) / N == tgas / p2d

// de / dt = edot

// Now we use our estimate of dT/de to get the estimated

// difference in the equilibrium

eqt2 = std::fmin(std::log(tgas[i-1]) + 0.1*dlogtem, t2[i-1]);

eqtdef = (eqt2 - t1[i-1])/(t2[i-1] - t1[i-1]);

eqk222 = my_rates->k22[indixe[i-1]-1] +

(my_rates->k22[indixe[i-1]+1-1] -my_rates->k22[indixe[i-1]-1])*eqtdef;

eqk132 = my_rates->k13[indixe[i-1]-1] +

(my_rates->k13[indixe[i-1]+1-1] -my_rates->k13[indixe[i-1]-1])*eqtdef;

heq2 = (-1. / (4.*eqk222)) * (eqk132-

std::sqrt(8.*eqk132*eqk222*

my_chemistry->HydrogenFractionByMass*d(i-1,j-1,k-1)+std::pow(eqk132,2.)));

eqt1 = std::fmax(std::log(tgas[i-1]) - 0.1*dlogtem, t1[i-1]);

eqtdef = (eqt1 - t1[i-1])/(t2[i-1] - t1[i-1]);

eqk221 = my_rates->k22[indixe[i-1]-1] +

(my_rates->k22[indixe[i-1]+1-1] -my_rates->k22[indixe[i-1]-1])*eqtdef;

eqk131 = my_rates->k13[indixe[i-1]-1] +

(my_rates->k13[indixe[i-1]+1-1] -my_rates->k13[indixe[i-1]-1])*eqtdef;

heq1 = (-1. / (4.*eqk221)) * (eqk131-

std::sqrt(8.*eqk131*eqk221*

my_chemistry->HydrogenFractionByMass*d(i-1,j-1,k-1)+std::pow(eqk131,2.)));

dheq = (std::fabs(heq2-heq1)/(std::exp(eqt2) - std::exp(eqt1)))

* (tgas[i-1]/p2d[i-1]) * edot[i-1];

heq = (-1. / (4.*k22[i-1])) * (k13[i-1]-

std::sqrt(8.*k13[i-1]*k22[i-1]*

my_chemistry->HydrogenFractionByMass*d(i-1,j-1,k-1)+std::pow(k13[i-1],2.)));

dtit[i-1] = std::fmin(dtit[i-1], 0.1*heq/dheq);

}

if (iter>10LL) {

dtit[i-1] = std::fmin(olddtit*1.5, dtit[i-1]);

}

} else if ((itmask_nr[i-1]!=MASK_FALSE) &&

(imp_eng[i-1]==0)) {

dtit[i-1] = grackle::impl::fmin(std::fabs(0.1*e(i-1,j-1,k-1)/edot[i-1]*d(i-1,j-1,k-1)),

(*dt)-ttot[i-1], 0.5*(*dt));

} else if ((itmask_nr[i-1]!=MASK_FALSE) &&

(imp_eng[i-1]==1)) {

dtit[i-1] = (*dt) - ttot[i-1];

} else {

dtit[i-1] = (*dt);

}

// Compute maximum timestep for cooling/heating

for (i = my_fields->grid_start[0] + 1; i<=(my_fields->grid_end[0] + 1); i++) {

if (itmask[i-1] != MASK_FALSE) {

// Set energy per unit volume of this cell based in the pressure

// (the gamma used here is the right one even for H2 since p2d

// is calculated with this gamma).

energy = std::fmax(p2d[i-1]/(my_chemistry->Gamma-1.), tiny8);

// If the temperature is at the bottom of the temperature look-up

// table and edot < 0, then shut off the cooling.

if (tgas[i-1] <= 1.01*my_chemistry->TemperatureStart &&

edot[i-1] < 0.)

{ edot[i-1] = tiny8; }

if (std::fabs(edot[i-1]) < tiny8) { edot[i-1] = tiny8; }

// Compute timestep for 10% change

dtit[i-1] = grackle::impl::fmin((double)(std::fabs(0.1*

energy/edot[i-1]) ),

(*dt)-ttot[i-1], dtit[i-1]);

if (dtit[i-1] != dtit[i-1]) {

OMP_PRAGMA_CRITICAL

{

eprintf("HUGE dtit :: %g %g %g %g %g %g %g\n",

energy,

edot [ i-1 ],

dtit [ i-1 ],

(*dt),

ttot [ i-1 ],

std::fabs ( 0.1 * energy / edot [ i-1 ] ),

(double) ( std::fabs ( 0.1 * energy / edot [ i-1 ] ) ));

}

// Update total and gas energy

if (my_chemistry->with_radiative_cooling == 1) {

for (i = my_fields->grid_start[0] + 1; i<=(my_fields->grid_end[0] + 1); i++) {

if (itmask[i-1] != MASK_FALSE) {

e(i-1,j-1,k-1) = e(i-1,j-1,k-1) +

(gr_float)(edot[i-1]/d(i-1,j-1,k-1)*dtit[i-1] );

}

if (my_chemistry->primordial_chemistry > 0) {

// Solve rate equations with one linearly implicit Gauss-Seidel

// sweep of a backward Euler method ---

FORTRAN_NAME(step_rate_g)(de.data(), HI.data(), HII.data(), HeI.data(), HeII.data(), HeIII.data(), d.data(),

HM.data(), H2I.data(), H2II.data(), DI.data(), DII.data(), HDI.data(), dtit.data(),

&my_fields->grid_dimension[0], &my_fields->grid_dimension[1], &my_fields->grid_dimension[2], &my_fields->grid_start[0], &my_fields->grid_end[0], &j, &k,

&my_chemistry->primordial_chemistry, &anydust,

k1.data(), k2.data(), k3.data(), k4.data(), k5.data(), k6.data(), k7.data(), k8.data(), k9.data(), k10.data(), k11.data(),

k12.data(), k13.data(), k14.data(), k15.data(), k16.data(), k17.data(), k18.data(), k19.data(), k22.data(),

&my_uvb_rates->k24, &my_uvb_rates->k25, &my_uvb_rates->k26, &my_uvb_rates->k27, &my_uvb_rates->k28, &my_uvb_rates->k29, &my_uvb_rates->k30,

k50.data(), k51.data(), k52.data(), k53.data(), k54.data(), k55.data(), k56.data(), k57.data(), k58.data(),

h2dust.data(), rhoH.data(),

k24shield.data(), k25shield.data(), k26shield.data(),

k28shield.data(), k29shield.data(), k30shield.data(), k31shield.data(),

HIp.data(), HIIp.data(), HeIp.data(), HeIIp.data(), HeIIIp.data(), dep.data(),

HMp.data(), H2Ip.data(), H2IIp.data(), DIp.data(), DIIp.data(), HDIp.data(),

dedot_prev.data(), HIdot_prev.data(),

&my_chemistry->use_radiative_transfer, &my_chemistry->radiative_transfer_hydrogen_only,

kphHI.data(), kphHeI.data(), kphHeII.data(),

itmask.data(), itmask_metal.data(),

DM.data(), HDII.data(), HeHII.data(), imetal, metal.data(),

&my_chemistry->metal_chemistry, &my_chemistry->dust_species, &my_chemistry->grain_growth, &my_chemistry->dust_sublimation,

CI.data(), CII.data(), CO.data(), CO2.data(),

OI.data(), OH.data(), H2O.data(), O2.data(),

SiI.data(), SiOI.data(), SiO2I.data(),

CH.data(), CH2.data(), COII.data(), OII.data(),

OHII.data(), H2OII.data(), H3OII.data(), O2II.data(),

Mg.data(), Al.data(), S.data(), Fe.data(),

SiM.data(), FeM.data(), Mg2SiO4.data(), MgSiO3.data(), Fe3O4.data(),

AC.data(), SiO2D.data(), MgO.data(), FeS.data(), Al2O3.data(),

reforg.data(), volorg.data(), H2Oice.data(),

k125.data(), k129.data(), k130.data(), k131.data(), k132.data(),

k133.data(), k134.data(), k135.data(), k136.data(), k137.data(),

k148.data(), k149.data(), k150.data(), k151.data(), k152.data(),

k153.data(),

kz15.data(), kz16.data(), kz17.data(), kz18.data(), kz19.data(),

kz20.data(), kz21.data(), kz22.data(), kz23.data(), kz24.data(),

kz25.data(), kz26.data(), kz27.data(), kz28.data(), kz29.data(),

kz30.data(), kz31.data(), kz32.data(), kz33.data(), kz34.data(),

kz35.data(), kz36.data(), kz37.data(), kz38.data(), kz39.data(),

kz40.data(), kz41.data(), kz42.data(), kz43.data(), kz44.data(),

kz45.data(), kz46.data(), kz47.data(), kz48.data(), kz49.data(),

kz50.data(), kz51.data(), kz52.data(), kz53.data(), kz54.data(),

DMp.data(), HDIIp.data(), HeHIIp.data(),

CIp.data(), CIIp.data(), COp.data(), CO2p.data(),

OIp.data(), OHp.data(), H2Op.data(), O2p.data(),

SiIp.data(), SiOIp.data(), SiO2Ip.data(),

CHp.data(), CH2p.data(), COIIp.data(), OIIp.data(),

OHIIp.data(), H2OIIp.data(), H3OIIp.data(), O2IIp.data(),

Mgp.data(), Alp.data(), Sp.data(), Fep.data(),

SiMp.data(), FeMp.data(), Mg2SiO4p.data(), MgSiO3p.data(), Fe3O4p.data(),

ACp.data(), SiO2Dp.data(), MgOp.data(), FeSp.data(), Al2O3p.data(),

reforgp.data(), volorgp.data(), H2Oicep.data(),

kdSiM.data(), kdFeM.data(), kdMg2SiO4.data(), kdMgSiO3.data(), kdFe3O4.data(),

kdAC.data(), kdSiO2D.data(), kdMgO.data(), kdFeS.data(), kdAl2O3.data(),

kdreforg.data(), kdvolorg.data(), kdH2Oice.data(),

&my_chemistry->radiative_transfer_HDI_dissociation, kdissHDI.data(), &my_chemistry->radiative_transfer_metal_ionization, kphCI.data(), kphOI.data(),

&my_chemistry->radiative_transfer_metal_dissociation, kdissCO.data(), kdissOH.data(), kdissH2O.data()

);

// Note 10/18/2024: the code from here to the comment "end Newton-Raphson scheme"

// should be put into its own function.

// Start Newton-Raphson scheme

for (i = my_fields->grid_start[0] + 1; i<=(my_fields->grid_end[0] + 1); i++) {

if (itmask_nr[i-1] != MASK_FALSE) {

// If density and temperature are low, update gas energy explicitly

if (my_chemistry->with_radiative_cooling == 1) {

if (imp_eng[i-1] == 0) {

e(i-1,j-1,k-1) = e(i-1,j-1,k-1) +

(gr_float)(edot[i-1]/d(i-1,j-1,k-1)*dtit[i-1] );

}

// initialize arrays

if (my_chemistry->primordial_chemistry > 0) { nsp = 6; }

if (my_chemistry->primordial_chemistry > 1) { nsp = nsp + 3; }

if (my_chemistry->primordial_chemistry > 2) { nsp = nsp + 3; }

if (my_chemistry->primordial_chemistry > 3) { nsp = nsp + 3; }

if (itmask_metal[i-1] != MASK_FALSE) {

if (my_chemistry->metal_chemistry == 1) {

nsp = nsp + 19;

if ( ( my_chemistry->grain_growth == 1 ) || ( my_chemistry->dust_sublimation == 1) ) {

if (my_chemistry->dust_species > 0) { nsp = nsp + 1; }

if (my_chemistry->dust_species > 1) { nsp = nsp + 3; }

}

if ( ( my_chemistry->grain_growth == 1 ) || ( my_chemistry->dust_sublimation == 1) ) {

if (my_chemistry->dust_species > 0) { nsp = nsp + 2; }

if (my_chemistry->dust_species > 1) { nsp = nsp + 8; }

if (my_chemistry->dust_species > 2) { nsp = nsp + 3; }

}

nsp = nsp + imp_eng[i-1];

idsp.reserve(nsp);

mtrx_data_.reserve(nsp * nsp);

mtrx = grackle::impl::View<double**>(mtrx_data_.data(), nsp, nsp);

vec.reserve(nsp);

if ( my_chemistry->primordial_chemistry > 0 ) {

dsp[ 1-1] = de(i-1,j-1,k-1);

dsp[ 2-1] = HI(i-1,j-1,k-1);

dsp[ 3-1] = HII(i-1,j-1,k-1);

dsp[ 4-1] = HeI(i-1,j-1,k-1);

dsp[ 5-1] = HeII(i-1,j-1,k-1);

dsp[ 6-1] = HeIII(i-1,j-1,k-1);

}

if ( my_chemistry->primordial_chemistry > 1 ) {

dsp[ 7-1] = HM(i-1,j-1,k-1);

dsp[ 8-1] = H2I(i-1,j-1,k-1);

dsp[ 9-1] = H2II(i-1,j-1,k-1);

}

if ( my_chemistry->primordial_chemistry > 2 ) {

dsp[10-1] = DI(i-1,j-1,k-1);

dsp[11-1] = DII(i-1,j-1,k-1);

dsp[12-1] = HDI(i-1,j-1,k-1);

}

if ( my_chemistry->primordial_chemistry > 3 ) {

dsp[13-1] = DM(i-1,j-1,k-1);

dsp[14-1] = HDII(i-1,j-1,k-1);

dsp[15-1] = HeHII(i-1,j-1,k-1);

}

if ( itmask_metal[i-1] != MASK_FALSE ) {

if ( my_chemistry->metal_chemistry == 1 ) {

dsp[16-1] = CI(i-1,j-1,k-1);

dsp[17-1] = CII(i-1,j-1,k-1);

dsp[18-1] = CO(i-1,j-1,k-1);

dsp[19-1] = CO2(i-1,j-1,k-1);

dsp[20-1] = OI(i-1,j-1,k-1);

dsp[21-1] = OH(i-1,j-1,k-1);

dsp[22-1] = H2O(i-1,j-1,k-1);

dsp[23-1] = O2(i-1,j-1,k-1);

dsp[24-1] = SiI(i-1,j-1,k-1);

dsp[25-1] = SiOI(i-1,j-1,k-1);

dsp[26-1] = SiO2I(i-1,j-1,k-1);

dsp[27-1] = CH(i-1,j-1,k-1);

dsp[28-1] = CH2(i-1,j-1,k-1);

dsp[29-1] = COII(i-1,j-1,k-1);

dsp[30-1] = OII(i-1,j-1,k-1);

dsp[31-1] = OHII(i-1,j-1,k-1);

dsp[32-1] = H2OII(i-1,j-1,k-1);

dsp[33-1] = H3OII(i-1,j-1,k-1);

dsp[34-1] = O2II(i-1,j-1,k-1);

if ( ( my_chemistry->grain_growth == 1 ) || ( my_chemistry->dust_sublimation == 1) ) {

if (my_chemistry->dust_species > 0) {

dsp[35-1] = Mg(i-1,j-1,k-1);

}

if (my_chemistry->dust_species > 1) {

dsp[36-1] = Al(i-1,j-1,k-1);

dsp[37-1] = S(i-1,j-1,k-1);

dsp[38-1] = Fe(i-1,j-1,k-1);

}

if ( ( my_chemistry->grain_growth == 1 ) || ( my_chemistry->dust_sublimation == 1) ) {

if (my_chemistry->dust_species > 0) {

dsp[39-1] = MgSiO3(i-1,j-1,k-1);

dsp[40-1] = AC(i-1,j-1,k-1);

}

if (my_chemistry->dust_species > 1) {

dsp[41-1] = SiM(i-1,j-1,k-1);

dsp[42-1] = FeM(i-1,j-1,k-1);

dsp[43-1] = Mg2SiO4(i-1,j-1,k-1);

dsp[44-1] = Fe3O4(i-1,j-1,k-1);

dsp[45-1] = SiO2D(i-1,j-1,k-1);

dsp[46-1] = MgO(i-1,j-1,k-1);

dsp[47-1] = FeS(i-1,j-1,k-1);

dsp[48-1] = Al2O3(i-1,j-1,k-1);

}

if (my_chemistry->dust_species > 2) {

dsp[49-1] = reforg(i-1,j-1,k-1);

dsp[50-1] = volorg(i-1,j-1,k-1);

dsp[51-1] = H2Oice(i-1,j-1,k-1);

}

dsp[i_eng-1] = e(i-1,j-1,k-1);

id = 0;

if (my_chemistry->primordial_chemistry > 0) {

for (isp = 1; isp<=(6); isp++) {

id = id + 1;

idsp[id-1] = isp;

}

if (my_chemistry->primordial_chemistry > 1) {

for (isp = 7; isp<=(9); isp++) {

id = id + 1;

idsp[id-1] = isp;

}

if (my_chemistry->primordial_chemistry > 2) {

for (isp = 10; isp<=(12); isp++) {

id = id + 1;

idsp[id-1] = isp;

}

if (my_chemistry->primordial_chemistry > 3) {

for (isp = 13; isp<=(15); isp++) {

id = id + 1;

idsp[id-1] = isp;

}

if (itmask_metal[i-1] != MASK_FALSE) {

if (my_chemistry->metal_chemistry == 1) {

for (isp = 16; isp<=(34); isp++) {

id = id + 1;

idsp[id-1] = isp;

}

if ( ( my_chemistry->grain_growth == 1 ) || ( my_chemistry->dust_sublimation == 1) ) {

if (my_chemistry->dust_species > 0) {

for (isp = 35; isp<=(35); isp++) {

id = id + 1;

idsp[id-1] = isp;

}

if (my_chemistry->dust_species > 1) {

for (isp = 36; isp<=(38); isp++) {

id = id + 1;

idsp[id-1] = isp;

}

if ( ( my_chemistry->grain_growth == 1 ) || ( my_chemistry->dust_sublimation == 1) ) {

if (my_chemistry->dust_species > 0) {

for (isp = 39; isp<=(40); isp++) {

id = id + 1;

idsp[id-1] = isp;

}

if (my_chemistry->dust_species > 1) {

for (isp = 41; isp<=(48); isp++) {

id = id + 1;

idsp[id-1] = isp;

}

if (my_chemistry->dust_species > 2) {

for (isp = 49; isp<=(51); isp++) {

id = id + 1;

idsp[id-1] = isp;

}

if ( imp_eng[i-1] ==1 ) {

id = id + 1;

idsp[id-1] = i_eng;

}

// Save arrays at ttot(i)

std::memcpy(dsp0.data(), dsp.data(), sizeof(double)*i_eng);

std::memset(ddsp.data(), 0, sizeof(double)*i_eng);

// Search for the timestep for which chemistry converges

ierror=1;

itr_time=0;

while ((ierror==1)) {

// If not converge, restore arrays at ttot(i)

std::memcpy(dsp.data(), dsp0.data(), sizeof(double)*i_eng);

std::memset(ddsp.data(), 0, sizeof(double)*i_eng);

// Iteration to solve ODEs// to get more accuracy

for (isp = 1; isp<=(nsp); isp++) {

vec[isp-1] = vec[isp-1]/d(i-1,j-1,k-1);

}

FORTRAN_NAME(gaussj_g)(&nsp, mtrx.data(), vec.data(), &ierror);

if(ierror == 1) {

goto label_9998;

}

// multiply with density again

for (isp = 1; isp<=(nsp); isp++) {

vec[isp-1] = vec[isp-1]*d(i-1,j-1,k-1);

}

for (isp = 1; isp<=(nsp); isp++) {

ddsp[idsp[isp-1]-1] = ddsp[idsp[isp-1]-1] + vec[isp-1];

dsp[idsp[isp-1]-1] = dsp[idsp[isp-1]-1] + vec[isp-1];

}

if (imp_eng[i-1] == 1) {

if( (my_chemistry->primordial_chemistry > 0) && (my_chemistry->with_radiative_cooling == 1) ) {

for (isp = 1; isp<=(nsp); isp++) {

if ( (dsp[idsp[isp-1]-1] != dsp[idsp[isp-1]-1])

|| (dsp[idsp[isp-1]-1] <= 0.) ) {

ierror = 1;

goto label_9997;

}

err_max = 0.e0;

for (isp = 1; isp<=(nsp); isp++) {

if(dsp[idsp[isp-1]-1] > tiny8) {

err = grackle::impl::dabs(vec[isp-1] / dsp[idsp[isp-1]-1]);

} else {

err = 0.e0;

}

if(err > err_max) {

err_max = err;

}

itr=itr+1;

}

label_9998:

label_9997:

label_9996:

// Check if the fractions are valid after an iteration

if( (my_chemistry->primordial_chemistry > 0) && (my_chemistry->with_radiative_cooling == 1) ) {

for (isp = 1; isp<=(nsp); isp++) {

if ( (dsp[idsp[isp-1]-1] != dsp[idsp[isp-1]-1])

|| (dsp[idsp[isp-1]-1] <= 0.) ) {

ierror = 1;

}

if(ierror == 1) {

dtit[i-1] = 0.5e0*dtit[i-1];

}

itr_time=itr_time+1;

}

if ( my_chemistry->primordial_chemistry > 0 ) {

de(i-1,j-1,k-1) = dsp[ 1-1];

HI(i-1,j-1,k-1) = dsp[ 2-1];

HII(i-1,j-1,k-1) = dsp[ 3-1];

HeI(i-1,j-1,k-1) = dsp[ 4-1];

HeII(i-1,j-1,k-1) = dsp[ 5-1];

HeIII(i-1,j-1,k-1) = dsp[ 6-1];

}

if ( my_chemistry->primordial_chemistry > 1 ) {

HM(i-1,j-1,k-1) = dsp[ 7-1];

H2I(i-1,j-1,k-1) = dsp[ 8-1];

H2II(i-1,j-1,k-1) = dsp[ 9-1];

}

if ( my_chemistry->primordial_chemistry > 2 ) {

DI(i-1,j-1,k-1) = dsp[10-1];

DII(i-1,j-1,k-1) = dsp[11-1];

HDI(i-1,j-1,k-1) = dsp[12-1];

}

if ( my_chemistry->primordial_chemistry > 3 ) {

DM(i-1,j-1,k-1) = dsp[13-1];

HDII(i-1,j-1,k-1) = dsp[14-1];

HeHII(i-1,j-1,k-1) = dsp[15-1];

}

if ( itmask_metal[i-1] != MASK_FALSE ) {

if ( my_chemistry->metal_chemistry == 1 ) {

CI(i-1,j-1,k-1) = dsp[16-1];

CII(i-1,j-1,k-1) = dsp[17-1];

CO(i-1,j-1,k-1) = dsp[18-1];

CO2(i-1,j-1,k-1) = dsp[19-1];

OI(i-1,j-1,k-1) = dsp[20-1];

OH(i-1,j-1,k-1) = dsp[21-1];

H2O(i-1,j-1,k-1) = dsp[22-1];

O2(i-1,j-1,k-1) = dsp[23-1];

SiI(i-1,j-1,k-1) = dsp[24-1];

SiOI(i-1,j-1,k-1) = dsp[25-1];

SiO2I(i-1,j-1,k-1) = dsp[26-1];

CH(i-1,j-1,k-1) = dsp[27-1];

CH2(i-1,j-1,k-1) = dsp[28-1];

COII(i-1,j-1,k-1) = dsp[29-1];

OII(i-1,j-1,k-1) = dsp[30-1];

OHII(i-1,j-1,k-1) = dsp[31-1];

H2OII(i-1,j-1,k-1) = dsp[32-1];

H3OII(i-1,j-1,k-1) = dsp[33-1];

O2II(i-1,j-1,k-1) = dsp[34-1];

if ( ( my_chemistry->grain_growth == 1 ) || ( my_chemistry->dust_sublimation == 1 ) ) {

if (my_chemistry->dust_species > 0) {

Mg(i-1,j-1,k-1) = dsp[35-1];

}

if (my_chemistry->dust_species > 1) {

Al(i-1,j-1,k-1) = dsp[36-1];

S(i-1,j-1,k-1) = dsp[37-1];

Fe(i-1,j-1,k-1) = dsp[38-1];

}

if ( ( my_chemistry->grain_growth == 1 ) || ( my_chemistry->dust_sublimation == 1 ) ) {

if (my_chemistry->dust_species > 0) {

MgSiO3(i-1,j-1,k-1) = dsp[39-1];

AC(i-1,j-1,k-1) = dsp[40-1];

}

if (my_chemistry->dust_species > 1) {

SiM(i-1,j-1,k-1) = dsp[41-1];

FeM(i-1,j-1,k-1) = dsp[42-1];

Mg2SiO4(i-1,j-1,k-1) = dsp[43-1];

Fe3O4(i-1,j-1,k-1) = dsp[44-1];

SiO2D(i-1,j-1,k-1) = dsp[45-1];

MgO(i-1,j-1,k-1) = dsp[46-1];

FeS(i-1,j-1,k-1) = dsp[47-1];

Al2O3(i-1,j-1,k-1) = dsp[48-1];

}

if (my_chemistry->dust_species > 2) {

reforg(i-1,j-1,k-1) = dsp[49-1];

volorg(i-1,j-1,k-1) = dsp[50-1];

H2Oice(i-1,j-1,k-1) = dsp[51-1];

}

e(i-1,j-1,k-1) = dsp[i_eng-1];

idsp.clear();

vec.clear();

mtrx = grackle::impl::View<double**>();

mtrx_data_.clear();

}

// return itmask

for (i = my_fields->grid_start[0] + 1; i<=(my_fields->grid_end[0] + 1); i++) {

itmask[i-1] = itmask_tmp[i-1];

}

// Add the timestep to the elapsed time for each cell and find

// minimum elapsed time step in this row

ttmin = huge8;

for (i = my_fields->grid_start[0] + 1; i<=(my_fields->grid_end[0] + 1); i++) {

ttot[i-1] = std::fmin(ttot[i-1] + dtit[i-1], (*dt));

if (std::fabs((*dt)-ttot[i-1]) <

tolerance*(*dt)) { itmask[i-1] = MASK_FALSE; }

if (ttot[i-1]<ttmin) { ttmin = ttot[i-1]; }

}

// If all cells are done (on this slice), then exit

if (std::fabs((*dt)-ttmin) < tolerance*(*dt)) { goto label_9999; }

// Next subcycle iteration

}

label_9999:

// Abort if iteration count exceeds maximum

if (iter > my_chemistry->max_iterations) {

OMP_PRAGMA_CRITICAL

{

printf("inside if statement solve rate cool: %d %d\n",

my_fields->grid_start[0],

my_fields->grid_end[0]);

eprintf("MULTI_COOL iter > %d at j,k = %d %d\n",

my_chemistry->max_iterations,

j,

k);

printf("FATAL error (2) in MULTI_COOL\n");

//_// PORT: write(0,'(" dt = ",1pe10.3," ttmin = ",1pe10.3)') dt, ttmin

//_// PORT: write(0,'((16(1pe8.1)))') (dtit(i),i=is+1,ie+1)

//_// PORT: write(0,'((16(1pe8.1)))') (ttot(i),i=is+1,ie+1)

//_// PORT: write(0,'((16(1pe8.1)))') (edot(i),i=is+1,ie+1)

//_// PORT: write(0,'((16(I3)))') (itmask(i),i=is+1,ie+1)

if (my_chemistry->exit_after_iterations_exceeded == 1) {

(*ierr) = 0;

}

// WARNING_MESSAGE

}

if (iter > my_chemistry->max_iterations/2) {

OMP_PRAGMA_CRITICAL

{

eprintf("MULTI_COOL iter,j,k = %d %d %d\n", iter, j, k);

}

// Next j,k

}

//_// PORT: #ifdef _OPENMP

//_// PORT: !$omp end parallel do

//_// PORT: #endif

// If an error has been produced, return now.

if ((*ierr) == 0) {

return;

}

// Convert densities back to comoving from proper

if (my_units->comoving_coordinates == 1) {

factor = (gr_float)(std::pow(my_units->a_value,3) );

FORTRAN_NAME(scale_fields_g)(

&my_chemistry->primordial_chemistry, imetal, &my_chemistry->use_dust_density_field, &my_chemistry->metal_chemistry,

&my_chemistry->metal_abundances, &my_chemistry->dust_species, &my_chemistry->multi_metals, &my_chemistry->grain_growth, &my_chemistry->dust_sublimation, &factor,

&my_fields->grid_start[0], &my_fields->grid_end[0], &my_fields->grid_start[1], &my_fields->grid_end[1], &my_fields->grid_start[2], &my_fields->grid_end[2], &my_fields->grid_dimension[0], &my_fields->grid_dimension[1], &my_fields->grid_dimension[2],

d.data(), de.data(), HI.data(), HII.data(), HeI.data(), HeII.data(), HeIII.data(),

HM.data(), H2I.data(), H2II.data(), DI.data(), DII.data(), HDI.data(), DM.data(), HDII.data(), HeHII.data(),

metal.data(), dust.data(),

CI.data(), CII.data(), CO.data(), CO2.data(),

OI.data(), OH.data(), H2O.data(), O2.data(),

SiI.data(), SiOI.data(), SiO2I.data(),

CH.data(), CH2.data(), COII.data(), OII.data(),

OHII.data(), H2OII.data(), H3OII.data(), O2II.data(),

Mg.data(), Al.data(), S.data(), Fe.data(),

SiM.data(), FeM.data(), Mg2SiO4.data(), MgSiO3.data(), Fe3O4.data(),

AC.data(), SiO2D.data(), MgO.data(), FeS.data(), Al2O3.data(),

reforg.data(), volorg.data(), H2Oice.data(),

metal_loc.data(), metal_C13.data(), metal_C20.data(), metal_C25.data(), metal_C30.data(),

metal_F13.data(), metal_F15.data(), metal_F50.data(), metal_F80.data(),

metal_P170.data(), metal_P200.data(), metal_Y19.data());

}

if (my_chemistry->primordial_chemistry > 0) {

// Correct the species to ensure consistency (i.e. type conservation)

#ifdef ABUNDANCE_CORRECTION

FORTRAN_NAME(make_consistent_g)(de.data(), HI.data(), HII.data(), HeI.data(), HeII.data(), HeIII.data(),

HM.data(), H2I.data(), H2II.data(), DI.data(), DII.data(), HDI.data(), metal.data(), dust.data(),

d.data(), &my_fields->grid_start[0], &my_fields->grid_end[0], &my_fields->grid_start[1], &my_fields->grid_end[1], &my_fields->grid_start[2], &my_fields->grid_end[2],

&my_fields->grid_dimension[0], &my_fields->grid_dimension[1], &my_fields->grid_dimension[2], &my_chemistry->primordial_chemistry, imetal, &my_chemistry->HydrogenFractionByMass, &my_chemistry->DeuteriumToHydrogenRatio,

&my_chemistry->use_dust_density_field, &my_chemistry->metal_chemistry, &my_chemistry->grain_growth, &dom,

DM.data(), HDII.data(), HeHII.data(),

CI.data(), CII.data(), CO.data(), CO2.data(),

OI.data(), OH.data(), H2O.data(), O2.data(),

SiI.data(), SiOI.data(), SiO2I.data(),

CH.data(), CH2.data(), COII.data(), OII.data(),

OHII.data(), H2OII.data(), H3OII.data(), O2II.data(),

Mg.data(), Al.data(), S.data(), Fe.data(),

SiM.data(), FeM.data(), Mg2SiO4.data(), MgSiO3.data(), Fe3O4.data(),

AC.data(), SiO2D.data(), MgO.data(), FeS.data(), Al2O3.data(),

reforg.data(), volorg.data(), H2Oice.data(),

&my_chemistry->multi_metals, &my_chemistry->metal_abundances, &my_chemistry->dust_species, &my_chemistry->use_multiple_dust_temperatures, &my_chemistry->dust_sublimation,

metal_loc.data(),

metal_C13.data(), metal_C20.data(), metal_C25.data(), metal_C30.data(),

metal_F13.data(), metal_F15.data(), metal_F50.data(), metal_F80.data(),

metal_P170.data(), metal_P200.data(), metal_Y19.data(),

&my_rates->SN0_N,

my_rates->SN0_XC, my_rates->SN0_XO, my_rates->SN0_XMg, my_rates->SN0_XAl, my_rates->SN0_XSi,

my_rates->SN0_XS, my_rates->SN0_XFe,

my_rates->SN0_fC, my_rates->SN0_fO, my_rates->SN0_fMg, my_rates->SN0_fAl, my_rates->SN0_fSi,

my_rates->SN0_fS, my_rates->SN0_fFe

);

FORTRAN_NAME(ceiling_species_g)(d.data(), de.data(), HI.data(), HII.data(), HeI.data(), HeII.data(), HeIII.data(),

HM.data(), H2I.data(), H2II.data(), DI.data(), DII.data(), HDI.data(), metal.data(), dust.data(),

&my_fields->grid_start[0], &my_fields->grid_end[0], &my_fields->grid_start[1], &my_fields->grid_end[1], &my_fields->grid_start[2], &my_fields->grid_end[2],

&my_fields->grid_dimension[0], &my_fields->grid_dimension[1], &my_fields->grid_dimension[2], &my_chemistry->primordial_chemistry, imetal, &my_chemistry->use_dust_density_field,

DM.data(), HDII.data(), HeHII.data(),

&my_chemistry->metal_abundances, &my_chemistry->metal_chemistry, &my_chemistry->dust_species, &my_chemistry->multi_metals,

&my_chemistry->grain_growth, &my_chemistry->dust_sublimation,

CI.data(), CII.data(), CO.data(), CO2.data(),

OI.data(), OH.data(), H2O.data(), O2.data(),

SiI.data(), SiOI.data(), SiO2I.data(),

CH.data(), CH2.data(), COII.data(), OII.data(),

OHII.data(), H2OII.data(), H3OII.data(), O2II.data(),

Mg.data(), Al.data(), S.data(), Fe.data(),

SiM.data(), FeM.data(), Mg2SiO4.data(), MgSiO3.data(), Fe3O4.data(),

AC.data(), SiO2D.data(), MgO.data(), FeS.data(), Al2O3.data(),

reforg.data(), volorg.data(), H2Oice.data(),

metal_loc.data(),

metal_C13.data(), metal_C20.data(), metal_C25.data(), metal_C30.data(),

metal_F13.data(), metal_F15.data(), metal_F50.data(), metal_F80.data(),

metal_P170.data(), metal_P200.data(), metal_Y19.data());

#endif

}

return;

}

} // extern "C"