From 24cfbc542d52ee405cd89cce0e5e045c7bc94f6e Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Tue, 22 Oct 2024 08:14:38 -0400 Subject: [PATCH 01/11] Standardize parameter declaration Previously, there were 2 different ways to declare parameters scattered throughout the codebase. This patch standardizes on a single approach that is easier to transcribe. --- src/clib/calc_grain_size_increment_1d.F | 8 +++---- src/clib/calc_tdust_1d_g.F | 29 +++++++++++-------------- src/clib/calc_tdust_3d_g.F | 3 +-- src/clib/calc_temp1d_cloudy_g.F | 7 +++--- src/clib/calc_temp_cloudy_g.F | 3 +-- src/clib/cool1d_multi_g.F | 9 ++++---- src/clib/lookup_cool_rates0d.F | 13 +++++------ src/clib/solve_rate_cool_g.F | 20 ++++++++--------- 8 files changed, 39 insertions(+), 53 deletions(-) diff --git a/src/clib/calc_grain_size_increment_1d.F b/src/clib/calc_grain_size_increment_1d.F index 160ed62d..76c922f3 100644 --- a/src/clib/calc_grain_size_increment_1d.F +++ b/src/clib/calc_grain_size_increment_1d.F @@ -457,9 +457,8 @@ subroutine calc_grain_size_increment_species_1d( real*8 dsp0, SN_sgsp, SN_kpsp real*8 SN_dsp0(SN0_N), SN_nsp0(SN0_N) real*8 drsp(in) - real*8 pi, mh - parameter (pi = pi_val) - parameter (mh = mass_h) + real*8, parameter :: pi = pi_val + real*8, parameter :: mh = mass_h ! debug real*8 SN_dsp(SN0_N), SN_msp(SN0_N), dsp1 integer iTd, iTd0 @@ -649,8 +648,7 @@ subroutine solve_cubic_equation(a, b, c, root) real*8 q, r, m real*8 th real*8 s,t - real*8 pi - parameter (pi = pi_val) + real*8, parameter :: pi = pi_val q = (a*a - 3.d0*b)/9.d0 r = (2.d0*a*a*a - 9.d0*a*b + 27.d0*c)/54.d0 diff --git a/src/clib/calc_tdust_1d_g.F b/src/clib/calc_tdust_1d_g.F index ba85091d..b6578456 100644 --- a/src/clib/calc_tdust_1d_g.F +++ b/src/clib/calc_tdust_1d_g.F @@ -65,18 +65,17 @@ subroutine calc_tdust_1d_g( ! Parameters integer idspecies - real*8 t_subl - parameter(t_subl = 1.5e3_DKIND) ! grain sublimation temperature - real*8 radf - parameter(radf = 4._DKIND * sigma_sb) - real*8 kgr1 - parameter(kgr1 = 4.0e-4_DKIND / 0.009387d0) + ! grain sublimation temperature + real*8, parameter :: t_subl = 1.5e3_DKIND + real*8, parameter :: radf = 4._DKIND * sigma_sb + real*8, parameter :: kgr1 = 4.0e-4_DKIND / 0.009387d0 !! should be normalized with local fgr. [GC20200701] - real*8 tol, bi_tol, minpert, gamma_isrf(in) - parameter(tol = 1.e-5_DKIND, bi_tol = 1.e-3_DKIND, - & minpert = 1.e-10_DKIND) - integer itmax, bi_itmax - parameter(itmax = 50, bi_itmax = 30) + real*8 gamma_isrf(in) + real*8, parameter :: tol = 1.e-5_DKIND + real*8, parameter :: bi_tol = 1.e-3_DKIND + real*8, parameter :: minpert = 1.e-10_DKIND + integer, parameter :: itmax = 50 + integer, parameter :: bi_itmax = 30 ! Locals @@ -384,9 +383,8 @@ subroutine calc_kappa_gr_g( ! Parameters - real*8 kgr1, kgr200 - parameter(kgr1 = 4.0e-4_DKIND / 0.009387d0 - & , kgr200 = 16.0_DKIND / 0.009387d0) + real*8, parameter :: kgr1 = 4.0e-4_DKIND / 0.009387d0 + real*8, parameter :: kgr200 = 16.0_DKIND / 0.009387d0 !! should be normalized with local fgr. [GC20200701] !! This value is valid only for Td < 50 K (Omukai 2000). @@ -498,8 +496,7 @@ subroutine calc_gr_balance_g( ! Parameters - real*8 radf - parameter(radf = 4._DKIND * sigma_sb) + real*8, parameter :: radf = 4._DKIND * sigma_sb ! Locals diff --git a/src/clib/calc_tdust_3d_g.F b/src/clib/calc_tdust_3d_g.F index 92b94e12..feaa0f15 100644 --- a/src/clib/calc_tdust_3d_g.F +++ b/src/clib/calc_tdust_3d_g.F @@ -187,8 +187,7 @@ subroutine calc_tdust_3d_g( ! Parameters - real*8 mh - parameter (mh = mass_h) + real*8, parameter :: mh = mass_h ! Locals diff --git a/src/clib/calc_temp1d_cloudy_g.F b/src/clib/calc_temp1d_cloudy_g.F index 0495d955..4aa4dd2a 100644 --- a/src/clib/calc_temp1d_cloudy_g.F +++ b/src/clib/calc_temp1d_cloudy_g.F @@ -82,10 +82,9 @@ subroutine calc_temp1d_cloudy_g(d, metal, e, rhoH, ! Parameters - integer ti_max - real*8 mu_metal - parameter (mu_metal = 16._DKIND) ! approx. mean molecular weight of metals - parameter (ti_max = 20) + ! approx. mean molecular weight of metals + real*8, parameter :: mu_metal = 16._DKIND + integer, parameter :: ti_max = 20 ! Locals diff --git a/src/clib/calc_temp_cloudy_g.F b/src/clib/calc_temp_cloudy_g.F index ac80f6f2..6f94bbae 100644 --- a/src/clib/calc_temp_cloudy_g.F +++ b/src/clib/calc_temp_cloudy_g.F @@ -86,8 +86,7 @@ subroutine calc_temp_cloudy_g(d, e, metal, temperature, ! Parameters - real*8 mh - parameter (mh = mass_h) + real*8, parameter :: mh = mass_h ! Locals diff --git a/src/clib/cool1d_multi_g.F b/src/clib/cool1d_multi_g.F index be4f9090..a4c7fa5a 100644 --- a/src/clib/cool1d_multi_g.F +++ b/src/clib/cool1d_multi_g.F @@ -208,11 +208,10 @@ subroutine cool1d_multi_g( ! Parameters - integer ti_max - real*8 mh, mu_metal - parameter (mh = mass_h) !DPC - parameter (mu_metal = 16._DKIND) ! approx. mean molecular weight of metals - parameter (ti_max = 20) + real*8, parameter :: mh = mass_h !DPC + ! approx. mean molecular weight of metals + real*8, parameter :: mu_metal = 16._DKIND + integer, parameter :: ti_max = 20 ! Locals diff --git a/src/clib/lookup_cool_rates0d.F b/src/clib/lookup_cool_rates0d.F index 39ad19fd..9db71086 100644 --- a/src/clib/lookup_cool_rates0d.F +++ b/src/clib/lookup_cool_rates0d.F @@ -249,21 +249,18 @@ subroutine lookup_cool_rates0d(output, dtit, ! Parameters - integer itmax - parameter (itmax = 10000) + integer, parameter :: itmax = 10000 #ifdef CONFIG_BFLOAT_4 - R_PREC tolerance - parameter (tolerance = 1.0e-05_RKIND) + R_PREC, parameter :: tolerance = 1.0e-05_RKIND #endif #ifdef CONFIG_BFLOAT_8 - R_PREC tolerance - parameter (tolerance = 1.0e-10_RKIND) + R_PREC, parameter :: tolerance = 1.0e-10_RKIND #endif - real*8 mh, pi - parameter (mh = mass_h, pi = pi_val) + real*8, parameter :: mh = mass_h + real*8, parameter :: pi = pi_val ! Locals diff --git a/src/clib/solve_rate_cool_g.F b/src/clib/solve_rate_cool_g.F index 6bc3b49c..7dacd865 100644 --- a/src/clib/solve_rate_cool_g.F +++ b/src/clib/solve_rate_cool_g.F @@ -444,17 +444,15 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de, ! Parameters #ifdef GRACKLE_FLOAT_4 - R_PREC tolerance - parameter (tolerance = 1.0e-05_RKIND) + R_PREC, parameter :: tolerance = 1.0e-05_RKIND #endif #ifdef GRACKLE_FLOAT_8 - R_PREC tolerance - parameter (tolerance = 1.0e-10_RKIND) + R_PREC, parameter :: tolerance = 1.0e-10_RKIND #endif - real*8 mh, pi - parameter (mh = mass_h, pi = pi_val) + real*8, parameter :: mh = mass_h + real*8, parameter :: pi = pi_val ! Locals @@ -2517,9 +2515,9 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, ! Parameters - real*8 everg, e24, e26 - parameter(everg = ev2erg, e24 = 13.6_DKIND, - & e26 = 24.6_DKIND) + real*8, parameter :: everg = ev2erg + real*8, parameter :: e24 = 13.6_DKIND + real*8, parameter :: e26 = 24.6_DKIND ! locals @@ -2596,8 +2594,8 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, real*8 d_Td(ndratec), d_Tg(nratec) integer idratec, iratec real*8 kd - real*8 fh, mh - parameter (mh = mass_h) !DPC + real*8 fh + real*8, parameter :: mh = mass_h !DPC ! locals for H2 self-shielding as WG+19 From 21fb999088c8555b3b2e54ca2ed41101b3f7c813 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Tue, 22 Oct 2024 09:01:09 -0400 Subject: [PATCH 02/11] standardized the format of precision-dependent tolerance --- src/clib/lookup_cool_rates0d.F | 6 ++---- src/clib/solve_rate_cool_g.F | 4 +--- 2 files changed, 3 insertions(+), 7 deletions(-) diff --git a/src/clib/lookup_cool_rates0d.F b/src/clib/lookup_cool_rates0d.F index 9db71086..4403da8e 100644 --- a/src/clib/lookup_cool_rates0d.F +++ b/src/clib/lookup_cool_rates0d.F @@ -251,11 +251,9 @@ subroutine lookup_cool_rates0d(output, dtit, integer, parameter :: itmax = 10000 -#ifdef CONFIG_BFLOAT_4 +#ifdef GRACKLE_FLOAT_4 R_PREC, parameter :: tolerance = 1.0e-05_RKIND -#endif - -#ifdef CONFIG_BFLOAT_8 +#else R_PREC, parameter :: tolerance = 1.0e-10_RKIND #endif diff --git a/src/clib/solve_rate_cool_g.F b/src/clib/solve_rate_cool_g.F index 7dacd865..0e859ee4 100644 --- a/src/clib/solve_rate_cool_g.F +++ b/src/clib/solve_rate_cool_g.F @@ -445,9 +445,7 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de, #ifdef GRACKLE_FLOAT_4 R_PREC, parameter :: tolerance = 1.0e-05_RKIND -#endif - -#ifdef GRACKLE_FLOAT_8 +#else R_PREC, parameter :: tolerance = 1.0e-10_RKIND #endif From 3f747da3aaa2741bb3224139cba4754d3c9398b7 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Sun, 27 Oct 2024 17:41:49 -0400 Subject: [PATCH 03/11] adding include-directive for the sake of consistency. --- src/clib/calc_grain_size_increment_1d.F | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/clib/calc_grain_size_increment_1d.F b/src/clib/calc_grain_size_increment_1d.F index 76c922f3..6bd59e26 100644 --- a/src/clib/calc_grain_size_increment_1d.F +++ b/src/clib/calc_grain_size_increment_1d.F @@ -434,6 +434,9 @@ subroutine calc_grain_size_increment_species_1d( implicit NONE + +#include "grackle_fortran_types.def" + ! input integer in, jn, kn, is, ie, j, k logical itmask(in) From ed0532807206cbbbe1b6a9c98862627fb542da57 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Tue, 29 Oct 2024 08:17:32 -0400 Subject: [PATCH 04/11] Convert most logicals to MASK_TYPE This patch converts nearly all usages of the ``logical`` type to the newly defined ``MASK_TYPE`` in the Fortran source files throughout Grackle. - The ``MASK_TYPE`` type is a custom datatype that just wraps a 32bit integer. - Because Fortran will not implicitly cast between logicals and integers we need to explicitly compare this value against ``MASK_TRUE`` and ``MASK_FALSE``. - For consistency with C semantics, if-statements that want to see if the mask is true always compare ``MASK_TYPE`` value is not equal to ``MASK_FALSE``. - The **only** ``logical`` variable that wasn't changed was the ``end_int`` variable that we pass into ``interpolate_3Dz_g``. We haven't touched this variable to avoid interfering with existing efforts to transcribe that function to C/C++ - This patch is extremely similar in spirit to GH PR grackle-project/grackle#227. But there a fewer minor distinctions: - I am actually more confident in the correctness of this patch (I used custom tools to automate the majority of these changes) - I converted more variables in this patch (that other PR **only** changed the ``itmask`` variable). I have explicitly confirmed that all tests continue to pass after the introduction of this PR --- src/clib/calc_grain_size_increment_1d.F | 8 +- src/clib/calc_tdust_1d_g.F | 44 +++--- src/clib/calc_tdust_3d_g.F | 10 +- src/clib/calc_temp1d_cloudy_g.F | 6 +- src/clib/calc_temp_cloudy_g.F | 4 +- src/clib/cool1d_cloudy_g.F | 4 +- src/clib/cool1d_cloudy_old_tables_g.F | 4 +- src/clib/cool1d_multi_g.F | 165 +++++++++++----------- src/clib/cool_multi_time_g.F | 4 +- src/clib/lookup_cool_rates0d.F | 34 +++-- src/clib/solve_rate_cool_g.F | 179 +++++++++++++----------- src/include/grackle_fortran_types.def | 6 + 12 files changed, 250 insertions(+), 218 deletions(-) diff --git a/src/clib/calc_grain_size_increment_1d.F b/src/clib/calc_grain_size_increment_1d.F index 6bd59e26..23b5f5b9 100644 --- a/src/clib/calc_grain_size_increment_1d.F +++ b/src/clib/calc_grain_size_increment_1d.F @@ -40,7 +40,7 @@ subroutine calc_grain_size_increment_1d( ! in integer in, jn, kn, is, ie, j, k - logical itmask(in) + MASK_TYPE itmask(in) integer immulti, imabund, idspecies, igrgr real*8 dom R_PREC d(in,jn,kn) @@ -369,7 +369,7 @@ subroutine calc_grain_size_increment_1d( endif do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then if (idspecies .gt. 0) then sgtot (i) = sgMgSiO3 (i) @@ -439,7 +439,7 @@ subroutine calc_grain_size_increment_species_1d( ! input integer in, jn, kn, is, ie, j, k - logical itmask(in) + MASK_TYPE itmask(in) integer igrgr integer iSN, nSN, SN0_N real*8 dom @@ -480,7 +480,7 @@ subroutine calc_grain_size_increment_species_1d( ! enddo do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then !!!!!!!!!!!!! !!!! if( dsp(i,j,k) .gt. 1.d-15*d(i,j,k) ) then !!!!!!!!!!!!! diff --git a/src/clib/calc_tdust_1d_g.F b/src/clib/calc_tdust_1d_g.F index b6578456..d0069cdb 100644 --- a/src/clib/calc_tdust_1d_g.F +++ b/src/clib/calc_tdust_1d_g.F @@ -60,7 +60,7 @@ subroutine calc_tdust_1d_g( ! Iteration mask - logical itmask(in) + MASK_TYPE itmask(in) ! Parameters @@ -89,7 +89,7 @@ subroutine calc_tdust_1d_g( & slope(in), tdplus(in), tdustnow(in), tdustold(in), & pert(in), & bi_t_mid(in), bi_t_high(in) - logical nm_itmask(in), bi_itmask(in) + MASK_TYPE nm_itmask(in), bi_itmask(in) !\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\///////////////////////////////// !======================================================================= @@ -105,7 +105,7 @@ subroutine calc_tdust_1d_g( Td_N(1) = gr_N(2) Td_Size = gr_N(2) do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then logalsp(:,i) = log10(alsp(:,i)) endif enddo @@ -119,7 +119,7 @@ subroutine calc_tdust_1d_g( ! Set local iteration mask and initial guess do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then gamma_isrf(i) = isrf(i) * gamma_isrfa(i) endif enddo @@ -127,17 +127,17 @@ subroutine calc_tdust_1d_g( do i = is+1, ie+1 nm_itmask(i) = itmask(i) bi_itmask(i) = itmask(i) - if ( nm_itmask(i) ) then + if ( nm_itmask(i) .ne. MASK_FALSE ) then if (trad .ge. tgas(i)) then tdustnow(i) = trad - nm_itmask(i) = .false. - bi_itmask(i) = .false. + nm_itmask(i) = MASK_FALSE + bi_itmask(i) = MASK_FALSE c_done = c_done + 1 nm_done = nm_done + 1 else if (tgas(i) .gt. t_subl) then ! Use bisection if T_gas > grain sublimation temperature. - nm_itmask(i) = .false. + nm_itmask(i) = MASK_FALSE nm_done = nm_done + 1 else tdustnow(i) = max(trad, @@ -158,7 +158,7 @@ subroutine calc_tdust_1d_g( ! Loop over slice do i = is+1, ie+1 - if ( nm_itmask(i) ) then + if ( nm_itmask(i) .ne. MASK_FALSE ) then tdplus(i) = max(1.e-3_DKIND, ((1._DKIND + pert(i)) & * tdustnow(i))) @@ -185,7 +185,7 @@ subroutine calc_tdust_1d_g( & gamma_isrf, nh, nm_itmask, solplus, in, is, ie) do i = is+1, ie+1 - if ( nm_itmask(i) ) then + if ( nm_itmask(i) .ne. MASK_FALSE ) then ! Use Newton's method to solve for Tdust @@ -203,13 +203,13 @@ subroutine calc_tdust_1d_g( ! If negative solution calculated, give up and wait for bisection step. if (tdustnow(i) .lt. trad) then - nm_itmask(i) = .false. + nm_itmask(i) = MASK_FALSE nm_done = nm_done + 1 ! Check for convergence of solution else if (abs(sol(i) / solplus(i)) .lt. tol) then - nm_itmask(i) = .false. + nm_itmask(i) = MASK_FALSE c_done = c_done + 1 - bi_itmask(i) = .false. + bi_itmask(i) = MASK_FALSE nm_done = nm_done + 1 endif @@ -234,7 +234,7 @@ subroutine calc_tdust_1d_g( ! If iteration count exceeded, try once more with bisection if (c_done .lt. c_total) then do i = is+1, ie+1 - if ( bi_itmask(i) ) then + if ( bi_itmask(i) .ne. MASK_FALSE ) then tdustnow(i) = trad ! bi_t_high(i) = tgas(i) bi_t_high(i) = 3e3_DKIND @@ -244,7 +244,7 @@ subroutine calc_tdust_1d_g( do iter = 1, bi_itmax do i = is+1, ie+1 - if ( bi_itmask(i) ) then + if ( bi_itmask(i) .ne. MASK_FALSE ) then bi_t_mid(i) = 0.5_DKIND * (tdustnow(i) + bi_t_high(i)) if (iter .eq. 1) then @@ -262,7 +262,7 @@ subroutine calc_tdust_1d_g( & gamma_isrf, nh, bi_itmask, sol, in, is, ie) do i = is+1, ie+1 - if ( bi_itmask(i) ) then + if ( bi_itmask(i) .ne. MASK_FALSE ) then if (sol(i) .gt. 0._DKIND) then tdustnow(i) = bi_t_mid(i) @@ -272,7 +272,7 @@ subroutine calc_tdust_1d_g( if ((abs(bi_t_high(i) - tdustnow(i)) / tdustnow(i)) & .le. bi_tol) then - bi_itmask(i) = .false. + bi_itmask(i) = MASK_FALSE c_done = c_done + 1 endif @@ -308,7 +308,7 @@ subroutine calc_tdust_1d_g( ! Copy values back to thrown slice do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then ! Check for bad solutions if (tdustnow(i) .lt. 0._DKIND) then @@ -379,7 +379,7 @@ subroutine calc_kappa_gr_g( ! Iteration mask - logical itmask(in) + MASK_TYPE itmask(in) ! Parameters @@ -407,7 +407,7 @@ subroutine calc_kappa_gr_g( !======================================================================= do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then if(idspecies.eq.0) then @@ -492,7 +492,7 @@ subroutine calc_gr_balance_g( ! Iteration mask - logical itmask(in) + MASK_TYPE itmask(in) ! Parameters @@ -510,7 +510,7 @@ subroutine calc_gr_balance_g( !======================================================================= do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then sol(i) = gamma_isrf(i) + radf * kgr(i) * & (trad4 - tdust(i)**4) + diff --git a/src/clib/calc_tdust_3d_g.F b/src/clib/calc_tdust_3d_g.F index feaa0f15..4637a406 100644 --- a/src/clib/calc_tdust_3d_g.F +++ b/src/clib/calc_tdust_3d_g.F @@ -217,7 +217,7 @@ subroutine calc_tdust_3d_g( & vibh(in) ! Iteration mask for multi_cool - logical itmask(in) + MASK_TYPE itmask(in) !\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\///////////////////////////////// !======================================================================= @@ -377,7 +377,7 @@ subroutine calc_tdust_3d_g( ! Set itmask to all true - itmask(i) = .true. + itmask(i) = MASK_TRUE enddo @@ -386,7 +386,7 @@ subroutine calc_tdust_3d_g( if (imetal .eq. 1) then do i = is+1, ie + 1 if (metal(i,j,k) .lt. 1.e-9_DKIND * d(i,j,k)) then - itmask(i) = .false. + itmask(i) = MASK_FALSE endif enddo endif @@ -430,7 +430,7 @@ subroutine calc_tdust_3d_g( endif do i = is+1, ie+1 - if(itmask(i)) then + if(itmask(i) .ne. MASK_FALSE) then ! Calculate metallicity if (imetal .eq. 1) then @@ -720,7 +720,7 @@ subroutine calc_tdust_3d_g( ! Copy slice values back to grid do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then if (itdmulti .eq. 0) then dust_temp(i,j,k) = tdust(i) else diff --git a/src/clib/calc_temp1d_cloudy_g.F b/src/clib/calc_temp1d_cloudy_g.F index 4aa4dd2a..47865dde 100644 --- a/src/clib/calc_temp1d_cloudy_g.F +++ b/src/clib/calc_temp1d_cloudy_g.F @@ -78,7 +78,7 @@ subroutine calc_temp1d_cloudy_g(d, metal, e, rhoH, ! Iteration mask - logical itmask(in) + MASK_TYPE itmask(in) ! Parameters @@ -147,14 +147,14 @@ subroutine calc_temp1d_cloudy_g(d, metal, e, rhoH, endif do i=is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then ! Calculate proper log(n_H) log_n_h(i) = log10(rhoH(i) * dom) endif enddo do i=is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then munew = 1._DKIND do ti = 1, ti_max muold = munew diff --git a/src/clib/calc_temp_cloudy_g.F b/src/clib/calc_temp_cloudy_g.F index 6f94bbae..5a047ae1 100644 --- a/src/clib/calc_temp_cloudy_g.F +++ b/src/clib/calc_temp_cloudy_g.F @@ -101,7 +101,7 @@ subroutine calc_temp_cloudy_g(d, e, metal, temperature, ! Iteration mask - logical itmask(in) + MASK_TYPE itmask(in) ! !\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\///////////////////////////////// !======================================================================= @@ -144,7 +144,7 @@ subroutine calc_temp_cloudy_g(d, e, metal, temperature, ! Initialize iteration mask to true for all cells. do i = is+1, ie+1 - itmask(i) = .true. + itmask(i) = MASK_TRUE rhoH(i) = fh * d(i,j,k) enddo diff --git a/src/clib/cool1d_cloudy_g.F b/src/clib/cool1d_cloudy_g.F index 6d77b318..8d7676cf 100644 --- a/src/clib/cool1d_cloudy_g.F +++ b/src/clib/cool1d_cloudy_g.F @@ -75,7 +75,7 @@ subroutine cool1d_cloudy_g(d, rhoH, metallicity, ! Iteration mask - logical itmask(in) + MASK_TYPE itmask(in) ! Parameters @@ -115,7 +115,7 @@ subroutine cool1d_cloudy_g(d, rhoH, metallicity, endif do i=is+1, ie+1 - if ( itmask(i) ) then + if (itmask(i) .ne. MASK_FALSE) then log10tem(i) = logtem(i) * inv_log10 diff --git a/src/clib/cool1d_cloudy_old_tables_g.F b/src/clib/cool1d_cloudy_old_tables_g.F index 2e95caea..e816d569 100644 --- a/src/clib/cool1d_cloudy_old_tables_g.F +++ b/src/clib/cool1d_cloudy_old_tables_g.F @@ -80,7 +80,7 @@ subroutine cool1D_cloudy_old_tables_g(d, de, rhoH, metallicity, ! Iteration mask - logical itmask(in) + MASK_TYPE itmask(in) ! Parameters @@ -124,7 +124,7 @@ subroutine cool1D_cloudy_old_tables_g(d, de, rhoH, metallicity, endif do i=is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then log10tem(i) = logtem(i) * inv_log10 diff --git a/src/clib/cool1d_multi_g.F b/src/clib/cool1d_multi_g.F index a4c7fa5a..c59bba07 100644 --- a/src/clib/cool1d_multi_g.F +++ b/src/clib/cool1d_multi_g.F @@ -358,8 +358,8 @@ subroutine cool1d_multi_g( & , gisrfreforg(in), gisrfvolorg(in), gisrfH2Oice(in) ! Iteration mask - logical itmask(in), anydust, interp - logical itmask_metal(in), itmask_tab(in) + MASK_TYPE itmask(in), anydust, interp + MASK_TYPE itmask_metal(in), itmask_tab(in) !!#define CALCULATE_TGAS_SELF_CONSISTENTLY #ifdef CALCULATE_TGAS_SELF_CONSISTENTLY integer iter_tgas @@ -375,13 +375,20 @@ subroutine cool1d_multi_g( ! Set flag for dust-related options - anydust = (idust .gt. 0) .or. (idustall .gt. 0) .or. - & (idustrec .gt. 0) + if ((idust .gt. 0) .or. (idustall .gt. 0) .or. + & (idustrec .gt. 0)) then + anydust = MASK_TRUE + else + anydust = MASK_FALSE + endif ! Set flag for needing interpolation variables - interp = (ispecies .gt. 0) .or. (idustall .gt. 0) - + if ((ispecies .gt. 0) .or. (idustall .gt. 0)) then + interp = MASK_TRUE + else + interp = MASK_FALSE + endif ! Set log values of start and end of lookup tables logtem0 = log(temstart) @@ -414,7 +421,7 @@ subroutine cool1d_multi_g( ! Initialize edot do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then edot(i) = 0._DKIND end if enddo @@ -422,7 +429,7 @@ subroutine cool1d_multi_g( ! Compute Pressure do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then p2d(i) = (gamma - 1._DKIND)*d(i,j,k)*e(i,j,k) end if enddo @@ -438,13 +445,13 @@ subroutine cool1d_multi_g( if (imetal .eq. 1) then do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then rhoH(i) = fh * (d(i,j,k) - metal(i,j,k)) endif enddo else do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then rhoH(i) = fh * d(i,j,k) endif enddo @@ -465,7 +472,7 @@ subroutine cool1d_multi_g( ! Compute mean molecular weight (and temperature) directly do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then mmw(i) = & (HeI(i,j,k) + HeII(i,j,k) + HeIII(i,j,k))/4._DKIND + & HI(i,j,k) + HII(i,j,k) + de(i,j,k) @@ -478,7 +485,7 @@ subroutine cool1d_multi_g( if (ispecies .gt. 1) then do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then mmw(i) = mmw(i) + & HM(i,j,k) + (H2I(i,j,k) + H2II(i,j,k))/2._DKIND rhoH(i) = rhoH(i) + H2I(i,j,k) + H2II(i,j,k) @@ -490,14 +497,14 @@ subroutine cool1d_multi_g( if (imetal .eq. 1) then do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then mmw(i) = mmw(i) + metal(i,j,k)/mu_metal end if enddo endif do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then tgas(i) = max(p2d(i)*utem/mmw(i), temstart) mmw(i) = d(i,j,k) / mmw(i) end if @@ -507,7 +514,7 @@ subroutine cool1d_multi_g( if (ispecies .gt. 1) then do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then nH2 = 0.5_DKIND*(H2I(i,j,k) + H2II(i,j,k)) nother = (HeI(i,j,k) + HeII(i,j,k) + & HeIII(i,j,k))/4._DKIND + @@ -552,19 +559,19 @@ subroutine cool1d_multi_g( if (iTfloor .eq. 1) then do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then if (tgas(i) .le. Tfloor_scalar) then edot(i) = tiny - itmask(i) = .false. + itmask(i) = MASK_FALSE endif endif enddo else if (iTfloor .eq. 2) then do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then if (tgas(i) .le. Tfloor(i,j,k)) then edot(i) = tiny - itmask(i) = .false. + itmask(i) = MASK_FALSE endif endif enddo @@ -574,20 +581,20 @@ subroutine cool1d_multi_g( if (imetal .eq. 1) then do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then metallicity(i) = metal(i,j,k) / d(i,j,k) / z_solar endif enddo else do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then metallicity(i) = tiny endif enddo endif do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then mynh(i) = rhoH(i) * dom end if enddo @@ -596,7 +603,7 @@ subroutine cool1d_multi_g( if (iter .eq. 1) then do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then tgasold(i) = tgas(i) end if enddo @@ -606,7 +613,7 @@ subroutine cool1d_multi_g( logdom = log10(dom) do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then logT(i) = log10(tgas(i)) if(icmbTfloor .eq. 1) & logTcmb(i) = log10(comp2) @@ -643,7 +650,7 @@ subroutine cool1d_multi_g( enddo do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then ! Compute log temperature and truncate if above/below table max/min @@ -656,9 +663,9 @@ subroutine cool1d_multi_g( ! Compute interpolation indices - if (interp) then + if (interp .ne. MASK_FALSE) then do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then ! Compute index into the table and precompute parts of linear interp @@ -677,7 +684,7 @@ subroutine cool1d_multi_g( if (ispecies .gt. 0) then do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then ! Lookup cooling values and do a linear temperature in log(T) @@ -712,7 +719,7 @@ subroutine cool1d_multi_g( ! Compute the cooling function do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then edot(i) = ( ! Collisional excitations @@ -768,7 +775,7 @@ subroutine cool1d_multi_g( if (ih2cr .eq. 3) then do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then lognhat = logH2I(i) - logdvdr(i) @@ -799,7 +806,7 @@ subroutine cool1d_multi_g( else if (ih2cr .eq. 2) then do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then gaHI(i) = gaHIa(indixe(i)) + tdef(i) & *(gaHIa(indixe(i)+1) - gaHIa(indixe(i))) gaH2(i) = gaH2a(indixe(i)) + tdef(i) @@ -820,7 +827,7 @@ subroutine cool1d_multi_g( enddo do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then #ifdef OPTICAL_DEPTH_FUDGE nH2 = 0.5_DKIND*H2I(i,j,k) nother = (HeI(i,j,k) + HeII(i,j,k) + @@ -858,7 +865,7 @@ subroutine cool1d_multi_g( else if (ih2cr .eq. 1) then do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then gpldl(i) = gpldla(indixe(i)) + tdef(i) & *(gpldla(indixe(i)+1) - gpldla(indixe(i))) gphdl(i) = gphdla(indixe(i)) + tdef(i) @@ -869,7 +876,7 @@ subroutine cool1d_multi_g( enddo do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then #define NO_OPTICAL_DEPTH_FUDGE #ifdef OPTICAL_DEPTH_FUDGE @@ -903,7 +910,7 @@ subroutine cool1d_multi_g( else if (ih2cr .eq. 0) then do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then hyd01k(i) = hyd01ka(indixe(i)) + tdef(i) & *(hyd01ka(indixe(i)+1)-hyd01ka(indixe(i))) h2k01(i) = h2k01a(indixe(i)) + tdef(i) @@ -920,7 +927,7 @@ subroutine cool1d_multi_g( enddo do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then qq = 1.2_DKIND*(HI(i,j,k)*dom)**0.77_DKIND + & (H2I(i,j,k)*dom/2._DKIND)**0.77_DKIND vibl = (HI(i,j,k)*hyd01k(i) + @@ -954,7 +961,7 @@ subroutine cool1d_multi_g( C Ripamonti & Abel 2003 if (iciecool.eq.1) then do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then c Only calculate if H2I(i) is a substantial fraction if (d(i,j,k)*dom.gt.1e10_DKIND) then ciefudge = 1._DKIND @@ -974,7 +981,7 @@ subroutine cool1d_multi_g( c CIE H2 cooling using Yoshida et al. (2006) else if (iciecool .eq. 2) then do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then cieY06(i) = cieY06a(indixe(i)) + tdef(i) & *(cieY06a(indixe(i)+1) - cieY06a(indixe(i))) LCIE(i) = - cieY06(i) * (H2I(i,j,k)/2.d0)**2 @@ -993,7 +1000,7 @@ subroutine cool1d_multi_g( if (ihdcr .eq. 1 ) then do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then lognhat = logHDI(i) - logdvdr(i) @@ -1024,7 +1031,7 @@ subroutine cool1d_multi_g( else if (ihdcr .eq. 0) then do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then c CMB cooling floor if (tgas(i) .gt. comp2) then hdlte(i) = hdltea(indixe(i)) + tdef(i) @@ -1039,7 +1046,7 @@ subroutine cool1d_multi_g( enddo do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then c old (incorrect) way: c hdlte1 = hdlte(i)/(HDI(i,j,k)*dom/2._DKIND) c hdlow1 = max(hdlow(i), tiny) @@ -1063,12 +1070,12 @@ subroutine cool1d_multi_g( if (metallicity(i) .ge. min_metallicity) then itmask_metal(i) = itmask(i) else - itmask_metal(i) = .false. + itmask_metal(i) = MASK_FALSE endif enddo else do i = is+1, ie+1 - itmask_metal(i) = .false. + itmask_metal(i) = MASK_FALSE enddo endif @@ -1112,16 +1119,16 @@ subroutine cool1d_multi_g( ! Calculate dust to gas ratio - if (anydust .or. (igammah .gt. 0)) then + if ((anydust .ne. MASK_FALSE) .or. (igammah .gt. 0)) then if (idustfield .gt. 0) then do i = is+1, ie+1 - if ( itmask_metal(i) ) then + if ( itmask_metal(i) .ne. MASK_FALSE ) then dust2gas(i) = dust(i,j,k) / d(i,j,k) endif enddo else do i = is+1, ie+1 - if ( itmask_metal(i) ) then + if ( itmask_metal(i) .ne. MASK_FALSE ) then dust2gas(i) = fgr * metallicity(i) endif enddo @@ -1130,25 +1137,25 @@ subroutine cool1d_multi_g( ! Calculate interstellar radiation field - if (anydust .or. (igammah .gt. 1)) then + if ((anydust .ne. MASK_FALSE) .or. (igammah .gt. 1)) then if (iisrffield .gt. 0) then do i = is+1, ie+1 - if ( itmask_metal(i) ) then + if ( itmask_metal(i) .ne. MASK_FALSE ) then myisrf(i) = isrf_habing(i,j,k) endif enddo else do i = is+1, ie+1 - if ( itmask_metal(i) ) then + if ( itmask_metal(i) .ne. MASK_FALSE ) then myisrf(i) = isrf endif enddo endif endif - if (anydust .or. (igammah .gt. 1)) then + if ((anydust .ne. MASK_FALSE) .or. (igammah .gt. 1)) then do i = is+1, ie+1 - if ( itmask_metal(i) ) then + if ( itmask_metal(i) .ne. MASK_FALSE ) then if (idspecies .eq. 0 ) then if (idustfield .gt. 0) then @@ -1201,12 +1208,12 @@ subroutine cool1d_multi_g( ! --- Gas to grain heat transfer --- - if (anydust) then + if (anydust .ne. MASK_FALSE) then ! Look up gas/grain heat transfer rates do i = is+1, ie+1 - if ( itmask_metal(i) ) then + if ( itmask_metal(i) .ne. MASK_FALSE ) then if(idspecies .eq. 0) then @@ -1379,7 +1386,7 @@ subroutine cool1d_multi_g( ! Calculate dust cooling rate do i = is+1, ie+1 - if ( itmask_metal(i) ) then + if ( itmask_metal(i) .ne. MASK_FALSE ) then if (idspecies .eq. 0) then @@ -1438,7 +1445,7 @@ subroutine cool1d_multi_g( if ( ipcont .eq. 1 ) then do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then !! primordial continuum opacity !! call interpolate_2D_g( @@ -1451,7 +1458,7 @@ subroutine cool1d_multi_g( else do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then alpha(i) = 0. endif enddo @@ -1461,9 +1468,9 @@ subroutine cool1d_multi_g( ! Add dust opacity. ! if (idspecies .eq. 0), dust opacity is overestimated at Td > 50 K ! We better not include dust opacity. - if ((anydust).and.(idspecies .gt. 0)) then + if ((anydust .ne. MASK_FALSE).and.(idspecies .gt. 0)) then do i = is+1, ie+1 - if ( itmask_metal(i) ) then + if ( itmask_metal(i) .ne. MASK_FALSE ) then if (itdmulti .eq. 0) then @@ -1502,7 +1509,7 @@ subroutine cool1d_multi_g( endif !! anydust do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then tau_con(i) = alpha(i) * lshield_con(i) endif enddo @@ -1514,7 +1521,7 @@ subroutine cool1d_multi_g( if (iradshield == 0) then ! no shielding do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then edot(i) = edot(i) + real(ipiht, DKIND)*( & piHI *HI (i,j,k) ! pi of HI & + piHeI *HeI (i,j,k)*0.25_DKIND ! pi of HeI @@ -1531,7 +1538,7 @@ subroutine cool1d_multi_g( ! do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then if (k24 .lt. tiny8) then fSShHI = 1._DKIND else @@ -1563,7 +1570,7 @@ subroutine cool1d_multi_g( ! do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then ! ! HI self shielding ratio ! @@ -1616,7 +1623,7 @@ subroutine cool1d_multi_g( ! do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then ! ! HI self shielding ratio ! @@ -1685,7 +1692,7 @@ subroutine cool1d_multi_g( ! Calculate electron density from mean molecular weight do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then myde(i) = 1 - mmw(i) * (3.0_DKIND * fh + 1.0_DKIND) / & 4.0_DKIND @@ -1706,7 +1713,7 @@ subroutine cool1d_multi_g( if (igammah .eq. 1) then do i = is + 1, ie + 1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then if ( tgas(i) > 2.d4 ) then gammaha_eff(i) = 0._DKIND else @@ -1719,7 +1726,7 @@ subroutine cool1d_multi_g( else if (igammah .eq. 2) then do i = is + 1, ie + 1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then if ( tgas(i) > 2.d4 ) then gammaha_eff(i) = 0._DKIND else @@ -1733,7 +1740,7 @@ subroutine cool1d_multi_g( else if (igammah .eq. 3) then do i = is + 1, ie + 1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then pe_X = myisrf(i) * dom_inv * sqrt(tgas(i)) / myde(i) pe_eps = & (4.9d-2 / @@ -1748,7 +1755,7 @@ subroutine cool1d_multi_g( if (igammah .gt. 0) then do i = is + 1, ie + 1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then edot(i) = edot(i) + gammaha_eff(i) * rhoH(i) * & dom_inv * dust2gas(i) / fgr endif @@ -1760,14 +1767,14 @@ subroutine cool1d_multi_g( if ((idustall .gt. 0) .or. (idustrec .gt. 0)) then do i = is + 1, ie + 1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then regr(i) = regra(indixe(i)) + tdef(i) & *(regra(indixe(i)+1) -regra(indixe(i))) endif enddo do i = is + 1, ie + 1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then grbeta = 0.74_DKIND / tgas(i)**0.068_DKIND edot(i) = edot(i) - & regr(i) * (myisrf(i)*dom_inv / myde(i))**grbeta * @@ -1780,7 +1787,7 @@ subroutine cool1d_multi_g( ! Compton cooling or heating and X-ray compton heating do i = is + 1, ie + 1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then edot(i) = edot(i) @@ -1799,7 +1806,7 @@ subroutine cool1d_multi_g( if (iradtrans .eq. 1) then do i = is + 1, ie + 1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then edot(i) = edot(i) + real(ipiht, DKIND) * photogamma(i,j,k) & / coolunit * HI(i,j,k) / dom @@ -1826,10 +1833,10 @@ subroutine cool1d_multi_g( ! Determine if the temperature is above the threshold to do tabulated cooling. do i = is+1, ie+1 itmask_tab(i) = itmask_metal(i) - if ( itmask_tab(i) ) then + if ( itmask_tab(i) .ne. MASK_FALSE ) then if (( tmcool .gt. 0.0d0 ) .and. & ( tgas(i) .lt. tmcool )) then - itmask_tab(i) = .false. + itmask_tab(i) = MASK_FALSE endif endif enddo @@ -1866,7 +1873,7 @@ subroutine cool1d_multi_g( ! C/O fine-structure cooling do i = is+1, ie+1 - if ( itmask_metal(i) ) then + if ( itmask_metal(i) .ne. MASK_FALSE ) then ! CI lognhat = logCI(i) - logdvdr(i) @@ -2025,7 +2032,7 @@ subroutine cool1d_multi_g( if (iVheat .eq. 1) then do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then edot(i) = edot(i) + Vheat(i,j,k) / coolunit / dom**2 end if enddo @@ -2035,7 +2042,7 @@ subroutine cool1d_multi_g( if (iMheat .eq. 1) then do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then edot(i) = edot(i) + Mheat(i,j,k) * d(i,j,k) * mh & / coolunit / dom end if @@ -2046,7 +2053,7 @@ subroutine cool1d_multi_g( ! Continuum opacity do i = is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then if ( tau_con(i) .gt. 1.d0 ) then if ( tau_con(i) .lt. 1.d2 ) then edot(i) = edot(i) * tau_con(i)**(-2.d0) @@ -2060,7 +2067,7 @@ subroutine cool1d_multi_g( ! Set tgasold do i=is+1, ie+1 - if ( itmask(i) ) then + if ( itmask(i) .ne. MASK_FALSE ) then tgasold(i) = tgas(i) end if enddo diff --git a/src/clib/cool_multi_time_g.F b/src/clib/cool_multi_time_g.F index a0338579..5e3e0a84 100644 --- a/src/clib/cool_multi_time_g.F +++ b/src/clib/cool_multi_time_g.F @@ -286,7 +286,7 @@ subroutine cool_multi_time_g( ! Iteration mask for multi_cool - logical itmask(in), itmask_metal(in) + MASK_TYPE itmask(in), itmask_metal(in) !\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\///////////////////////////////// !======================================================================= @@ -347,7 +347,7 @@ subroutine cool_multi_time_g( j = mod(t,dj) + js+1 do i = is+1, ie+1 - itmask(i) = .true. + itmask(i) = MASK_TRUE end do ! Compute the cooling rate diff --git a/src/clib/lookup_cool_rates0d.F b/src/clib/lookup_cool_rates0d.F index 4403da8e..11e1f7f9 100644 --- a/src/clib/lookup_cool_rates0d.F +++ b/src/clib/lookup_cool_rates0d.F @@ -141,7 +141,7 @@ subroutine lookup_cool_rates0d(output, dtit, & iH2shield, iradshield, & iradtrans, irt_honly, iisrffield & ,imchem, igrgr, ipcont - logical itmask, itmask_metal, anydust + MASK_TYPE itmask, itmask_metal, anydust real*8 dx_cgs, aye, temstart, temend, gamma, & utim, uxyz, uaye, urho, utem, fh, z_solar, @@ -432,7 +432,7 @@ subroutine lookup_cool_rates0d(output, dtit, HDII = dsp(14) HeHII = dsp(15) endif - if ( itmask_metal ) then + if ( itmask_metal .ne. MASK_FALSE ) then if (imchem .eq. 1) then CI = dsp(16) CII = dsp(17) @@ -909,7 +909,8 @@ subroutine lookup_cool_rates0d(output, dtit, scoef = scoef & + kdissHDI * HDI /3.0_DKIND endif - if ((imchem .eq. 1) .and. (itmask_metal )) then + if ((imchem .eq. 1) .and. + & (itmask_metal .ne. MASK_FALSE)) then if (idissZ .gt. 0) then scoef = scoef & + kdissOH * OH /17.0_DKIND @@ -918,8 +919,8 @@ subroutine lookup_cool_rates0d(output, dtit, endif endif - if (anydust) then - if(itmask_metal ) then + if (anydust .ne. MASK_FALSE) then + if(itmask_metal .ne. MASK_FALSE ) then acoef = acoef + 2._DKIND * h2dust * rhoH endif endif @@ -949,7 +950,8 @@ subroutine lookup_cool_rates0d(output, dtit, & + k152 * HeHII / 5._DKIND endif - if ( (imchem .eq. 1) .and. (itmask_metal ) ) then + if ((imchem .eq. 1) .and. + & (itmask_metal .ne. MASK_FALSE)) then scoef = scoef & + kz20 * CI * H2I / 24._DKIND & + kz21 * OI * H2I / 32._DKIND @@ -1033,7 +1035,8 @@ subroutine lookup_cool_rates0d(output, dtit, & + k149 * HeI / 4._DKIND endif - if ( (imchem .eq. 1) .and. (itmask_metal ) ) then + if ((imchem .eq. 1) .and. + & (itmask_metal .ne. MASK_FALSE)) then scoef = scoef & + kz39 * OII * HI / 16._DKIND & + kz43 * COII * HI / 28._DKIND @@ -1065,7 +1068,8 @@ subroutine lookup_cool_rates0d(output, dtit, if ( (iradtrans .eq. 1) .and. (irt_honly .eq. 1) ) & scoef = scoef + kphHI * HI if (iradtrans .eq. 1) then - if ((imchem .eq. 1) .and. (itmask_metal )) then + if ((imchem .eq. 1) .and. + & (itmask_metal .ne. MASK_FALSE)) then if (iionZ .gt. 0) then scoef = scoef & + kphCI * CI /12.0_DKIND @@ -1100,7 +1104,8 @@ subroutine lookup_cool_rates0d(output, dtit, & + k153 * HeHII / 5._DKIND endif - if ( (imchem .eq. 1) .and. (itmask_metal ) ) then + if ((imchem .eq. 1) .and. + & (itmask_metal .ne. MASK_FALSE)) then scoef = scoef acoef = acoef & + kz44 * CII / 12._DKIND @@ -1124,8 +1129,8 @@ subroutine lookup_cool_rates0d(output, dtit, & + k12 *de ) & + k29shield + k31shield - if (anydust) then - if(itmask_metal ) then + if (anydust .ne. MASK_FALSE) then + if(itmask_metal .ne. MASK_FALSE ) then scoef = scoef + 2._DKIND * h2dust * & HI * rhoH endif @@ -1141,7 +1146,8 @@ subroutine lookup_cool_rates0d(output, dtit, & + k54 * DI / 2._DKIND endif - if ( (imchem .eq. 1) .and. (itmask_metal ) ) then + if ((imchem .eq. 1) .and. + & (itmask_metal .ne. MASK_FALSE)) then scoef = scoef + 2._DKIND * ( 0._DKIND & + kz15 * HI * CH / 13._DKIND & + kz16 * HI * CH2 / 14._DKIND @@ -1380,7 +1386,7 @@ subroutine lookup_cool_rates0d(output, dtit, ! if (imchem .eq. 1) then - if (itmask_metal ) then + if (itmask_metal .ne. MASK_FALSE ) then C***** CI ********** scoef = 0._DKIND + 12._DKIND * ( 0._DKIND @@ -1820,7 +1826,7 @@ subroutine lookup_cool_rates0d(output, dtit, ! if ( ( igrgr .eq. 1 ) .or. ( idsub .eq. 1) ) then - if (itmask_metal ) then + if (itmask_metal .ne. MASK_FALSE ) then if (idspecies .gt. 0) then C***** MgSiO3 ********** diff --git a/src/clib/solve_rate_cool_g.F b/src/clib/solve_rate_cool_g.F index 0e859ee4..478d1df7 100644 --- a/src/clib/solve_rate_cool_g.F +++ b/src/clib/solve_rate_cool_g.F @@ -541,9 +541,9 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de, ! Iteration mask - logical itmask(in), anydust - logical itmask_tmp(in), itmask_nr(in) - logical itmask_metal(in) + MASK_TYPE itmask(in), anydust + MASK_TYPE itmask_tmp(in), itmask_nr(in) + MASK_TYPE itmask_metal(in) integer itr, imp_eng(in), itr_time integer nsp, isp, jsp, id real*8 dspj, err, err_max @@ -565,7 +565,11 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de, ! Set flag for dust-related options - anydust = (idust .gt. 0) .or. (idustall .gt. 0) + if ((idust .gt. 0) .or. (idustall .gt. 0)) then + anydust = MASK_TRUE + else + anydust = MASK_FALSE + endif ! ignore metal chemistry/cooling below this metallicity min_metallicity = 1.d-9 / z_solar @@ -725,7 +729,7 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de, ! Initialize iteration mask to true for all cells. do i = is+1, ie+1 - itmask(i) = .true. + itmask(i) = MASK_TRUE enddo ! If we are using coupled radiation with intermediate stepping, @@ -735,9 +739,9 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de, if (iradcoupled .eq. 1 .and. iradstep .eq. 1) then do i = is+1, ie+1 if (kphHI(i,j,k) .gt. 0) then - itmask(i) = .true. + itmask(i) = MASK_TRUE else - itmask(i) = .false. + itmask(i) = MASK_FALSE endif enddo endif @@ -746,9 +750,9 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de, if (iradcoupled .eq. 1 .and. iradstep .eq. 0) then do i = is+1, ie + 1 if (kphHI(i,j,k) .gt. 0) then - itmask(i) = .false. + itmask(i) = MASK_FALSE else - itmask(i) = .true. + itmask(i) = MASK_TRUE endif enddo endif @@ -771,7 +775,7 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de, do iter = 1, itmax do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then dtit(i) = huge8 endif enddo @@ -1008,7 +1012,7 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de, itmask_tmp = itmask itmask_nr = itmask do i = is+1, ie+1 - if ( itmask_tmp(i) ) then + if ( itmask_tmp(i) .ne. MASK_FALSE ) then if ( ( (imetal .eq. 0) & .and. (ddom(i) .lt. 1.d8) ) @@ -1017,16 +1021,16 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de, & .and. (ddom(i) .lt. 1.d8) ) & .or. ( (metallicity(i) .gt. min_metallicity) & .and. (ddom(i) .lt. 1.d6) ) ) ) ) then - itmask_nr(i) = .false. + itmask_nr(i) = MASK_FALSE else - itmask(i) = .false. + itmask(i) = MASK_FALSE endif endif enddo do i = is+1, ie+1 - if (itmask_nr(i)) then + if (itmask_nr(i) .ne. MASK_FALSE) then if ( (icool .eq. 1) .and. (ispecies .gt. 1) .and. & ((ddom(i) .gt. 1.d7) & .and.(tgas(i) .gt. 1650.d0)) ) then @@ -1040,7 +1044,7 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de, ! Find timestep that keeps relative chemical changes below 10% do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then ! Bound from below to prevent numerical errors if (abs(dedot(i)) .lt. tiny8) @@ -1122,11 +1126,13 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de, dtit(i) = min(olddtit*1.5_DKIND, dtit(i)) endif - else if ((itmask_nr(i)).and.(imp_eng(i).eq.0)) then + else if ((itmask_nr(i).ne.MASK_FALSE).and. + & (imp_eng(i).eq.0)) then dtit(i) = min(abs(0.1_DKIND*e(i,j,k)/edot(i)*d(i,j,k)), & dt-ttot(i), 0.5_DKIND*dt) - else if ((itmask_nr(i)).and.(imp_eng(i).eq.1)) then + else if ((itmask_nr(i).ne.MASK_FALSE).and. + & (imp_eng(i).eq.1)) then dtit(i) = dt - ttot(i) else ! itmask @@ -1139,7 +1145,7 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de, ! Compute maximum timestep for cooling/heating do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then ! Set energy per unit volume of this cell based in the pressure ! (the gamma used here is the right one even for H2 since p2d ! is calculated with this gamma). @@ -1179,7 +1185,7 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de, if (icool .eq. 1) then do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then e(i,j,k) = e(i,j,k) + & real(edot(i)/d(i,j,k)*dtit(i), RKIND) @@ -1253,7 +1259,7 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de, ! should be put into its own function. ! Start Newton-Raphson scheme do i = is+1, ie+1 - if (itmask_nr(i)) then + if (itmask_nr(i) .ne. MASK_FALSE) then ! If density and temperature are low, update gas energy explicitly @@ -1269,7 +1275,7 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de, if (ispecies .gt. 1) nsp = nsp + 3 if (ispecies .gt. 2) nsp = nsp + 3 if (ispecies .gt. 3) nsp = nsp + 3 - if (itmask_metal(i)) then + if (itmask_metal(i) .ne. MASK_FALSE) then if (imchem .eq. 1) then nsp = nsp + 19 if ( ( igrgr .eq. 1 ) .or. ( idsub .eq. 1) ) then @@ -1311,7 +1317,7 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de, dsp(14) = HDII(i,j,k) dsp(15) = HeHII(i,j,k) endif - if ( itmask_metal(i) ) then + if ( itmask_metal(i) .ne. MASK_FALSE ) then if ( imchem .eq. 1 ) then dsp(16) = CI(i,j,k) dsp(17) = CII(i,j,k) @@ -1392,7 +1398,7 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de, idsp(id) = isp enddo endif - if (itmask_metal(i)) then + if (itmask_metal(i) .ne. MASK_FALSE) then if (imchem .eq. 1) then do isp = 16, 34 id = id + 1 @@ -1850,7 +1856,7 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de, HDII(i,j,k) = dsp(14) HeHII(i,j,k) = dsp(15) endif - if ( itmask_metal(i) ) then + if ( itmask_metal(i) .ne. MASK_FALSE ) then if ( imchem .eq. 1 ) then CI(i,j,k) = dsp(16) CII(i,j,k) = dsp(17) @@ -1928,7 +1934,7 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de, do i = is+1, ie+1 ttot(i) = min(ttot(i) + dtit(i), dt) if (abs(dt-ttot(i)) .lt. - & tolerance*dt) itmask(i) = .false. + & tolerance*dt) itmask(i) = MASK_FALSE if (ttot(i).lt.ttmin) ttmin = ttot(i) enddo @@ -1955,7 +1961,7 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de, write(0,'((16(1pe8.1)))') (dtit(i),i=is+1,ie+1) write(0,'((16(1pe8.1)))') (ttot(i),i=is+1,ie+1) write(0,'((16(1pe8.1)))') (edot(i),i=is+1,ie+1) - write(0,'((16(l3)))') (itmask(i),i=is+1,ie+1) + write(0,'((16(I3)))') (itmask(i),i=is+1,ie+1) if (exititmax .eq. 1) then ierr = 0 @@ -2393,7 +2399,7 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, real*8 temstart, temend, tgas1d(in), mmw(in), dom, & dtemstart, dtemend real*8 coolunit, tbase1, xbase1, dx_cgs, c_ljeans - logical itmask(in), anydust + MASK_TYPE itmask(in), anydust ! Chemistry rates as a function of temperature @@ -2600,7 +2606,7 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, real*8 tgas_touse, ngas_touse, aWG2019 real*8 nSSh, nratio - logical itmask_metal(in) + MASK_TYPE itmask_metal(in) ! debug integer item @@ -2612,7 +2618,7 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, dlogtem = (log(temend) - log(temstart))/real(nratec-1, DKIND) do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then ! Compute temp-centered temperature (and log) ! logtem(i) = log(0.5_DKIND*(tgas(i)+tgasold(i))) @@ -2653,7 +2659,7 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, if (ispecies .gt. 1) then do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then k7(i) = k7a(indixe(i)) + & (k7a(indixe(i)+1) -k7a(indixe(i)))*tdef(i) k8(i) = k8a(indixe(i)) + @@ -2697,7 +2703,7 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, do n1 = 1, 14 do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then k13dd(i,n1) = k13dda(indixe(i),n1) + & (k13dda(indixe(i)+1,n1) - & k13dda(indixe(i) ,n1) )*tdef(i) @@ -2711,7 +2717,7 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, if (ispecies .gt. 2) then do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then k50(i) = k50a(indixe(i)) + & (k50a(indixe(i)+1) -k50a(indixe(i)))*tdef(i) k51(i) = k51a(indixe(i)) + @@ -2734,7 +2740,7 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, if (ispecies .gt. 3) then do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then k125(i) = k125a(indixe(i)) + & (k125a(indixe(i)+1) -k125a(indixe(i)))*tdef(i) k129(i) = k129a(indixe(i)) + @@ -2775,7 +2781,7 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, if (imchem .eq. 1) then do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then kz15(i) = kz15a(indixe(i)) + & (kz15a(indixe(i)+1) -kz15a(indixe(i)))*tdef(i) kz16(i) = kz16a(indixe(i)) + @@ -2862,7 +2868,7 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, ! Compute grain size increment - if ( (anydust) .and. (idspecies .gt. 0) ) then + if ( (anydust .ne. MASK_FALSE) .and. (idspecies .gt. 0) ) then call calc_grain_size_increment_1d( & immulti, imabund, idspecies, igrgr, itmask_metal @@ -2900,7 +2906,7 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, ! Look-up for H2 formation on dust - if (anydust) then + if (anydust .ne. MASK_FALSE) then d_logtem0 = log(dtemstart) d_logtem9 = log(dtemend) @@ -2910,7 +2916,7 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, if(idspecies .eq. 0) then do i = is+1, ie+1 - if (itmask_metal(i)) then + if (itmask_metal(i) .ne. MASK_FALSE) then ! Assume dust melting at T > 1500 K @@ -2970,7 +2976,7 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, enddo do i = is+1, ie+1 - if (itmask_metal(i)) then + if (itmask_metal(i) .ne. MASK_FALSE) then if (itdmulti .eq. 0) then @@ -3108,7 +3114,7 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, ! Compute grain growth rate do i = is+1, ie+1 - if (itmask_metal(i)) then + if (itmask_metal(i) .ne. MASK_FALSE) then if (igrgr .eq. 1) then @@ -3315,7 +3321,7 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, ! Include approximate self-shielding factors if requested do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then k24shield(i) = k24 k25shield(i) = k25 k26shield(i) = k26 @@ -3334,13 +3340,13 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, ! your hydro code, add kdissH2I later if (iradtrans == 0 .or. iuseH2shield == 1) then do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then k31shield(i) = k31 endif enddo else do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then k31shield(i) = k31 + kdissH2I(i,j,k) endif enddo @@ -3349,7 +3355,7 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, if (iH2shield .gt. 0) then do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then ! Calculate a Sobolev-like length assuming a 3D grid. if (iH2shield == 1) then @@ -3420,7 +3426,7 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, if (iradtrans == 1 .and. iuseH2shield == 1) then C write(*,*) 'kdissH2I included' do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then k31shield(i) = k31shield(i) + kdissH2I(i,j,k) endif enddo @@ -3429,7 +3435,7 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, ! Custom H2 shielding if (iH2shieldcustom .gt. 0) then do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then k31shield(i) = f_shield_custom(i,j,k) * k31shield(i) endif enddo @@ -3439,7 +3445,7 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, if (iradshield > 0) then ! Compute shielding factors do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then ! Compute shielding factor for H nSSh = 6.73e-3_DKIND * @@ -3494,7 +3500,7 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, ! using same scaling. (rate k29) ! do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then if (k24 .lt. tiny8) then k24shield(i) = 0._DKIND @@ -3528,7 +3534,7 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, ! do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then if (k24 .lt. tiny8) then k24shield(i) = 0._DKIND @@ -3577,7 +3583,7 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, ! in HI and HeI, but ignoring HeII heating entirely ! do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then if (k24 .lt. tiny8) then k24shield(i) = 0._DKIND @@ -3628,7 +3634,7 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, ! (see calc_rate.src) do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then x = max(HII(i,j,k)/(HI(i,j,k)+HII(i,j,k)), 1.0e-4_DKIND) factor = 0.3908_DKIND*(1._DKIND - & x**0.4092_DKIND)**1.7592_DKIND @@ -3655,7 +3661,7 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k, #ifdef USE_DENSITY_DEPENDENT_H2_DISSOCIATION_RATE if (ispecies .gt. 1 .and. ithreebody .eq. 0) then do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then nh = min(HI(i,j,k)*dom, 1.0e9_DKIND) k13(i) = tiny8 if (tgas1d(i) .ge. 500._DKIND .and. @@ -3732,7 +3738,7 @@ subroutine rate_timestep_g( & iradtrans, irt_honly real*8 dedot(in), HIdot(in), dom real*8 edot(in) - logical itmask(in), anydust + MASK_TYPE itmask(in), anydust ! Density fields @@ -3781,7 +3787,7 @@ subroutine rate_timestep_g( if (ispecies .eq. 1) then do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then ! Compute the electron density rate-of-change dedot(i) = @@ -3813,7 +3819,7 @@ subroutine rate_timestep_g( ! Include molecular hydrogen rates for HIdot do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then HIdot(i) = & - k1(i) *de(i,j,k) *HI(i,j,k) & - k7(i) *de(i,j,k) *HI(i,j,k) @@ -3837,7 +3843,7 @@ subroutine rate_timestep_g( ! Add H2 formation on dust grains - if (anydust) then + if (anydust .ne. MASK_FALSE) then if (metal(i,j,k) .gt. 1.e-9_DKIND * d(i,j,k)) then HIdot(i) = HIdot(i) & - 2._DKIND * h2dust(i) * rhoH(i) * HI(i,j,k) @@ -3899,7 +3905,7 @@ subroutine rate_timestep_g( !! endif - if (anydust) then + if (anydust .ne. MASK_FALSE) then if (metal(i,j,k) .gt. 1.e-9_DKIND * d(i,j,k)) then H2delta(i) = H2delta(i) + & h2dust(i) * HI(i,j,k) * rhoH(i) * @@ -3966,7 +3972,7 @@ subroutine rate_timestep_g( if (iradtrans .eq. 1) then if (irt_honly .eq. 0) then do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then HIdot(i) = HIdot(i) - kphHI(i,j,k)*HI(i,j,k) dedot(i) = dedot(i) + kphHI(i,j,k)*HI(i,j,k) & + kphHeI(i,j,k) * HeI(i,j,k) / 4._DKIND @@ -3975,7 +3981,7 @@ subroutine rate_timestep_g( enddo else do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then HIdot(i) = HIdot(i) - kphHI(i,j,k)*HI(i,j,k) dedot(i) = dedot(i) + kphHI(i,j,k)*HI(i,j,k) endif @@ -3983,7 +3989,7 @@ subroutine rate_timestep_g( endif if ((ispecies .gt. 2).and.(idissHDI .gt. 0)) then do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then HIdot(i) = HIdot(i) & + kdissHDI(i,j,k) * HDI(i,j,k)/3.0_DKIND endif @@ -3991,7 +3997,7 @@ subroutine rate_timestep_g( endif if ((imchem .gt. 0).and.(iionZ .gt. 0)) then do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then dedot(i) = dedot(i) & + kphCI(i,j,k) * CI(i,j,k)/12.0_DKIND & + kphOI(i,j,k) * OI(i,j,k)/16.0_DKIND @@ -4000,7 +4006,7 @@ subroutine rate_timestep_g( endif if ((imchem .gt. 0).and.(idissZ .gt. 0)) then do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then HIdot(i) = HIdot(i) & + kdissOH (i,j,k) * OH(i,j,k) /17.0_DKIND & + kdissH2O(i,j,k) * H2O(i,j,k)/18.0_DKIND @@ -4086,7 +4092,7 @@ subroutine step_rate_g(de, HI, HII, HeI, HeII, HeIII, d, integer ispecies, in, jn, kn, is, ie, j, k, & iradtrans, irt_honly real*8 dtit(in), dedot_prev(in), HIdot_prev(in) - logical itmask(in), itmask_metal(in), anydust + MASK_TYPE itmask(in), itmask_metal(in), anydust ! Density fields @@ -4181,7 +4187,7 @@ subroutine step_rate_g(de, HI, HII, HeI, HeII, HeIII, d, if (ispecies .eq. 1) then do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then ! 1) HI @@ -4266,7 +4272,7 @@ subroutine step_rate_g(de, HI, HII, HeI, HeII, HeIII, d, ! --- (B) Do helium chemistry in any case: (for all ispecies values) --- do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then ! 4) HeI @@ -4330,7 +4336,7 @@ subroutine step_rate_g(de, HI, HII, HeI, HeII, HeIII, d, ! First, do HI/HII with molecular hydrogen terms do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then ! 1) HI ! @@ -4361,7 +4367,8 @@ subroutine step_rate_g(de, HI, HII, HeI, HeII, HeIII, d, scoef = scoef & + kdissHDI(i,j,k) * HDI(i,j,k)/3.0_DKIND endif - if ((imchem .eq. 1) .and. (itmask_metal(i))) then + if ( (imchem .eq. 1) .and. + & (itmask_metal(i) .ne. MASK_FALSE) ) then if (idissZ .gt. 0) then scoef = scoef & + kdissOH (i,j,k) * OH(i,j,k) /17.0_DKIND @@ -4370,8 +4377,8 @@ subroutine step_rate_g(de, HI, HII, HeI, HeII, HeIII, d, endif endif - if (anydust) then - if(itmask_metal(i)) then + if (anydust .ne. MASK_FALSE) then + if(itmask_metal(i) .ne. MASK_FALSE) then acoef = acoef + 2._DKIND * h2dust(i) * rhoH(i) endif endif @@ -4401,7 +4408,8 @@ subroutine step_rate_g(de, HI, HII, HeI, HeII, HeIII, d, & + k152(i) * HeHII(i,j,k) / 5._DKIND endif - if ( (imchem .eq. 1) .and. (itmask_metal(i)) ) then + if ( (imchem .eq. 1) .and. + & (itmask_metal(i) .ne. MASK_FALSE) ) then scoef = scoef & + kz20(i) * CI(i,j,k) * H2I(i,j,k) / 24._DKIND & + kz21(i) * OI(i,j,k) * H2I(i,j,k) / 32._DKIND @@ -4485,7 +4493,8 @@ subroutine step_rate_g(de, HI, HII, HeI, HeII, HeIII, d, & + k149(i) * HeI(i,j,k) / 4._DKIND endif - if ( (imchem .eq. 1) .and. (itmask_metal(i)) ) then + if ( (imchem .eq. 1) .and. + & (itmask_metal(i) .ne. MASK_FALSE) ) then scoef = scoef & + kz39(i) * OII(i,j,k) * HI(i,j,k) / 16._DKIND & + kz43(i) * COII(i,j,k) * HI(i,j,k) / 28._DKIND @@ -4517,7 +4526,8 @@ subroutine step_rate_g(de, HI, HII, HeI, HeII, HeIII, d, if ( (iradtrans .eq. 1) .and. (irt_honly .eq. 1) ) & scoef = scoef + kphHI(i,j,k) * HIp(i) if (iradtrans .eq. 1) then - if ((imchem .eq. 1) .and. (itmask_metal(i))) then + if ( (imchem .eq. 1) .and. + & (itmask_metal(i) .ne. MASK_FALSE) ) then if (iionZ .gt. 0) then scoef = scoef & + kphCI(i,j,k) * CI(i,j,k)/12.0_DKIND @@ -4552,7 +4562,8 @@ subroutine step_rate_g(de, HI, HII, HeI, HeII, HeIII, d, & + k153(i) * HeHII(i,j,k) / 5._DKIND endif - if ( (imchem .eq. 1) .and. (itmask_metal(i)) ) then + if ( (imchem .eq. 1) .and. + & (itmask_metal(i) .ne. MASK_FALSE) ) then scoef = scoef acoef = acoef & + kz44(i) * CII(i,j,k) / 12._DKIND @@ -4576,8 +4587,8 @@ subroutine step_rate_g(de, HI, HII, HeI, HeII, HeIII, d, & + k12(i)*de(i,j,k) ) & + k29shield(i) + k31shield(i) - if (anydust) then - if(itmask_metal(i)) then + if (anydust .ne. MASK_FALSE) then + if(itmask_metal(i) .ne. MASK_FALSE) then scoef = scoef + 2._DKIND * h2dust(i) * & HI(i,j,k) * rhoH(i) endif @@ -4593,7 +4604,8 @@ subroutine step_rate_g(de, HI, HII, HeI, HeII, HeIII, d, & + k54(i) * DI (i,j,k) / 2._DKIND endif #endif - if ( (imchem .eq. 1) .and. (itmask_metal(i)) ) then + if ( (imchem .eq. 1) .and. + & (itmask_metal(i) .ne. MASK_FALSE) ) then scoef = scoef + 2._DKIND * ( 0._DKIND & + kz15(i) * HI(i,j,k) * CH(i,j,k) / 13._DKIND & + kz16(i) * HI(i,j,k) * CH2(i,j,k) / 14._DKIND @@ -4689,7 +4701,7 @@ subroutine step_rate_g(de, HI, HII, HeI, HeII, HeIII, d, ! if (ispecies .gt. 2) then do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then ! ! 1) DI ! @@ -4780,7 +4792,7 @@ subroutine step_rate_g(de, HI, HII, HeI, HeII, HeIII, d, ! if (ispecies .gt. 3) then do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then ! ! 1) DM ! @@ -4832,7 +4844,7 @@ subroutine step_rate_g(de, HI, HII, HeI, HeII, HeIII, d, ! if (imchem .eq. 1) then do i = is+1, ie+1 - if (itmask_metal(i)) then + if (itmask_metal(i) .ne. MASK_FALSE) then C***** CI ********** scoef = 0._DKIND + 12._DKIND * ( 0._DKIND @@ -5272,7 +5284,7 @@ subroutine step_rate_g(de, HI, HII, HeI, HeII, HeIII, d, ! if ( ( igrgr .eq. 1 ) .or. ( idsub .eq. 1) ) then do i = is+1, ie+1 - if (itmask_metal(i)) then + if (itmask_metal(i) .ne. MASK_FALSE) then if (idspecies .gt. 0) then C***** MgSiO3 ********** @@ -5417,7 +5429,7 @@ subroutine step_rate_g(de, HI, HII, HeI, HeII, HeIII, d, ! --- (E) Set densities from 1D temps to 3D fields --- do i = is+1, ie+1 - if (itmask(i)) then + if (itmask(i) .ne. MASK_FALSE) then HIdot_prev(i) = abs(HI(i,j,k)-HIp(i)) / & max(real(dtit(i), DKIND), tiny8) HI(i,j,k) = max(real(HIp(i), RKIND), tiny) @@ -5441,7 +5453,8 @@ subroutine step_rate_g(de, HI, HII, HeI, HeII, HeIII, d, if (ispecies .gt. 3) & de(i,j,k) = de(i,j,k) - DM(i,j,k)/2._RKIND & + HDII(i,j,k)/3._RKIND + HeHII(i,j,k)/5._RKIND - if ( (imchem .eq. 1) .and. (itmask_metal(i)) ) + if ( (imchem .eq. 1) .and. + & (itmask_metal(i) .ne. MASK_FALSE) ) & de(i,j,k) = de(i,j,k) & + CII(i,j,k)/12._RKIND + COII(i,j,k)/28._RKIND & + OII(i,j,k)/16._RKIND + OHII(i,j,k)/17._RKIND @@ -5469,7 +5482,8 @@ subroutine step_rate_g(de, HI, HII, HeI, HeII, HeIII, d, HeHII(i,j,k) = max(real(HeHIIp(i), RKIND), tiny) endif - if ( (imchem .eq. 1) .and. (itmask_metal(i)) ) then + if ( (imchem .eq. 1) .and. + & (itmask_metal(i) .ne. MASK_FALSE) ) then CI(i,j,k) = max(real(CIp(i) , RKIND), tiny) CII(i,j,k) = max(real(CIIp(i) , RKIND), tiny) CO(i,j,k) = max(real(COp(i) , RKIND), tiny) @@ -5647,7 +5661,6 @@ subroutine make_consistent_g(de, HI, HII, HeI, HeII, HeIII, & , Sig(in), Sg(in), Feg(in) real*8 Cd(in), Od(in), Mgd(in), Ald(in) & , Sid(in), Sd(in), Fed(in) -! logical itmask_metal(in) ! Loop over all zones diff --git a/src/include/grackle_fortran_types.def b/src/include/grackle_fortran_types.def index 0dcff7b8..bf352c56 100644 --- a/src/include/grackle_fortran_types.def +++ b/src/include/grackle_fortran_types.def @@ -14,6 +14,12 @@ #include "grackle_float.h" +#define MASK_TYPE integer*4 + integer, parameter :: MASK_KIND=4 +#define MASK_TRUE 1_MASK_KIND +#define MASK_FALSE 0_MASK_KIND + + #ifdef GRACKLE_FLOAT_4 #define tiny 1.e-20 #define huge 1.e+20 From 2d24754b0541254764ad6d0f97c2854ce2954df4 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Sun, 27 Oct 2024 17:42:56 -0400 Subject: [PATCH 05/11] minor refactoring of cool1d_multi_g so that it is easier to transcribe --- src/clib/cool1d_multi_g.F | 46 +++++++++++++++++++-------------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/src/clib/cool1d_multi_g.F b/src/clib/cool1d_multi_g.F index c59bba07..b5d02638 100644 --- a/src/clib/cool1d_multi_g.F +++ b/src/clib/cool1d_multi_g.F @@ -361,10 +361,10 @@ subroutine cool1d_multi_g( MASK_TYPE itmask(in), anydust, interp MASK_TYPE itmask_metal(in), itmask_tab(in) !!#define CALCULATE_TGAS_SELF_CONSISTENTLY -#ifdef CALCULATE_TGAS_SELF_CONSISTENTLY +!#ifdef CALCULATE_TGAS_SELF_CONSISTENTLY integer iter_tgas real*8 tgas_err, tgas0 -#endif /* NOT important */ +!#endif /* NOT important */ ! debug real*8 edotunit integer i_max @@ -519,36 +519,36 @@ subroutine cool1d_multi_g( nother = (HeI(i,j,k) + HeII(i,j,k) + & HeIII(i,j,k))/4._DKIND + & HI(i,j,k) + HII(i,j,k) + de(i,j,k) -#ifdef CALCULATE_TGAS_SELF_CONSISTENTLY + iter_tgas = 0 tgas_err = huge8 do while ((iter_tgas .lt. 100) & .and.(tgas_err .gt. 1.d-3)) - tgas0 = tgas(i) -#endif - if (nH2/nother .gt. 1.0e-3_DKIND) then - x = 6100._DKIND/tgas(i) ! not quite self-consistent - if (x .gt. 10._DKIND) then - gamma2 = 0.5_DKIND*5._DKIND + tgas0 = tgas(i) + if (nH2/nother .gt. 1.0e-3_DKIND) then + x = 6100._DKIND/tgas(i) ! not quite self-consistent + if (x .gt. 10._DKIND) then + gamma2 = 0.5_DKIND*5._DKIND + else + gamma2 = 0.5_DKIND*(5._DKIND + 2._DKIND*x**2 + & * exp(x)/(exp(x)-1)**2) + endif else - gamma2 = 0.5_DKIND*(5._DKIND + 2._DKIND*x**2 * - & exp(x)/(exp(x)-1)**2) + gamma2 = 2.5_DKIND endif - else - gamma2 = 2.5_DKIND - endif - gamma2 = 1._DKIND + (nH2 + nother)/ - & (nH2*gamma2 + nother/(gamma-1._DKIND)) + gamma2 = 1._DKIND + (nH2 + nother)/ + & (nH2*gamma2 + nother/(gamma-1._DKIND)) #ifdef CALCULATE_TGAS_SELF_CONSISTENTLY - tgas(i) = max((gamma2 - 1._DKIND)*mmw(i)*e(i,j,k)*utem - & , temstart) - tgas_err = dabs(tgas0 - tgas(i)) / tgas0 - iter_tgas = iter_tgas + 1 - end do + tgas(i) = max((gamma2 - 1._DKIND)*mmw(i)*e(i,j,k)* + & utem, temstart) + tgas_err = dabs(tgas0 - tgas(i)) / tgas0 + iter_tgas = iter_tgas + 1 #else - tgas(i) = tgas(i) * (gamma2 - 1._DKIND)/ - & (gamma - 1._DKIND) + tgas(i) = tgas(i) * (gamma2 - 1._DKIND)/ + & (gamma - 1._DKIND) + iter_tgas = 101 #endif + end do end if enddo endif From 0de37a139af12d1b189fefca248a95f369736170 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Sun, 27 Oct 2024 17:43:17 -0400 Subject: [PATCH 06/11] fixing a type declaration. --- src/clib/solve_rate_cool_g.F | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/clib/solve_rate_cool_g.F b/src/clib/solve_rate_cool_g.F index 478d1df7..be3d3730 100644 --- a/src/clib/solve_rate_cool_g.F +++ b/src/clib/solve_rate_cool_g.F @@ -4166,7 +4166,7 @@ subroutine step_rate_g(de, HI, HII, HeI, HeII, HeIII, d, & , OHIIp(in) , H2OIIp(in) , H3OIIp(in) & , O2IIp(in) , Mgp(in) , Alp(in) & , Sp(in) , Fep(in) - real*8 SiMp(in), FeMp(in), Mg2SiO4p(in) + R_PREC SiMp(in), FeMp(in), Mg2SiO4p(in) & , MgSiO3p(in), Fe3O4p(in), ACp(in) & , SiO2Dp(in), MgOp(in), FeSp(in) & , Al2O3p(in) From ecc61f2d6f765007d4ea63f3193d8522b8f8cb6f Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Wed, 4 Dec 2024 09:52:06 -0500 Subject: [PATCH 07/11] reformatted some lines into a more-friendly format for transcription. --- src/clib/solve_rate_cool_g.F | 57 ++++++++++++++++++++++++------------ 1 file changed, 38 insertions(+), 19 deletions(-) diff --git a/src/clib/solve_rate_cool_g.F b/src/clib/solve_rate_cool_g.F index be3d3730..b725fdf7 100644 --- a/src/clib/solve_rate_cool_g.F +++ b/src/clib/solve_rate_cool_g.F @@ -5782,27 +5782,46 @@ subroutine make_consistent_g(de, HI, HII, HeI, HeII, HeIII, C enddo nSN = 12 - SN_i( 1) = 1; SN_metal(:, 1) = metal_loc(:,j,k) - SN_i( 2) = 2; SN_metal(:, 2) = metal_C13(:,j,k) - SN_i( 3) = 3; SN_metal(:, 3) = metal_C20(:,j,k) - SN_i( 4) = 4; SN_metal(:, 4) = metal_C25(:,j,k) - SN_i( 5) = 5; SN_metal(:, 5) = metal_C30(:,j,k) - SN_i( 6) = 6; SN_metal(:, 6) = metal_F13(:,j,k) - SN_i( 7) = 7; SN_metal(:, 7) = metal_F15(:,j,k) - SN_i( 8) = 8; SN_metal(:, 8) = metal_F50(:,j,k) - SN_i( 9) = 9; SN_metal(:, 9) = metal_F80(:,j,k) - SN_i(10) =10; SN_metal(:,10) = metal_P170(:,j,k) - SN_i(11) =11; SN_metal(:,11) = metal_P200(:,j,k) - SN_i(12) =12; SN_metal(:,12) = metal_Y19(:,j,k) + SN_i( 1) = 1 + SN_metal(:, 1) = metal_loc(:,j,k) + SN_i( 2) = 2 + SN_metal(:, 2) = metal_C13(:,j,k) + SN_i( 3) = 3 + SN_metal(:, 3) = metal_C20(:,j,k) + SN_i( 4) = 4 + SN_metal(:, 4) = metal_C25(:,j,k) + SN_i( 5) = 5 + SN_metal(:, 5) = metal_C30(:,j,k) + SN_i( 6) = 6 + SN_metal(:, 6) = metal_F13(:,j,k) + SN_i( 7) = 7 + SN_metal(:, 7) = metal_F15(:,j,k) + SN_i( 8) = 8 + SN_metal(:, 8) = metal_F50(:,j,k) + SN_i( 9) = 9 + SN_metal(:, 9) = metal_F80(:,j,k) + SN_i(10) =10 + SN_metal(:,10) = metal_P170(:,j,k) + SN_i(11) =11 + SN_metal(:,11) = metal_P200(:,j,k) + SN_i(12) =12 + SN_metal(:,12) = metal_Y19(:,j,k) do i = is+1, ie+1 - Ct(i) = 0._DKIND; Cg(i) = 0._DKIND - Ot(i) = 0._DKIND; Og(i) = 0._DKIND - Mgt(i) = 0._DKIND; Mgg(i) = 0._DKIND - Alt(i) = 0._DKIND; Alg(i) = 0._DKIND - Sit(i) = 0._DKIND; Sig(i) = 0._DKIND - St(i) = 0._DKIND; Sg(i) = 0._DKIND - Fet(i) = 0._DKIND; Feg(i) = 0._DKIND + Ct(i) = 0._DKIND + Cg(i) = 0._DKIND + Ot(i) = 0._DKIND + Og(i) = 0._DKIND + Mgt(i) = 0._DKIND + Mgg(i) = 0._DKIND + Alt(i) = 0._DKIND + Alg(i) = 0._DKIND + Sit(i) = 0._DKIND + Sig(i) = 0._DKIND + St(i) = 0._DKIND + Sg(i) = 0._DKIND + Fet(i) = 0._DKIND + Feg(i) = 0._DKIND do iSN = 1, nSN iSN0 = SN_i(iSN) From db6d0dd0eb0227221c8c4403b2cc34779318f6c0 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Wed, 4 Dec 2024 09:58:22 -0500 Subject: [PATCH 08/11] Removed unused argument from lookup_cool_rates0d --- src/clib/lookup_cool_rates0d.F | 4 ++-- src/clib/solve_rate_cool_g.F | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/clib/lookup_cool_rates0d.F b/src/clib/lookup_cool_rates0d.F index 11e1f7f9..53bd6e07 100644 --- a/src/clib/lookup_cool_rates0d.F +++ b/src/clib/lookup_cool_rates0d.F @@ -3,7 +3,7 @@ ! calculate rates - subroutine lookup_cool_rates0d(output, dtit, + subroutine lookup_cool_rates0d(dtit, & d, u, v, w, & nsp, dsp, dspdot, nratec, & iexpand, ispecies, imetal, imcool, @@ -133,7 +133,7 @@ subroutine lookup_cool_rates0d(output, dtit, ! General Arguments - integer output, nratec, nsp, + integer nratec, nsp, & iexpand, ih2co, ipiht, ispecies, imetal, idim, & imcool, idust, idustall, idustfield, idustrec, & igammah, ih2optical, iciecool, ih2cr, ihdcr, ithreebody, diff --git a/src/clib/solve_rate_cool_g.F b/src/clib/solve_rate_cool_g.F index b725fdf7..26ddfd1f 100644 --- a/src/clib/solve_rate_cool_g.F +++ b/src/clib/solve_rate_cool_g.F @@ -1471,7 +1471,7 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de, go to 9996 endif - call lookup_cool_rates0d(itr, dtit(i), + call lookup_cool_rates0d(dtit(i), & d(i,j,k), u(i,j,k), v(i,j,k), w(i,j,k), & nsp, dsp, dspdot, nratec, & iexpand, ispecies, imetal, imcool, @@ -1607,7 +1607,7 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de, endif enddo - call lookup_cool_rates0d(1, dtit(i), + call lookup_cool_rates0d(dtit(i), & d(i,j,k), u(i,j,k), v(i,j,k), w(i,j,k), & nsp, dsp1, dspdot1, nratec, & iexpand, ispecies, imetal, imcool, From a0f244cfbec870549f4d825a988d2644ffc08167 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Wed, 4 Dec 2024 10:04:03 -0500 Subject: [PATCH 09/11] Define dummy_iter_arg within cool_multi_time_g The `cool1d_multi_g` subroutine expects an argument `iter`. Previously, the `cool_multi_time_g` subroutine would pass a value of 1 directly to this argument. To simplify the process of transcription, we now store the value inside of a local variable called `dummy_iter_arg` (if we didn't do this, we would need to insert logic into our transcription routine to inject a custom variable to hold the value of 1 and then pass a pointer to that value into `cool1d_multi_g`). --- src/clib/cool_multi_time_g.F | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/clib/cool_multi_time_g.F b/src/clib/cool_multi_time_g.F index 5e3e0a84..a39a6fa1 100644 --- a/src/clib/cool_multi_time_g.F +++ b/src/clib/cool_multi_time_g.F @@ -283,6 +283,7 @@ subroutine cool_multi_time_g( & brem(in), cieco(in), & hyd01k(in), h2k01(in), vibh(in), roth(in), rotl(in), & gpldl(in), gphdl(in), hdlte(in), hdlow(in) + integer dummy_iter_arg ! Iteration mask for multi_cool @@ -351,13 +352,15 @@ subroutine cool_multi_time_g( end do ! Compute the cooling rate + dummy_iter_arg=1 call cool1d_multi_g( & d, e, u, v, w, de, HI, HII, HeI, HeII, HeIII, & in, jn, kn, nratec, & iexpand, ispecies, imetal, imcool, & idust, idustall, idustfield, idustrec, - & idim, is, ie, j, k, ih2co, ipiht, 1, igammah, + & idim, is, ie, j, k, ih2co, ipiht, + & dummy_iter_arg, igammah, & aye, temstart, temend, z_solar, fgr, & utem, uxyz, uaye, urho, utim, & gamma, fh, From 2916da1f66503305326eff4219bcfa544f3ab56c Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Wed, 4 Dec 2024 10:22:36 -0500 Subject: [PATCH 10/11] update readme for this special development branch --- README.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/README.md b/README.md index 3953979a..f3b347e9 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,9 @@ +> [!IMPORTANT] +> This is a special version of the [brittonsmith:gen2024](https://github.com/brittonsmith/grackle/tree/gen2024) branch (i.e. the branch of changes proposed for merging in the grackle-project/grackle#177 Pull Request). +> +> This branch includes additional changes that are needed to simplify the transcription process to C++. There are pending PRs to merge all of these changes into the gen2024 branch ([see this list of PRs](https://github.com/brittonsmith/grackle/pulls/mabruzzo)) + + # Grackle [![Users' Mailing List](https://img.shields.io/badge/Users-List-lightgrey.svg)](https://groups.google.com/forum/#!forum/grackle-cooling-users) From 234a6769a81444119337330c17b5decf4ff1af2c Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Wed, 4 Dec 2024 10:31:11 -0500 Subject: [PATCH 11/11] Minor tweaks to improve transcription I made some minor tweaks to where we pass in the pointers to the i-dimension of grid_dimensions, grid_start, and grid_end. This is a purely superficial change, but will allow automated transcription to produce better code. --- src/clib/calculate_cooling_time.c | 6 +++--- src/clib/calculate_dust_temperature.c | 6 +++--- src/clib/calculate_temperature.c | 6 +++--- src/clib/solve_chemistry.c | 6 +++--- 4 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/clib/calculate_cooling_time.c b/src/clib/calculate_cooling_time.c index 029fe898..4dbe9573 100644 --- a/src/clib/calculate_cooling_time.c +++ b/src/clib/calculate_cooling_time.c @@ -224,7 +224,7 @@ int local_calculate_cooling_time(chemistry_data *my_chemistry, my_fields->HeII_density, my_fields->HeIII_density, cooling_time, - my_fields->grid_dimension, + my_fields->grid_dimension+0, my_fields->grid_dimension+1, my_fields->grid_dimension+2, &my_chemistry->NumberOfTemperatureBins, @@ -237,10 +237,10 @@ int local_calculate_cooling_time(chemistry_data *my_chemistry, &my_chemistry->use_dust_density_field, &my_chemistry->dust_recombination_cooling, &(my_fields->grid_rank), - my_fields->grid_start, + my_fields->grid_start+0, my_fields->grid_start+1, my_fields->grid_start+2, - my_fields->grid_end, + my_fields->grid_end+0, my_fields->grid_end+1, my_fields->grid_end+2, &my_chemistry->ih2co, diff --git a/src/clib/calculate_dust_temperature.c b/src/clib/calculate_dust_temperature.c index 0030a5bf..d9029346 100644 --- a/src/clib/calculate_dust_temperature.c +++ b/src/clib/calculate_dust_temperature.c @@ -135,17 +135,17 @@ int local_calculate_dust_temperature(chemistry_data *my_chemistry, my_fields->HM_density, my_fields->H2I_density, my_fields->H2II_density, - my_fields->grid_dimension, + my_fields->grid_dimension+0, my_fields->grid_dimension+1, my_fields->grid_dimension+2, &my_chemistry->NumberOfTemperatureBins, &my_units->comoving_coordinates, &my_chemistry->primordial_chemistry, &(my_fields->grid_rank), - my_fields->grid_start, + my_fields->grid_start+0, my_fields->grid_start+1, my_fields->grid_start+2, - my_fields->grid_end, + my_fields->grid_end+0, my_fields->grid_end+1, my_fields->grid_end+2, &my_units->a_value, diff --git a/src/clib/calculate_temperature.c b/src/clib/calculate_temperature.c index 4466db60..c684f608 100644 --- a/src/clib/calculate_temperature.c +++ b/src/clib/calculate_temperature.c @@ -189,15 +189,15 @@ int local_calculate_temperature_table(chemistry_data *my_chemistry, my_fields->internal_energy, my_fields->metal_density, temperature, - my_fields->grid_dimension, + my_fields->grid_dimension+0, my_fields->grid_dimension+1, my_fields->grid_dimension+2, &my_units->comoving_coordinates, &metal_field_present, - my_fields->grid_start, + my_fields->grid_start+0, my_fields->grid_start+1, my_fields->grid_start+2, - my_fields->grid_end, + my_fields->grid_end+0, my_fields->grid_end+1, my_fields->grid_end+2, &my_units->a_value, diff --git a/src/clib/solve_chemistry.c b/src/clib/solve_chemistry.c index 9dcb4add..e5b305f8 100644 --- a/src/clib/solve_chemistry.c +++ b/src/clib/solve_chemistry.c @@ -257,7 +257,7 @@ int local_solve_chemistry(chemistry_data *my_chemistry, my_fields->HeI_density, my_fields->HeII_density, my_fields->HeIII_density, - my_fields->grid_dimension, + my_fields->grid_dimension+0, my_fields->grid_dimension+1, my_fields->grid_dimension+2, &my_chemistry->NumberOfTemperatureBins, @@ -269,10 +269,10 @@ int local_solve_chemistry(chemistry_data *my_chemistry, &my_chemistry->dust_chemistry, &my_chemistry->use_dust_density_field, &(my_fields->grid_rank), - my_fields->grid_start, + my_fields->grid_start+0, my_fields->grid_start+1, my_fields->grid_start+2, - my_fields->grid_end, + my_fields->grid_end+0, my_fields->grid_end+1, my_fields->grid_end+2, &my_chemistry->ih2co,