Implemented Sunny Wong rotation brunt change

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

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

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) &

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

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

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
 

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