diff --git a/star/private/brunt.f90 b/star/private/brunt.f90 index 19cc1527e..2e68fefae 100644 --- a/star/private/brunt.f90 +++ b/star/private/brunt.f90 @@ -56,13 +56,14 @@ subroutine do_brunt_B(s,nzlo,nzhi,ierr) ierr = 0 nz = s% nz - - if (.not. s% calculate_Brunt_B) then - call set_nan(s% brunt_B(1:nz)) - call set_nan(s% unsmoothed_brunt_B(1:nz)) - return + + if (s% fill_arrays_with_nans) then + call set_nan(s% brunt_B(nzlo:nzhi)) + call set_nan(s% unsmoothed_brunt_B(nzlo:nzhi)) end if + if (.not. s% calculate_Brunt_B) return + if (s% use_other_brunt) then call s% other_brunt(s% id, ierr) if (ierr /= 0) then @@ -88,6 +89,7 @@ subroutine do_brunt_B(s,nzlo,nzhi,ierr) end do allocate(smoothing_array(nz)) + if(s% fill_arrays_with_nans) call fill_with_nans(smoothing_array) call smooth_brunt_B(smoothing_array) if (s% use_other_brunt_smoothing) then call s% other_brunt_smoothing(s% id, ierr) @@ -132,15 +134,20 @@ subroutine do_brunt_N2(s,nzlo,nzhi,ierr) ierr = 0 nz = s% nz - if (.not. (s% calculate_Brunt_B .and. s% calculate_Brunt_N2)) then - call set_nan(s% brunt_N2(1:nz)) - call set_nan(s% brunt_N2_composition_term(1:nz)) - return + if (s% fill_arrays_with_NaNs) then + call set_nan(s% brunt_N2(nzlo:nzhi)) + call set_nan(s% brunt_N2_composition_term(nzlo:nzhi)) end if + if (.not. (s% calculate_Brunt_B .and. s% calculate_Brunt_N2)) return + allocate(rho_P_chiT_chiRho(nz), rho_P_chiT_chiRho_face(nz)) + if(s% fill_arrays_with_NaNs) then + call fill_with_nans(rho_P_chiT_chiRho) + call fill_with_nans(rho_P_chiT_chiRho_face) + end if - do k=1,nz + do k=nzlo,nzhi rho_P_chiT_chiRho(k) = (s% rho(k)/s% Peos(k))*(s% chiT(k)/s% chiRho(k)) ! correct for difference between gravitational mass density and baryonic mass density (rho) if (s% use_mass_corrections) then @@ -156,7 +163,7 @@ subroutine do_brunt_N2(s,nzlo,nzhi,ierr) return end if - do k=1,nz ! clip B and calculate N^2 from B + do k=nzlo,nzhi ! clip B and calculate N^2 from B if (abs(s% brunt_B(k)) < s% min_magnitude_brunt_B .or. & s% gradT(k) == 0 .or. is_bad(s% gradT_sub_grada(k))) then s% brunt_B(k) = 0 @@ -178,7 +185,7 @@ subroutine do_brunt_N2(s,nzlo,nzhi,ierr) end do if (s% brunt_N2_coefficient /= 1d0) then - do k=1,nz + do k=nzlo,nzhi s% brunt_N2(k) = s% brunt_N2_coefficient*s% brunt_N2(k) s% brunt_N2_composition_term(k) = & s% brunt_N2_coefficient*s% brunt_N2_composition_term(k) @@ -210,6 +217,12 @@ subroutine do_brunt_B_MHM_form(s, nzlo, nzhi, ierr) species = s% species allocate(T_face(nz), rho_face(nz), chiT_face(nz), chiRho_face(nz)) + if(s% fill_arrays_with_nans) then + call fill_with_nans(T_face) + call fill_with_nans(rho_face) + call fill_with_nans(chiT_face) + call fill_with_nans(chiRho_face) + end if call get_face_values(s, s% chiT, chiT_face, ierr) if (ierr /= 0) return diff --git a/star/private/photo_in.f90 b/star/private/photo_in.f90 index 25b393b01..73c8e04a0 100644 --- a/star/private/photo_in.f90 +++ b/star/private/photo_in.f90 @@ -147,8 +147,8 @@ subroutine read_star_photo(s, fname, ierr) read(iounit, iostat=ierr) & s% dq(1:nz), s% xa(:,1:nz), s% xh(:,1:nz), & s% omega(1:nz), s% j_rot(1:nz), s% mlt_vc(1:nz), s% conv_vel(1:nz), & - s% D_ST_start(1:nz), s% nu_ST_start(1:nz), & ! needed for ST time smoothing - s% have_ST_start_info + s% D_ST_start(1:nz), s% nu_ST_start(1:nz), & ! needed for ST time smoothing + s% have_ST_start_info, s% mstar_old ! needed for ST time smoothing call read_part_number(iounit) if (failed('rsp_num_periods')) return diff --git a/star/private/photo_out.f90 b/star/private/photo_out.f90 index 9eed9c50c..a9d487e29 100644 --- a/star/private/photo_out.f90 +++ b/star/private/photo_out.f90 @@ -84,8 +84,8 @@ subroutine output_star_photo(s,iounit,ierr) write(iounit) & s% dq(1:nz), s% xa(:,1:nz), s% xh(:,1:nz), & s% omega(1:nz), s% j_rot(1:nz), s% mlt_vc(1:nz), s% conv_vel(1:nz), & - s% D_ST_start(1:nz), s% nu_ST_start(1:nz), & ! needed for ST time smoothing - s% have_ST_start_info + s% D_ST_start(1:nz), s% nu_ST_start(1:nz), & ! needed for ST time smoothing + s% have_ST_start_info, s% mstar_old ! needed for ST time smoothing call write_part_number(iounit) write(iounit) & diff --git a/star/private/rotation_mix_info.f90 b/star/private/rotation_mix_info.f90 index 527be285e..bae903b1c 100644 --- a/star/private/rotation_mix_info.f90 +++ b/star/private/rotation_mix_info.f90 @@ -354,6 +354,18 @@ subroutine mix_D_omega allocate(sig(nz), rhs(nz), d(nz), du(nz), dl(nz), bp(nz), vp(nz), xp(nz), x(nz)) + if(s% fill_arrays_with_NaNs) then + call fill_with_NaNs(sig) + call fill_with_NaNs(rhs) + call fill_with_NaNs(d) + call fill_with_NaNs(du) + call fill_with_NaNs(dl) + call fill_with_NaNs(bp) + call fill_with_NaNs(vp) + call fill_with_NaNs(xp) + call fill_with_NaNs(x) + end if + rate = min(s% D_omega_mixing_rate, 1d0/dt) do k=2,nz-1 if (s% D_omega(k) == 0 .or. s% D_omega(k+1) == 0) then @@ -479,6 +491,62 @@ subroutine setup(ierr) smooth_work(nz,num_instabilities), saved(nz,num_instabilities), & unstable(num_instabilities,nz)) + if(s% fill_arrays_with_NaNs) then + call fill_with_NaNs(csound) + call fill_with_NaNs(rho) + call fill_with_NaNs(T) + call fill_with_NaNs(P) + call fill_with_NaNs(cp) + call fill_with_NaNs(cv) + call fill_with_NaNs(chiRho) + call fill_with_NaNs(abar) + call fill_with_NaNs(zbar) + call fill_with_NaNs(opacity) + call fill_with_NaNs(kap_cond) + call fill_with_NaNs(gamma1) + call fill_with_NaNs(mu_alt) + call fill_with_NaNs(omega) + call fill_with_NaNs(cell_dr) + call fill_with_NaNs(eps_nuc) + call fill_with_NaNs(enu) + call fill_with_NaNs(L_neu) + call fill_with_NaNs(gradT_sub_grada) + call fill_with_NaNs(delta) + call fill_with_NaNs(scale_height) + call fill_with_NaNs(dRho) + call fill_with_NaNs(dr) + call fill_with_NaNs(dPressure) + call fill_with_NaNs(domega) + call fill_with_NaNs(d_mu) + call fill_with_NaNs(d_j_rot) + call fill_with_NaNs(dRho_dr) + call fill_with_NaNs(dRho_dr_ad) + call fill_with_NaNs(dr2omega) + call fill_with_NaNs(H_T) + call fill_with_NaNs(domega_dlnR) + call fill_with_NaNs(Hj) + call fill_with_NaNs(dlnR_domega) + call fill_with_NaNs(t_dyn) + call fill_with_NaNs(t_kh) + call fill_with_NaNs(Ri_mu) + call fill_with_NaNs(Ri_T) + call fill_with_NaNs(ve0) + call fill_with_NaNs(ve_mu) + call fill_with_NaNs(v_ssi) + call fill_with_NaNs(h_ssi) + call fill_with_NaNs(Ris_1) + call fill_with_NaNs(Ris_2) + call fill_with_NaNs(v_es) + call fill_with_NaNs(H_es) + call fill_with_NaNs(v_gsf) + call fill_with_NaNs(H_gsf) + call fill_with_NaNs(N2) + call fill_with_NaNs(N2_mu) + call fill_with_NaNs_2d(smooth_work) + call fill_with_NaNs_2d(saved) + end if + + ! interpolate by mass to get values at cell boundaries enu00 = s% eps_nuc_neu_total(1) + s% non_nuc_neu(1) enu(1) = enu00 @@ -589,10 +657,8 @@ subroutine setup(ierr) Hj(1) = Hj(2); Hj(nz) = Hj(nz-1) do i = 1, nz - dlnRho_dlnP = s% grad_density(i) - dlnT_dlnP = s% grad_temperature(i) - N2(i) = -grav(i)*(1/gamma1(i) - dlnRho_dlnP)/scale_height(i) - N2_mu(i) = -grav(i)/scale_height(i)*(1/chiRho(i) - delta(i)*dlnT_dlnP - dlnRho_dlnP) + N2(i) = s%brunt_N2(i) + N2_mu(i) = s%brunt_N2_composition_term(i) end do do k=1,nz diff --git a/star/test_suite/do1_test_source b/star/test_suite/do1_test_source index 52d1bc469..3ad101721 100755 --- a/star/test_suite/do1_test_source +++ b/star/test_suite/do1_test_source @@ -4,6 +4,10 @@ # Slow cases do_one ppisn "Successful test: evolved 100 days past first relax" "final.mod" x150 + + +return +do_one ppisn "Successful test: evolved 100 days past first relax" "final.mod" x150 do_one 1M_pre_ms_to_wd "stop because log_surface_luminosity <= log_L_lower_limit" "final.mod" auto do_one ccsn_IIp "shock has reached target location 1" "shock_part5.mod" auto do_one 1M_thermohaline "all values are within tolerances" "final.mod" auto diff --git a/star/test_suite/ppisn/inlist_ppisn b/star/test_suite/ppisn/inlist_ppisn index 20554ac7d..10750d6ac 100644 --- a/star/test_suite/ppisn/inlist_ppisn +++ b/star/test_suite/ppisn/inlist_ppisn @@ -64,7 +64,6 @@ &controls - fill_arrays_with_nans = .false. report_solver_progress = .false. max_resid_jump_limit = 1d99 diff --git a/star_data/public/star_data_def.inc b/star_data/public/star_data_def.inc index d0b7e1740..8cc626462 100644 --- a/star_data/public/star_data_def.inc +++ b/star_data/public/star_data_def.inc @@ -1,7 +1,7 @@ character(len=24) :: version_number ! mesa version from file $MESA_DIR/data/version_number - integer, parameter :: star_def_version = 17 + integer, parameter :: star_def_version = 18 integer, parameter :: nz_alloc_extra = 200