! ***********************************************************************
!
!   Copyright (C) 2012  Bill Paxton
!
!   this file is part of mesa.
!
!   mesa is free software; you can redistribute it and/or modify
!   it under the terms of the gnu general library public license as published
!   by the free software foundation; either version 2 of the license, or
!   (at your option) any later version.
!
!   mesa is distributed in the hope that it will be useful, 
!   but without any warranty; without even the implied warranty of
!   merchantability or fitness for a particular purpose.  see the
!   gnu library general public license for more details.
!
!   you should have received a copy of the gnu library general public license
!   along with this software; if not, write to the free software
!   foundation, inc., 59 temple place, suite 330, boston, ma 02111-1307 usa
!
! *********************************************************************** 
      module run_binary_extras 

      use star_lib
      use star_def
      use const_def
      use math_lib
      use binary_def
      
      implicit none
      
      include "binary_test_suite_extras_def.inc"

      contains

      include "binary_test_suite_extras.inc"
      
      subroutine extras_binary_controls(binary_id, ierr)
         integer :: binary_id
         integer, intent(out) :: ierr
         type (binary_info), pointer :: b
         ierr = 0

         call binary_ptr(binary_id, b, ierr)
         if (ierr /= 0) then
            write(*,*) 'failed in binary_ptr'
            return
         end if

         ! Set these function pointers to point to the functions you wish to use in
         ! your run_binary_extras. Any which are not set, default to a null_ version
         ! which does nothing.
         b% how_many_extra_binary_history_header_items => how_many_extra_binary_history_header_items
         b% data_for_extra_binary_history_header_items => data_for_extra_binary_history_header_items
         b% how_many_extra_binary_history_columns => how_many_extra_binary_history_columns
         b% data_for_extra_binary_history_columns => data_for_extra_binary_history_columns

         b% extras_binary_startup=> extras_binary_startup
         b% extras_binary_start_step=> extras_binary_start_step
         b% extras_binary_check_model=> extras_binary_check_model
         b% extras_binary_finish_step => extras_binary_finish_step
         b% extras_binary_after_evolve=> extras_binary_after_evolve

         ! CARB prescription
         b% other_jdot_mb => jdot_mb_routine

         ! Once you have set the function pointers you want, 
         ! then uncomment this (or set it in your star_job inlist)
         ! to disable the printed warning message,
         b% warn_binary_extra =.false.
         
      end subroutine extras_binary_controls


      integer function how_many_extra_binary_history_header_items(binary_id)
         use binary_def, only: binary_info
         integer, intent(in) :: binary_id
         how_many_extra_binary_history_header_items = 0
      end function how_many_extra_binary_history_header_items

      subroutine data_for_extra_binary_history_header_items( &
           binary_id, n, names, vals, ierr)
         type (binary_info), pointer :: b
         integer, intent(in) :: binary_id, n
         character (len=maxlen_binary_history_column_name) :: names(n)
         real(dp) :: vals(n)
         integer, intent(out) :: ierr
         ierr = 0
         call binary_ptr(binary_id, b, ierr)
         if (ierr /= 0) then
            write(*,*) 'failed in binary_ptr'
            return
         end if
      end subroutine data_for_extra_binary_history_header_items


      integer function how_many_extra_binary_history_columns(binary_id)
         use binary_def, only: binary_info
         integer, intent(in) :: binary_id
         how_many_extra_binary_history_columns = 18
      end function how_many_extra_binary_history_columns



      subroutine data_for_extra_binary_history_columns(binary_id, n, names, vals, ierr)
         use const_def, only: dp
         type (binary_info), pointer :: b
         type (star_info), pointer :: s
         integer, intent(in) :: binary_id
         integer, intent(in) :: n
         character (len=maxlen_binary_history_column_name) :: names(n)
         real(dp) :: vals(n)
         integer, intent(out) :: ierr
         integer i,j
         real(dp) :: R_conv_base, R_core, M_conv_env, M_rad_env, M_core, M_env
         logical conv


         ierr = 0
         call binary_ptr(binary_id, b, ierr)
         if (ierr /= 0) then
            write(*,*) 'failed in binary_ptr'
            return
         end if

         s => b% s_donor



         R_conv_base = 0d0
         R_core = 0d0
         M_core = 0d0
         M_conv_env = 0d0
         M_rad_env = 0d0

         conv = .false.


         do i = 1, s% nz-1

           if(.not.conv.and.s% mlt_mixing_type(i)==1)then ! Find outermost convective cell

            conv = .true.

            do j = i, s% nz-1

             if(conv.and.s%mlt_mixing_type(j)/=1)then ! find base of convective envelope
              conv = .false.
              R_conv_base = s% rmid(j)
              M_conv_env = s% m(1) - s% m(j)
             end if

             if(s% X(j)<1d-2)exit! Exit at core surface

            end do

            R_core = s% rmid(j)
            M_core = s% m(j)
            M_rad_env = s% m(1) - M_conv_env - M_core

            exit

           end if ! END if(.not.conv.and.s% mlt_mixing_type(i)==1)then 

         end do ! END do i = 1, s% nz-1



            
         names(1) = "turnover_time"
         vals(1) = s% xtra(1)

         names(2) = "mag_field"
         vals(2) = s% xtra(2)

         names(3) = "conv_env_r"
         vals(3) = s% xtra(3)

         names(4) = "conv_env_m"
         vals(4) = s% xtra(4)

         names(5) = "sonic_cross_time"
         vals(5) = s% xtra(5)

         names(6) = "turnover_time"
         vals(6) = s% xtra(6)

         names(7) = "mass_transfer_beta"
         vals(7) = s% xtra(7)

         names(8) = "mdotSTABLE"
         vals(8) = s% xtra(8)

         names(9) = "mdotRG"
         vals(9) = s% xtra(9)



         ! ==========  STELLAR STRUCTURE ==============

         names(10) = 'R_conv_base'
         vals(10) =  R_conv_base/Rsun

         names(11) = 'R_core'
         vals(11) =  R_core/Rsun

         names(12) = 'M_conv_env'
         vals(12) =  M_conv_env/Msun

         names(13) = 'M_rad_env'
         vals(13) =  M_rad_env/Msun

         names(14) = 'M_core'
         vals(14) =  M_core/Msun

         names(15) = 'M_star'
         vals(15) =  s% m(1)/Msun

         names(16) = 'R_star'
         vals(16) =  s% r(1)/Rsun

         names(17) = "qconv_env"
         vals(17) = s% xtra(10)

         names(18) = "jdot_scale"
         vals(18) = s% xtra(11)



         
      end subroutine data_for_extra_binary_history_columns
      
      
      integer function extras_binary_startup(binary_id,restart,ierr)
         type (binary_info), pointer :: b
         type (star_info), pointer :: s
         integer, intent(in) :: binary_id
         integer, intent(out) :: ierr
         logical, intent(in) :: restart

         call binary_ptr(binary_id, b, ierr)

         if (ierr /= 0) then ! failure in  binary_ptr
            return
         end if

         call test_suite_startup(b, restart, ierr)


         s => b% s_donor

         ! MESSAGE
         write(*,*) ""
         write(*,*) ""
         write(*,*) "*************************************"
         write(*,*) "*************************************"
         write(*,*) " --- STARTING SIMULATION ---"
         write(*,*) "*************************************"
         write(*,*) "*************************************"
         write(*,*) ""
         write(*,*) ""


         ! Saving star mass ( in Msun )
         s% xtra(12) = b% m(b% d_i)/Msun 

         ! Saving mass transfer rate ( in Msun/yr )
         s% xtra(13) = (abs(b% mtransfer_rate)/Msun)*secyer


         extras_binary_startup = keep_going

      end function extras_binary_startup
      
      integer function extras_binary_start_step(binary_id,ierr)
         type (binary_info), pointer :: b
         type (star_info), pointer :: s
         integer, intent(in) :: binary_id
         integer, intent(out) :: ierr

         extras_binary_start_step = keep_going
         call binary_ptr(binary_id, b, ierr)
         if (ierr /= 0) then ! failure in  binary_ptr
            return
         end if

      
      end function  extras_binary_start_step
      
      !Return either keep_going, retry or terminate
      integer function extras_binary_check_model(binary_id)
         type (binary_info), pointer :: b
         integer, intent(in) :: binary_id
         integer :: ierr
         call binary_ptr(binary_id, b, ierr)
         if (ierr /= 0) then ! failure in  binary_ptr
            return
         end if  
         extras_binary_check_model = keep_going
        
      end function extras_binary_check_model

      
      
      ! returns either keep_going or terminate.
      ! note: cannot request retry; extras_check_model can do that.
      integer function extras_binary_finish_step(binary_id)
         type (binary_info), pointer :: b
         type (star_info), pointer :: s
         integer, intent(in) :: binary_id
         integer :: ierr
         integer :: check
         real(dp) :: massWD,mdotWD,mdotRG,mdotSTABLE

         call binary_ptr(binary_id, b, ierr)
         if (ierr /= 0) then ! failure in  binary_ptr
            return
         end if  


         s => b% s_donor



         ! ===============
         ! Checking whether the WD mass should increase or not 
         ! (it depends on the mass transfer rate)
         ! ===============


         massWD = b% m(b% a_i) ! WD mass in g

         mdotWD = b% mtransfer_rate
         mdotWD = (abs(mdotWD)/Msun)*secyer ! mass transfer rate in Msun/yr


         ! Criterion from Wolf et al. (2013, ApJ, 777, 136)

         mdotSTABLE = -8.97878d0 + 3.36623d0*(massWD/Msun) -1.10131d0*(massWD/Msun)*(massWD/Msun)
         mdotSTABLE = pow(10.d0,mdotSTABLE)

         mdotRG = -8.21704 + 2.89808*(massWD/Msun) -1.05277*(massWD/Msun)*(massWD/Msun)
         mdotRG = pow(10.d0,mdotRG)



        
         ! In the dwarf-nova-like regime, we assume that the mass is entirely lost due to nova eruptions
         if (mdotWD < mdotSTABLE) then
           b% mass_transfer_beta = 1.d0
         endif


         ! In the stable regime, we assume that the mass is entirely accreted due to stable hydrogen burning
         if ( (mdotWD > mdotSTABLE) .and. (mdotWD < mdotRG) ) then
           b% mass_transfer_beta = 0.d0
         endif


         ! In the red-giant-like regime, we assume that the maximum accretion is mdotRG, 
         ! i.e. the remaining mass will be eventually lost due to fast isotropic winds
         if (mdotWD > mdotRG) then
           b% mass_transfer_beta = 1.d0 - (mdotRG / mdotWD)
         endif



         ! From this moment on, the hydrogen envelope is gone
         if ( (s% star_mass - s% he_core_mass) .lt. 1d-12) then

           s% atm_option = 'table'
           s% atm_table = 'tau_1'

           s% always_use_dedt_form_of_energy_eqn = .true.
           s% convergence_ignore_equL_residuals = .false.
           s% okay_to_reduce_gradT_excess = .false.

           s% tol_correction_norm = 3d-4
           s% tol_max_correction = 3d-2

         end if


         s% xtra(7) = b% mass_transfer_beta
         s% xtra(8) = log10(mdotSTABLE)
         s% xtra(9) = log10(mdotRG)




         ! =================================================
         !    SAVING PROFILES AT KEY EVOLUTIONARY STAGES
         ! =================================================

         if ( ( s% he_core_mass .gt. 1d-2 ) .and. &
              ( s% x_logical_ctrl(1).eqv..false. ) ) then 
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "     BEGIN: helium core"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileBEGHe'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(1) = .true.
          end if  

           
         if ( ( s% center_h1 .lt. 1d-4 ) .and. &
              ( s% x_logical_ctrl(2).eqv..false. ) ) then 
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "END: main sequence life"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileENDMS'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(2) = .true.
          end if  


         if ( ( mdotWD .gt. 1d-15 ) .and. ( s% xtra(13) .lt. 1d-15 ) .and. &
              ( s% x_logical_ctrl(16).eqv..false. ) ) then
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "BEGIN: roche lobe overflow"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileBEGMT'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(16) = .true.
         end if


         if ( ( mdotWD .lt. 1d-15 ) .and. ( s% xtra(13) .gt. 1d-15 ) .and. &
              ( s% x_logical_ctrl(14).eqv..false. ) ) then
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "BEGIN: detached phase"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileBEGDP'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(14) = .true.
         end if

         if ( ( mdotWD .gt. 1d-15 ) .and. ( s% xtra(13) .lt. 1d-15 ) .and. &
              ( s% x_logical_ctrl(14).eqv..true. ) .and. &
              ( s% x_logical_ctrl(15).eqv..false. ) ) then
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "END: detached phase"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileENDDP'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(15) = .true.
         end if

         
         if ( ( (s% star_mass - s% he_core_mass) .lt. 1d-12) .and. &
              ( s% x_logical_ctrl(3).eqv..false. ) ) then
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "END: hydrogen envelope is gone"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileENDENV'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(3) = .true.
         end if


         ! --------------
         !    MASSES
         ! --------------

         if ( (s% xtra(12).gt.1.00) .and. ((b% m(b% d_i)/Msun).lt.1.00) .and. &
              ( s% x_logical_ctrl(17).eqv..false. ) ) then
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "DONOR MASS: 1.00 Msun"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileM100'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(17) = .true.
         end if

         if ( (s% xtra(12).gt.0.90) .and. ((b% m(b% d_i)/Msun).lt.0.90) .and. &
              ( s% x_logical_ctrl(18).eqv..false. ) ) then
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "DONOR MASS: 0.90 Msun"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileM090'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(18) = .true.
         end if

         if ( (s% xtra(12).gt.0.80) .and. ((b% m(b% d_i)/Msun).lt.0.80) .and. &
              ( s% x_logical_ctrl(19).eqv..false. ) ) then
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "DONOR MASS: 0.80 Msun"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileM080'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(19) = .true.
         end if

         if ( (s% xtra(12).gt.0.70) .and. ((b% m(b% d_i)/Msun).lt.0.70) .and. &
              ( s% x_logical_ctrl(20).eqv..false. ) ) then
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "DONOR MASS: 0.70 Msun"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileM070'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(20) = .true.
         end if

         if ( (s% xtra(12).gt.0.60) .and. ((b% m(b% d_i)/Msun).lt.0.60) .and. &
              ( s% x_logical_ctrl(21).eqv..false. ) ) then
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "DONOR MASS: 0.60 Msun"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileM060'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(21) = .true.
         end if

         if ( (s% xtra(12).gt.0.50) .and. ((b% m(b% d_i)/Msun).lt.0.50) .and. &
              ( s% x_logical_ctrl(22).eqv..false. ) ) then
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "DONOR MASS: 0.50 Msun"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileM050'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(22) = .true.
         end if

         if ( (s% xtra(12).gt.0.40) .and. ((b% m(b% d_i)/Msun).lt.0.40) .and. &
              ( s% x_logical_ctrl(23).eqv..false. ) ) then
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "DONOR MASS: 0.40 Msun"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileM040'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(23) = .true.
         end if

         if ( (s% xtra(12).gt.0.30) .and. ((b% m(b% d_i)/Msun).lt.0.30) .and. &
              ( s% x_logical_ctrl(24).eqv..false. ) ) then
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "DONOR MASS: 0.30 Msun"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileM030'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(24) = .true.
         end if

         if ( (s% xtra(12).gt.0.20) .and. ((b% m(b% d_i)/Msun).lt.0.20) .and. &
              ( s% x_logical_ctrl(25).eqv..false. ) ) then
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "DONOR MASS: 0.20 Msun"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileM020'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(25) = .true.
         end if

         if ( (s% xtra(12).gt.0.15) .and. ((b% m(b% d_i)/Msun).lt.0.15) .and. &
              ( s% x_logical_ctrl(26).eqv..false. ) ) then
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "DONOR MASS: 0.15 Msun"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileM015'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(26) = .true.
         end if






         if ( (s% xtra(12).gt.0.10) .and. ((b% m(b% d_i)/Msun).lt.0.10) .and. &
              ( s% x_logical_ctrl(4).eqv..false. ) ) then
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "DONOR MASS: 0.10 Msun"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileM010'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(4) = .true.
         end if
         
         if ( (s% xtra(12).gt.0.09) .and. ((b% m(b% d_i)/Msun).lt.0.09) .and. &
              ( s% x_logical_ctrl(5).eqv..false. ) ) then
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "DONOR MASS: 0.09 Msun"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileM009'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(5) = .true.
         end if

         if ( (s% xtra(12).gt.0.08) .and. ((b% m(b% d_i)/Msun).lt.0.08) .and. &
              ( s% x_logical_ctrl(6).eqv..false. ) ) then
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "DONOR MASS: 0.08 Msun"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileM008'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(6) = .true.
         end if

         if ( (s% xtra(12).gt.0.07) .and. ((b% m(b% d_i)/Msun).lt.0.07) .and. &
              ( s% x_logical_ctrl(7).eqv..false. ) ) then
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "DONOR MASS: 0.07 Msun"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileM007'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(7) = .true.
         end if

         if ( (s% xtra(12).gt.0.06) .and. ((b% m(b% d_i)/Msun).lt.0.06) .and. &
              ( s% x_logical_ctrl(8).eqv..false. ) ) then
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "DONOR MASS: 0.06 Msun"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileM006'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(8) = .true.
         end if

         if ( (s% xtra(12).gt.0.05) .and. ((b% m(b% d_i)/Msun).lt.0.05) .and. &
              ( s% x_logical_ctrl(9).eqv..false. ) ) then
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "DONOR MASS: 0.05 Msun"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileM005'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(9) = .true.
         end if

         if ( (s% xtra(12).gt.0.04) .and. ((b% m(b% d_i)/Msun).lt.0.04) .and. &
              ( s% x_logical_ctrl(10).eqv..false. ) ) then
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "DONOR MASS: 0.04 Msun"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileM004'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(10) = .true.
         end if

         if ( (s% xtra(12).gt.0.03) .and. ((b% m(b% d_i)/Msun).lt.0.03) .and. &
              ( s% x_logical_ctrl(11).eqv..false. ) ) then
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "DONOR MASS: 0.03 Msun"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileM003'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(11) = .true.
         end if


         if ( (s% xtra(12).gt.0.02) .and. ((b% m(b% d_i)/Msun).lt.0.02) .and. &
              ( s% x_logical_ctrl(12).eqv..false. ) ) then
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "DONOR MASS: 0.02 Msun"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileM002'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(12) = .true.
         end if


         if ( (s% xtra(12).gt.0.01) .and. ((b% m(b% d_i)/Msun).lt.0.01) .and. &
              ( s% x_logical_ctrl(13).eqv..false. ) ) then
            write(*,*) ""
            write(*,*) ""
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) "DONOR MASS: 0.01 Msun"
            write(*,*) "*************************************"
            write(*,*) "*************************************"
            write(*,*) ""
            write(*,*) ""
            s% profile_data_prefix = 'profileM001'
            s% need_to_update_history_now = .true.
            s% need_to_save_profiles_now = .true.
            s% x_logical_ctrl(13) = .true.
         end if



         ! Saving star mass ( in Msun )
         s% xtra(12) = b% m(b% d_i)/Msun 

         ! Saving mass transfer rate ( in Msun/yr )
         s% xtra(13) = (abs(b% mtransfer_rate)/Msun)*secyer


         check = 0

         ! Checking if donor has some helium in its core at the onset of mass transfer
         if ( (b% rl_relative_gap(b% d_i) + 1.d0) > 0.99 ) then
          if ( (s% star_mass - s% he_core_mass) .gt. 0.5 ) then
           if ( s% he_core_mass < 0.01 ) then
             check = 1
             write(*,*) "stopping: Roche lobe factor (",(b% rl_relative_gap(b% d_i) + 1.d0), & 
                        ") is greater than 0.99 and there is no helium in the core"
           end if
          end if
         end if


         if ( abs(b% mtransfer_rate)/Msun*secyer > 1d-4 ) then
             check = 1
             write(*,*) "stopping: mass transfer rate (",abs(b% mtransfer_rate)/Msun*secyer, & 
                        ") is greater than 1d-4"
         end if


         if (check.eq.1) then
             extras_binary_finish_step = terminate
         else
             extras_binary_finish_step = keep_going
         end if 

         
      end function extras_binary_finish_step
      
      subroutine extras_binary_after_evolve(binary_id, ierr)
         type (binary_info), pointer :: b
         integer, intent(in) :: binary_id
         integer, intent(out) :: ierr
         call binary_ptr(binary_id, b, ierr)
         if (ierr /= 0) then ! failure in  binary_ptr
            return
         end if      
         
 

         call test_suite_after_evolve(b, ierr)

      end subroutine extras_binary_after_evolve     


      subroutine check_jdot_mb_conditions(b, s, apply_jdot_mb, qconv_env)
         type (binary_info), pointer :: b
         type (star_info), pointer :: s
         logical, intent(out) :: apply_jdot_mb
         real(dp), intent(out) :: qconv_env
         
         real(dp) :: qrad_core
         integer :: i, k, id

         include 'formats.inc'

         ! calculate how much of inner region is convective
         qrad_core = 0d0

         do k = s% nz, 1, -1
            if (s% q(k) > b% jdot_mb_qlim_for_check_rad_core .and. & 
               (qrad_core == 0d0 .or. s% mixing_type(k) /= convective_mixing)) exit
            if (s% mixing_type(k) == convective_mixing) &
               qrad_core = qrad_core + s% dq(k)
         end do

         ! calculate how much of the envelope is convective
         qconv_env = 0d0

         do k = 1, s% nz
            if (s% q(k) < b% jdot_mb_qlim_for_check_conv_env .and. &
               (qconv_env == 0d0 .or. s% mixing_type(k) /= convective_mixing)) exit
            if (s% mixing_type(k) == convective_mixing) &
               qconv_env = qconv_env + s% dq(k)
         end do

         apply_jdot_mb = .true.

         if (qconv_env < b% jdot_mb_min_qconv_env) then
            apply_jdot_mb = .false.
            return
         end if

         if (qconv_env > b% jdot_mb_max_qconv_env) then
            apply_jdot_mb = .false.
            return
         end if

         if (qrad_core > b% jdot_mb_max_qrad_core) then
            apply_jdot_mb = .false.
            return
         end if
         
      end subroutine check_jdot_mb_conditions





      subroutine jdot_mb_routine(binary_id, ierr)
         integer, intent(in) :: binary_id
         integer, intent(out) :: ierr
         type (binary_info), pointer :: b
         logical :: apply_mdot_mb
         logical :: apply_jdot_mb
         real(dp) :: rsun4,two_pi_div_p3, qconv_env, jdot_scale
