Fixing fpe's in binary (#486)

Work on fpe's in binary as well as a few other places in star.

Author: rjfarmer

binary/private/binary_evolve.f90

@@ -188,6 +188,7 @@ subroutine binarydata_init(b, doing_restart)
             b% CE_lambda2 = 0d0
             b% CE_Ebind1 = 0d0
             b% CE_Ebind2 = 0d0
+            b% mtransfer_rate = 0
 
             b% num_tries = 0
 
@@ -337,6 +338,11 @@ integer function binary_evolve_step(b)
             b% eccentricity = b% eccentricity + get_edot(b) *b% time_step*secyer
             if (b% eccentricity < b% min_eccentricity) b% eccentricity = b% min_eccentricity
             if (b% eccentricity > b% max_eccentricity) b% eccentricity = b% max_eccentricity
+         else
+            b% edot_tidal = 0d0
+            b% edot_enhance = 0d0
+            b% extra_edot = 0d0
+            b% edot = 0d0
          end if
 
          !use new eccentricity to calculate new time coordinate
@@ -504,17 +510,19 @@ integer function binary_finish_step(b)
 
          binary_finish_step = keep_going
          ! update change factor in case mtransfer_rate has changed
-         if(b% mtransfer_rate_old /= b% mtransfer_rate .and. &
-             b% mtransfer_rate /= 0 .and. b% mtransfer_rate_old /= 0) then
-            if(b% mtransfer_rate < b% mtransfer_rate_old) then
-               b% change_factor = b% change_factor*(1d0-b% implicit_lambda) + b% implicit_lambda* &
-                  (1+b% change_factor_fraction*(b% mtransfer_rate/b% mtransfer_rate_old-1))
-            else
-               b% change_factor = b% change_factor*(1d0-b% implicit_lambda) + b% implicit_lambda* &
-                   (1+b% change_factor_fraction*(b% mtransfer_rate_old/b% mtransfer_rate-1))
+         if(.not. b% doing_first_model_of_run) then
+            if(b% mtransfer_rate_old /= b% mtransfer_rate .and. &
+               b% mtransfer_rate /= 0 .and. b% mtransfer_rate_old /= 0) then
+               if(b% mtransfer_rate < b% mtransfer_rate_old) then
+                  b% change_factor = b% change_factor*(1d0-b% implicit_lambda) + b% implicit_lambda* &
+                     (1+b% change_factor_fraction*(b% mtransfer_rate/b% mtransfer_rate_old-1))
+               else
+                  b% change_factor = b% change_factor*(1d0-b% implicit_lambda) + b% implicit_lambda* &
+                     (1+b% change_factor_fraction*(b% mtransfer_rate_old/b% mtransfer_rate-1))
+               end if
+               if(b% change_factor > b% max_change_factor) b% change_factor = b% max_change_factor
+               if(b% change_factor < b% min_change_factor) b% change_factor = b% min_change_factor
             end if
-            if(b% change_factor > b% max_change_factor) b% change_factor = b% max_change_factor
-            if(b% change_factor < b% min_change_factor) b% change_factor = b% min_change_factor
          end if
 
          ! store all variables into "old"

binary/private/binary_timestep.f90

@@ -126,39 +126,46 @@ integer function binary_pick_next_timestep(b)
 
          if (b% max_timestep < 0) b% max_timestep = b% s_donor% dt
 
-
          b% env(b% d_i) = s% star_mass - s% he_core_mass
          if (b% point_mass_i == 0) then
             b% env(b% a_i) = b% s_accretor% star_mass - b% s_accretor% he_core_mass
          end if
 
-         if (b% env_old(b% d_i) /= 0) then
-            env_change = b% env(b% d_i) - b% env_old(b% d_i)
+         if(.not. b% doing_first_model_of_run) then
+            if (b% env_old(b% d_i) /= 0) then
+               env_change = b% env(b% d_i) - b% env_old(b% d_i)
+            else
+               env_change = 0
+            end if
+            
+            if (b% rl_relative_gap_old(b% d_i) /= 0) then
+               rel_gap_change = b% rl_relative_gap_old(b% d_i) - b% rl_relative_gap(b% d_i)
+            else
+               rel_gap_change = 0
+            end if
+            
+            if (b% angular_momentum_j_old /= 0) then
+               j_change = b% angular_momentum_j - b% angular_momentum_j_old
+            else
+               j_change = 0
+            end if
+            
+            if (b% separation_old /= 0) then
+               sep_change = b% separation - b% separation_old
+            else
+               sep_change = 0
+            end if
+            if (b% eccentricity_old /= 0) then
+               e_change = b% eccentricity - b% eccentricity_old
+            else
+               e_change = 0
+            end if
          else
             env_change = 0
