API tools faq
paste
mabruzzo

preserve-args solve_rate_cool_g-cpp.C

mabruzzo
Dec 4th, 2024
19
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 103.36 KB | None | 0 0
raw download clone embed print report
  1. #include <cstdio>
  2. #include <vector>
  3.  
  4. #include "grackle.h"
  5. #include "fortran_func_decls.h"
  6. #include "utils-cpp.hpp"
  7. #include "solve_rate_cool_g-cpp.h"
  8.  
  9. //_// TODO: ADD ANY OTHER INCLUDE DIRECTIVES
  10.  
  11. extern "C" {
  12. void solve_rate_cool_g(
  13.   int* icool, gr_float* d_data_, gr_float* e_data_, gr_float* u_data_,
  14.   gr_float* v_data_, gr_float* w_data_, gr_float* de_data_, gr_float* HI_data_,
  15.   gr_float* HII_data_, gr_float* HeI_data_, gr_float* HeII_data_,
  16.   gr_float* HeIII_data_, int* in, int* jn, int* kn, int* nratec, int* iexpand,
  17.   int* ispecies, int* imetal, int* imcool, int* idust, int* idustall,
  18.   int* idustfield, int* idim, int* is, int* js, int* ks, int* ie, int* je,
  19.   int* ke, int* ih2co, int* ipiht, int* idustrec, int* igammah, double* dx,
  20.   double* dt, double* aye, double* temstart, double* temend, double* utem,
  21.   double* uxyz, double* uaye, double* urho, double* utim, double* gamma,
  22.   double* fh, double* dtoh, double* z_solar, double* fgr, double* k1a,
  23.   double* k2a, double* k3a, double* k4a, double* k5a, double* k6a, double* k7a,
  24.   double* k8a, double* k9a, double* k10a, double* k11a, double* k12a,
  25.   double* k13a, double* k13dda, double* k14a, double* k15a, double* k16a,
  26.   double* k17a, double* k18a, double* k19a, double* k22a, double* k24,
  27.   double* k25, double* k26, double* k27, double* k28, double* k29, double* k30,
  28.   double* k31, double* k50a, double* k51a, double* k52a, double* k53a,
  29.   double* k54a, double* k55a, double* k56a, double* k57a, double* k58a,
  30.   int* ndratec, double* dtemstart, double* dtemend, double* h2dusta,
  31.   double* ncrna, double* ncrd1a, double* ncrd2a, double* ceHIa, double* ceHeIa,
  32.   double* ceHeIIa, double* ciHIa, double* ciHeIa, double* ciHeISa,
  33.   double* ciHeIIa, double* reHIIa, double* reHeII1a, double* reHeII2a,
  34.   double* reHeIIIa, double* brema, double* compa, double* gammaha, double* isrf,
  35.   double* regra, double* gamma_isrfa, double* comp_xraya, double* comp_temp,
  36.   double* piHI, double* piHeI, double* piHeII, gr_float* HM_data_,
  37.   gr_float* H2I_data_, gr_float* H2II_data_, gr_float* DI_data_,
  38.   gr_float* DII_data_, gr_float* HDI_data_, gr_float* metal_data_,
  39.   gr_float* dust_data_, double* hyd01ka, double* h2k01a, double* vibha,
  40.   double* rotha, double* rotla, double* gpldla, double* gphdla, double* hdltea,
  41.   double* hdlowa, double* gaHIa, double* gaH2a, double* gaHea, double* gaHpa,
  42.   double* gaela, double* h2ltea, double* gasgra, int* iH2shield,
  43.   int* iradshield, double* avgsighi, double* avgsighei, double* avgsigheii,
  44.   int* iradtrans, int* iradcoupled, int* iradstep, int* irt_honly,
  45.   gr_float* kphHI_data_, gr_float* kphHeI_data_, gr_float* kphHeII_data_,
  46.   gr_float* kdissH2I_data_, gr_float* photogamma_data_,
  47.   gr_float* xH2shield_data_, int* ierr, int* ih2optical, int* iciecool,
  48.   int* ithreebody, int* ih2cr, int* ihdcr, double* ciecoa, int* icmbTfloor,
  49.   int* iClHeat, double* clEleFra, long long* priGridRank, long long* priGridDim,
  50.   double* priPar1, double* priPar2, double* priPar3, double* priPar4,
  51.   double* priPar5, long long* priDataSize, double* priCooling,
  52.   double* priHeating, double* priMMW, long long* metGridRank,
  53.   long long* metGridDim, double* metPar1, double* metPar2, double* metPar3,
  54.   double* metPar4, double* metPar5, long long* metDataSize, double* metCooling,
  55.   double* metHeating, int* clnew, int* iVheat, int* iMheat,
  56.   gr_float* Vheat_data_, gr_float* Mheat_data_, int* iTfloor,
  57.   double* Tfloor_scalar, gr_float* Tfloor_data_, int* imchem, int* igrgr,
  58.   int* ipcont, double* tmcool, gr_float* DM_data_, gr_float* HDII_data_,
  59.   gr_float* HeHII_data_, gr_float* CI_data_, gr_float* CII_data_,
  60.   gr_float* CO_data_, gr_float* CO2_data_, gr_float* OI_data_,
  61.   gr_float* OH_data_, gr_float* H2O_data_, gr_float* O2_data_,
  62.   gr_float* SiI_data_, gr_float* SiOI_data_, gr_float* SiO2I_data_,
  63.   gr_float* CH_data_, gr_float* CH2_data_, gr_float* COII_data_,
  64.   gr_float* OII_data_, gr_float* OHII_data_, gr_float* H2OII_data_,
  65.   gr_float* H3OII_data_, gr_float* O2II_data_, gr_float* Mg_data_,
  66.   gr_float* Al_data_, gr_float* S_data_, gr_float* Fe_data_,
  67.   gr_float* SiM_data_, gr_float* FeM_data_, gr_float* Mg2SiO4_data_,
  68.   gr_float* MgSiO3_data_, gr_float* Fe3O4_data_, gr_float* AC_data_,
  69.   gr_float* SiO2D_data_, gr_float* MgO_data_, gr_float* FeS_data_,
  70.   gr_float* Al2O3_data_, gr_float* reforg_data_, gr_float* volorg_data_,
  71.   gr_float* H2Oice_data_, double* k125a, double* k129a, double* k130a,
  72.   double* k131a, double* k132a, double* k133a, double* k134a, double* k135a,
  73.   double* k136a, double* k137a, double* k148a, double* k149a, double* k150a,
  74.   double* k151a, double* k152a, double* k153a, double* kz15a, double* kz16a,
  75.   double* kz17a, double* kz18a, double* kz19a, double* kz20a, double* kz21a,
  76.   double* kz22a, double* kz23a, double* kz24a, double* kz25a, double* kz26a,
  77.   double* kz27a, double* kz28a, double* kz29a, double* kz30a, double* kz31a,
  78.   double* kz32a, double* kz33a, double* kz34a, double* kz35a, double* kz36a,
  79.   double* kz37a, double* kz38a, double* kz39a, double* kz40a, double* kz41a,
  80.   double* kz42a, double* kz43a, double* kz44a, double* kz45a, double* kz46a,
  81.   double* kz47a, double* kz48a, double* kz49a, double* kz50a, double* kz51a,
  82.   double* kz52a, double* kz53a, double* kz54a, double* cieY06a, int* LH2_N,
  83.   int* LH2_Size, double* LH2_D, double* LH2_T, double* LH2_H, double* LH2_dD,
  84.   double* LH2_dT, double* LH2_dH, double* LH2_L, int* LHD_N, int* LHD_Size,
  85.   double* LHD_D, double* LHD_T, double* LHD_H, double* LHD_dD, double* LHD_dT,
  86.   double* LHD_dH, double* LHD_L, int* LCI_N, int* LCI_Size, double* LCI_D,
  87.   double* LCI_T, double* LCI_H, double* LCI_dD, double* LCI_dT, double* LCI_dH,
  88.   double* LCI_L, int* LCII_N, int* LCII_Size, double* LCII_D, double* LCII_T,
  89.   double* LCII_H, double* LCII_dD, double* LCII_dT, double* LCII_dH,
  90.   double* LCII_L, int* LOI_N, int* LOI_Size, double* LOI_D, double* LOI_T,
  91.   double* LOI_H, double* LOI_dD, double* LOI_dT, double* LOI_dH, double* LOI_L,
  92.   int* LCO_N, int* LCO_Size, double* LCO_D, double* LCO_T, double* LCO_H,
  93.   double* LCO_dD, double* LCO_dT, double* LCO_dH, double* LCO_L, int* LOH_N,
  94.   int* LOH_Size, double* LOH_D, double* LOH_T, double* LOH_H, double* LOH_dD,
  95.   double* LOH_dT, double* LOH_dH, double* LOH_L, int* LH2O_N, int* LH2O_Size,
  96.   double* LH2O_D, double* LH2O_T, double* LH2O_H, double* LH2O_dD,
  97.   double* LH2O_dT, double* LH2O_dH, double* LH2O_L, int* alphap_N,
  98.   int* alphap_Size, double* alphap_D, double* alphap_T, double* alphap_dD,
  99.   double* alphap_dT, double* alphap_Data, int* immulti, int* imabund,
  100.   int* idspecies, int* itdmulti, int* idsub, gr_float* metal_loc_data_,
  101.   gr_float* metal_C13_data_, gr_float* metal_C20_data_,
  102.   gr_float* metal_C25_data_, gr_float* metal_C30_data_,
  103.   gr_float* metal_F13_data_, gr_float* metal_F15_data_,
  104.   gr_float* metal_F50_data_, gr_float* metal_F80_data_,
  105.   gr_float* metal_P170_data_, gr_float* metal_P200_data_,
  106.   gr_float* metal_Y19_data_, int* SN0_N, double* SN0_XC, double* SN0_XO,
  107.   double* SN0_XMg, double* SN0_XAl, double* SN0_XSi, double* SN0_XS,
  108.   double* SN0_XFe, double* SN0_fC, double* SN0_fO, double* SN0_fMg,
  109.   double* SN0_fAl, double* SN0_fSi, double* SN0_fS, double* SN0_fFe,
  110.   double* SN0_fSiM, double* SN0_fFeM, double* SN0_fMg2SiO4, double* SN0_fMgSiO3,
  111.   double* SN0_fFe3O4, double* SN0_fAC, double* SN0_fSiO2D, double* SN0_fMgO,
  112.   double* SN0_fFeS, double* SN0_fAl2O3, double* SN0_freforg,
  113.   double* SN0_fvolorg, double* SN0_fH2Oice, double* SN0_r0SiM,
  114.   double* SN0_r0FeM, double* SN0_r0Mg2SiO4, double* SN0_r0MgSiO3,
  115.   double* SN0_r0Fe3O4, double* SN0_r0AC, double* SN0_r0SiO2D, double* SN0_r0MgO,
  116.   double* SN0_r0FeS, double* SN0_r0Al2O3, double* SN0_r0reforg,
  117.   double* SN0_r0volorg, double* SN0_r0H2Oice, int* gr_N, int* gr_Size,
  118.   double* gr_dT, double* gr_Td, double* SN0_kpSiM, double* SN0_kpFeM,
  119.   double* SN0_kpMg2SiO4, double* SN0_kpMgSiO3, double* SN0_kpFe3O4,
  120.   double* SN0_kpAC, double* SN0_kpSiO2D, double* SN0_kpMgO, double* SN0_kpFeS,
  121.   double* SN0_kpAl2O3, double* SN0_kpreforg, double* SN0_kpvolorg,
  122.   double* SN0_kpH2Oice, double* h2dustSa, double* h2dustCa, double* gasgr2a,
  123.   double* gamma_isrf2a, double* grogra, int* idissHDI, gr_float* kdissHDI_data_,
  124.   int* iionZ, gr_float* kphCI_data_, gr_float* kphOI_data_, int* idissZ,
  125.   gr_float* kdissCO_data_, gr_float* kdissOH_data_, gr_float* kdissH2O_data_,
  126.   int* iuseH2shield, int* iisrffield, gr_float* isrf_habing,
  127.   int* iH2shieldcustom, gr_float* f_shield_custom, int* itmax, int* exititmax
  128. )
  129. {
  130.  
  131.  
  132.   // SOLVE MULTI-SPECIES RATE EQUATIONS AND RADIATIVE COOLING
  133.  
  134.   // written by: Yu Zhang, Peter Anninos and Tom Abel
  135.   // date:
  136.   // modified1:  January, 1996 by Greg Bryan; converted to KRONOS
  137.   // modified2:  October, 1996 by GB; adapted to AMR
  138.   // modified3:  May,     1999 by GB and Tom Abel, 3bodyH2, solver, HD
  139.   // modified4:  June,    2005 by GB to solve rate & cool at same time
  140.   // modified5:  April,   2009 by JHW to include radiative transfer
  141.   // modified6:  September, 2009 by BDS to include cloudy cooling
  142.  
  143.   // PURPOSE:
  144.   //   Solve the multi-species rate and cool equations.
  145.  
  146.   // INPUTS:
  147.   //   icool    - flag to update energy from radiative cooling
  148.   //   in,jn,kn - dimensions of 3D fields
  149.  
  150.   //   d        - total density field
  151.   //   de       - electron density field
  152.   //   HI,HII   - H density fields (neutral & ionized)
  153.   //   HeI/II/III - He density fields
  154.   //   DI/II    - D density fields (neutral & ionized)
  155.   //   HDI      - neutral HD molecule density field
  156.   //   HM       - H- density field
  157.   //   H2I      - H_2 (molecular H) density field
  158.   //   H2II     - H_2+ density field
  159.   //   metal    - metal density field
  160.   //   dust     - dust density field
  161.   //   kph*     - photoionization fields
  162.   //   gamma*   - photoheating fields
  163.   //   f_shield_custom - custom H2 shielding factor
  164.  
  165.   //   is,ie    - start and end indices of active region (zero based)
  166.   //   iexpand  - comoving coordinates flag (0 = off, 1 = on)
  167.   //   idim     - dimensionality (rank) of problem
  168.   //   ispecies - chemistry module (1 - H/He only, 2 - molecular H, 3 - D)
  169.   //   imetal   - flag if metal field is active (0 = no, 1 = yes)
  170.   //   imcool   - flag if there is metal cooling
  171.   //   idust    - flag for H2 formation on dust grains
  172.   //   idustall - flag for dust (0 - none, 1 - heating/cooling + H2 form.)
  173.   //   idustfield - flag if a dust density field is present
  174.   //   iisrffield - flag if a field for the interstellar radiation field is present
  175.   //   ih2co    - flag to include H2 cooling (1 = on, 0 = off)
  176.   //   ipiht    - flag to include photoionization heating (1 = on, 0 = off)
  177.   //   idustrec - flag to include dust recombination cooling (1 = on, -1 = off)
  178.   //   iH2shield - flag for approximate self-shielding of H2 (Wolcott-Green+ 2011)
  179.   //   iradshield - flag for approximate self-shielding of UV background
  180.   //   avgsighi   - spectrum averaged ionization crs for HI for use with shielding
  181.   //   avgsighei  - spectrum averaged ionization crs for HeI for use with shielding
  182.   //   avgsigheii - spectrum averaged ionization crs for HeII for use with shielding
  183.   //   iradtrans - flag to include radiative transfer (1 = on, 0 = off)
  184.   //   iradcoupled - flag to indicate coupled radiative transfer
  185.   //   iradstep  - flag to indicate intermediate coupled radiative transfer timestep
  186.   //   irt_honly - flag to indicate applying RT ionization and heating to HI only
  187.   //   iH2shieldcustom - flag to indicate a custom H2 shielding factor is provided
  188.  
  189.   //     fh       - Hydrogen mass fraction (typically 0.76)
  190.   //     dtoh     - Deuterium to H mass ratio
  191.   //     z_solar  - Solar metal mass fraction
  192.   //     fgr      - the local dust to gas ratio (by mass)
  193.   //     dt       - timestep to integrate over
  194.   //     aye      - expansion factor (in code units)
  195.  
  196.   //     utim     - time units (i.e. code units to CGS conversion factor)
  197.   //     uaye     - expansion factor conversion factor (uaye = 1/(1+zinit))
  198.   //     urho     - density units
  199.   //     uxyz     - length units
  200.   //     utem     - temperature(-like) units
  201.  
  202.   //     temstart, temend - start and end of temperature range for rate table
  203.   //     nratec   - dimensions of chemical rate arrays (functions of temperature)
  204.   //     dtemstart, dtemend - start and end of dust temperature range
  205.   //     ndratec  - extra dimension for H2 formation on dust rate (dust temperature)
  206.  
  207.   //     icmbTfloor - flag to include temperature floor from cmb
  208.   //     iClHeat    - flag to include cloudy heating
  209.   //     priGridRank - rank of cloudy primordial cooling data grid
  210.   //     priGridDim  - array containing dimensions of cloudy primordial data
  211.   //     priPar1, priPar2, priPar3 - arrays containing primordial grid parameter values
  212.   //     priDataSize - total size of flattened 1D primordial cooling data array
  213.   //     priCooling  - primordial cooling data
  214.   //     priHeating  - primordial heating data
  215.   //     priMMW      - primordial mmw data
  216.   //     metGridRank - rank of cloudy metal cooling data grid
  217.   //     metGridDim  - array containing dimensions of cloudy metal data
  218.   //     metPar1, metPar2, metPar3 - arrays containing metal grid parameter values
  219.   //     metDataSize - total size of flattened 1D metal cooling data array
  220.   //     metCooling  - metal cooling data
  221.   //     metHeating  - metal heating data
  222.   //     iVheat      - flag for using volumetric heating rate
  223.   //     iMheat      - flag for using specific heating rate
  224.   //     Vheat       - array of volumetric heating rates
  225.   //     Mheat       - array of specific heating rates
  226.   //     iTfloor     - flag for using temperature floor field
  227.   //     Tfloor_scalar - scalar temperature floor value
  228.   //     Tfloor      - array of temperature floor values
  229.   //     itmax       - maximum allowed sub-cycle iterations
  230.   //     exititmax   - flag to exit if max iterations exceeded
  231.  
  232.   //   OUTPUTS:
  233.   //     update chemical rate densities (HI, HII, etc)
  234.  
  235.   //   PARAMETERS:
  236.   //     mh      - H mass in cgs units
  237.  
  238.   // -----------------------------------------------------------------------
  239.  
  240.  
  241.   // General Arguments
  242.  
  243.   int ierror;
  244.  
  245.   // -- removed line (previously just declared arg types) --
  246.  
  247.   // Density, energy and velocity fields fields
  248.  
  249.   grackle::impl::View<gr_float***> de(de_data_, (*in), (*jn), (*kn));
  250.   grackle::impl::View<gr_float***> HI(HI_data_, (*in), (*jn), (*kn));
  251.   grackle::impl::View<gr_float***> HII(HII_data_, (*in), (*jn), (*kn));
  252.   grackle::impl::View<gr_float***> HeI(HeI_data_, (*in), (*jn), (*kn));
  253.   grackle::impl::View<gr_float***> HeII(HeII_data_, (*in), (*jn), (*kn));
  254.   grackle::impl::View<gr_float***> HeIII(HeIII_data_, (*in), (*jn), (*kn));
  255.   grackle::impl::View<gr_float***> HM(HM_data_, (*in), (*jn), (*kn));
  256.   grackle::impl::View<gr_float***> H2I(H2I_data_, (*in), (*jn), (*kn));
  257.   grackle::impl::View<gr_float***> H2II(H2II_data_, (*in), (*jn), (*kn));
  258.   grackle::impl::View<gr_float***> DI(DI_data_, (*in), (*jn), (*kn));
  259.   grackle::impl::View<gr_float***> DII(DII_data_, (*in), (*jn), (*kn));
  260.   grackle::impl::View<gr_float***> HDI(HDI_data_, (*in), (*jn), (*kn));
  261.   grackle::impl::View<gr_float***> d(d_data_, (*in), (*jn), (*kn));
  262.   grackle::impl::View<gr_float***> e(e_data_, (*in), (*jn), (*kn));
  263.   grackle::impl::View<gr_float***> u(u_data_, (*in), (*jn), (*kn));
  264.   grackle::impl::View<gr_float***> v(v_data_, (*in), (*jn), (*kn));
  265.   grackle::impl::View<gr_float***> w(w_data_, (*in), (*jn), (*kn));
  266.   grackle::impl::View<gr_float***> metal(metal_data_, (*in), (*jn), (*kn));
  267.   grackle::impl::View<gr_float***> dust(dust_data_, (*in), (*jn), (*kn));
  268.   grackle::impl::View<gr_float***> Vheat(Vheat_data_, (*in), (*jn), (*kn));
  269.   grackle::impl::View<gr_float***> Mheat(Mheat_data_, (*in), (*jn), (*kn));
  270.   grackle::impl::View<gr_float***> Tfloor(Tfloor_data_, (*in), (*jn), (*kn));
  271.   grackle::impl::View<gr_float***> DM(DM_data_, (*in), (*jn), (*kn));
  272.   grackle::impl::View<gr_float***> HDII(HDII_data_, (*in), (*jn), (*kn));
  273.   grackle::impl::View<gr_float***> HeHII(HeHII_data_, (*in), (*jn), (*kn));
  274.   grackle::impl::View<gr_float***> CI(CI_data_, (*in), (*jn), (*kn));
  275.   grackle::impl::View<gr_float***> CII(CII_data_, (*in), (*jn), (*kn));
  276.   grackle::impl::View<gr_float***> CO(CO_data_, (*in), (*jn), (*kn));
  277.   grackle::impl::View<gr_float***> CO2(CO2_data_, (*in), (*jn), (*kn));
  278.   grackle::impl::View<gr_float***> OI(OI_data_, (*in), (*jn), (*kn));
  279.   grackle::impl::View<gr_float***> OH(OH_data_, (*in), (*jn), (*kn));
  280.   grackle::impl::View<gr_float***> H2O(H2O_data_, (*in), (*jn), (*kn));
  281.   grackle::impl::View<gr_float***> O2(O2_data_, (*in), (*jn), (*kn));
  282.   grackle::impl::View<gr_float***> SiI(SiI_data_, (*in), (*jn), (*kn));
  283.   grackle::impl::View<gr_float***> SiOI(SiOI_data_, (*in), (*jn), (*kn));
  284.   grackle::impl::View<gr_float***> SiO2I(SiO2I_data_, (*in), (*jn), (*kn));
  285.   grackle::impl::View<gr_float***> CH(CH_data_, (*in), (*jn), (*kn));
  286.   grackle::impl::View<gr_float***> CH2(CH2_data_, (*in), (*jn), (*kn));
  287.   grackle::impl::View<gr_float***> COII(COII_data_, (*in), (*jn), (*kn));
  288.   grackle::impl::View<gr_float***> OII(OII_data_, (*in), (*jn), (*kn));
  289.   grackle::impl::View<gr_float***> OHII(OHII_data_, (*in), (*jn), (*kn));
  290.   grackle::impl::View<gr_float***> H2OII(H2OII_data_, (*in), (*jn), (*kn));
  291.   grackle::impl::View<gr_float***> H3OII(H3OII_data_, (*in), (*jn), (*kn));
  292.   grackle::impl::View<gr_float***> O2II(O2II_data_, (*in), (*jn), (*kn));
  293.   grackle::impl::View<gr_float***> Mg(Mg_data_, (*in), (*jn), (*kn));
  294.   grackle::impl::View<gr_float***> Al(Al_data_, (*in), (*jn), (*kn));
  295.   grackle::impl::View<gr_float***> S(S_data_, (*in), (*jn), (*kn));
  296.   grackle::impl::View<gr_float***> Fe(Fe_data_, (*in), (*jn), (*kn));
  297.   grackle::impl::View<gr_float***> SiM(SiM_data_, (*in), (*jn), (*kn));
  298.   grackle::impl::View<gr_float***> FeM(FeM_data_, (*in), (*jn), (*kn));
  299.   grackle::impl::View<gr_float***> Mg2SiO4(Mg2SiO4_data_, (*in), (*jn), (*kn));
  300.   grackle::impl::View<gr_float***> MgSiO3(MgSiO3_data_, (*in), (*jn), (*kn));
  301.   grackle::impl::View<gr_float***> Fe3O4(Fe3O4_data_, (*in), (*jn), (*kn));
  302.   grackle::impl::View<gr_float***> AC(AC_data_, (*in), (*jn), (*kn));
  303.   grackle::impl::View<gr_float***> SiO2D(SiO2D_data_, (*in), (*jn), (*kn));
  304.   grackle::impl::View<gr_float***> MgO(MgO_data_, (*in), (*jn), (*kn));
  305.   grackle::impl::View<gr_float***> FeS(FeS_data_, (*in), (*jn), (*kn));
  306.   grackle::impl::View<gr_float***> Al2O3(Al2O3_data_, (*in), (*jn), (*kn));
  307.   grackle::impl::View<gr_float***> reforg(reforg_data_, (*in), (*jn), (*kn));
  308.   grackle::impl::View<gr_float***> volorg(volorg_data_, (*in), (*jn), (*kn));
  309.   grackle::impl::View<gr_float***> H2Oice(H2Oice_data_, (*in), (*jn), (*kn));
  310.   grackle::impl::View<gr_float***> metal_loc(metal_loc_data_, (*in), (*jn), (*kn));
  311.   grackle::impl::View<gr_float***> metal_C13(metal_C13_data_, (*in), (*jn), (*kn));
  312.   grackle::impl::View<gr_float***> metal_C20(metal_C20_data_, (*in), (*jn), (*kn));
  313.   grackle::impl::View<gr_float***> metal_C25(metal_C25_data_, (*in), (*jn), (*kn));
  314.   grackle::impl::View<gr_float***> metal_C30(metal_C30_data_, (*in), (*jn), (*kn));
  315.   grackle::impl::View<gr_float***> metal_F13(metal_F13_data_, (*in), (*jn), (*kn));
  316.   grackle::impl::View<gr_float***> metal_F15(metal_F15_data_, (*in), (*jn), (*kn));
  317.   grackle::impl::View<gr_float***> metal_F50(metal_F50_data_, (*in), (*jn), (*kn));
  318.   grackle::impl::View<gr_float***> metal_F80(metal_F80_data_, (*in), (*jn), (*kn));
  319.   grackle::impl::View<gr_float***> metal_P170(metal_P170_data_, (*in), (*jn), (*kn));
  320.   grackle::impl::View<gr_float***> metal_P200(metal_P200_data_, (*in), (*jn), (*kn));
  321.   grackle::impl::View<gr_float***> metal_Y19(metal_Y19_data_, (*in), (*jn), (*kn));
  322.  
  323.   // Radiative transfer fields
  324.  
  325.   grackle::impl::View<gr_float***> kphHI(kphHI_data_, (*in), (*jn), (*kn));
  326.   grackle::impl::View<gr_float***> kphHeI(kphHeI_data_, (*in), (*jn), (*kn));
  327.   grackle::impl::View<gr_float***> kphHeII(kphHeII_data_, (*in), (*jn), (*kn));
  328.   grackle::impl::View<gr_float***> kdissH2I(kdissH2I_data_, (*in), (*jn), (*kn));
  329.   grackle::impl::View<gr_float***> photogamma(photogamma_data_, (*in), (*jn), (*kn));
  330.  
  331.   // -- removed line (previously just declared arg types) --
  332.   grackle::impl::View<gr_float***> kdissHDI(kdissHDI_data_, (*in), (*jn), (*kn));
  333.   grackle::impl::View<gr_float***> kphCI(kphCI_data_, (*in), (*jn), (*kn));
  334.   grackle::impl::View<gr_float***> kphOI(kphOI_data_, (*in), (*jn), (*kn));
  335.   grackle::impl::View<gr_float***> kdissCO(kdissCO_data_, (*in), (*jn), (*kn));
  336.   grackle::impl::View<gr_float***> kdissOH(kdissOH_data_, (*in), (*jn), (*kn));
  337.   grackle::impl::View<gr_float***> kdissH2O(kdissH2O_data_, (*in), (*jn), (*kn));
  338.   // -- removed line (previously just declared arg types) --
  339.  
  340.   // H2 self-shielding length-scale field
  341.  
  342.   grackle::impl::View<gr_float***> xH2shield(xH2shield_data_, (*in), (*jn), (*kn));
  343.  
  344.   // Interstellar radiation field for dust heating
  345.  
  346.  
  347.   // Custom H2 shielding factor
  348.  
  349.  
  350.   // Cooling tables (coolings rates as a function of temperature)
  351.  
  352.   // -- removed line (previously just declared arg types) --
  353.   // -- removed line (previously just declared arg types) --
  354.   // -- removed line (previously just declared arg types) --
  355.   // -- removed line (previously just declared arg types) --
  356.   // -- removed line (previously just declared arg types) --
  357.   // -- removed line (previously just declared arg types) --
  358.   // -- removed line (previously just declared arg types) --
  359.   // -- removed line (previously just declared arg types) --
  360.   // -- removed line (previously just declared arg types) --
  361.   // -- removed line (previously just declared arg types) --
  362.   // -- removed line (previously just declared arg types) --
  363.   // -- removed line (previously just declared arg types) --
  364.   // -- removed line (previously just declared arg types) --
  365.   // -- removed line (previously just declared arg types) --
  366.   // -- removed line (previously just declared arg types) --
  367.   // -- removed line (previously just declared arg types) --
  368.   // -- removed line (previously just declared arg types) --
  369.   // -- removed line (previously just declared arg types) --
  370.   // -- removed line (previously just declared arg types) --
  371.   // -- removed line (previously just declared arg types) --
  372.   // -- removed line (previously just declared arg types) --
  373.   // -- removed line (previously just declared arg types) --
  374.   // -- removed line (previously just declared arg types) --
  375.   // -- removed line (previously just declared arg types) --
  376.   // -- removed line (previously just declared arg types) --
  377.   // opacity table
  378.   // -- removed line (previously just declared arg types) --
  379.   // -- removed line (previously just declared arg types) --
  380.   // -- removed line (previously just declared arg types) --
  381.  
  382.   // -- removed line (previously just declared arg types) --
  383.  
  384.   // Chemistry tables (rates as a function of temperature)
  385.  
  386.   // -- removed line (previously just declared arg types) --
  387.  
  388.   // Cloudy cooling data
  389.  
  390.   // -- removed line (previously just declared arg types) --
  391.   // -- removed line (previously just declared arg types) --
  392.   // -- removed line (previously just declared arg types) --
  393.  
  394.   // Parameters
  395.  
  396. #ifdef GRACKLE_FLOAT_4
  397.   const gr_float tolerance = (gr_float)(1.0e-05);
  398. #else
  399.   const gr_float tolerance = (gr_float)(1.0e-10);
  400. #endif
  401.  
  402.   const double mh = mass_h;
  403.   const double pi = pi_val;
  404.  
  405.   // Locals
  406.  
  407.   int i, j, k, iter;
  408.   int t, dj, dk;
  409.   double ttmin, dom, energy, comp1, comp2;
  410.   double coolunit, dbase1, tbase1, xbase1, chunit, uvel;
  411.   double heq1, heq2, eqk221, eqk222, eqk131, eqk132, eqt1, eqt2, eqtdef, dheq, heq, dlogtem, dx_cgs, c_ljeans, min_metallicity;
  412.   gr_float factor;
  413.  
  414.   // row temporaries
  415.  
  416.   std::vector<long long> indixe((*in));
  417.   double olddtit;
  418.   std::vector<double> t1((*in));
  419.   std::vector<double> t2((*in));
  420.   std::vector<double> logtem((*in));
  421.   std::vector<double> tdef((*in));
  422.   std::vector<double> dtit((*in));
  423.   std::vector<double> ttot((*in));
  424.   std::vector<double> p2d((*in));
  425.   std::vector<double> tgas((*in));
  426.   std::vector<double> tgasold((*in));
  427.   std::vector<double> tdust((*in));
  428.   std::vector<double> metallicity((*in));
  429.   std::vector<double> dust2gas((*in));
  430.   std::vector<double> rhoH((*in));
  431.   std::vector<double> mmw((*in));
  432.   std::vector<double> mynh((*in));
  433.   std::vector<double> myde((*in));
  434.   std::vector<double> gammaha_eff((*in));
  435.   std::vector<double> gasgr_tdust((*in));
  436.   std::vector<double> regr((*in));
  437.   std::vector<double> ddom((*in));
  438.  
  439.   // Rate equation row temporaries
  440.  
  441.   std::vector<double> HIp((*in));
  442.   std::vector<double> HIIp((*in));
  443.   std::vector<double> HeIp((*in));
  444.   std::vector<double> HeIIp((*in));
  445.   std::vector<double> HeIIIp((*in));
  446.   std::vector<double> HMp((*in));
  447.   std::vector<double> H2Ip((*in));
  448.   std::vector<double> H2IIp((*in));
  449.   std::vector<double> dep((*in));
  450.   std::vector<double> dedot((*in));
  451.   std::vector<double> HIdot((*in));
  452.   std::vector<double> dedot_prev((*in));
  453.   std::vector<double> DIp((*in));
  454.   std::vector<double> DIIp((*in));
  455.   std::vector<double> HDIp((*in));
  456.   std::vector<double> HIdot_prev((*in));
  457.   std::vector<double> k24shield((*in));
  458.   std::vector<double> k25shield((*in));
  459.   std::vector<double> k26shield((*in));
  460.   std::vector<double> k28shield((*in));
  461.   std::vector<double> k29shield((*in));
  462.   std::vector<double> k30shield((*in));
  463.   std::vector<double> k31shield((*in));
  464.   std::vector<double> k1((*in));
  465.   std::vector<double> k2((*in));
  466.   std::vector<double> k3((*in));
  467.   std::vector<double> k4((*in));
  468.   std::vector<double> k5((*in));
  469.   std::vector<double> k6((*in));
  470.   std::vector<double> k7((*in));
  471.   std::vector<double> k8((*in));
  472.   std::vector<double> k9((*in));
  473.   std::vector<double> k10((*in));
  474.   std::vector<double> k11((*in));
  475.   std::vector<double> k12((*in));
  476.   std::vector<double> k13((*in));
  477.   std::vector<double> k14((*in));
  478.   std::vector<double> k15((*in));
  479.   std::vector<double> k16((*in));
  480.   std::vector<double> k17((*in));
  481.   std::vector<double> k18((*in));
  482.   std::vector<double> k19((*in));
  483.   std::vector<double> k22((*in));
  484.   std::vector<double> k50((*in));
  485.   std::vector<double> k51((*in));
  486.   std::vector<double> k52((*in));
  487.   std::vector<double> k53((*in));
  488.   std::vector<double> k54((*in));
  489.   std::vector<double> k55((*in));
  490.   std::vector<double> k56((*in));
  491.   std::vector<double> k57((*in));
  492.   std::vector<double> k58((*in));
  493.   std::vector<double> k13dd((*in) * 14);
  494.   std::vector<double> h2dust((*in));
  495.   std::vector<double> ncrn((*in));
  496.   std::vector<double> ncrd1((*in));
  497.   std::vector<double> ncrd2((*in));
  498.   std::vector<double> DMp((*in));
  499.   std::vector<double> HDIIp((*in));
  500.   std::vector<double> HeHIIp((*in));
  501.   std::vector<double> CIp((*in));
  502.   std::vector<double> CIIp((*in));
  503.   std::vector<double> COp((*in));
  504.   std::vector<double> CO2p((*in));
  505.   std::vector<double> OIp((*in));
  506.   std::vector<double> OHp((*in));
  507.   std::vector<double> H2Op((*in));
  508.   std::vector<double> O2p((*in));
  509.   std::vector<double> SiIp((*in));
  510.   std::vector<double> SiOIp((*in));
  511.   std::vector<double> SiO2Ip((*in));
  512.   std::vector<double> CHp((*in));
  513.   std::vector<double> CH2p((*in));
  514.   std::vector<double> COIIp((*in));
  515.   std::vector<double> OIIp((*in));
  516.   std::vector<double> OHIIp((*in));
  517.   std::vector<double> H2OIIp((*in));
  518.   std::vector<double> H3OIIp((*in));
  519.   std::vector<double> O2IIp((*in));
  520.   std::vector<double> Mgp((*in));
  521.   std::vector<double> Alp((*in));
  522.   std::vector<double> Sp((*in));
  523.   std::vector<double> Fep((*in));
  524.   std::vector<gr_float> SiMp((*in));
  525.   std::vector<gr_float> FeMp((*in));
  526.   std::vector<gr_float> Mg2SiO4p((*in));
  527.   std::vector<gr_float> MgSiO3p((*in));
  528.   std::vector<gr_float> Fe3O4p((*in));
  529.   std::vector<gr_float> ACp((*in));
  530.   std::vector<gr_float> SiO2Dp((*in));
  531.   std::vector<gr_float> MgOp((*in));
  532.   std::vector<gr_float> FeSp((*in));
  533.   std::vector<gr_float> Al2O3p((*in));
  534.   std::vector<gr_float> reforgp((*in));
  535.   std::vector<gr_float> volorgp((*in));
  536.   std::vector<gr_float> H2Oicep((*in));
  537.  
  538.   std::vector<double> k125((*in));
  539.   std::vector<double> k129((*in));
  540.   std::vector<double> k130((*in));
  541.   std::vector<double> k131((*in));
  542.   std::vector<double> k132((*in));
  543.   std::vector<double> k133((*in));
  544.   std::vector<double> k134((*in));
  545.   std::vector<double> k135((*in));
  546.   std::vector<double> k136((*in));
  547.   std::vector<double> k137((*in));
  548.   std::vector<double> k148((*in));
  549.   std::vector<double> k149((*in));
  550.   std::vector<double> k150((*in));
  551.   std::vector<double> k151((*in));
  552.   std::vector<double> k152((*in));
  553.   std::vector<double> k153((*in));
  554.   std::vector<double> kz15((*in));
  555.   std::vector<double> kz16((*in));
  556.   std::vector<double> kz17((*in));
  557.   std::vector<double> kz18((*in));
  558.   std::vector<double> kz19((*in));
  559.   std::vector<double> kz20((*in));
  560.   std::vector<double> kz21((*in));
  561.   std::vector<double> kz22((*in));
  562.   std::vector<double> kz23((*in));
  563.   std::vector<double> kz24((*in));
  564.   std::vector<double> kz25((*in));
  565.   std::vector<double> kz26((*in));
  566.   std::vector<double> kz27((*in));
  567.   std::vector<double> kz28((*in));
  568.   std::vector<double> kz29((*in));
  569.   std::vector<double> kz30((*in));
  570.   std::vector<double> kz31((*in));
  571.   std::vector<double> kz32((*in));
  572.   std::vector<double> kz33((*in));
  573.   std::vector<double> kz34((*in));
  574.   std::vector<double> kz35((*in));
  575.   std::vector<double> kz36((*in));
  576.   std::vector<double> kz37((*in));
  577.   std::vector<double> kz38((*in));
  578.   std::vector<double> kz39((*in));
  579.   std::vector<double> kz40((*in));
  580.   std::vector<double> kz41((*in));
  581.   std::vector<double> kz42((*in));
  582.   std::vector<double> kz43((*in));
  583.   std::vector<double> kz44((*in));
  584.   std::vector<double> kz45((*in));
  585.   std::vector<double> kz46((*in));
  586.   std::vector<double> kz47((*in));
  587.   std::vector<double> kz48((*in));
  588.   std::vector<double> kz49((*in));
  589.   std::vector<double> kz50((*in));
  590.   std::vector<double> kz51((*in));
  591.   std::vector<double> kz52((*in));
  592.   std::vector<double> kz53((*in));
  593.   std::vector<double> kz54((*in));
  594.   // -- removed line (previously just declared arg types) --
  595.   std::vector<double> kdSiM((*in));
  596.   std::vector<double> kdFeM((*in));
  597.   std::vector<double> kdMg2SiO4((*in));
  598.   std::vector<double> kdMgSiO3((*in));
  599.   std::vector<double> kdFe3O4((*in));
  600.   std::vector<double> kdAC((*in));
  601.   std::vector<double> kdSiO2D((*in));
  602.   std::vector<double> kdMgO((*in));
  603.   std::vector<double> kdFeS((*in));
  604.   std::vector<double> kdAl2O3((*in));
  605.   std::vector<double> kdreforg((*in));
  606.   std::vector<double> kdvolorg((*in));
  607.   std::vector<double> kdH2Oice((*in));
  608.   // grain temperature
  609.   std::vector<double> tSiM((*in));
  610.   std::vector<double> tFeM((*in));
  611.   std::vector<double> tMg2SiO4((*in));
  612.   std::vector<double> tMgSiO3((*in));
  613.   std::vector<double> tFe3O4((*in));
  614.   std::vector<double> tAC((*in));
  615.   std::vector<double> tSiO2D((*in));
  616.   std::vector<double> tMgO((*in));
  617.   std::vector<double> tFeS((*in));
  618.   std::vector<double> tAl2O3((*in));
  619.   std::vector<double> treforg((*in));
  620.   std::vector<double> tvolorg((*in));
  621.   std::vector<double> tH2Oice((*in));
  622.  
  623.   // Cooling/heating row locals
  624.  
  625.   std::vector<double> ceHI((*in));
  626.   std::vector<double> ceHeI((*in));
  627.   std::vector<double> ceHeII((*in));
  628.   std::vector<double> ciHI((*in));
  629.   std::vector<double> ciHeI((*in));
  630.   std::vector<double> ciHeIS((*in));
  631.   std::vector<double> ciHeII((*in));
  632.   std::vector<double> reHII((*in));
  633.   std::vector<double> reHeII1((*in));
  634.   std::vector<double> reHeII2((*in));
  635.   std::vector<double> reHeIII((*in));
  636.   std::vector<double> brem((*in));
  637.   std::vector<double> edot((*in));
  638.   std::vector<double> hyd01k((*in));
  639.   std::vector<double> h2k01((*in));
  640.   std::vector<double> vibh((*in));
  641.   std::vector<double> roth((*in));
  642.   std::vector<double> rotl((*in));
  643.   std::vector<double> gpldl((*in));
  644.   std::vector<double> gphdl((*in));
  645.   std::vector<double> hdlte((*in));
  646.   std::vector<double> hdlow((*in));
  647.   std::vector<double> cieco((*in));
  648.  
  649.   // Iteration mask
  650.  
  651.   gr_mask_type anydust;
  652.   std::vector<gr_mask_type> itmask((*in));
  653.   std::vector<gr_mask_type> itmask_tmp((*in));
  654.   std::vector<gr_mask_type> itmask_nr((*in));
  655.   std::vector<gr_mask_type> itmask_metal((*in));
  656.   int itr, itr_time;
  657.   std::vector<int> imp_eng((*in));
  658.   int nsp, isp, jsp, id;
  659.   double dspj, err, err_max;
  660.   const int i_eng = 52;
  661.   std::vector<double> dsp(i_eng);
  662.   std::vector<double> dsp0(i_eng);
  663.   std::vector<double> dsp1(i_eng);
  664.   std::vector<double> dspdot(i_eng);
  665.   std::vector<double> dspdot1(i_eng);
  666.   std::vector<double> ddsp(i_eng);
  667.   std::vector<double> der_data_(i_eng * i_eng);
  668.   grackle::impl::View<double**> der(der_data_.data(), i_eng, i_eng);
  669.   std::vector<int> idsp;
  670.   std::vector<double> mtrx_data_;
  671.   grackle::impl::View<double**> mtrx;
  672.   std::vector<double> vec;
  673.   const double eps = 1.e-4;
  674.  
  675.  
  676.   // \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\/////////////////////////////////
  677.   // =======================================================================
  678.  
  679.   // Set error indicator
  680.  
  681.   (*ierr) = 1;
  682.  
  683.   // Set flag for dust-related options
  684.  
  685.   if (((*idust) > 0)  ||  ((*idustall) > 0))  {
  686.     anydust = MASK_TRUE;
  687.   } else {
  688.     anydust = MASK_FALSE;
  689.   }
  690.  
  691.   // ignore metal chemistry/cooling below this metallicity
  692.   min_metallicity = 1.e-9 / (*z_solar);
  693.      
  694.   // Set units
  695.  
  696.   dom      = (*urho)*(std::pow((*aye),3))/mh;
  697.   tbase1   = (*utim);
  698.   xbase1   = (*uxyz)/((*aye)*(*uaye));    // uxyz is [x]*a      = [x]*[a]*a'        '
  699.   dbase1   = (*urho)*std::pow(((*aye)*(*uaye)),3); // urho is [dens]/a^3 = [dens]/([a]*a')^3 '
  700.   coolunit = (std::pow((*uaye),5) * std::pow(xbase1,2) * std::pow(mh,2)) / (std::pow(tbase1,3) * dbase1);
  701.   uvel = ((*uxyz)/(*aye)) / (*utim);
  702.   // chunit = (1.60218e-12_DKIND)/(2._DKIND*uvel*uvel*mh)   ! 1 eV per H2 formed
  703.   chunit = (1.60218e-12)/(uvel*uvel*mh);            // 1 eV per REACTION (Feb 2020, Gen Chiaki)
  704.  
  705.   dx_cgs = (*dx) * xbase1;
  706.   c_ljeans = std::sqrt(((*gamma) * pi * kboltz) /
  707.        (GravConst * mh * dbase1));
  708.  
  709.   dlogtem = (std::log((*temend)) - std::log((*temstart)))/(double)((*nratec)-1 );
  710.  
  711.   // We better make consistent at first GC202002
  712.  
  713.   if ((*ispecies) > 0)  {
  714.  
  715. #define ABUNDANCE_CORRECTION
  716. #ifdef ABUNDANCE_CORRECTION
  717.      FORTRAN_NAME(make_consistent_g)(de.data(), HI.data(), HII.data(), HeI.data(), HeII.data(), HeIII.data(),
  718.                          HM.data(), H2I.data(), H2II.data(), DI.data(), DII.data(), HDI.data(), metal.data(), dust.data(),
  719.                          d.data(), is, ie, js, je, ks, ke,
  720.                          in, jn, kn, ispecies, imetal, fh, dtoh,
  721.                         idustfield, imchem, igrgr, &dom,
  722.                         DM.data(), HDII.data(), HeHII.data(),
  723.                         CI.data(), CII.data(), CO.data(), CO2.data(),
  724.                         OI.data(), OH.data(), H2O.data(), O2.data(),
  725.                         SiI.data(), SiOI.data(), SiO2I.data(),
  726.                         CH.data(), CH2.data(), COII.data(), OII.data(),
  727.                         OHII.data(), H2OII.data(), H3OII.data(), O2II.data(),
  728.                         Mg.data(), Al.data(), S.data(), Fe.data(),
  729.                         SiM.data(), FeM.data(), Mg2SiO4.data(), MgSiO3.data(), Fe3O4.data(),
  730.                         AC.data(), SiO2D.data(), MgO.data(), FeS.data(), Al2O3.data(),
  731.                         reforg.data(), volorg.data(), H2Oice.data(),
  732.                         immulti, imabund, idspecies, itdmulti, idsub,
  733.                         metal_loc.data(),
  734.                         metal_C13.data(), metal_C20.data(), metal_C25.data(), metal_C30.data(),
  735.                         metal_F13.data(), metal_F15.data(), metal_F50.data(), metal_F80.data(),
  736.                         metal_P170.data(), metal_P200.data(), metal_Y19.data(),
  737.                         SN0_N,
  738.                         SN0_XC, SN0_XO, SN0_XMg, SN0_XAl, SN0_XSi,
  739.                         SN0_XS, SN0_XFe,
  740.                         SN0_fC, SN0_fO, SN0_fMg, SN0_fAl, SN0_fSi,
  741.                         SN0_fS, SN0_fFe
  742.                           );
  743. #endif
  744.  
  745.   }
  746.  
  747.   // Convert densities from comoving to proper
  748.  
  749.   if ((*iexpand) == 1)  {
  750.  
  751.     factor = (gr_float)(std::pow((*aye),(-3)) );
  752.  
  753.      FORTRAN_NAME(scale_fields_g)(
  754.       ispecies, imetal, idustfield, imchem,
  755.       imabund, idspecies, immulti, igrgr, idsub, &factor,
  756.       is, ie, js, je, ks, ke, in, jn, kn,
  757.       d.data(), de.data(), HI.data(), HII.data(), HeI.data(), HeII.data(), HeIII.data(),
  758.       HM.data(), H2I.data(), H2II.data(), DI.data(), DII.data(), HDI.data(), DM.data(), HDII.data(), HeHII.data(),
  759.       metal.data(), dust.data(),
  760.       CI.data(), CII.data(), CO.data(), CO2.data(),
  761.       OI.data(), OH.data(), H2O.data(), O2.data(),
  762.       SiI.data(), SiOI.data(), SiO2I.data(),
  763.       CH.data(), CH2.data(), COII.data(), OII.data(),
  764.       OHII.data(), H2OII.data(), H3OII.data(), O2II.data(),
  765.       Mg.data(), Al.data(), S.data(), Fe.data(),
  766.       SiM.data(), FeM.data(), Mg2SiO4.data(), MgSiO3.data(), Fe3O4.data(),
  767.       AC.data(), SiO2D.data(), MgO.data(), FeS.data(), Al2O3.data(),
  768.       reforg.data(), volorg.data(), H2Oice.data(),
  769.       metal_loc.data(), metal_C13.data(), metal_C20.data(), metal_C25.data(), metal_C30.data(),
  770.       metal_F13.data(), metal_F15.data(), metal_F50.data(), metal_F80.data(),
  771.       metal_P170.data(), metal_P200.data(), metal_Y19.data());
  772.  
  773.   }
  774.  
  775. #ifdef ABUNDANCE_CORRECTION
  776.    FORTRAN_NAME(ceiling_species_g)(d.data(), de.data(), HI.data(), HII.data(), HeI.data(), HeII.data(), HeIII.data(),
  777.                        HM.data(), H2I.data(), H2II.data(), DI.data(), DII.data(), HDI.data(), metal.data(), dust.data(),
  778.                        is, ie, js, je, ks, ke,
  779.                        in, jn, kn, ispecies, imetal, idustfield,
  780.                       DM.data(), HDII.data(), HeHII.data(),
  781.                       imabund, imchem, idspecies, immulti,
  782.                       igrgr, idsub,
  783.                       CI.data(), CII.data(), CO.data(), CO2.data(),
  784.                       OI.data(), OH.data(), H2O.data(), O2.data(),
  785.                       SiI.data(), SiOI.data(), SiO2I.data(),
  786.                       CH.data(), CH2.data(), COII.data(), OII.data(),
  787.                       OHII.data(), H2OII.data(), H3OII.data(), O2II.data(),
  788.                       Mg.data(), Al.data(), S.data(), Fe.data(),
  789.                       SiM.data(), FeM.data(), Mg2SiO4.data(), MgSiO3.data(), Fe3O4.data(),
  790.                       AC.data(), SiO2D.data(), MgO.data(), FeS.data(), Al2O3.data(),
  791.                       reforg.data(), volorg.data(), H2Oice.data(),
  792.                       metal_loc.data(),
  793.                       metal_C13.data(), metal_C20.data(), metal_C25.data(), metal_C30.data(),
  794.                       metal_F13.data(), metal_F15.data(), metal_F50.data(), metal_F80.data(),
  795.                       metal_P170.data(), metal_P200.data(), metal_Y19.data());
  796. #endif
  797.  
  798.  
  799.   // Loop over zones, and do an entire i-column in one go
  800.   dk = (*ke) - (*ks) + 1;
  801.   dj = (*je) - (*js) + 1;
  802.  
  803.   // parallelize the k and j loops with OpenMP
  804.   // flat j and k loops for better parallelism
  805.   //_// PORT: #ifdef _OPENMP
  806.   //_// PORT: ! ierr is declared as shared and should be modified with atomic operation
  807.   //_// PORT: !$omp parallel do schedule(runtime) private(
  808.   //_// PORT: !$omp&   i, j, k, iter,
  809.   //_// PORT: !$omp&   ttmin, energy, comp1, comp2,
  810.   //_// PORT: !$omp&   heq1, heq2, eqk221, eqk222, eqk131, eqk132,
  811.   //_// PORT: !$omp&   eqt1, eqt2, eqtdef, dheq, heq,
  812.   //_// PORT: !$omp&   indixe,
  813.   //_// PORT: !$omp&   t1, t2, logtem, tdef,
  814.   //_// PORT: !$omp&   dtit, ttot, p2d, tgas, tgasold,
  815.   //_// PORT: !$omp&   tdust, metallicity, dust2gas, rhoH, mmw,
  816.   //_// PORT: !$omp&   mynh, myde, gammaha_eff, gasgr_tdust, regr, ddom,
  817.   //_// PORT: !$omp&   olddtit,
  818.   //_// PORT: !$omp&   HIp, HIIp, HeIp, HeIIp, HeIIIp,
  819.   //_// PORT: !$omp&   HMp, H2Ip, H2IIp,
  820.   //_// PORT: !$omp&   dep, dedot,HIdot, dedot_prev,
  821.   //_// PORT: !$omp&   DIp, DIIp, HDIp, HIdot_prev,
  822.   //_// PORT: !$omp&   k24shield, k25shield, k26shield,
  823.   //_// PORT: !$omp&   k28shield, k29shield, k30shield,
  824.   //_// PORT: !$omp&   k31shield,
  825.   //_// PORT: !$omp&   k1 , k2 , k3 , k4 , k5,
  826.   //_// PORT: !$omp&   k6 , k7 , k8 , k9 , k10,
  827.   //_// PORT: !$omp&   k11, k12, k13, k14, k15,
  828.   //_// PORT: !$omp&   k16, k17, k18, k19, k22,
  829.   //_// PORT: !$omp&   k50, k51, k52, k53, k54,
  830.   //_// PORT: !$omp&   k55, k56, k57, k58, k13dd, h2dust,
  831.   //_// PORT: !$omp&   ncrn, ncrd1, ncrd2,
  832.   //_// PORT: !$omp&   ceHI, ceHeI, ceHeII,
  833.   //_// PORT: !$omp&   ciHI, ciHeI, ciHeIS, ciHeII,
  834.   //_// PORT: !$omp&   reHII, reHeII1, reHeII2, reHeIII,
  835.   //_// PORT: !$omp&   brem, edot,
  836.   //_// PORT: !$omp&   hyd01k, h2k01, vibh, roth, rotl,
  837.   //_// PORT: !$omp&   gpldl, gphdl, hdlte, hdlow, cieco,
  838.   //_// PORT: !$omp&   itmask, itmask_metal )
  839.   //_// PORT: #endif
  840.   for (t = 0; t<=(dk * dj - 1); t++) {
  841.     k = t/dj      + (*ks)+1;
  842.     j = grackle::impl::mod(t,dj) + (*js)+1;
  843.  
  844.     // tolerance = 1.0e-06_DKIND * dt
  845.  
  846.     // Initialize iteration mask to true for all cells.
  847.  
  848.     for (i = (*is) + 1; i<=((*ie) + 1); i++) {
  849.       itmask[i-1] = MASK_TRUE;
  850.     }
  851.  
  852.     // If we are using coupled radiation with intermediate stepping,
  853.     // set iteration mask to include only cells with radiation in the
  854.     // intermediate coupled chemistry / energy step
  855.     if ((*iradtrans) == 1)  {
  856.       if ((*iradcoupled) == 1  &&  (*iradstep) == 1)  {
  857.         for (i = (*is) + 1; i<=((*ie) + 1); i++) {
  858.           if (kphHI(i-1,j-1,k-1) > 0)  {
  859.             itmask[i-1] = MASK_TRUE;
  860.           } else {
  861.             itmask[i-1] = MASK_FALSE;
  862.           }
  863.         }
  864.       }
  865.  
  866.       // Normal rate solver, but don't double count cells with radiation
  867.       if ((*iradcoupled) == 1  &&  (*iradstep) == 0)  {
  868.         for (i = (*is) + 1; i<=((*ie) + 1); i++) {
  869.           if (kphHI(i-1,j-1,k-1) > 0)  {
  870.             itmask[i-1] = MASK_FALSE;
  871.           } else {
  872.             itmask[i-1] = MASK_TRUE;
  873.           }
  874.         }
  875.       }
  876.     }
  877.  
  878.     // Set time elapsed to zero for each cell in 1D section
  879.  
  880.     for (i = (*is) + 1; i<=((*ie) + 1); i++) {
  881.       ttot[i-1] = 0.;
  882.     }
  883.  
  884.     // A useful slice variable since we do this a lot
  885.  
  886.     for (i = (*is) + 1; i<=((*ie) + 1); i++) {
  887.       ddom[i-1] = d(i-1,j-1,k-1) * dom;
  888.     }
  889.  
  890.     // ------------------ Loop over subcycles ----------------
  891.  
  892.     for (iter = 1; iter<=((*itmax)); iter++) {
  893.  
  894.       for (i = (*is) + 1; i<=((*ie) + 1); i++) {
  895.         if (itmask[i-1] != MASK_FALSE)  {
  896.           dtit[i-1] = huge8;
  897.         }
  898.       }
  899.  
  900.       // Compute the cooling rate, tgas, tdust, and metallicity for this row
  901.  
  902.        FORTRAN_NAME(cool1d_multi_g)(
  903.                 d.data(), e.data(), u.data(), v.data(), w.data(), de.data(), HI.data(), HII.data(), HeI.data(), HeII.data(), HeIII.data(),
  904.                 in, jn, kn, nratec,
  905.                 iexpand, ispecies, imetal, imcool,
  906.                 idust, idustall, idustfield, idustrec,
  907.                 idim, is, ie, &j, &k, ih2co, ipiht, &iter, igammah,
  908.                 aye, temstart, temend, z_solar, fgr,
  909.                 utem, uxyz, uaye, urho, utim,
  910.                 gamma, fh,
  911.                 ceHIa, ceHeIa, ceHeIIa, ciHIa, ciHeIa,
  912.                 ciHeISa, ciHeIIa, reHIIa, reHeII1a,
  913.                 reHeII2a, reHeIIIa, brema, compa, gammaha,
  914.                 isrf, regra, gamma_isrfa, comp_xraya, comp_temp,
  915.                 piHI, piHeI, piHeII, &comp1, &comp2,
  916.                 HM.data(), H2I.data(), H2II.data(), DI.data(), DII.data(), HDI.data(), metal.data(), dust.data(),
  917.                 hyd01ka, h2k01a, vibha, rotha, rotla,
  918.                 hyd01k.data(), h2k01.data(), vibh.data(), roth.data(), rotl.data(),
  919.                 gpldla, gphdla, gpldl.data(), gphdl.data(),
  920.                 hdltea, hdlowa, hdlte.data(), hdlow.data(),
  921.                 gaHIa, gaH2a, gaHea, gaHpa, gaela,
  922.                 h2ltea, gasgra,
  923.                 ceHI.data(), ceHeI.data(), ceHeII.data(), ciHI.data(), ciHeI.data(), ciHeIS.data(), ciHeII.data(),
  924.                 reHII.data(), reHeII1.data(), reHeII2.data(), reHeIII.data(), brem.data(),
  925.                 indixe.data(), t1.data(), t2.data(), logtem.data(), tdef.data(), edot.data(),
  926.                 tgas.data(), tgasold.data(), mmw.data(), p2d.data(), tdust.data(), metallicity.data(),
  927.                 dust2gas.data(), rhoH.data(), mynh.data(), myde.data(),
  928.                 gammaha_eff.data(), gasgr_tdust.data(), regr.data(),
  929.                 iradshield, avgsighi, avgsighei, avgsigheii,
  930.                 k24, k26,
  931.                 iradtrans, photogamma.data(),
  932.                 ih2optical, iciecool, ih2cr, ihdcr, ciecoa, cieco.data(),
  933.                 icmbTfloor, iClHeat, clEleFra,
  934.                 priGridRank, priGridDim,
  935.                 priPar1, priPar2, priPar3, priPar4, priPar5,
  936.                 priDataSize, priCooling, priHeating, priMMW,
  937.                 metGridRank, metGridDim,
  938.                 metPar1, metPar2, metPar3, metPar4, metPar5,
  939.                 metDataSize, metCooling, metHeating, clnew,
  940.                 iVheat, iMheat, Vheat.data(), Mheat.data(),
  941.                 iTfloor, Tfloor_scalar, Tfloor.data(),
  942.                 iisrffield, isrf_habing, itmask.data(), itmask_metal.data(),
  943.                imchem, igrgr, ipcont, tmcool,
  944.                DM.data(), HDII.data(), HeHII.data(),
  945.                CI.data(), CII.data(), CO.data(), CO2.data(),
  946.                OI.data(), OH.data(), H2O.data(), O2.data(),
  947.                SiI.data(), SiOI.data(), SiO2I.data(),
  948.                CH.data(), CH2.data(), COII.data(), OII.data(),
  949.                OHII.data(), H2OII.data(), H3OII.data(), O2II.data(),
  950.                Mg.data(), Al.data(), S.data(), Fe.data(),
  951.                SiM.data(), FeM.data(), Mg2SiO4.data(), MgSiO3.data(), Fe3O4.data(),
  952.                AC.data(), SiO2D.data(), MgO.data(), FeS.data(), Al2O3.data(),
  953.                reforg.data(), volorg.data(), H2Oice.data(),
  954.                cieY06a,
  955.                LH2_N, LH2_Size,
  956.                LH2_D, LH2_T, LH2_H,
  957.                LH2_dD, LH2_dT, LH2_dH, LH2_L,
  958.                LHD_N, LHD_Size,
  959.                LHD_D, LHD_T, LHD_H,
  960.                LHD_dD, LHD_dT, LHD_dH, LHD_L,
  961.                LCI_N, LCI_Size,
  962.                LCI_D, LCI_T, LCI_H,
  963.                LCI_dD, LCI_dT, LCI_dH, LCI_L,
  964.                LCII_N, LCII_Size,
  965.                LCII_D, LCII_T, LCII_H,
  966.                LCII_dD, LCII_dT, LCII_dH, LCII_L,
  967.                LOI_N, LOI_Size,
  968.                LOI_D, LOI_T, LOI_H,
  969.                LOI_dD, LOI_dT, LOI_dH, LOI_L,
  970.                LCO_N, LCO_Size,
  971.                LCO_D, LCO_T, LCO_H,
  972.                LCO_dD, LCO_dT, LCO_dH, LCO_L,
  973.                LOH_N, LOH_Size,
  974.                LOH_D, LOH_T, LOH_H,
  975.                LOH_dD, LOH_dT, LOH_dH, LOH_L,
  976.                LH2O_N, LH2O_Size,
  977.                LH2O_D, LH2O_T, LH2O_H,
  978.                LH2O_dD, LH2O_dT, LH2O_dH, LH2O_L,
  979.                alphap_N, alphap_Size,
  980.                alphap_D, alphap_T, alphap_dD, alphap_dT,
  981.                alphap_Data,
  982.                immulti, imabund, idspecies, itdmulti, idsub,
  983.                metal_loc.data(),
  984.                metal_C13.data(), metal_C20.data(), metal_C25.data(), metal_C30.data(),
  985.                metal_F13.data(), metal_F15.data(), metal_F50.data(), metal_F80.data(),
  986.                metal_P170.data(), metal_P200.data(), metal_Y19.data(),
  987.                SN0_N,
  988.                SN0_fSiM, SN0_fFeM, SN0_fMg2SiO4, SN0_fMgSiO3,
  989.                SN0_fFe3O4, SN0_fAC, SN0_fSiO2D, SN0_fMgO,
  990.                SN0_fFeS, SN0_fAl2O3,
  991.                SN0_freforg, SN0_fvolorg, SN0_fH2Oice,
  992.                SN0_r0SiM, SN0_r0FeM, SN0_r0Mg2SiO4, SN0_r0MgSiO3,
  993.                SN0_r0Fe3O4, SN0_r0AC, SN0_r0SiO2D, SN0_r0MgO,
  994.                SN0_r0FeS, SN0_r0Al2O3,
  995.                SN0_r0reforg, SN0_r0volorg, SN0_r0H2Oice,
  996.                gr_N, gr_Size, gr_dT, gr_Td,
  997.                SN0_kpSiM, SN0_kpFeM, SN0_kpMg2SiO4, SN0_kpMgSiO3,
  998.                SN0_kpFe3O4, SN0_kpAC, SN0_kpSiO2D, SN0_kpMgO,
  999.                SN0_kpFeS, SN0_kpAl2O3,
  1000.                SN0_kpreforg, SN0_kpvolorg, SN0_kpH2Oice,
  1001.                tSiM.data(), tFeM.data(), tMg2SiO4.data(), tMgSiO3.data(), tFe3O4.data(),
  1002.                tAC.data(), tSiO2D.data(), tMgO.data(), tFeS.data(), tAl2O3.data(),
  1003.                treforg.data(), tvolorg.data(), tH2Oice.data(),
  1004.                gasgr2a, gamma_isrf2a
  1005.           );
  1006.  
  1007.       if ((*ispecies) > 0)  {
  1008.  
  1009.         // Look-up rates as a function of temperature for 1D set of zones
  1010.         //  (maybe should add itmask to this call)
  1011.  
  1012.          FORTRAN_NAME(lookup_cool_rates1d_g)(temstart, temend, nratec, &j, &k,
  1013.                  is, ie, ithreebody,
  1014.                  in, jn, kn, ispecies, &anydust,
  1015.                  iH2shield, iradshield,
  1016.                  tgas.data(), mmw.data(), d.data(), HI.data(), HII.data(), HeI.data(), HeII.data(), HeIII.data(),
  1017.                  HM.data(), H2I.data(), H2II.data(), DI.data(), DII.data(), HDI.data(),
  1018.                  tdust.data(), dust2gas.data(),
  1019.                  k1a, k2a, k3a, k4a, k5a, k6a, k7a, k8a, k9a, k10a,
  1020.                  k11a, k12a, k13a, k13dda, k14a, k15a, k16a,
  1021.                  k17a, k18a, k19a, k22a,
  1022.                  k50a, k51a, k52a, k53a, k54a, k55a, k56a,
  1023.                  k57a, k58a, ndratec, dtemstart, dtemend, h2dusta,
  1024.                  ncrna, ncrd1a, ncrd2a,
  1025.                  avgsighi, avgsighei, avgsigheii, piHI, piHeI,
  1026.                  k1.data(), k2.data(), k3.data(), k4.data(), k5.data(), k6.data(), k7.data(), k8.data(), k9.data(), k10.data(),
  1027.                  k11.data(), k12.data(), k13.data(), k14.data(), k15.data(), k16.data(), k17.data(), k18.data(),
  1028.                  k19.data(), k22.data(), k24, k25, k26, k28, k29, k30, k31,
  1029.                  k50.data(), k51.data(), k52.data(), k53.data(), k54.data(), k55.data(), k56.data(), k57.data(),
  1030.                  k58.data(), k13dd.data(), k24shield.data(), k25shield.data(), k26shield.data(),
  1031.                  k28shield.data(), k29shield.data(), k30shield.data(),
  1032.                  k31shield.data(), h2dust.data(), ncrn.data(), ncrd1.data(), ncrd2.data(),
  1033.                  t1.data(), t2.data(), tdef.data(), logtem.data(), indixe.data(),
  1034.                  &dom, &coolunit, &tbase1, &xbase1, &dx_cgs, &c_ljeans,
  1035.                  iradtrans, kdissH2I.data(), xH2shield.data(), itmask.data(),
  1036.                  itmask_metal.data(),
  1037.                  fh, metal.data(),
  1038.                  DM.data(), HDII.data(), HeHII.data(), imetal, imchem, igrgr,
  1039.                  CI.data(), CII.data(), CO.data(), CO2.data(),
  1040.                  OI.data(), OH.data(), H2O.data(), O2.data(),
  1041.                  SiI.data(), SiOI.data(), SiO2I.data(),
  1042.                  CH.data(), CH2.data(), COII.data(), OII.data(),
  1043.                  OHII.data(), H2OII.data(), H3OII.data(), O2II.data(),
  1044.                  Mg.data(), Al.data(), S.data(), Fe.data(),
  1045.                  SiM.data(), FeM.data(), Mg2SiO4.data(), MgSiO3.data(), Fe3O4.data(),
  1046.                  AC.data(), SiO2D.data(), MgO.data(), FeS.data(), Al2O3.data(),
  1047.                  reforg.data(), volorg.data(), H2Oice.data(),
  1048.                  k125a, k129a, k130a, k131a, k132a,
  1049.                  k133a, k134a, k135a, k136a, k137a,
  1050.                  k148a, k149a, k150a, k151a, k152a,
  1051.                  k153a,
  1052.                  kz15a, kz16a, kz17a, kz18a, kz19a,
  1053.                  kz20a, kz21a, kz22a, kz23a, kz24a,
  1054.                  kz25a, kz26a, kz27a, kz28a, kz29a,
  1055.                  kz30a, kz31a, kz32a, kz33a, kz34a,
  1056.                  kz35a, kz36a, kz37a, kz38a, kz39a,
  1057.                  kz40a, kz41a, kz42a, kz43a, kz44a,
  1058.                  kz45a, kz46a, kz47a, kz48a, kz49a,
  1059.                  kz50a, kz51a, kz52a, kz53a, kz54a,
  1060.                  k125.data(),  k129.data(),  k130.data(),  k131.data(),  k132.data(),
  1061.                  k133.data(),  k134.data(),  k135.data(),  k136.data(),  k137.data(),
  1062.                  k148.data(),  k149.data(),  k150.data(),  k151.data(),  k152.data(),
  1063.                  k153.data(),
  1064.                  kz15.data(),  kz16.data(),  kz17.data(),  kz18.data(),  kz19.data(),
  1065.                  kz20.data(),  kz21.data(),  kz22.data(),  kz23.data(),  kz24.data(),
  1066.                  kz25.data(),  kz26.data(),  kz27.data(),  kz28.data(),  kz29.data(),
  1067.                  kz30.data(),  kz31.data(),  kz32.data(),  kz33.data(),  kz34.data(),
  1068.                  kz35.data(),  kz36.data(),  kz37.data(),  kz38.data(),  kz39.data(),
  1069.                  kz40.data(),  kz41.data(),  kz42.data(),  kz43.data(),  kz44.data(),
  1070.                  kz45.data(),  kz46.data(),  kz47.data(),  kz48.data(),  kz49.data(),
  1071.                  kz50.data(),  kz51.data(),  kz52.data(),  kz53.data(),  kz54.data(),
  1072.                  immulti, imabund, idspecies, itdmulti, idsub,
  1073.                  metal_loc.data(),
  1074.                  metal_C13.data(), metal_C20.data(), metal_C25.data(), metal_C30.data(),
  1075.                  metal_F13.data(), metal_F15.data(), metal_F50.data(), metal_F80.data(),
  1076.                  metal_P170.data(), metal_P200.data(), metal_Y19.data(),
  1077.                  SN0_N,
  1078.                  SN0_fSiM, SN0_fFeM, SN0_fMg2SiO4, SN0_fMgSiO3,
  1079.                  SN0_fFe3O4, SN0_fAC, SN0_fSiO2D, SN0_fMgO,
  1080.                  SN0_fFeS, SN0_fAl2O3,
  1081.                  SN0_freforg, SN0_fvolorg, SN0_fH2Oice,
  1082.                  SN0_r0SiM, SN0_r0FeM, SN0_r0Mg2SiO4, SN0_r0MgSiO3,
  1083.                  SN0_r0Fe3O4, SN0_r0AC, SN0_r0SiO2D, SN0_r0MgO,
  1084.                  SN0_r0FeS, SN0_r0Al2O3,
  1085.                  SN0_r0reforg, SN0_r0volorg, SN0_r0H2Oice,
  1086.                  gr_N, gr_Size, gr_dT, gr_Td,
  1087.                  SN0_kpSiM, SN0_kpFeM, SN0_kpMg2SiO4, SN0_kpMgSiO3,
  1088.                  SN0_kpFe3O4, SN0_kpAC, SN0_kpSiO2D, SN0_kpMgO,
  1089.                  SN0_kpFeS, SN0_kpAl2O3,
  1090.                  SN0_kpreforg, SN0_kpvolorg, SN0_kpH2Oice,
  1091.                  h2dustSa, h2dustCa, rhoH.data(), grogra, dt,
  1092.                  kdSiM.data(), kdFeM.data(), kdMg2SiO4.data(),
  1093.                  kdMgSiO3.data(), kdFe3O4.data(), kdAC.data(), kdSiO2D.data(), kdMgO.data(), kdFeS.data(),
  1094.                  kdAl2O3.data(), kdreforg.data(), kdvolorg.data(), kdH2Oice.data(),
  1095.                  tSiM.data(), tFeM.data(), tMg2SiO4.data(), tMgSiO3.data(), tFe3O4.data(),
  1096.                  tAC.data(), tSiO2D.data(), tMgO.data(), tFeS.data(), tAl2O3.data(),
  1097.                  treforg.data(), tvolorg.data(), tH2Oice.data(), iuseH2shield
  1098.         );
  1099.  
  1100.         // Compute dedot and HIdot, the rates of change of de and HI
  1101.         //   (should add itmask to this call)
  1102.  
  1103.          FORTRAN_NAME(rate_timestep_g)(
  1104.                        dedot.data(), HIdot.data(), ispecies, &anydust,
  1105.                        de.data(), HI.data(), HII.data(), HeI.data(), HeII.data(), HeIII.data(), d.data(),
  1106.                        HM.data(), H2I.data(), H2II.data(),
  1107.                        in, jn, kn, is, ie, &j, &k,
  1108.                        k1.data(), k2.data(), k3.data(), k4.data(), k5.data(), k6.data(), k7.data(), k8.data(), k9.data(), k10.data(), k11.data(),
  1109.                        k12.data(), k13.data(), k14.data(), k15.data(), k16.data(), k17.data(), k18.data(), k19.data(), k22.data(),
  1110.                        k24, k25, k26, k27, k28, k29, k30,
  1111.                        k50.data(), k51.data(), k52.data(), k53.data(), k54.data(), k55.data(), k56.data(), k57.data(), k58.data(),
  1112.                        h2dust.data(), ncrn.data(), ncrd1.data(), ncrd2.data(), rhoH.data(),
  1113.                        k24shield.data(), k25shield.data(), k26shield.data(),
  1114.                        k28shield.data(), k29shield.data(), k30shield.data(), k31shield.data(),
  1115.                        iradtrans, irt_honly,
  1116.                        kphHI.data(), kphHeI.data(), kphHeII.data(),
  1117.                        itmask.data(), edot.data(), &chunit, &dom, metal.data(),
  1118.                       HDI.data(), imchem, CI.data(), OI.data(), OH.data(), CO.data(), H2O.data(),
  1119.                       idissHDI, kdissHDI.data(), iionZ, kphCI.data(), kphOI.data(),
  1120.                       idissZ, kdissCO.data(), kdissOH.data(), kdissH2O.data()
  1121.                             );
  1122.  
  1123.         // move itmask temporary array
  1124.         // then split cells with low densities
  1125.         //    => Gauss-Seidel scheme
  1126.         //              and with high densities
  1127.         //    => Newton-Raphson scheme
  1128.  
  1129.         std::memcpy(itmask_tmp.data(), itmask.data(), sizeof(gr_mask_type)*(*in));
  1130.         std::memcpy(itmask_nr.data(), itmask.data(), sizeof(gr_mask_type)*(*in));
  1131.         for (i = (*is) + 1; i<=((*ie) + 1); i++) {
  1132.           if ( itmask_tmp[i-1] != MASK_FALSE )  {
  1133.  
  1134.             if ( ( ((*imetal) == 0)
  1135.               &&  (ddom[i-1] < 1.e8) )
  1136.              ||  ( ((*imetal) == 1)
  1137.               &&  ( ( (metallicity[i-1] <= min_metallicity)
  1138.                   &&  (ddom[i-1] < 1.e8) )
  1139.                  ||  ( (metallicity[i-1] > min_metallicity)
  1140.                   &&  (ddom[i-1] < 1.e6) ) ) ) )  {
  1141.               itmask_nr[i-1] = MASK_FALSE;
  1142.             } else {
  1143.               itmask[i-1] = MASK_FALSE;
  1144.             }
  1145.  
  1146.           }
  1147.         }
  1148.  
  1149.         for (i = (*is) + 1; i<=((*ie) + 1); i++) {
  1150.           if (itmask_nr[i-1] != MASK_FALSE)  {
  1151.             if ( ((*icool) == 1)  &&  ((*ispecies) > 1)  &&
  1152.                  ((ddom[i-1] > 1.e7)
  1153.              && (tgas[i-1] > 1650.e0)) )  {
  1154.               imp_eng[i-1] = 1;
  1155.             } else {
  1156.               imp_eng[i-1] = 0;
  1157.             }
  1158.           }
  1159.         }
  1160.  
  1161.         // Find timestep that keeps relative chemical changes below 10%
  1162.  
  1163.         for (i = (*is) + 1; i<=((*ie) + 1); i++) {
  1164.           if (itmask[i-1] != MASK_FALSE)  {
  1165.             // Bound from below to prevent numerical errors
  1166.                
  1167.             if (std::fabs(dedot[i-1]) < tiny8)
  1168.                        { dedot[i-1] = std::fmin(tiny,de(i-1,j-1,k-1)); }
  1169.             if (std::fabs(HIdot[i-1]) < tiny8)
  1170.                        { HIdot[i-1] = std::fmin(tiny,HI(i-1,j-1,k-1)); }
  1171.  
  1172.             // If the net rate is almost perfectly balanced then set
  1173.             //     it to zero (since it is zero to available precision)
  1174.  
  1175.             if (std::fmin(std::fabs(k1[i-1]* de(i-1,j-1,k-1)*HI(i-1,j-1,k-1)),
  1176.                     std::fabs(k2[i-1]*HII(i-1,j-1,k-1)*de(i-1,j-1,k-1)))/
  1177.                 std::fmax(std::fabs(dedot[i-1]),std::fabs(HIdot[i-1])) >
  1178.                  1.0e6)  {
  1179.               dedot[i-1] = tiny8;
  1180.               HIdot[i-1] = tiny8;
  1181.             }
  1182.  
  1183.             // If the iteration count is high then take the smaller of
  1184.             //   the calculated dedot and last time step's actual dedot.
  1185.             //   This is intended to get around the problem of a low
  1186.             //   electron or HI fraction which is in equilibrium with high
  1187.             //   individual terms (which all nearly cancel).
  1188.  
  1189.             if (iter > 50)  {
  1190.               dedot[i-1] = std::fmin(std::fabs(dedot[i-1]), std::fabs(dedot_prev[i-1]));
  1191.               HIdot[i-1] = std::fmin(std::fabs(HIdot[i-1]), std::fabs(HIdot_prev[i-1]));
  1192.             }
  1193.  
  1194.             // compute minimum rate timestep
  1195.  
  1196.             olddtit = dtit[i-1];
  1197.             dtit[i-1] = grackle::impl::fmin(std::fabs(0.1*de(i-1,j-1,k-1)/dedot[i-1]),
  1198.                  std::fabs(0.1*HI(i-1,j-1,k-1)/HIdot[i-1]),
  1199.                  (*dt)-ttot[i-1], 0.5*(*dt));
  1200.  
  1201.             if (ddom[i-1] > 1.e8  &&
  1202.                  edot[i-1] > 0.       &&
  1203.                 (*ispecies) > 1)  {
  1204.               // Equilibrium value for H is:
  1205.               // H = (-1._DKIND / (4*k22)) * (k13 - sqrt(8 k13 k22 rho + k13^2))
  1206.               // We now want this to change by 10% or less, but we're only
  1207.               // differentiating by dT.  We have de/dt.  We need dT/de.
  1208.               // T = (g-1)*p2d*utem/N; tgas == (g-1)(p2d*utem/N)
  1209.               // dH_eq / dt = (dH_eq/dT) * (dT/de) * (de/dt)
  1210.               // dH_eq / dT (see above; we can calculate the derivative here)
  1211.               // dT / de = utem * (gamma - 1._DKIND) / N == tgas / p2d
  1212.               // de / dt = edot
  1213.               // Now we use our estimate of dT/de to get the estimated
  1214.               // difference in the equilibrium
  1215.               eqt2 = std::fmin(std::log(tgas[i-1]) + 0.1*dlogtem, t2[i-1]);
  1216.               eqtdef = (eqt2 - t1[i-1])/(t2[i-1] - t1[i-1]);
  1217.               eqk222 = k22a[indixe[i-1]-1] +
  1218.                 (k22a[indixe[i-1]+1-1] -k22a[indixe[i-1]-1])*eqtdef;
  1219.               eqk132 = k13a[indixe[i-1]-1] +
  1220.                 (k13a[indixe[i-1]+1-1] -k13a[indixe[i-1]-1])*eqtdef;
  1221.               heq2 = (-1. / (4.*eqk222)) * (eqk132-
  1222.                    std::sqrt(8.*eqk132*eqk222*
  1223.                    (*fh)*d(i-1,j-1,k-1)+std::pow(eqk132,2.)));
  1224.  
  1225.               eqt1 = std::fmax(std::log(tgas[i-1]) - 0.1*dlogtem, t1[i-1]);
  1226.               eqtdef = (eqt1 - t1[i-1])/(t2[i-1] - t1[i-1]);
  1227.               eqk221 = k22a[indixe[i-1]-1] +
  1228.                 (k22a[indixe[i-1]+1-1] -k22a[indixe[i-1]-1])*eqtdef;
  1229.               eqk131 = k13a[indixe[i-1]-1] +
  1230.                 (k13a[indixe[i-1]+1-1] -k13a[indixe[i-1]-1])*eqtdef;
  1231.               heq1 = (-1. / (4.*eqk221)) * (eqk131-
  1232.                    std::sqrt(8.*eqk131*eqk221*
  1233.                    (*fh)*d(i-1,j-1,k-1)+std::pow(eqk131,2.)));
  1234.  
  1235.               dheq = (std::fabs(heq2-heq1)/(std::exp(eqt2) - std::exp(eqt1)))
  1236.                    * (tgas[i-1]/p2d[i-1]) * edot[i-1];
  1237.               heq = (-1. / (4.*k22[i-1])) * (k13[i-1]-
  1238.                    std::sqrt(8.*k13[i-1]*k22[i-1]*
  1239.                    (*fh)*d(i-1,j-1,k-1)+std::pow(k13[i-1],2.)));
  1240.               dtit[i-1] = std::fmin(dtit[i-1], 0.1*heq/dheq);
  1241.             }
  1242.             if (iter>10LL)  {
  1243.               dtit[i-1] = std::fmin(olddtit*1.5, dtit[i-1]);
  1244.             }
  1245.  
  1246.           } else if ((itmask_nr[i-1]!=MASK_FALSE) &&
  1247.                    (imp_eng[i-1]==0))  {
  1248.             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)),
  1249.                  (*dt)-ttot[i-1], 0.5*(*dt));
  1250.  
  1251.           } else if ((itmask_nr[i-1]!=MASK_FALSE) &&
  1252.                    (imp_eng[i-1]==1))  {
  1253.             dtit[i-1] = (*dt) - ttot[i-1];
  1254.  
  1255.           } else {
  1256.             dtit[i-1] = (*dt);
  1257.           }
  1258.         }
  1259.  
  1260.       }
  1261.  
  1262.       // Compute maximum timestep for cooling/heating
  1263.  
  1264.       for (i = (*is) + 1; i<=((*ie) + 1); i++) {
  1265.         if (itmask[i-1] != MASK_FALSE)  {
  1266.           // Set energy per unit volume of this cell based in the pressure
  1267.           // (the gamma used here is the right one even for H2 since p2d
  1268.           //  is calculated with this gamma).
  1269.  
  1270.           energy = std::fmax(p2d[i-1]/((*gamma)-1.), tiny8);
  1271.  
  1272.           // If the temperature is at the bottom of the temperature look-up
  1273.           // table and edot < 0, then shut off the cooling.
  1274.  
  1275.           if (tgas[i-1] <= 1.01*(*temstart)  &&
  1276.                edot[i-1] < 0.)
  1277.                { edot[i-1] = tiny8; }
  1278.           if (std::fabs(edot[i-1]) < tiny8) { edot[i-1] = tiny8; }
  1279.  
  1280.           // Compute timestep for 10% change
  1281.  
  1282.           dtit[i-1] = grackle::impl::fmin((double)(std::fabs(0.1*
  1283.             energy/edot[i-1]) ),
  1284.             (*dt)-ttot[i-1], dtit[i-1]);
  1285.  
  1286.           if (dtit[i-1] != dtit[i-1])  {
  1287.             OMP_PRAGMA_CRITICAL
  1288.             {
  1289.               eprintf("HUGE dtit ::  %g %g %g %g %g %g %g\n",
  1290.                       energy,
  1291.                       edot [ i-1 ],
  1292.                       dtit [ i-1 ],
  1293.                       (*dt),
  1294.                       ttot [ i-1 ],
  1295.                       std::fabs ( 0.1 * energy / edot [ i-1 ] ),
  1296.                       (double) ( std::fabs ( 0.1 * energy / edot [ i-1 ] ) ));
  1297.             }
  1298.           }
  1299.  
  1300.         }
  1301.       }
  1302.  
  1303.       // Update total and gas energy
  1304.  
  1305.       if ((*icool) == 1)  {
  1306.         for (i = (*is) + 1; i<=((*ie) + 1); i++) {
  1307.           if (itmask[i-1] != MASK_FALSE)  {
  1308.             e(i-1,j-1,k-1)  = e(i-1,j-1,k-1) +
  1309.                     (gr_float)(edot[i-1]/d(i-1,j-1,k-1)*dtit[i-1] );
  1310.  
  1311.           }
  1312.         }
  1313.       }
  1314.  
  1315.       if ((*ispecies) > 0)  {
  1316.  
  1317.         // Solve rate equations with one linearly implicit Gauss-Seidel
  1318.         // sweep of a backward Euler method ---
  1319.  
  1320.          FORTRAN_NAME(step_rate_g)(de.data(), HI.data(), HII.data(), HeI.data(), HeII.data(), HeIII.data(), d.data(),
  1321.                        HM.data(), H2I.data(), H2II.data(), DI.data(), DII.data(), HDI.data(), dtit.data(),
  1322.                        in, jn, kn, is, ie, &j, &k,
  1323.                        ispecies, &anydust,
  1324.                        k1.data(), k2.data(), k3.data(), k4.data(), k5.data(), k6.data(), k7.data(), k8.data(), k9.data(), k10.data(), k11.data(),
  1325.                        k12.data(), k13.data(), k14.data(), k15.data(), k16.data(), k17.data(), k18.data(), k19.data(), k22.data(),
  1326.                        k24, k25, k26, k27, k28, k29, k30,
  1327.                        k50.data(), k51.data(), k52.data(), k53.data(), k54.data(), k55.data(), k56.data(), k57.data(), k58.data(),
  1328.                        h2dust.data(), rhoH.data(),
  1329.                        k24shield.data(), k25shield.data(), k26shield.data(),
  1330.                        k28shield.data(), k29shield.data(), k30shield.data(), k31shield.data(),
  1331.                        HIp.data(), HIIp.data(), HeIp.data(), HeIIp.data(), HeIIIp.data(), dep.data(),
  1332.                        HMp.data(), H2Ip.data(), H2IIp.data(), DIp.data(), DIIp.data(), HDIp.data(),
  1333.                        dedot_prev.data(), HIdot_prev.data(),
  1334.                        iradtrans, irt_honly,
  1335.                        kphHI.data(), kphHeI.data(), kphHeII.data(),
  1336.                        itmask.data(), itmask_metal.data(),
  1337.                       DM.data(), HDII.data(), HeHII.data(), imetal, metal.data(),
  1338.                       imchem, idspecies, igrgr, idsub,
  1339.                       CI.data(), CII.data(), CO.data(), CO2.data(),
  1340.                       OI.data(), OH.data(), H2O.data(), O2.data(),
  1341.                       SiI.data(), SiOI.data(), SiO2I.data(),
  1342.                       CH.data(), CH2.data(), COII.data(), OII.data(),
  1343.                       OHII.data(), H2OII.data(), H3OII.data(), O2II.data(),
  1344.                       Mg.data(), Al.data(), S.data(), Fe.data(),
  1345.                       SiM.data(), FeM.data(), Mg2SiO4.data(), MgSiO3.data(), Fe3O4.data(),
  1346.                       AC.data(), SiO2D.data(), MgO.data(), FeS.data(), Al2O3.data(),
  1347.                       reforg.data(), volorg.data(), H2Oice.data(),
  1348.                       k125.data(), k129.data(), k130.data(), k131.data(), k132.data(),
  1349.                       k133.data(), k134.data(), k135.data(), k136.data(), k137.data(),
  1350.                       k148.data(), k149.data(), k150.data(), k151.data(), k152.data(),
  1351.                       k153.data(),
  1352.                       kz15.data(),  kz16.data(),  kz17.data(),  kz18.data(),  kz19.data(),
  1353.                       kz20.data(),  kz21.data(),  kz22.data(),  kz23.data(),  kz24.data(),
  1354.                       kz25.data(),  kz26.data(),  kz27.data(),  kz28.data(),  kz29.data(),
  1355.                       kz30.data(),  kz31.data(),  kz32.data(),  kz33.data(),  kz34.data(),
  1356.                       kz35.data(),  kz36.data(),  kz37.data(),  kz38.data(),  kz39.data(),
  1357.                       kz40.data(),  kz41.data(),  kz42.data(),  kz43.data(),  kz44.data(),
  1358.                       kz45.data(),  kz46.data(),  kz47.data(),  kz48.data(),  kz49.data(),
  1359.                       kz50.data(),  kz51.data(),  kz52.data(),  kz53.data(),  kz54.data(),
  1360.                       DMp.data(), HDIIp.data(), HeHIIp.data(),
  1361.                       CIp.data(), CIIp.data(), COp.data(), CO2p.data(),
  1362.                       OIp.data(), OHp.data(), H2Op.data(), O2p.data(),
  1363.                       SiIp.data(), SiOIp.data(), SiO2Ip.data(),
  1364.                       CHp.data(), CH2p.data(), COIIp.data(), OIIp.data(),
  1365.                       OHIIp.data(), H2OIIp.data(), H3OIIp.data(), O2IIp.data(),
  1366.                       Mgp.data(), Alp.data(), Sp.data(), Fep.data(),
  1367.                       SiMp.data(), FeMp.data(), Mg2SiO4p.data(), MgSiO3p.data(), Fe3O4p.data(),
  1368.                       ACp.data(), SiO2Dp.data(), MgOp.data(), FeSp.data(), Al2O3p.data(),
  1369.                       reforgp.data(), volorgp.data(), H2Oicep.data(),
  1370.                       kdSiM.data(), kdFeM.data(), kdMg2SiO4.data(), kdMgSiO3.data(), kdFe3O4.data(),
  1371.                       kdAC.data(), kdSiO2D.data(), kdMgO.data(), kdFeS.data(), kdAl2O3.data(),
  1372.                       kdreforg.data(), kdvolorg.data(), kdH2Oice.data(),
  1373.                       idissHDI, kdissHDI.data(), iionZ, kphCI.data(), kphOI.data(),
  1374.                       idissZ, kdissCO.data(), kdissOH.data(), kdissH2O.data()
  1375.              );
  1376.  
  1377.         // Note 10/18/2024: the code from here to the comment "end Newton-Raphson scheme"
  1378.         // should be put into its own function.
  1379.         //       Start Newton-Raphson scheme
  1380.         for (i = (*is) + 1; i<=((*ie) + 1); i++) {
  1381.           if (itmask_nr[i-1] != MASK_FALSE)  {
  1382.  
  1383.             // If density and temperature are low, update gas energy explicitly
  1384.  
  1385.             if ((*icool) == 1)  {
  1386.               if (imp_eng[i-1] == 0)  {
  1387.                 e(i-1,j-1,k-1)  = e(i-1,j-1,k-1) +
  1388.                      (gr_float)(edot[i-1]/d(i-1,j-1,k-1)*dtit[i-1] );
  1389.               }
  1390.             }
  1391.  
  1392.             // initialize arrays
  1393.             if ((*ispecies) > 0) { nsp = 6; }
  1394.             if ((*ispecies) > 1) { nsp = nsp + 3; }
  1395.             if ((*ispecies) > 2) { nsp = nsp + 3; }
  1396.             if ((*ispecies) > 3) { nsp = nsp + 3; }
  1397.             if (itmask_metal[i-1] != MASK_FALSE)  {
  1398.               if ((*imchem) == 1)  {
  1399.                 nsp = nsp + 19;
  1400.                 if ( ( (*igrgr) == 1 )  ||  ( (*idsub) == 1) )  {
  1401.                   if ((*idspecies) > 0) { nsp = nsp + 1; }
  1402.                   if ((*idspecies) > 1) { nsp = nsp + 3; }
  1403.                 }
  1404.               }
  1405.               if ( ( (*igrgr) == 1 )  ||  ( (*idsub) == 1) )  {
  1406.                 if ((*idspecies) > 0) { nsp = nsp + 2; }
  1407.                 if ((*idspecies) > 1) { nsp = nsp + 8; }
  1408.                 if ((*idspecies) > 2) { nsp = nsp + 3; }
  1409.               }
  1410.             }
  1411.             nsp = nsp + imp_eng[i-1];
  1412.             idsp.reserve(nsp);
  1413.             mtrx_data_.reserve(nsp * nsp);
  1414.             mtrx = grackle::impl::View<double**>(mtrx_data_.data(), nsp, nsp);
  1415.             vec.reserve(nsp);
  1416.  
  1417.             if ( (*ispecies) > 0 )  {
  1418.               dsp[ 1-1] = de(i-1,j-1,k-1);
  1419.               dsp[ 2-1] = HI(i-1,j-1,k-1);
  1420.               dsp[ 3-1] = HII(i-1,j-1,k-1);
  1421.               dsp[ 4-1] = HeI(i-1,j-1,k-1);
  1422.               dsp[ 5-1] = HeII(i-1,j-1,k-1);
  1423.               dsp[ 6-1] = HeIII(i-1,j-1,k-1);
  1424.             }
  1425.             if ( (*ispecies) > 1 )  {
  1426.               dsp[ 7-1] = HM(i-1,j-1,k-1);
  1427.               dsp[ 8-1] = H2I(i-1,j-1,k-1);
  1428.               dsp[ 9-1] = H2II(i-1,j-1,k-1);
  1429.             }
  1430.             if ( (*ispecies) > 2 )  {
  1431.               dsp[10-1] = DI(i-1,j-1,k-1);
  1432.               dsp[11-1] = DII(i-1,j-1,k-1);
  1433.               dsp[12-1] = HDI(i-1,j-1,k-1);
  1434.             }
  1435.             if ( (*ispecies) > 3 )  {
  1436.               dsp[13-1] = DM(i-1,j-1,k-1);
  1437.               dsp[14-1] = HDII(i-1,j-1,k-1);
  1438.               dsp[15-1] = HeHII(i-1,j-1,k-1);
  1439.             }
  1440.             if ( itmask_metal[i-1] != MASK_FALSE )  {
  1441.               if ( (*imchem) == 1 )  {
  1442.                 dsp[16-1] = CI(i-1,j-1,k-1);
  1443.                 dsp[17-1] = CII(i-1,j-1,k-1);
  1444.                 dsp[18-1] = CO(i-1,j-1,k-1);
  1445.                 dsp[19-1] = CO2(i-1,j-1,k-1);
  1446.                 dsp[20-1] = OI(i-1,j-1,k-1);
  1447.                 dsp[21-1] = OH(i-1,j-1,k-1);
  1448.                 dsp[22-1] = H2O(i-1,j-1,k-1);
  1449.                 dsp[23-1] = O2(i-1,j-1,k-1);
  1450.                 dsp[24-1] = SiI(i-1,j-1,k-1);
  1451.                 dsp[25-1] = SiOI(i-1,j-1,k-1);
  1452.                 dsp[26-1] = SiO2I(i-1,j-1,k-1);
  1453.                 dsp[27-1] = CH(i-1,j-1,k-1);
  1454.                 dsp[28-1] = CH2(i-1,j-1,k-1);
  1455.                 dsp[29-1] = COII(i-1,j-1,k-1);
  1456.                 dsp[30-1] = OII(i-1,j-1,k-1);
  1457.                 dsp[31-1] = OHII(i-1,j-1,k-1);
  1458.                 dsp[32-1] = H2OII(i-1,j-1,k-1);
  1459.                 dsp[33-1] = H3OII(i-1,j-1,k-1);
  1460.                 dsp[34-1] = O2II(i-1,j-1,k-1);
  1461.                 if ( ( (*igrgr) == 1 )  ||  ( (*idsub) == 1) )  {
  1462.                   if ((*idspecies) > 0)  {
  1463.                     dsp[35-1] = Mg(i-1,j-1,k-1);
  1464.                   }
  1465.                   if ((*idspecies) > 1)  {
  1466.                     dsp[36-1] = Al(i-1,j-1,k-1);
  1467.                     dsp[37-1] = S(i-1,j-1,k-1);
  1468.                     dsp[38-1] = Fe(i-1,j-1,k-1);
  1469.                   }
  1470.                 }
  1471.               }
  1472.               if ( ( (*igrgr) == 1 )  ||  ( (*idsub) == 1) )  {
  1473.                 if ((*idspecies) > 0)  {
  1474.                   dsp[39-1] = MgSiO3(i-1,j-1,k-1);
  1475.                   dsp[40-1] = AC(i-1,j-1,k-1);
  1476.                 }
  1477.                 if ((*idspecies) > 1)  {
  1478.                   dsp[41-1] = SiM(i-1,j-1,k-1);
  1479.                   dsp[42-1] = FeM(i-1,j-1,k-1);
  1480.                   dsp[43-1] = Mg2SiO4(i-1,j-1,k-1);
  1481.                   dsp[44-1] = Fe3O4(i-1,j-1,k-1);
  1482.                   dsp[45-1] = SiO2D(i-1,j-1,k-1);
  1483.                   dsp[46-1] = MgO(i-1,j-1,k-1);
  1484.                   dsp[47-1] = FeS(i-1,j-1,k-1);
  1485.                   dsp[48-1] = Al2O3(i-1,j-1,k-1);
  1486.                 }
  1487.                 if ((*idspecies) > 2)  {
  1488.                   dsp[49-1] = reforg(i-1,j-1,k-1);
  1489.                   dsp[50-1] = volorg(i-1,j-1,k-1);
  1490.                   dsp[51-1] = H2Oice(i-1,j-1,k-1);
  1491.                 }
  1492.               }
  1493.             }
  1494.             dsp[i_eng-1] = e(i-1,j-1,k-1);
  1495.  
  1496.             id = 0;
  1497.             if ((*ispecies) > 0)  {
  1498.               for (isp = 1; isp<=(6); isp++) {
  1499.                 id = id + 1;
  1500.                 idsp[id-1] = isp;
  1501.               }
  1502.             }
  1503.             if ((*ispecies) > 1)  {
  1504.               for (isp = 7; isp<=(9); isp++) {
  1505.                 id = id + 1;
  1506.                 idsp[id-1] = isp;
  1507.               }
  1508.             }
  1509.             if ((*ispecies) > 2)  {
  1510.               for (isp = 10; isp<=(12); isp++) {
  1511.                 id = id + 1;
  1512.                 idsp[id-1] = isp;
  1513.               }
  1514.             }
  1515.             if ((*ispecies) > 3)  {
  1516.               for (isp = 13; isp<=(15); isp++) {
  1517.                 id = id + 1;
  1518.                 idsp[id-1] = isp;
  1519.               }
  1520.             }
  1521.             if (itmask_metal[i-1] != MASK_FALSE)  {
  1522.               if ((*imchem) == 1)  {
  1523.                 for (isp = 16; isp<=(34); isp++) {
  1524.                   id = id + 1;
  1525.                   idsp[id-1] = isp;
  1526.                 }
  1527.                 if ( ( (*igrgr) == 1 )  ||  ( (*idsub) == 1) )  {
  1528.                   if ((*idspecies) > 0)  {
  1529.                     for (isp = 35; isp<=(35); isp++) {
  1530.                       id = id + 1;
  1531.                       idsp[id-1] = isp;
  1532.                     }
  1533.                   }
  1534.                   if ((*idspecies) > 1)  {
  1535.                     for (isp = 36; isp<=(38); isp++) {
  1536.                       id = id + 1;
  1537.                       idsp[id-1] = isp;
  1538.                     }
  1539.                   }
  1540.                 }
  1541.               }
  1542.               if ( ( (*igrgr) == 1 )  ||  ( (*idsub) == 1) )  {
  1543.                 if ((*idspecies) > 0)  {
  1544.                   for (isp = 39; isp<=(40); isp++) {
  1545.                     id = id + 1;
  1546.                     idsp[id-1] = isp;
  1547.                   }
  1548.                 }
  1549.                 if ((*idspecies) > 1)  {
  1550.                   for (isp = 41; isp<=(48); isp++) {
  1551.                     id = id + 1;
  1552.                     idsp[id-1] = isp;
  1553.                   }
  1554.                 }
  1555.                 if ((*idspecies) > 2)  {
  1556.                   for (isp = 49; isp<=(51); isp++) {
  1557.                     id = id + 1;
  1558.                     idsp[id-1] = isp;
  1559.                   }
  1560.                 }
  1561.               }
  1562.             }
  1563.             if ( imp_eng[i-1] ==1 )  {
  1564.               id = id + 1;
  1565.               idsp[id-1] = i_eng;
  1566.             }
  1567.  
  1568.             // Save arrays at ttot(i)
  1569.  
  1570.             std::memcpy(dsp0.data(), dsp.data(), sizeof(double)*i_eng);
  1571.             std::memset(ddsp.data(), 0, sizeof(double)*i_eng);
  1572.  
  1573.             // Search for the timestep for which chemistry converges
  1574.  
  1575.             ierror=1;
  1576.             itr_time=0;
  1577.             while ((ierror==1)) {
  1578.  
  1579.               // If not converge, restore arrays at ttot(i)
  1580.  
  1581.               std::memcpy(dsp.data(), dsp0.data(), sizeof(double)*i_eng);
  1582.               std::memset(ddsp.data(), 0, sizeof(double)*i_eng);
  1583.  
  1584.               // Iteration to solve ODEs
  1585.  
  1586.               err_max=1.e2;
  1587.               itr=0;
  1588.               while (err_max>1.e-8) {
  1589.                 if(itr>=20)  {
  1590.                   ierror = 1;
  1591.                   goto label_9996;
  1592.                 }
  1593.  
  1594.                  FORTRAN_NAME(lookup_cool_rates0d)(&dtit[i-1],
  1595.                 &d(i-1,j-1,k-1), &u(i-1,j-1,k-1), &v(i-1,j-1,k-1), &w(i-1,j-1,k-1),
  1596.                 &nsp, dsp.data(), dspdot.data(), nratec,
  1597.                 iexpand, ispecies, imetal, imcool,
  1598.                 idust, idustall, idustfield, idustrec,
  1599.                 ih2co, ipiht, &iter, igammah,
  1600.                 aye, temstart, temend, z_solar, fgr,
  1601.                 utem, uxyz, uaye, urho, utim,
  1602.                 gamma, fh,
  1603.                 ceHIa, ceHeIa, ceHeIIa, ciHIa, ciHeIa,
  1604.                 ciHeISa, ciHeIIa, reHIIa, reHeII1a,
  1605.                 reHeII2a, reHeIIIa, brema, compa, gammaha,
  1606.                 isrf, regra, gamma_isrfa, comp_xraya, comp_temp,
  1607.                 piHI, piHeI, piHeII, &comp1, &comp2,
  1608.                 &metal(i-1,j-1,k-1), &dust(i-1,j-1,k-1),
  1609.                 hyd01ka, h2k01a, vibha, rotha, rotla,
  1610.                 &hyd01k[i-1], &h2k01[i-1], &vibh[i-1], &roth[i-1], &rotl[i-1],
  1611.                 gpldla, gphdla, &gpldl[i-1], &gphdl[i-1],
  1612.                 hdltea, hdlowa, &hdlte[i-1], &hdlow[i-1],
  1613.                 gaHIa, gaH2a, gaHea,
  1614.                 gaHpa, gaela,
  1615.                 h2ltea, gasgra,
  1616.                 &ceHI[i-1], &ceHeI[i-1], &ceHeII[i-1], &ciHI[i-1], &ciHeI[i-1],
  1617.                 &ciHeIS[i-1], &ciHeII[i-1],
  1618.                 &reHII[i-1], &reHeII1[i-1], &reHeII2[i-1], &reHeIII[i-1], &brem[i-1],
  1619.                 &indixe[i-1], &t1[i-1], &t2[i-1], &logtem[i-1], &tdef[i-1], &edot[i-1],
  1620.                 &tgas[i-1], &tgasold[i-1], &mmw[i-1], &p2d[i-1], &tdust[i-1], &metallicity[i-1],
  1621.                 &dust2gas[i-1], &rhoH[i-1], &mynh[i-1], &myde[i-1],
  1622.                 &gammaha_eff[i-1], &gasgr_tdust[i-1], &regr[i-1],
  1623.                 iradshield, avgsighi, avgsighei, avgsigheii,
  1624.                 iradtrans, irt_honly,
  1625.                 ih2optical, iciecool, ih2cr, ihdcr, ciecoa, &cieco[i-1],
  1626.                 icmbTfloor, iClHeat, clEleFra,
  1627.                 priGridRank, priGridDim,
  1628.                 priPar1, priPar2, priPar3, priPar4, priPar5,
  1629.                 priDataSize, priCooling, priHeating, priMMW,
  1630.                 metGridRank, metGridDim,
  1631.                 metPar1, metPar2, metPar3, metPar4, metPar5,
  1632.                 metDataSize, metCooling, metHeating, clnew,
  1633.                 iVheat, iMheat, &Vheat(i-1,j-1,k-1), &Mheat(i-1,j-1,k-1),
  1634.                 iTfloor, Tfloor_scalar, &Tfloor(i-1,j-1,k-1),
  1635.                 iisrffield, isrf_habing,
  1636.                 ithreebody, &anydust, iH2shield,
  1637.                 k1a, k2a, k3a, k4a, k5a, k6a, k7a, k8a, k9a, k10a,
  1638.                 k11a, k12a, k13a, k13dda, k14a, k15a, k16a,
  1639.                 k17a, k18a, k19a, k22a,
  1640.                 k24, k25, k26, k27, k28, k29, k30, k31,
  1641.                 k50a, k51a, k52a, k53a, k54a, k55a, k56a,
  1642.                 k57a, k58a, ndratec, dtemstart, dtemend, h2dusta,
  1643.                 ncrna, ncrd1a, ncrd2a,
  1644.                 &h2dust[i-1], &ncrn[i-1], &ncrd1[i-1], &ncrd2[i-1],
  1645.                 &dom, &coolunit, &tbase1, &xbase1, &dx_cgs, &c_ljeans,
  1646.                 &kphHI(i-1,j-1,k-1), &kphHeI(i-1,j-1,k-1), &kphHeII(i-1,j-1,k-1), &kdissH2I(i-1,j-1,k-1),
  1647.                 &photogamma(i-1,j-1,k-1), &xH2shield(i-1,j-1,k-1), &chunit, &itmask_nr[i-1],
  1648.                 &itmask_metal[i-1],
  1649.                  imchem, igrgr, ipcont, tmcool,
  1650.                  k125a, k129a, k130a, k131a, k132a,
  1651.                  k133a, k134a, k135a, k136a, k137a,
  1652.                  k148a, k149a, k150a, k151a, k152a,
  1653.                  k153a,
  1654.                  kz15a, kz16a, kz17a, kz18a, kz19a,
  1655.                  kz20a, kz21a, kz22a, kz23a, kz24a,
  1656.                  kz25a, kz26a, kz27a, kz28a, kz29a,
  1657.                  kz30a, kz31a, kz32a, kz33a, kz34a,
  1658.                  kz35a, kz36a, kz37a, kz38a, kz39a,
  1659.                  kz40a, kz41a, kz42a, kz43a, kz44a,
  1660.                  kz45a, kz46a, kz47a, kz48a, kz49a,
  1661.                  kz50a, kz51a, kz52a, kz53a, kz54a,
  1662.                  cieY06a,
  1663.                  LH2_N, LH2_Size,
  1664.                  LH2_D, LH2_T, LH2_H,
  1665.                  LH2_dD, LH2_dT, LH2_dH, LH2_L,
  1666.                  LHD_N, LHD_Size,
  1667.                  LHD_D, LHD_T, LHD_H,
  1668.                  LHD_dD, LHD_dT, LHD_dH, LHD_L,
  1669.                  LCI_N, LCI_Size,
  1670.                  LCI_D, LCI_T, LCI_H,
  1671.                  LCI_dD, LCI_dT, LCI_dH, LCI_L,
  1672.                  LCII_N, LCII_Size,
  1673.                  LCII_D, LCII_T, LCII_H,
  1674.                  LCII_dD, LCII_dT, LCII_dH, LCII_L,
  1675.                  LOI_N, LOI_Size,
  1676.                  LOI_D, LOI_T, LOI_H,
  1677.                  LOI_dD, LOI_dT, LOI_dH, LOI_L,
  1678.                  LCO_N, LCO_Size,
  1679.                  LCO_D, LCO_T, LCO_H,
  1680.                  LCO_dD, LCO_dT, LCO_dH, LCO_L,
  1681.                  LOH_N, LOH_Size,
  1682.                  LOH_D, LOH_T, LOH_H,
  1683.                  LOH_dD, LOH_dT, LOH_dH, LOH_L,
  1684.                  LH2O_N, LH2O_Size,
  1685.                  LH2O_D, LH2O_T, LH2O_H,
  1686.                  LH2O_dD, LH2O_dT, LH2O_dH, LH2O_L,
  1687.                  alphap_N, alphap_Size,
  1688.                  alphap_D, alphap_T, alphap_dD, alphap_dT,
  1689.                  alphap_Data,
  1690.                  immulti, imabund, idspecies, itdmulti, idsub,
  1691.                  &metal_loc(i-1,j-1,k-1), &metal_C13(i-1,j-1,k-1), &metal_C20(i-1,j-1,k-1),
  1692.                  &metal_C25(i-1,j-1,k-1), &metal_C30(i-1,j-1,k-1), &metal_F13(i-1,j-1,k-1),
  1693.                  &metal_F15(i-1,j-1,k-1), &metal_F50(i-1,j-1,k-1), &metal_F80(i-1,j-1,k-1),
  1694.                  &metal_P170(i-1,j-1,k-1), &metal_P200(i-1,j-1,k-1), &metal_Y19(i-1,j-1,k-1),
  1695.                  SN0_N,
  1696.                  SN0_fSiM, SN0_fFeM, SN0_fMg2SiO4, SN0_fMgSiO3,
  1697.                  SN0_fFe3O4, SN0_fAC, SN0_fSiO2D, SN0_fMgO,
  1698.                  SN0_fFeS, SN0_fAl2O3,
  1699.                  SN0_freforg, SN0_fvolorg, SN0_fH2Oice,
  1700.                  SN0_r0SiM, SN0_r0FeM, SN0_r0Mg2SiO4, SN0_r0MgSiO3,
  1701.                  SN0_r0Fe3O4, SN0_r0AC, SN0_r0SiO2D, SN0_r0MgO,
  1702.                  SN0_r0FeS, SN0_r0Al2O3,
  1703.                  SN0_r0reforg, SN0_r0volorg, SN0_r0H2Oice,
  1704.                  gr_N, gr_Size, gr_dT, gr_Td,
  1705.                  SN0_kpSiM, SN0_kpFeM, SN0_kpMg2SiO4, SN0_kpMgSiO3,
  1706.                  SN0_kpFe3O4, SN0_kpAC, SN0_kpSiO2D, SN0_kpMgO,
  1707.                  SN0_kpFeS, SN0_kpAl2O3,
  1708.                  SN0_kpreforg, SN0_kpvolorg, SN0_kpH2Oice,
  1709.                  h2dustSa, h2dustCa, grogra,
  1710.                  &tSiM[i-1], &tFeM[i-1], &tMg2SiO4[i-1], &tMgSiO3[i-1], &tFe3O4[i-1],
  1711.                  &tAC[i-1], &tSiO2D[i-1], &tMgO[i-1], &tFeS[i-1], &tAl2O3[i-1],
  1712.                  &treforg[i-1], &tvolorg[i-1], &tH2Oice[i-1],
  1713.                  gasgr2a, gamma_isrf2a,
  1714.                  &imp_eng[i-1],
  1715.                  idissHDI, &kdissHDI(i-1,j-1,k-1), iionZ, &kphCI(i-1,j-1,k-1), &kphOI(i-1,j-1,k-1),
  1716.                  idissZ, &kdissCO(i-1,j-1,k-1), &kdissOH(i-1,j-1,k-1), &kdissH2O(i-1,j-1,k-1),
  1717.                  iuseH2shield
  1718.                 );
  1719.  
  1720.                 for (jsp = 1; jsp<=(nsp); jsp++) {
  1721.                   dspj = eps * dsp[idsp[jsp-1]-1];
  1722.                   for (isp = 1; isp<=(nsp); isp++) {
  1723.                     if(isp == jsp)  {
  1724.                       dsp1[idsp[isp-1]-1] = dsp[idsp[isp-1]-1] + dspj;
  1725.                     } else {
  1726.                       dsp1[idsp[isp-1]-1] = dsp[idsp[isp-1]-1];
  1727.                     }
  1728.                   }
  1729.  
  1730.                    FORTRAN_NAME(lookup_cool_rates0d)(&dtit[i-1],
  1731.                   &d(i-1,j-1,k-1), &u(i-1,j-1,k-1), &v(i-1,j-1,k-1), &w(i-1,j-1,k-1),
  1732.                   &nsp, dsp1.data(), dspdot1.data(), nratec,
  1733.                   iexpand, ispecies, imetal, imcool,
  1734.                   idust, idustall, idustfield, idustrec,
  1735.                   ih2co, ipiht, &iter, igammah,
  1736.                   aye, temstart, temend, z_solar, fgr,
  1737.                   utem, uxyz, uaye, urho, utim,
  1738.                   gamma, fh,
  1739.                   ceHIa, ceHeIa, ceHeIIa, ciHIa, ciHeIa,
  1740.                   ciHeISa, ciHeIIa, reHIIa, reHeII1a,
  1741.                   reHeII2a, reHeIIIa, brema, compa, gammaha,
  1742.                   isrf, regra, gamma_isrfa, comp_xraya, comp_temp,
  1743.                   piHI, piHeI, piHeII, &comp1, &comp2,
  1744.                   &metal(i-1,j-1,k-1), &dust(i-1,j-1,k-1),
  1745.                   hyd01ka, h2k01a, vibha, rotha, rotla,
  1746.                   &hyd01k[i-1], &h2k01[i-1], &vibh[i-1], &roth[i-1], &rotl[i-1],
  1747.                   gpldla, gphdla, &gpldl[i-1], &gphdl[i-1],
  1748.                   hdltea, hdlowa, &hdlte[i-1], &hdlow[i-1],
  1749.                   gaHIa, gaH2a, gaHea,
  1750.                   gaHpa, gaela,
  1751.                   h2ltea, gasgra,
  1752.                   &ceHI[i-1], &ceHeI[i-1], &ceHeII[i-1], &ciHI[i-1], &ciHeI[i-1],
  1753.                   &ciHeIS[i-1], &ciHeII[i-1],
  1754.                   &reHII[i-1], &reHeII1[i-1], &reHeII2[i-1], &reHeIII[i-1], &brem[i-1],
  1755.                   &indixe[i-1], &t1[i-1], &t2[i-1], &logtem[i-1], &tdef[i-1], &edot[i-1],
  1756.                   &tgas[i-1], &tgasold[i-1], &mmw[i-1], &p2d[i-1], &tdust[i-1], &metallicity[i-1],
  1757.                   &dust2gas[i-1], &rhoH[i-1], &mynh[i-1], &myde[i-1],
  1758.                   &gammaha_eff[i-1], &gasgr_tdust[i-1], &regr[i-1],
  1759.                   iradshield, avgsighi, avgsighei, avgsigheii,
  1760.                   iradtrans, irt_honly,
  1761.                   ih2optical, iciecool, ih2cr, ihdcr, ciecoa, &cieco[i-1],
  1762.                   icmbTfloor, iClHeat, clEleFra,
  1763.                   priGridRank, priGridDim,
  1764.                   priPar1, priPar2, priPar3, priPar4, priPar5,
  1765.                   priDataSize, priCooling, priHeating, priMMW,
  1766.                   metGridRank, metGridDim,
  1767.                   metPar1, metPar2, metPar3, metPar4, metPar5,
  1768.                   metDataSize, metCooling, metHeating, clnew,
  1769.                   iVheat, iMheat, &Vheat(i-1,j-1,k-1), &Mheat(i-1,j-1,k-1),
  1770.                   iTfloor, Tfloor_scalar, &Tfloor(i-1,j-1,k-1),
  1771.                   iisrffield, isrf_habing,
  1772.                   ithreebody, &anydust, iH2shield,
  1773.                   k1a, k2a, k3a, k4a, k5a, k6a, k7a, k8a, k9a, k10a,
  1774.                   k11a, k12a, k13a, k13dda, k14a, k15a, k16a,
  1775.                   k17a, k18a, k19a, k22a,
  1776.                   k24, k25, k26, k27, k28, k29, k30, k31,
  1777.                   k50a, k51a, k52a, k53a, k54a, k55a, k56a,
  1778.                   k57a, k58a, ndratec, dtemstart, dtemend, h2dusta,
  1779.                   ncrna, ncrd1a, ncrd2a,
  1780.                   &h2dust[i-1], &ncrn[i-1], &ncrd1[i-1], &ncrd2[i-1],
  1781.                   &dom, &coolunit, &tbase1, &xbase1, &dx_cgs, &c_ljeans,
  1782.                   &kphHI(i-1,j-1,k-1), &kphHeI(i-1,j-1,k-1), &kphHeII(i-1,j-1,k-1), &kdissH2I(i-1,j-1,k-1),
  1783.                   &photogamma(i-1,j-1,k-1), &xH2shield(i-1,j-1,k-1), &chunit, &itmask_nr[i-1],
  1784.                   &itmask_metal[i-1],
  1785.                    imchem, igrgr, ipcont, tmcool,
  1786.                    k125a, k129a, k130a, k131a, k132a,
  1787.                    k133a, k134a, k135a, k136a, k137a,
  1788.                    k148a, k149a, k150a, k151a, k152a,
  1789.                    k153a,
  1790.                    kz15a, kz16a, kz17a, kz18a, kz19a,
  1791.                    kz20a, kz21a, kz22a, kz23a, kz24a,
  1792.                    kz25a, kz26a, kz27a, kz28a, kz29a,
  1793.                    kz30a, kz31a, kz32a, kz33a, kz34a,
  1794.                    kz35a, kz36a, kz37a, kz38a, kz39a,
  1795.                    kz40a, kz41a, kz42a, kz43a, kz44a,
  1796.                    kz45a, kz46a, kz47a, kz48a, kz49a,
  1797.                    kz50a, kz51a, kz52a, kz53a, kz54a,
  1798.                    cieY06a,
  1799.                    LH2_N, LH2_Size,
  1800.                    LH2_D, LH2_T, LH2_H,
  1801.                    LH2_dD, LH2_dT, LH2_dH, LH2_L,
  1802.                    LHD_N, LHD_Size,
  1803.                    LHD_D, LHD_T, LHD_H,
  1804.                    LHD_dD, LHD_dT, LHD_dH, LHD_L,
  1805.                    LCI_N, LCI_Size,
  1806.                    LCI_D, LCI_T, LCI_H,
  1807.                    LCI_dD, LCI_dT, LCI_dH, LCI_L,
  1808.                    LCII_N, LCII_Size,
  1809.                    LCII_D, LCII_T, LCII_H,
  1810.                    LCII_dD, LCII_dT, LCII_dH, LCII_L,
  1811.                    LOI_N, LOI_Size,
  1812.                    LOI_D, LOI_T, LOI_H,
  1813.                    LOI_dD, LOI_dT, LOI_dH, LOI_L,
  1814.                    LCO_N, LCO_Size,
  1815.                    LCO_D, LCO_T, LCO_H,
  1816.                    LCO_dD, LCO_dT, LCO_dH, LCO_L,
  1817.                    LOH_N, LOH_Size,
  1818.                    LOH_D, LOH_T, LOH_H,
  1819.                    LOH_dD, LOH_dT, LOH_dH, LOH_L,
  1820.                    LH2O_N, LH2O_Size,
  1821.                    LH2O_D, LH2O_T, LH2O_H,
  1822.                    LH2O_dD, LH2O_dT, LH2O_dH, LH2O_L,
  1823.                    alphap_N, alphap_Size,
  1824.                    alphap_D, alphap_T, alphap_dD, alphap_dT,
  1825.                    alphap_Data,
  1826.                    immulti, imabund, idspecies, itdmulti, idsub,
  1827.                    &metal_loc(i-1,j-1,k-1), &metal_C13(i-1,j-1,k-1), &metal_C20(i-1,j-1,k-1),
  1828.                    &metal_C25(i-1,j-1,k-1), &metal_C30(i-1,j-1,k-1), &metal_F13(i-1,j-1,k-1),
  1829.                    &metal_F15(i-1,j-1,k-1), &metal_F50(i-1,j-1,k-1), &metal_F80(i-1,j-1,k-1),
  1830.                    &metal_P170(i-1,j-1,k-1), &metal_P200(i-1,j-1,k-1), &metal_Y19(i-1,j-1,k-1),
  1831.                    SN0_N,
  1832.                    SN0_fSiM, SN0_fFeM, SN0_fMg2SiO4, SN0_fMgSiO3,
  1833.                    SN0_fFe3O4, SN0_fAC, SN0_fSiO2D, SN0_fMgO,
  1834.                    SN0_fFeS, SN0_fAl2O3,
  1835.                    SN0_freforg, SN0_fvolorg, SN0_fH2Oice,
  1836.                    SN0_r0SiM, SN0_r0FeM, SN0_r0Mg2SiO4, SN0_r0MgSiO3,
  1837.                    SN0_r0Fe3O4, SN0_r0AC, SN0_r0SiO2D, SN0_r0MgO,
  1838.                    SN0_r0FeS, SN0_r0Al2O3,
  1839.                    SN0_r0reforg, SN0_r0volorg, SN0_r0H2Oice,
  1840.                    gr_N, gr_Size, gr_dT, gr_Td,
  1841.                    SN0_kpSiM, SN0_kpFeM, SN0_kpMg2SiO4, SN0_kpMgSiO3,
  1842.                    SN0_kpFe3O4, SN0_kpAC, SN0_kpSiO2D, SN0_kpMgO,
  1843.                    SN0_kpFeS, SN0_kpAl2O3,
  1844.                    SN0_kpreforg, SN0_kpvolorg, SN0_kpH2Oice,
  1845.                    h2dustSa, h2dustCa, grogra,
  1846.                    &tSiM[i-1], &tFeM[i-1], &tMg2SiO4[i-1], &tMgSiO3[i-1], &tFe3O4[i-1],
  1847.                    &tAC[i-1], &tSiO2D[i-1], &tMgO[i-1], &tFeS[i-1], &tAl2O3[i-1],
  1848.                    &treforg[i-1], &tvolorg[i-1], &tH2Oice[i-1],
  1849.                    gasgr2a, gamma_isrf2a,
  1850.                    &imp_eng[i-1],
  1851.                    idissHDI, &kdissHDI(i-1,j-1,k-1), iionZ, &kphCI(i-1,j-1,k-1), &kphOI(i-1,j-1,k-1),
  1852.                    idissZ, &kdissCO(i-1,j-1,k-1), &kdissOH(i-1,j-1,k-1), &kdissH2O(i-1,j-1,k-1),
  1853.                    iuseH2shield
  1854.                   );
  1855.  
  1856.                   for (isp = 1; isp<=(nsp); isp++) {
  1857.                     if ( (dsp[idsp[isp-1]-1]==0.e0)
  1858.                      &&  (dspdot1[idsp[isp-1]-1]
  1859.                      ==  dspdot[idsp[isp-1]-1]) )  {
  1860.                       der(idsp[isp-1]-1,idsp[jsp-1]-1) = 0.e0;
  1861.                     } else {
  1862.                       der(idsp[isp-1]-1,idsp[jsp-1]-1) =
  1863.                          (dspdot1[idsp[isp-1]-1]
  1864.                          - dspdot[idsp[isp-1]-1]) / dspj;
  1865.                     }
  1866.                   }
  1867.  
  1868.                 }
  1869.  
  1870.                 for (isp = 1; isp<=(nsp); isp++) {
  1871.                   for (jsp = 1; jsp<=(nsp); jsp++) {
  1872.                     if(isp == jsp)  {
  1873.                       mtrx(isp-1,jsp-1) = 1.e0 - dtit[i-1]
  1874.                          * der(idsp[isp-1]-1,idsp[jsp-1]-1);
  1875.                     } else {
  1876.                       mtrx(isp-1,jsp-1) =      - dtit[i-1]
  1877.                          * der(idsp[isp-1]-1,idsp[jsp-1]-1);
  1878.                     }
  1879.                   }
  1880.                 }
  1881.  
  1882.                 for (isp = 1; isp<=(nsp); isp++) {
  1883.                   vec[isp-1] = dspdot[idsp[isp-1]-1] * dtit[i-1]
  1884.                          - ddsp[idsp[isp-1]-1];
  1885.                 }
  1886.  
  1887.                 // to get more accuracy
  1888.                 for (isp = 1; isp<=(nsp); isp++) {
  1889.                   vec[isp-1] = vec[isp-1]/d(i-1,j-1,k-1);
  1890.                 }
  1891.  
  1892.                  FORTRAN_NAME(gaussj_g)(&nsp, mtrx.data(), vec.data(), &ierror);
  1893.                 if(ierror == 1)  {
  1894.                   goto label_9998;
  1895.                 }
  1896.  
  1897.                 // multiply with density again
  1898.                 for (isp = 1; isp<=(nsp); isp++) {
  1899.                   vec[isp-1] = vec[isp-1]*d(i-1,j-1,k-1);
  1900.                 }
  1901.  
  1902.                 for (isp = 1; isp<=(nsp); isp++) {
  1903.                   ddsp[idsp[isp-1]-1] = ddsp[idsp[isp-1]-1] + vec[isp-1];
  1904.                   dsp[idsp[isp-1]-1]  = dsp[idsp[isp-1]-1]  + vec[isp-1];
  1905.                 }
  1906.  
  1907.                 if (imp_eng[i-1] == 1)  {
  1908.                   if( ((*ispecies) > 0)  &&  ((*icool) == 1) )  {
  1909.                     for (isp = 1; isp<=(nsp); isp++) {
  1910.                       if ( (dsp[idsp[isp-1]-1] != dsp[idsp[isp-1]-1])
  1911.                        ||  (dsp[idsp[isp-1]-1] <= 0.) )  {
  1912.                         ierror = 1;
  1913.                         goto label_9997;
  1914.                       }
  1915.                     }
  1916.                   }
  1917.                 }
  1918.  
  1919.                 err_max = 0.e0;
  1920.                 for (isp = 1; isp<=(nsp); isp++) {
  1921.                   if(dsp[idsp[isp-1]-1] > tiny8)  {
  1922.                     err = grackle::impl::dabs(vec[isp-1] / dsp[idsp[isp-1]-1]);
  1923.                   } else {
  1924.                     err = 0.e0;
  1925.                   }
  1926.                   if(err > err_max)  {
  1927.                     err_max = err;
  1928.                   }
  1929.                 }
  1930.  
  1931.                 itr=itr+1;
  1932.               }
  1933.  
  1934. label_9998:
  1935. label_9997:
  1936. label_9996:
  1937.  
  1938.  
  1939.               // Check if the fractions are valid after an iteration
  1940.  
  1941.               if( ((*ispecies) > 0)  &&  ((*icool) == 1) )  {
  1942.                 for (isp = 1; isp<=(nsp); isp++) {
  1943.                   if ( (dsp[idsp[isp-1]-1] != dsp[idsp[isp-1]-1])
  1944.                    ||  (dsp[idsp[isp-1]-1] <= 0.) )  {
  1945.                     ierror = 1;
  1946.                   }
  1947.                 }
  1948.               }
  1949.               if(ierror == 1)  {
  1950.                 dtit[i-1] = 0.5e0*dtit[i-1];
  1951.               }
  1952.  
  1953.               itr_time=itr_time+1;
  1954.             }
  1955.  
  1956.             if ( (*ispecies) > 0 )  {
  1957.               de(i-1,j-1,k-1)      = dsp[ 1-1];
  1958.               HI(i-1,j-1,k-1)      = dsp[ 2-1];
  1959.               HII(i-1,j-1,k-1)     = dsp[ 3-1];
  1960.               HeI(i-1,j-1,k-1)     = dsp[ 4-1];
  1961.               HeII(i-1,j-1,k-1)    = dsp[ 5-1];
  1962.               HeIII(i-1,j-1,k-1)   = dsp[ 6-1];
  1963.             }
  1964.             if ( (*ispecies) > 1 )  {
  1965.               HM(i-1,j-1,k-1)      = dsp[ 7-1];
  1966.               H2I(i-1,j-1,k-1)     = dsp[ 8-1];
  1967.               H2II(i-1,j-1,k-1)    = dsp[ 9-1];
  1968.             }
  1969.             if ( (*ispecies) > 2 )  {
  1970.               DI(i-1,j-1,k-1)      = dsp[10-1];
  1971.               DII(i-1,j-1,k-1)     = dsp[11-1];
  1972.               HDI(i-1,j-1,k-1)     = dsp[12-1];
  1973.             }
  1974.             if ( (*ispecies) > 3 )  {
  1975.               DM(i-1,j-1,k-1)      = dsp[13-1];
  1976.               HDII(i-1,j-1,k-1)    = dsp[14-1];
  1977.               HeHII(i-1,j-1,k-1)   = dsp[15-1];
  1978.             }
  1979.             if ( itmask_metal[i-1] != MASK_FALSE )  {
  1980.               if ( (*imchem) == 1 )  {
  1981.                 CI(i-1,j-1,k-1)      = dsp[16-1];
  1982.                 CII(i-1,j-1,k-1)     = dsp[17-1];
  1983.                 CO(i-1,j-1,k-1)      = dsp[18-1];
  1984.                 CO2(i-1,j-1,k-1)     = dsp[19-1];
  1985.                 OI(i-1,j-1,k-1)      = dsp[20-1];
  1986.                 OH(i-1,j-1,k-1)      = dsp[21-1];
  1987.                 H2O(i-1,j-1,k-1)     = dsp[22-1];
  1988.                 O2(i-1,j-1,k-1)      = dsp[23-1];
  1989.                 SiI(i-1,j-1,k-1)     = dsp[24-1];
  1990.                 SiOI(i-1,j-1,k-1)    = dsp[25-1];
  1991.                 SiO2I(i-1,j-1,k-1)   = dsp[26-1];
  1992.                 CH(i-1,j-1,k-1)      = dsp[27-1];
  1993.                 CH2(i-1,j-1,k-1)     = dsp[28-1];
  1994.                 COII(i-1,j-1,k-1)    = dsp[29-1];
  1995.                 OII(i-1,j-1,k-1)     = dsp[30-1];
  1996.                 OHII(i-1,j-1,k-1)    = dsp[31-1];
  1997.                 H2OII(i-1,j-1,k-1)   = dsp[32-1];
  1998.                 H3OII(i-1,j-1,k-1)   = dsp[33-1];
  1999.                 O2II(i-1,j-1,k-1)    = dsp[34-1];
  2000.                 if ( ( (*igrgr) == 1 )  ||  ( (*idsub) == 1 ) )  {
  2001.                   if ((*idspecies) > 0)  {
  2002.                     Mg(i-1,j-1,k-1)      = dsp[35-1];
  2003.                   }
  2004.                   if ((*idspecies) > 1)  {
  2005.                     Al(i-1,j-1,k-1)      = dsp[36-1];
  2006.                     S(i-1,j-1,k-1)       = dsp[37-1];
  2007.                     Fe(i-1,j-1,k-1)      = dsp[38-1];
  2008.                   }
  2009.                 }
  2010.               }
  2011.               if ( ( (*igrgr) == 1 )  ||  ( (*idsub) == 1 ) )  {
  2012.                 if ((*idspecies) > 0)  {
  2013.                   MgSiO3(i-1,j-1,k-1)  = dsp[39-1];
  2014.                   AC(i-1,j-1,k-1)      = dsp[40-1];
  2015.                 }
  2016.                 if ((*idspecies) > 1)  {
  2017.                   SiM(i-1,j-1,k-1)     = dsp[41-1];
  2018.                   FeM(i-1,j-1,k-1)     = dsp[42-1];
  2019.                   Mg2SiO4(i-1,j-1,k-1) = dsp[43-1];
  2020.                   Fe3O4(i-1,j-1,k-1)   = dsp[44-1];
  2021.                   SiO2D(i-1,j-1,k-1)   = dsp[45-1];
  2022.                   MgO(i-1,j-1,k-1)     = dsp[46-1];
  2023.                   FeS(i-1,j-1,k-1)     = dsp[47-1];
  2024.                   Al2O3(i-1,j-1,k-1)   = dsp[48-1];
  2025.                 }
  2026.                 if ((*idspecies) > 2)  {
  2027.                   reforg(i-1,j-1,k-1)  = dsp[49-1];
  2028.                   volorg(i-1,j-1,k-1)  = dsp[50-1];
  2029.                   H2Oice(i-1,j-1,k-1)  = dsp[51-1];
  2030.                 }
  2031.               }
  2032.             }
  2033.  
  2034.             e(i-1,j-1,k-1)     = dsp[i_eng-1];
  2035.  
  2036.             idsp.clear();
  2037.             vec.clear();
  2038.             mtrx = grackle::impl::View<double**>();
  2039.             mtrx_data_.clear();
  2040.  
  2041.           }
  2042.         }
  2043.  
  2044.       }
  2045.  
  2046.       // return itmask
  2047.       for (i = (*is) + 1; i<=((*ie) + 1); i++) {
  2048.         itmask[i-1] = itmask_tmp[i-1];
  2049.       }
  2050.  
  2051.       // Add the timestep to the elapsed time for each cell and find
  2052.       //  minimum elapsed time step in this row
  2053.  
  2054.       ttmin = huge8;
  2055.       for (i = (*is) + 1; i<=((*ie) + 1); i++) {
  2056.         ttot[i-1] = std::fmin(ttot[i-1] + dtit[i-1], (*dt));
  2057.         if (std::fabs((*dt)-ttot[i-1]) <
  2058.              tolerance*(*dt)) { itmask[i-1] = MASK_FALSE; }
  2059.         if (ttot[i-1]<ttmin) { ttmin = ttot[i-1]; }
  2060.       }
  2061.  
  2062.       // If all cells are done (on this slice), then exit
  2063.  
  2064.       if (std::fabs((*dt)-ttmin) < tolerance*(*dt)) { goto label_9999; }
  2065.  
  2066.       // Next subcycle iteration
  2067.  
  2068.     }
  2069.  
  2070. label_9999:
  2071.  
  2072.     // Abort if iteration count exceeds maximum
  2073.  
  2074.     if (iter > (*itmax))  {
  2075.       OMP_PRAGMA_CRITICAL
  2076.       {
  2077.         printf("inside if statement solve rate cool: %d %d\n", (*is), (*ie));
  2078.         eprintf("MULTI_COOL iter >  %d  at j,k = %d %d\n", (*itmax), j, k);
  2079.         printf("FATAL error (2) in MULTI_COOL\n");
  2080.   //_// PORT:             write(0,'(" dt = ",1pe10.3," ttmin = ",1pe10.3)') dt, ttmin
  2081.   //_// PORT:             write(0,'((16(1pe8.1)))') (dtit(i),i=is+1,ie+1)
  2082.   //_// PORT:             write(0,'((16(1pe8.1)))') (ttot(i),i=is+1,ie+1)
  2083.   //_// PORT:             write(0,'((16(1pe8.1)))') (edot(i),i=is+1,ie+1)
  2084.   //_// PORT:             write(0,'((16(I3)))') (itmask(i),i=is+1,ie+1)
  2085.  
  2086.         if ((*exititmax) == 1)  {
  2087.           (*ierr) = 0;
  2088.         }
  2089.       }
  2090.       // WARNING_MESSAGE
  2091.     }
  2092.  
  2093.     if (iter > (*itmax)/2)  {
  2094.       OMP_PRAGMA_CRITICAL
  2095.       {
  2096.         eprintf("MULTI_COOL iter,j,k = %d %d %d\n", iter, j, k);
  2097.       }
  2098.     }
  2099.    
  2100.     // Next j,k
  2101.   }
  2102.   //_// PORT: #ifdef _OPENMP
  2103.   //_// PORT: !$omp end parallel do
  2104.   //_// PORT: #endif
  2105.  
  2106.   // If an error has been produced, return now.
  2107.  
  2108.   if ((*ierr) == 0)  {
  2109.     return;
  2110.   }
  2111.  
  2112.   // Convert densities back to comoving from proper
  2113.  
  2114.   if ((*iexpand) == 1)  {
  2115.  
  2116.     factor = (gr_float)(std::pow((*aye),3) );
  2117.  
  2118.      FORTRAN_NAME(scale_fields_g)(
  2119.       ispecies, imetal, idustfield, imchem,
  2120.       imabund, idspecies, immulti, igrgr, idsub, &factor,
  2121.       is, ie, js, je, ks, ke, in, jn, kn,
  2122.       d.data(), de.data(), HI.data(), HII.data(), HeI.data(), HeII.data(), HeIII.data(),
  2123.       HM.data(), H2I.data(), H2II.data(), DI.data(), DII.data(), HDI.data(), DM.data(), HDII.data(), HeHII.data(),
  2124.       metal.data(), dust.data(),
  2125.       CI.data(), CII.data(), CO.data(), CO2.data(),
  2126.       OI.data(), OH.data(), H2O.data(), O2.data(),
  2127.       SiI.data(), SiOI.data(), SiO2I.data(),
  2128.       CH.data(), CH2.data(), COII.data(), OII.data(),
  2129.       OHII.data(), H2OII.data(), H3OII.data(), O2II.data(),
  2130.       Mg.data(), Al.data(), S.data(), Fe.data(),
  2131.       SiM.data(), FeM.data(), Mg2SiO4.data(), MgSiO3.data(), Fe3O4.data(),
  2132.       AC.data(), SiO2D.data(), MgO.data(), FeS.data(), Al2O3.data(),
  2133.       reforg.data(), volorg.data(), H2Oice.data(),
  2134.       metal_loc.data(), metal_C13.data(), metal_C20.data(), metal_C25.data(), metal_C30.data(),
  2135.       metal_F13.data(), metal_F15.data(), metal_F50.data(), metal_F80.data(),
  2136.       metal_P170.data(), metal_P200.data(), metal_Y19.data());
  2137.  
  2138.   }
  2139.  
  2140.   if ((*ispecies) > 0)  {
  2141.  
  2142.     // Correct the species to ensure consistency (i.e. type conservation)
  2143.  
  2144. #ifdef ABUNDANCE_CORRECTION
  2145.      FORTRAN_NAME(make_consistent_g)(de.data(), HI.data(), HII.data(), HeI.data(), HeII.data(), HeIII.data(),
  2146.                          HM.data(), H2I.data(), H2II.data(), DI.data(), DII.data(), HDI.data(), metal.data(), dust.data(),
  2147.                          d.data(), is, ie, js, je, ks, ke,
  2148.                          in, jn, kn, ispecies, imetal, fh, dtoh,
  2149.                         idustfield, imchem, igrgr, &dom,
  2150.                         DM.data(), HDII.data(), HeHII.data(),
  2151.                         CI.data(), CII.data(), CO.data(), CO2.data(),
  2152.                         OI.data(), OH.data(), H2O.data(), O2.data(),
  2153.                         SiI.data(), SiOI.data(), SiO2I.data(),
  2154.                         CH.data(), CH2.data(), COII.data(), OII.data(),
  2155.                         OHII.data(), H2OII.data(), H3OII.data(), O2II.data(),
  2156.                         Mg.data(), Al.data(), S.data(), Fe.data(),
  2157.                         SiM.data(), FeM.data(), Mg2SiO4.data(), MgSiO3.data(), Fe3O4.data(),
  2158.                         AC.data(), SiO2D.data(), MgO.data(), FeS.data(), Al2O3.data(),
  2159.                         reforg.data(), volorg.data(), H2Oice.data(),
  2160.                         immulti, imabund, idspecies, itdmulti, idsub,
  2161.                         metal_loc.data(),
  2162.                         metal_C13.data(), metal_C20.data(), metal_C25.data(), metal_C30.data(),
  2163.                         metal_F13.data(), metal_F15.data(), metal_F50.data(), metal_F80.data(),
  2164.                         metal_P170.data(), metal_P200.data(), metal_Y19.data(),
  2165.                         SN0_N,
  2166.                         SN0_XC, SN0_XO, SN0_XMg, SN0_XAl, SN0_XSi,
  2167.                         SN0_XS, SN0_XFe,
  2168.                         SN0_fC, SN0_fO, SN0_fMg, SN0_fAl, SN0_fSi,
  2169.                         SN0_fS, SN0_fFe
  2170.                           );
  2171.  
  2172.      FORTRAN_NAME(ceiling_species_g)(d.data(), de.data(), HI.data(), HII.data(), HeI.data(), HeII.data(), HeIII.data(),
  2173.                          HM.data(), H2I.data(), H2II.data(), DI.data(), DII.data(), HDI.data(), metal.data(), dust.data(),
  2174.                          is, ie, js, je, ks, ke,
  2175.                          in, jn, kn, ispecies, imetal, idustfield,
  2176.                         DM.data(), HDII.data(), HeHII.data(),
  2177.                         imabund, imchem, idspecies, immulti,
  2178.                         igrgr, idsub,
  2179.                         CI.data(), CII.data(), CO.data(), CO2.data(),
  2180.                         OI.data(), OH.data(), H2O.data(), O2.data(),
  2181.                         SiI.data(), SiOI.data(), SiO2I.data(),
  2182.                         CH.data(), CH2.data(), COII.data(), OII.data(),
  2183.                         OHII.data(), H2OII.data(), H3OII.data(), O2II.data(),
  2184.                         Mg.data(), Al.data(), S.data(), Fe.data(),
  2185.                         SiM.data(), FeM.data(), Mg2SiO4.data(), MgSiO3.data(), Fe3O4.data(),
  2186.                         AC.data(), SiO2D.data(), MgO.data(), FeS.data(), Al2O3.data(),
  2187.                         reforg.data(), volorg.data(), H2Oice.data(),
  2188.                         metal_loc.data(),
  2189.                         metal_C13.data(), metal_C20.data(), metal_C25.data(), metal_C30.data(),
  2190.                         metal_F13.data(), metal_F15.data(), metal_F50.data(), metal_F80.data(),
  2191.                         metal_P170.data(), metal_P200.data(), metal_Y19.data());
  2192. #endif
  2193.  
  2194.   }
  2195.  
  2196.   return;
  2197. }
  2198.  
  2199. }  // extern "C"
  2200.  
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!