! ===========================
! NEW VARIABLES
! ===========================
         type (star_info), pointer :: s
         integer :: k, nz
         real(dp) :: bs4, clight5, cgrav3
         real(dp) :: turnover_time, tt_temp, tt_temp_scaled, tt_old, tt_diff
         real(dp) :: vel, vel_ratio, vel_diff, upper_lim, lower_lim, scaled_vel
         real(dp) :: eps_nuc_lim, eps_nuc
         real(dp) :: dr, tau_lim, delta_mag_chk
         real(dp) :: two_pi_div_p2, K2
         real(dp) :: tt_ratio, tt4
         real(dp) :: rot_ratio, rot4
         real(dp) :: rad4
         real(dp) :: v_esc2, v_mod2
         real(dp) :: alfven_no_R, R_alfven
         real(dp) :: jdot_mb_old, jdot_mb_new
         real(dp) :: conv_env_r, conv_env_m, sonic_cross_time, mag_field
         logical :: conv_env_found
! ===========================
         ierr = 0
         call binary_ptr(binary_id, b, ierr)
         if (ierr /= 0) then
            write(*,*) 'failed in binary_ptr'
            return
         end if


! ::::::::::::::::::::::::::::::::::::::::
! ::::::::::::::::::::::::::::::::::::::::
!
!         This routine was adapted from 
!       Van & Ivanova (2019, ApJ, 886, L31)
!    
!       https://zenodo.org/record/3647683
!
! ::::::::::::::::::::::::::::::::::::::::
! ::::::::::::::::::::::::::::::::::::::::