-         end if
-         
-         if (b% rl_relative_gap_old(b% d_i) /= 0) then
-            rel_gap_change = b% rl_relative_gap_old(b% d_i) - b% rl_relative_gap(b% d_i)
-         else
             rel_gap_change = 0
-         end if
-         
-         if (b% angular_momentum_j_old /= 0) then
-            j_change = b% angular_momentum_j - b% angular_momentum_j_old
-         else
             j_change = 0
-         end if
-         
-         if (b% separation_old /= 0) then
-            sep_change = b% separation - b% separation_old
-         else
             sep_change = 0
-         end if
-         if (b% eccentricity_old /= 0) then
-             e_change = b% eccentricity - b% eccentricity_old
-         else
-             e_change = 0
+            e_change = 0
          end if
 
          ! get limits for dt based on relative changes

star/private/history.f90

@@ -1950,10 +1950,12 @@ subroutine history_getval(&
             end do
 
          case(h_i_rot_total)
-            val = 0d0
-            do k = 1, s% nz
-               val = val + s% dm_bar(k) * s%i_rot(k)% val
-            end do
+            if(s% rotation_flag) then
+               val = 0d0
+               do k = 1, s% nz
+                  val = val + s% dm_bar(k) * s%i_rot(k)% val
+               end do
+            end if
          case(h_surf_avg_j_rot)
             val = if_rot(s% j_rot_avg_surf)
          case(h_surf_avg_omega)
@@ -2728,9 +2730,9 @@ subroutine history_getval(&
             int_val = s% k_below_const_q
             is_int_val = .true.
          case (h_q_below_const_q)
-            val = s% q(s% k_below_const_q)
+            if(s% k_below_const_q>0) val = s% q(s% k_below_const_q)
          case (h_logxq_below_const_q)
-            val = safe_log10(sum(s% dq(1:s% k_below_const_q - 1)))
+            if(s% k_below_const_q>0) val = safe_log10(sum(s% dq(1:s% k_below_const_q - 1)))
 
          case (h_k_const_mass)
             int_val = s% k_const_mass

star/test_suite/make_env/src/run_star_extras.f90

@@ -358,6 +358,7 @@ subroutine create_env(id, s, ierr)
 
          s% model_number = 0
          s% star_age = 0
+         s% crystal_core_boundary_mass = -1d0 
 
          write(*,2) 'start.mod', nz
          call star_write_model(id, 'start.mod', ierr)
@@ -479,6 +480,7 @@ subroutine do1_cell(ierr) ! uses r(k), T_m1, P_m1, logRho_m1
                s% lnT(k) = logT*ln10
                s% rho(k) = exp10(logRho)
                s% T(k) = T_00
+               s% rho_start(k) = s% rho(k)
                call star_do_eos_for_cell(s% id, k, ierr)
                if (ierr /= 0) then
                   write(*, *) 'Call star_do_eos_for_cell failed', k
@@ -498,6 +500,7 @@ subroutine do1_cell(ierr) ! uses r(k), T_m1, P_m1, logRho_m1
                s% csound_start(k) = s% csound(k)
                s% mlt_gradT_fraction = -1d0
                s% adjust_mlt_gradT_fraction(k) = -1d0
+               s% gradL_composition_term(k) = 0d0
 
                ! skipping use_other_alpha_mlt and other_gradr_factor
                s% alpha_mlt(k) = s% mixing_length_alpha

star/test_suite/starspots/src/run_star_extras.f90

@@ -165,7 +165,7 @@ subroutine YREC_spots_other_mlt_results(id, k, MLT_option, &  ! NOTE: k=0 is a v
 
         !------------------------------
          !if (s% star_age >= 10d0) then
-         if (.not. s% doing_relax) then
+         if (.not. s% doing_relax .and. .not. s% doing_first_model_of_run) then
             xspot_of_r = (P - PB_i)/P 
             gradr_spot = gradr/( fspot*pow4(xspot_of_r) + 1d0 - fspot)
          else