! ===========================
! INITIALIZE VARIABLES AND SET CONSTANTS
! ===========================

         s => b% s_donor
         nz = s% nz
         vel_ratio = s% x_ctrl(1)
         tau_lim = s% x_ctrl(2)
         tt_old = s% x_ctrl(3)

         conv_env_found = .false.

         sonic_cross_time = 0.d0
         mag_field = 0.d0
         conv_env_r = 0.d0
         conv_env_m = 0.d0
         turnover_time = 0.d0
         tt_temp = 0.d0
         tt_temp_scaled = 0.d0

         eps_nuc_lim = 1.0d-2
         vel_diff = 0.d0
         scaled_vel = 0.d0

         b% jdot_mb = 0.d0
         rsun4 = pow4(rsun)

         two_pi_div_p3 = (2.0d0*pi/b% period)*(2.0d0*pi/b% period)*(2.0d0*pi/b% period)
         two_pi_div_p2 = (2.0*pi/b% period)*(2.0*pi/b% period)

         K2 = 0.07 * 0.07

! ===========================
! INITIAL TURNOVER TIME CALCULATION
! ===========================

         do k = nz, 1, -1 ! beginning of do loop to calculate convective turnover time

                eps_nuc = s% eps_nuc(k)
                ! check if the cell we are looping through satisfies our convection criteria
                if ((s% gradr(k) .gt. s% grada(k)) .and. (eps_nuc .lt. eps_nuc_lim)) then
                    ! toggle the boolean to begin integration
                    conv_env_found = .true.
                end if

                ! only enter this portion if the convective boolean is true
                ! this loop will go from the innermost cell that is convective to 
                ! the surface. This is to try and smooth through any numeric issues
                ! with convective zones appearing and disappearing in MESA.
                if (conv_env_found) then

                    ! loop to calculate the size of the cell, the innermost cell
                    ! needs special consideration as it is above the core
                    if (k .lt. s% nz) then
                        dr = (s% r(k) - s% r(k + 1))
                    else
                        dr = (s% r(k) - s% R_center)
                    end if
                    
                    ! determine the convective velocity inside each given cell
                    if (s% mixing_type(k) == convective_mixing) then

                        ! need to ensure that the convective velocity is within
                        ! our defined limits, if they are outside of these limits
                        ! set them to be the max/min value allowed.
                        vel = s% conv_vel(k)
                        lower_lim = vel_ratio * s% csound(k)
                        upper_lim = 1.0 * s% csound(k)

                        if (vel .lt. lower_lim) then
                            vel = lower_lim
                        else if (vel .gt. upper_lim) then
                            vel = upper_lim
                        end if

                        ! Calculating convective envelope mass and radius
                        conv_env_r = conv_env_r + dr
                        conv_env_m = conv_env_m + s% dm(k)

                    
                    ! if the cell isnt defined by MESA to be convective take the
                    ! convective velocity to be equal to sound speed
                    else
                        vel = s% csound(k)
                    end if

                    ! Final check involving the opacity of the given cell. If the 
                    ! cell isn't near the surface (low tau) then include it in our integration
                    if (s% tau(k) .gt. tau_lim) then
                        sonic_cross_time = sonic_cross_time + (dr / s% csound(k))
                        tt_temp = tt_temp + (dr / vel)
                    end if
                end if

         end do ! end of do loop to calculate convective turnover time


         conv_env_m = conv_env_m/Msun
         conv_env_r = conv_env_r/Rsun

         ! reset the boolean just in case
         conv_env_found = .false.


! ===========================
! TURNOVER TIME CHECK, THIS IS TO TRY AND AVOID LARGE CHANGES
! ===========================

         ! simply set the turnover time to the internal variable calculated above
         turnover_time = tt_temp

         ! Empirical Relation
         !turnover_time = 1.16 - 1.49*(log10(b% m(b% d_i)/Msun))-&
         !             0.54*(log10(b% m(b% d_i)/Msun))**2
         !turnover_time = (10**(turnover_time))*(24.0*60.0*60.0)

         if (s% model_number .gt. 1) then
                ! calculate the variables used to check if our system is rapidly evolving
                tt_diff = abs(tt_old - tt_temp) / tt_old
                delta_mag_chk = s% dt / tt_old

                ! write (*,*) "tt_diff = ", tt_diff
                ! write (*,*) "delta_mag = ", delta_mag_chk
                ! write (*,*) "turnover_time = ", turnover_time
                ! write (*,*) "tt_old = ", tt_old

                ! check if timesteps are very small or if the relative change is very large
                if (tt_diff .gt. delta_mag_chk) then 
                    write (*,*) "large change, adjusting accordingly"
                    turnover_time = tt_old + (tt_temp - tt_old) * min((s% dt / tt_old), 0.5)
                    mag_field = (turnover_time / 2.8d6) * (2073600. / b% period)
 
                end if ! end of timestep/relative change check
         end if

         ! remember the current values to be used as comparison in the next step
         tt_old = turnover_time
         s% x_ctrl(3) = tt_old


! ===========================
! PRE-REQUISITIES FOR MAGNETIC BRAKING CALCULATION
! ===========================

         tt_ratio = turnover_time / 2.8d6
         tt4 = pow4(tt_ratio)

         rot_ratio = (2073600. / b% period )
         rot4 = pow4(rot_ratio)

         rad4 = pow4(b% r(b% d_i))

         v_esc2 = 2.0 * standard_cgrav * b% m(b% d_i) / b% r(b% d_i)
         v_mod2 = v_esc2 + 2.0 * two_pi_div_p2 * b% r(b% d_i) * b% r(b% d_i) / K2 

         alfven_no_R = rad4 * rot4 * tt4 / (b% mdot_system_wind(b% d_i) * &
                       b% mdot_system_wind(b% d_i)) * (1.0 / v_mod2)

         R_alfven = b% r(b% d_i) * alfven_no_R**(1.d0/3.d0)


! ===========================
! MAGNETIC BRAKING CALCULATION
! ===========================

         call check_jdot_mb_conditions(b, b% s_donor, apply_jdot_mb, qconv_env)

         if (apply_jdot_mb .or. b% keep_mb_on) then
            jdot_scale = 1d0
            if (b% jdot_mb_scale_for_low_qconv_env) then
               !scale jdot for tiny convective envelope, from Podsiadlowski et al. (2002)
               !The Astrophysical Journal, Volume 565, Issue 2, pp. 1107-1133
               if (qconv_env > b% jdot_mb_mass_frac_for_scale) then
                  jdot_scale = 1d0
               else
                  jdot_scale = exp(-b% jdot_mb_mass_frac_for_scale/max(1d-99,qconv_env)+1)
               end if
            end if


            ! ------------  RVJ ----------------
            !  Rappaport, Verbunt, and Joss, 1983, APJ, 275, 713
            if(b% magnetic_braking_gamma .ge. 0)then
              b% jdot_mb = -3.8d-30*b% m(b% d_i)*rsun4* &         
                           pow(min(b% r(b% d_i),b% rl(b% d_i))/rsun,b% magnetic_braking_gamma)* &
                           two_pi_div_p3*jdot_scale
            end if


            ! ------------  CARB ----------------
            !  Van and Ivanova, 2019, ApJ, 886, L31
            if(b% magnetic_braking_gamma .lt. 0)then
              b% jdot_mb = (2.0/3.0) * (2.0*pi/b% period) * b% mdot_system_wind(b% d_i) * &
                           R_alfven * R_alfven * jdot_scale
            end if



            b% using_jdot_mb(b% d_i) = .true.

            if ((apply_jdot_mb .or. b% keep_mb_on) .and. .not. b% using_jdot_mb_old(b% d_i)) then
               write(*,*) 'turn on magnetic braking for star ', b% d_i
            end if

         else if (.not. (apply_jdot_mb .or. b% keep_mb_on) .and. b% using_jdot_mb_old(b% d_i)) then
            ! required mdot for the implicit scheme may drop drastically,
            ! so its neccesary to increase change factor to avoid implicit 
            ! scheme from getting stuck
            b% change_factor = b% max_change_factor
            b% using_jdot_mb(b% d_i) = .false.
            write(*,*) 'turn off magnetic braking for star ', b% d_i
         end if

         if (b% point_mass_i == 0 .and. b% include_accretor_mb) then
            call check_jdot_mb_conditions(b, b% s_accretor, apply_jdot_mb, qconv_env)
            if (apply_jdot_mb .or. b% keep_mb_on) then
               jdot_scale = 1d0
               if (b% jdot_mb_scale_for_low_qconv_env) then
                  !scale jdot for tiny convective envelope, from Podsiadlowski et al. (2002)
                  !The Astrophysical Journal, Volume 565, Issue 2, pp. 1107-1133
                  if (qconv_env > b% jdot_mb_mass_frac_for_scale) then
                     jdot_scale = 1d0
                  else
                     jdot_scale = exp(-b% jdot_mb_mass_frac_for_scale/max(1d-99,qconv_env)+1)
                  end if
               end if
               

               ! ------------  RVJ ----------------
               !  Rappaport, Verbunt, and Joss, 1983, APJ, 275, 713
               if(b% magnetic_braking_gamma .ge. 0)then
                 b% jdot_mb = b% jdot_mb - &
                           3.8d-30*b% m(b% a_i)*rsun4* &
                           pow(min(b% r(b% a_i),b% rl(b% a_i))/rsun,b% magnetic_braking_gamma)* &
                           two_pi_div_p3*jdot_scale
               end if


               ! ------------  CARB ----------------
               !  Van and Ivanova, 2019, ApJ, 886, L31
               if(b% magnetic_braking_gamma .lt. 0)then
                 b% jdot_mb = b% jdot_mb + &
                              (2.0/3.0) * (2.0*pi/b% period) * b% mdot_system_wind(b% d_i) * &
                           R_alfven * R_alfven * jdot_scale
               end if


               b% using_jdot_mb(b% a_i) = .true.
               if ((apply_jdot_mb .or. b% keep_mb_on) .and. .not. b% using_jdot_mb_old(b% a_i)) then
                  write(*,*) 'turn on magnetic braking for star ', b% a_i
               end if
            else if (.not. (apply_jdot_mb .or. b% keep_mb_on) .and. b% using_jdot_mb_old(b% a_i)) then
               ! required mdot for the implicit scheme may drop drastically,
               ! so its neccesary to increase change factor to avoid implicit 
               ! scheme from getting stuck
               b% change_factor = b% max_change_factor
               b% using_jdot_mb(b% a_i) = .false.
               write(*,*) 'turn off magnetic braking for star ', b% a_i
            end if
         else if (b% model_twins_flag .and. b% include_accretor_mb) then
            b% jdot_mb = 2*b% jdot_mb
            b% using_jdot_mb(b% a_i) = .true.
         end if


         s% xtra(1) = turnover_time
         s% xtra(2) = mag_field
         s% xtra(3) = conv_env_r
         s% xtra(4) = conv_env_m
         s% xtra(5) = sonic_cross_time
         s% xtra(6) = 1.16 - 1.49*(log10(b% m(b% d_i)/Msun))-&
                      0.54*(log10(b% m(b% d_i)/Msun))**2
         s% xtra(6) = (10**(s% xtra(6)))*(24.0*60.0*60.0)

         s% xtra(10) = qconv_env
         s% xtra(11) = jdot_scale


      end subroutine jdot_mb_routine

      
      end module run_binary_extras
