! ***********************************************************************
!
!   Copyright (C) 2010  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_star_extras

      use star_lib
      use star_def
      use const_def
      use crlibm_lib
      use utils_lib
      use chem_def !to be able to call isotopes in a network independent fashion 
                   !only assuming the isotopes used ARE in the network           
      
      implicit none
      
      integer :: time0, time1, clock_rate, redo_count
      
      contains
      
      
      subroutine extras_controls(id, ierr)
         integer, intent(in) :: id
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         
         s% extras_startup => extras_startup
         s% extras_start_step => extras_start_step
         s% extras_check_model => extras_check_model
         s% extras_finish_step => extras_finish_step
         s% extras_after_evolve => extras_after_evolve
         s% how_many_extra_history_columns => how_many_extra_history_columns
         s% data_for_extra_history_columns => data_for_extra_history_columns
         s% how_many_extra_profile_columns => how_many_extra_profile_columns
         s% data_for_extra_profile_columns => data_for_extra_profile_columns  
         ! this is optional
         ! s% other_wind => brott_wind !!Mathieu commented
         ! remove eps_grav from the very surface, improves a bit stability
         !s% use_other_eps_grav = .true.
         !s% other_eps_grav => my_eps_grav

         !decorator for pgplot
         s% Profile_Panels1_pgstar_decorator => gamma4_3_decorator
         s% Profile_Panels2_pgstar_decorator => vesc_decorator
      end subroutine extras_controls

      
      subroutine vesc_decorator(id, xmin, xmax, ymin, ymax, plot_num, ierr)
         integer, intent(in) :: id
         !Not doubles
         real,intent(in) :: xmin, xmax, ymin, ymax 
         integer, intent(in) :: plot_num
         integer, intent(out) :: ierr
         real :: v_esc
         integer :: k
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
        
         if ((plot_num==1) .and. (s% lxtra1)) then
            do k=1, s% nz
               v_esc = s% x_ctrl(1)*sqrt(2.d0*s% cgrav(k)*s% m(k)/(s% r(k)))
               if (s% xh(s% i_u, k) >= v_esc) then
                  call pgslw(10)
                  call pgsci(clr_IndianRed)
                  call pgline(2,(/xmin,xmax/),(/v_esc,v_esc/))
                  call pgline(2,(/s%m(k),s% m(k)/),(/ymin,ymax/))
                  exit 
               end if
            end do
         end if

      end subroutine vesc_decorator


      subroutine gamma4_3_decorator(id, xmin, xmax, ymin, ymax, plot_num, ierr)
         integer, intent(in) :: id
         !Not doubles
         real,intent(in) :: xmin, xmax, ymin, ymax 
         integer, intent(in) :: plot_num
         integer, intent(out) :: ierr
         real :: gamma1_crit
         ierr = 0

         gamma1_crit = 4.0/3.0
         if (plot_num==2) then
            if ((ymin <= gamma1_crit) .and. (ymax >= gamma1_crit)) then
               call pgslw(10)
               call pgsci(clr_IndianRed)
               call pgline(2,(/xmin,xmax/),(/gamma1_crit,gamma1_crit/))
            end if
         end if

      end subroutine gamma4_3_decorator

      !! already commented by Pablo
      ! subroutine my_eps_grav(id, k, dt, ierr)
      !    integer, intent(in) :: id, k
      !    real(dp), intent(in) :: dt
      !    integer, intent(out) :: ierr
      !    type (star_info), pointer :: s
      !    ierr = 0
      !    call star_ptr(id, s, ierr)
      !    if (ierr /= 0) return
      !    if (s% q(k) < 0.9999) then
      !       return
      !    end if

      !    s% eps_grav(k) = 0d0
      !    s% d_eps_grav_dlndm1(k) = 0d0
      !    s% d_eps_grav_dlnd00(k) = 0d0
      !    s% d_eps_grav_dlndp1(k) = 0d0
      !    s% d_eps_grav_dlnTm1(k) = 0d0
      !    s% d_eps_grav_dlnT00(k) = 0d0
      !    s% d_eps_grav_dlnTp1(k) = 0d0
      !    s% d_eps_grav_dlnPgasm1_const_T(k) = 0d0
      !    s% d_eps_grav_dlnPgas00_const_T(k) = 0d0
      !    s% d_eps_grav_dlnPgasp1_const_T(k) = 0d0
      !    s% d_eps_grav_dlnTm1_const_Pgas(k) = 0d0
      !    s% d_eps_grav_dlnT00_const_Pgas(k) = 0d0
      !    s% d_eps_grav_dlnTp1_const_Pgas(k) = 0d0
      !    s% d_eps_grav_dlnR00(k) = 0d0
      !    s% d_eps_grav_dlnRp1(k) = 0d0
      !    s% d_eps_grav_dL00(k) = 0d0
      !    s% d_eps_grav_dLp1(k) = 0d0
      !    s% d_eps_grav_dv00(k) = 0d0
      !    s% d_eps_grav_dvp1(k) = 0d0
      !    s% eps_grav_dm_term_const_q(k) = 0d0
      !    s% eps_grav_dt_term_const_q(k) = 0d0
         
      ! end subroutine my_eps_grav


      !! commented by Mathieu
      ! subroutine brott_wind(id, Lsurf, Msurf, Rsurf, Tsurf, w, ierr)
      !    use star_def
      !    integer, intent(in) :: id
      !    real(dp), intent(in) :: Lsurf, Msurf, Rsurf, Tsurf ! surface values (cgs)
      !    ! NOTE: surface is outermost cell. not necessarily at photosphere.
      !    ! NOTE: don't assume that vars are set at this point.
      !    ! so if you want values other than those given as args,
      !    ! you should use values from s% xh(:,:) and s% xa(:,:) only.
      !    ! rather than things like s% Teff or s% lnT(:) which have not been set yet.
      !    real(dp), intent(out) :: w ! wind in units of Msun/year (value is >= 0)
      !    integer, intent(out) :: ierr

      !    integer :: h1, he4
      !    real(dp) :: Xs, Ys, Z_div_Z_solar, Teff_jump, alfa, L1, M1, R1, T1, &
      !       vink_wind, nieu_wind, hamann_wind
      !    type (star_info), pointer :: s
      !    ierr = 0
      !    call star_ptr(id, s, ierr)
      !    if (ierr /= 0) return

      !    L1 = Lsurf
      !    M1 = Msurf
      !    R1 = Rsurf
      !    T1 = Tsurf

      !    h1 = s% net_iso(ih1)
      !    he4 = s% net_iso(ihe4)
      !    Xs = s% xa(h1,1)
      !    Ys = s% xa(he4,1)
      !    ! Z=0.017 is Z from Grevesse et al. 1996
      !    Z_div_Z_solar = s% Zbase/0.02d0
      !    ! use Vink et al 2001, eqns 14 and 15 to set "jump" temperature
      !    Teff_jump = 1d3*(61.2d0 + 2.59d0*(-13.636d0 + 0.889d0*log10_cr(Z_div_Z_solar)))

      !    vink_wind = 0d0
      !    nieu_wind = 0d0
      !    hamann_wind = 0d0
      !    w = 0

      !    call eval_Vink_wind(vink_wind)
      !    call eval_Nieuwenhuijzen_wind(nieu_wind)
      !    call eval_Hamann_wind(hamann_wind)

      !    ! use 1/10 hamann
      !    hamann_wind = hamann_wind/10d0

      !    if (T1 < Teff_jump) then
      !       ! low T wind
      !       w = max(vink_wind, nieu_wind)
      !    else
      !       ! high T wind
      !       alfa = 0d0
      !       if (Xs > 0.7d0) then
      !          alfa = 1d0
      !       else if (Xs > 0.4d0 .and. Xs < 0.7d0) then
      !          alfa = (Xs - 0.4d0)/0.3d0
      !       end if
      !       w = alfa * vink_wind + (1d0-alfa) * hamann_wind
      !    end if

      !    ierr = 0

      !    contains

      !    subroutine eval_Vink_wind(w)
      !       real(dp), intent(inout) :: w
      !       real(dp) :: alfa, w1, w2, logMdot, dT, vinf_div_vesc

      !       ! alfa = 1 for hot side, = 0 for cool side
      !       if (T1 > 27500d0) then
      !          alfa = 1
      !       else if (T1 < 22500d0) then
      !          alfa = 0
      !       else
      !          dT = 100d0
      !          if (T1 > Teff_jump + dT) then
      !             alfa = 1
      !          else if (T1 < Teff_jump - dT) then
      !             alfa = 0
      !          else
      !             alfa = (T1 - (Teff_jump - dT)) / (2*dT)
      !          end if
      !       end if
            
      !       if (alfa > 0) then ! eval hot side wind (eqn 24)
      !          vinf_div_vesc = 2.6d0 ! this is the hot side galactic value
      !          vinf_div_vesc = vinf_div_vesc*pow_cr(Z_div_Z_solar,0.13d0) ! corrected for Z
      !          logMdot = &
      !             - 6.697d0 &
      !             + 2.194d0*log10_cr(L1/Lsun/1d5) &
      !             - 1.313d0*log10_cr(M1/Msun/30) &
      !             - 1.226d0*log10_cr(vinf_div_vesc/2d0) &
      !             + 0.933d0*log10_cr(T1/4d4) &
      !             - 10.92d0*pow2(log10_cr(T1/4d4)) &
      !             + 0.85d0*log10_cr(Z_div_Z_solar)
      !          w1 = exp10_cr(logMdot)
      !       else
      !          w1 = 0
      !       end if
            
      !       if (alfa < 1) then ! eval cool side wind (eqn 25)
      !          vinf_div_vesc = 1.3d0 ! this is the cool side galactic value
      !          vinf_div_vesc = vinf_div_vesc*pow_cr(Z_div_Z_solar,0.13d0) ! corrected for Z
      !          logMdot = &
      !             - 6.688d0 &
      !             + 2.210d0*log10_cr(L1/Lsun/1d5) &
      !             - 1.339d0*log10_cr(M1/Msun/30) &
      !             - 1.601d0*log10_cr(vinf_div_vesc/2d0) &
      !             + 1.07d0*log10_cr(T1/2d4) &
      !             + 0.85d0*log10_cr(Z_div_Z_solar)
      !          w2 = exp10_cr(logMdot)
      !       else
      !          w2 = 0
      !       end if
            
      !       w = alfa*w1 + (1 - alfa)*w2
            
      !    end subroutine eval_Vink_wind

      !    subroutine eval_Nieuwenhuijzen_wind(w)
      !       ! Nieuwenhuijzen, H.; de Jager, C. 1990, A&A, 231, 134 (eqn 2)
      !       real(dp), intent(out) :: w
      !       real(dp) :: log10w
      !       include 'formats'
      !       log10w = -14.02d0 &
      !                +1.24d0*log10_cr(L1/Lsun) &
      !                +0.16d0*log10_cr(M1/Msun) &
      !                +0.81d0*log10_cr(R1/Rsun) &
      !                +0.85d0*log10_cr(Z_div_Z_solar)
      !       w = exp10_cr(log10w)
      !    end subroutine eval_Nieuwenhuijzen_wind

      !    subroutine eval_Hamann_wind(w)
      !       ! Hamann, W.-R.; Koesterke, L.; Wessolowski, U. 1995, A&A, 299, 151
      !       real(dp), intent(out) :: w
      !       real(dp) :: log10w
      !       include 'formats'
      !       log10w = -11.95d0 &
      !                +1.5d0*log10_cr(L1/Lsun) &
      !                -2.85d0*Xs &
      !                + 0.85d0*log10_cr(Z_div_Z_solar)
      !       w = exp10_cr(log10w)
      !    end subroutine eval_Hamann_wind

      ! end subroutine brott_wind

      integer function extras_startup(id, restart, ierr)
         integer, intent(in) :: id
         logical, intent(in) :: restart
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         include 'formats'
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         extras_startup = 0
         call system_clock(time0,clock_rate)
         if (.not. restart) then
            call alloc_extra_info(s)
            ! lxtra1 is true while hydro is on
            s% lxtra1 = .false.
            ! lxtra2 is true after hydro is turned on for the first time
            s% lxtra2 = .false.
            ! lxtra3 is true after He depletion
            s% lxtra3 = .false.
            ! xtra 1 contains the time at which a pulse starts
            s% xtra1 = -1d0
            ! xtra2 contains the estimate for the dynamical time when that happens
            s% xtra2 = -1d0
            ! xtra29 stores surface pressure from previous step
            s% xtra29 = -1d0
            !! to number pulse output files
            s% ixtra1 = 0 
            ! ixtra2 counts the steps during which the conditions to relax a model
            ! are met
            s% ixtra2 = 0
         else ! it is a restart
            call unpack_extra_info(s)
            if (s% lxtra1) then
               call star_read_controls(id, 'inlist_hydro_on', ierr)
	       if(ierr/=0) call mesa_error(__FILE__,__LINE__)

               call star_set_v_flag(id, .false., ierr)
	       if(ierr/=0) call mesa_error(__FILE__,__LINE__)
               call star_set_u_flag(id, .true., ierr)
	       if(ierr/=0) call mesa_error(__FILE__,__LINE__)
               call star_set_RTI_flag(id, .true., ierr)
	       if(ierr/=0) call mesa_error(__FILE__,__LINE__)
            else if (s% lxtra2) then
               call star_read_controls(id, 'inlist_hydro_off', ierr)
	       if(ierr/=0) call mesa_error(__FILE__,__LINE__)
            end if
            if (s% lxtra2) then
               call star_read_controls(id, 'inlist_after_first_pulse', ierr)
	       if(ierr/=0) call mesa_error(__FILE__,__LINE__)
            end if
         end if
      end function extras_startup
      
      
      subroutine extras_after_evolve(id, id_extra, ierr)
         integer, intent(in) :: id, id_extra
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         double precision :: dt
         character (len=strlen) :: test
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         call system_clock(time1,clock_rate)
         dt = dble(time1 - time0) / clock_rate / 60
         write(*,'(/,a50,f12.2,99i10/)') 'runtime (minutes), retries, backups, steps', &
            dt, s% num_retries, s% num_backups, s% model_number
         ierr = 0
      end subroutine extras_after_evolve
      

      ! returns either keep_going, retry, backup, or terminate.
      integer function extras_check_model(id, id_extra)
         integer, intent(in) :: id, id_extra
         integer :: ierr, k
         real(dp) :: conv_fraction
         type (star_info), pointer :: s
         include 'formats'
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         extras_check_model = keep_going         

         ! write(*,*) "Change in surface pressure", log10_cr(abs(s% P(1)/s% xtra29)), log10_cr(abs(s% P(1)))
         if(s% xtra29 > 0d0 .and. log10_cr(abs(s% P(1)/s% xtra29)) > 0.2) then
            extras_check_model = retry
            write(*,*) "retry due to large change in surface pressure"
            return
         end if

      end function extras_check_model



      !!! MATHIEU'S OUTPUT
      integer function how_many_extra_history_columns(id, id_extra)
         integer, intent(in) :: id, id_extra
         integer :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         how_many_extra_history_columns = 10
      end function how_many_extra_history_columns
      
      
      subroutine data_for_extra_history_columns(id, id_extra, n, names, vals, ierr)
         integer, intent(in) :: id, id_extra, n
         character (len=maxlen_history_column_name) :: names(n)
         real(dp) :: vals(n)
         integer, intent(out) :: ierr
         real(dp) :: norm, norm2, avg_P, dt, tmp
         integer :: k
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         call system_clock(time1,clock_rate)
         dt = dble(time1 - time0) / clock_rate / 60

         names(1) = "mass_lowgamma"
         vals(1) = 0 !initialization                                                                                                                         
         names(2) = "avg_delta_gamma" ! \int \Gamma_1(r)P(r) dV / \int P(r) dV N.B: dV = 4pir^2dr in spherical symmetry 
         vals(2) = 0
         names(3) = "sigma_ad" ! see Stothers 1999 MNRAS 205, 365-372
         vals(3) = 0
         names(4) = 'runtime_minutes'
         vals(4) = dt
         names(5) = "center_s"
         vals(5) = 0
         names(6) = "kh_timescale_maxL_Lneu"
         vals(6) = 0
         names(7) = "r_95percentile_mass"
         vals(7) = 0
         names(8) = "lum_95percentile_mass"
         vals(8) = 0

         norm2 = 0
         avg_P = 0
         norm = 0  ! normalization    
         
         names(9) = 'vel_flag'
         if(s%u_flag) then
            vals(9) = 1.0
         else
            vals(9) = 0.0
         end if

	names(10) = 'pulse_num'
	vals(10) = s% ixtra1
         
         tmp=0d0
         k = 1
         do while (k .lt. s%nz)
            vals(2) = vals(2) + (s% gamma1(k)-4.d0/3.d0)*s% P(k)*s% dm(k)/s% rho(k)
            norm = norm + s%P(k)*s%dm(k)/s%rho(k)
            norm2 = norm2 + (1.d0/3.0d0)*s%r(k)*s%r(k)*s%r(k)*s%r(k)*s%rho(k)
            avg_P = avg_P + s%P(k)*s%r(k)*s%r(k)
            if (s% gamma1(k) < 4d0/3d0) then
               vals(1) = vals(1)+ s%dm(k)/Msun
            end if
            tmp = tmp + (abs(s% eps_nuc_neu_total(k)) + abs(s% non_nuc_neu(k)))*(s% dq(k)) ! to calculate the total neutrino luminosity in cgs units
            if (abs(s%q(k)-0.95d0) < abs(s%q(k+1)-0.95d0)) then
               ! if the present cell is closer to the 95th percentile by mass than the next one
               vals(7) = s% r(k)
               vals(8) = s% L(k)
            end if
            k = k+1
         end do
         ! normalize avg_gamma                                                                                                                               
         vals(2) = vals(2)/norm
         ! evaluate the frequency of the adiabatic pulsations (star unstable if negative)                                                                    
         vals(3) = ((3d0*vals(6)-4d0)*avg_P)/norm2
         vals(5) = s%entropy(s%nz)
         vals(6)= (standard_cgrav*s%star_mass*Msun*s%star_mass*Msun)/(2d0*s%r(1)*max(s%L(1), tmp)) !in cgs units
         vals(6) = vals(6)/(1.0d0*secyer) ! in years

       end subroutine data_for_extra_history_columns


      
      integer function how_many_extra_profile_columns(id, id_extra)
         use star_def, only: star_info
         integer, intent(in) :: id, id_extra
         integer :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         how_many_extra_profile_columns = 0
      end function how_many_extra_profile_columns
      
      
      subroutine data_for_extra_profile_columns(id, id_extra, n, nz, names, vals, ierr)
         use star_def, only: star_info, maxlen_profile_column_name
         use const_def, only: dp
         integer, intent(in) :: id, id_extra, n, nz
         character (len=maxlen_profile_column_name) :: names(n)
         real(dp) :: vals(nz,n)
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         integer :: k
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
      end subroutine data_for_extra_profile_columns
      

      integer function extras_start_step(id, id_extra)
         use utils_lib, only:alloc_iounit, free_iounit
         use crlibm_lib
         integer, intent(in) :: id, id_extra
         integer :: ierr
         type (star_info), pointer :: s
         integer :: k, k0, num_pts, species, model_number
         real(dp) :: v_esc, time, gamma1_integral, integral_norm, tdyn, &
            max_center_cell_dq, avg_v_div_vesc, energy_removed_layers, dt_next, dt
         real(dp), pointer :: xq(:), xa(:,:), entropy(:), conv_vel(:), xq_new(:)
         character (len=256) :: pulse_fname, photo_fname
         integer :: unit, pulse_num, line_num
         real(dp) :: core_mass
         logical :: just_did_relax
         logical :: exists
         include 'formats'
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return

         !check for helium depletion
         if(.not. s% lxtra3 .and. s% center_he4 < 1d-4) then
            s% lxtra3 = .true. 
         end if

         if (s% m(1)/Msun >= 45d0) then ! activate this only at high masses --  threshold might need to be checked
            ! this is to be able to compute the He core burning phase at the high mass end
            ! turn off at end of convective core phase (He free by then) (cf. the if below this one)
            if (.not. s%lxtra3 .and. & !he is depleted
                 (s% mlt_mixing_type(s% nz) == convective_mixing .or. &
                 s% mlt_mixing_type(s% nz) == softened_convective_mixing .or. &
                 s% mlt_mixing_type(s% nz) == overshoot_mixing) ) & ! the center is not convective or undershooting
                 then
               ! write(*,*) "core is convective!"
               ! -- predictive mixing: better treatment of convective core growth
               s% predictive_mix(1) = .true.
               s% predictive_zone_type(1) = 'burn_He'
               s% predictive_zone_loc(1) = 'any'
               s% predictive_bdy_loc(1) = 'any'
               
               s% predictive_mix(2) = .true.
               s% predictive_zone_type(2) = 'burn_H'
               s% predictive_zone_loc(2) = 'any'
               s% predictive_bdy_loc(2) = 'any'
               
               s% predictive_mix(3) = .true.
               s% predictive_zone_type(3) = 'nonburn'
               s% predictive_zone_loc(3) = 'any'
               s% predictive_bdy_loc(3) = 'any'
               
               s% predictive_mix(4) = .true.
               s% predictive_zone_type(4) = 'burn_Z' ! although mixing might not be instantaneous anymore at this stage...
               s% predictive_zone_loc(4) = 'any'
               s% predictive_bdy_loc(4) = 'any'
               
               ! smooth entropy gradients at edges of convective regions (will make the entropy relax fail)
               s% smooth_gradT = .true.
               s% gradT_smooth_low = -0.01d0
               s% gradT_smooth_high = 0.02d0
               
               ! MLT++ (see MESA paper 2) -- very similar to defaults
               s% gradT_excess_f1 = 1d-4
               s% gradT_excess_f2 = 1d-2
               s% gradT_excess_age_fraction = 0.9d0
               
               s% gradT_excess_lambda1 = 1.0
               s% gradT_excess_beta1 = 0.35
               
               s% gradT_excess_lambda2 = 0.5
               s% gradT_excess_beta2 = 0.25
               
               s% gradT_excess_dlambda = 0.1
               s% gradT_excess_dbeta = 0.1
               
               ! outer boundary: this should be the defaults and have no effect, but who knows... 
               s% which_atm_option = 'simple_photosphere' 
               s% Pextra_factor = -1
               
               ! this includes the change in mean molecular weight in the energy equation
               s% include_dmu_dt_in_eps_grav = .true.
               
               s% min_dt_for_increases_in_convection_velocity = 1d-2
               
               s% delta_lgT_cntr_limit = 0.0025
               s% delta_lgT_cntr_hard_limit = 0.005
               s% dt_div_min_dr_div_cs_limit = 10
               
               ! make the model crash if things go bad during He core burning...save computing time
               s% min_timestep_limit = 10 ! (seconds)
               
               
            else ! core is convective or still burning He
               ! write(*,*) "core is NOT convective anymore!"
               ! reset to defaults the settings above
               s% predictive_mix(1) = .false.
               s% predictive_zone_type(1) = ''
               s% predictive_zone_loc(1) = ''
               s% predictive_bdy_loc(1) = ''
               
               s% predictive_mix(2) = .false.
               s% predictive_zone_type(2) = ''
               s% predictive_zone_loc(2) = ''
               s% predictive_bdy_loc(2) = ''
               
               s% predictive_mix(3) = .false.
               s% predictive_zone_type(3) = ''
               s% predictive_zone_loc(3) = ''
               s% predictive_bdy_loc(3) = ''
               
               s% predictive_mix(4) = .false.
               s% predictive_zone_type(4) = '' ! although mixing might not be instantaneous anymore at this stage...
               s% predictive_zone_loc(4) = ''
               s% predictive_bdy_loc(4) = ''
               
               
               s% smooth_gradT = .false.
               s% gradT_smooth_low = -0.005d0
               s% gradT_smooth_high = 0.01d0
               
               s% gradT_excess_f1 = 1d-4
               s% gradT_excess_f2 = 1d-3
               s% gradT_excess_age_fraction = 0.9d0
               
               s% gradT_excess_lambda1 = 1.0
               s% gradT_excess_beta1 = 0.35
               
               s% gradT_excess_lambda2 = 0.5
               s% gradT_excess_beta2 = 0.25
               
               s% gradT_excess_dlambda = 0.1
               s% gradT_excess_dbeta = 0.1
               
               s% include_dmu_dt_in_eps_grav = .false.
               
               s% min_dt_for_increases_in_convection_velocity = -1
               
               s% delta_lgT_cntr_limit = 0.01d0
               s% delta_lgT_cntr_hard_limit = -1
               s% dt_div_min_dr_div_cs_limit = -1
               
               
               s% min_timestep_limit = 1d-12 ! (seconds)
               
            end if
         end if

         if(s% lxtra3) then
            s% delta_lgT_cntr_limit = s% x_ctrl(12)
            s% delta_lgT_cntr_hard_limit = 2*s% x_ctrl(12)
            ! write(*,*) "check dt limit", s% delta_lgT_cntr_limit
         end if


         s% job% pgstar_flag = .true.
         just_did_relax = .false.
         ! get point where v<vesc
         do k = s% nz, 1, -1  !going outward from the center
            v_esc = s% x_ctrl(1)*sqrt(2*s% cgrav(k)*s% m(k)/(s% r(k)))
            if (s% u(k) > v_esc) then
               exit
            end if
         end do  !now k is the innermost cell faster than the local v_esc

         if (k > s%x_integer_ctrl(13)  .and. s% xtra1 > 0d0) then

            avg_v_div_vesc = 0d0
            energy_removed_layers = 0d0
            do k0 = 1, k
               avg_v_div_vesc = avg_v_div_vesc + s% dm(k0)*s% u(k0)/sqrt(2*s% cgrav(k0)*s% m(k0)/(s% r(k0)))
               energy_removed_layers = energy_removed_layers + &
                  0.5d0*s% dm(k0)*s% u(k0)*s% u(k0) - s% dm(k0)*s% cgrav(k0)*s% m(k0)/s% r(k0)
            end do 
            avg_v_div_vesc = avg_v_div_vesc/(s% m(1) - s% m(k))

            ! compute gamma integral in what will remain of the star
            integral_norm = 0d0
            gamma1_integral = 0d0
            do k0=k,s% nz
               integral_norm = integral_norm + s% P(k0)*s% dm(k0)/s% rho(k0)
               gamma1_integral = gamma1_integral + &
                  (s% gamma1(k0)-4.d0/3.d0)*s% P(k0)*s% dm(k0)/s% rho(k0)
            end do
            gamma1_integral = gamma1_integral/max(1d-99,integral_norm)

            do k0 = k, s% nz  ! go downwards toward the center from k
               if (s% q(k0) < s% x_ctrl(7)*s% q(k)) then
                  exit
               end if
            end do ! now k0 is the mass location 0.98* mass faster than v_esc (0.98 = x_ctrl(7), check inlist_project)
 
            write(*,*) 'Layers above', s% x_ctrl(1), 'times vesc will be removed', k0, s% x_ctrl(7), k, s% q(k)
            write(*,*) 'checking for conditions inside q=', s% x_ctrl(7), 'of material that will remain'
            write(*,*) 'check time left', s% xtra1 + s% xtra2*s% x_ctrl(6) - s% time
            write(*,*) 'max vel inside fraction of cutoff', maxval(abs(s% u(k0:s% nz)))/1e5
            write(*,*) 'max v/cs inside fraction of cutoff', maxval(abs(s% u(k0:s% nz)/s% csound(k0:s% nz)))
            write(*,*) 'average v/vesc outside cutoff', avg_v_div_vesc
            write(*,*) 'Kinetic plus potential energy outside cutoff', energy_removed_layers
            write(*,*) 'mass inside cutoff', k, s% m(k)/Msun
            write(*,*) 'relax counter', s% ixtra2
            ! If enough time has passed after a pulse, then turn off Riemann hydro
            ! also, wait for the inner 98% of what would remain of the star to be below a threshold
            ! in v/cs, and below a threshold in v
            ! Verify also that the conditions to turn on Riemann hydro are not still satisfied
            if (s% lxtra1 .and. s% xtra1 > 0d0 .and. s% time > s% xtra1 + s% xtra2*s% x_ctrl(6)&
               .and. maxval(abs(s% u(k0:s% nz)))/1e5 < s% x_ctrl(8) &
               .and. maxval(abs(s% u(k0:s% nz)/s% csound(k0:s% nz))) < s% x_ctrl(9) &
               .and. log10_cr(s% T(s% nz)) < s% x_ctrl(10) &
               .and. (.not. (log10_cr(s% T(s% nz)) > s% x_ctrl(2) &
               .and. gamma1_integral < s% x_ctrl(3)) .or. (maxval(abs(s% u(k0:s% nz)))/1e5 < s% x_ctrl(14) &
               .and. maxval(abs(s% u(k0:s% nz)/s% csound(k0:s% nz))) < s% x_ctrl(15))) &
               .and. avg_v_div_vesc > s% x_ctrl(13)&
               .and. energy_removed_layers > 0d0 &
               .and. (log10_cr(s% power_neutrinos) < s% x_ctrl(16)) .or. (maxval(abs(s% u(k0:s% nz)))/1e5 < s% x_ctrl(14) &
               .and. maxval(abs(s% u(k0:s% nz)/s% csound(k0:s% nz))) < s% x_ctrl(15))) then
               s% ixtra2 = s% ixtra2 + 1
            else
               s% ixtra2 = 0
            end if
            if(s% ixtra2 >= s% x_integer_ctrl(1)) then
               ! here "k" is the index of the last cell to be taken. From 1 to k it's the ejecta.
               write(*,*) "saving info on the ejected material!"
               unit=alloc_iounit(ierr)
               if (ierr /= 0) then
                  print *, "error in alloc_unit",unit, ierr
                  stop
               end if
               pulse_num = s% ixtra1 + 1
               s% ixtra1 = pulse_num !update for next pulse...
               write(pulse_fname, '(a,i0.3,a)') 'pulse', pulse_num, '.data'
               pulse_fname = 'LOGS/'//trim(pulse_fname)
               print *, pulse_fname
               open(unit=unit, file=trim(pulse_fname), action='write', iostat=ierr)
               if (ierr /= 0) then
                  write(*, *) 'Failed to open pulse file ', trim(pulse_fname)
                  stop
               end if
               ! header
               write(unit, *) 'star_age', char(9), s% star_age
               write(unit, *) 'mass_removed', char(9), (s% m(1) - s% m(k))/Msun
               write(unit, *) 'num_zones_removed',char(9), k 
               write(unit, *) char(9), 'zone', char(9), 'velocity',char(9), 'dm',char(9), &
                    'logR', char(9), 'log_Rho',char(9), 'log_T',char(9), &
                    !     'logPgas',char(9),
                    'energy',char(9), 'entropy',char(9), 'opacity',char(9), &
                     'h1',char(9),'he3',char(9),'he4',char(9),'c12',char(9),'n14',char(9),'o16',char(9),'ne20',char(9), & !! this should be improved to be net independent
                     'mg24',char(9),'si28',char(9),'s32',char(9),'ar36',char(9),'ca40',char(9),'ti44',char(9), &
                     'cr48',char(9),'cr56',char(9),'fe52',char(9),'fe54',char(9),'fe56',char(9),'co56',char(9),'ni56'
               do line_num=1, k
                  write(unit, *)  line_num, char(9), & 
                       s% xh(s% i_u, line_num),char(9), s% dm(line_num), char(9), &
                       log10_cr(exp_cr(s% xh(s% i_lnR, line_num))), char(9), &
                       log10_cr(exp_cr(s% xh(s% i_lnd, line_num))), char(9), &
                       log10_cr(exp_cr(s% xh(s% i_lnT, line_num))), char(9), &
                       ! log_cr(EXP(s% xh(s% i_lnPgas, line_num))), char(9), &
                       s% energy(line_num),  char(9), s% entropy(line_num), char(9), s% opacity(line_num), char(9), & 
                       s%xa(s%net_iso(ih1), line_num)+s%xa(s%net_iso(iprot), line_num), char(9), & ! sum iprot and h1 from approx nets
                       s%xa(s%net_iso(ihe3), line_num), char(9), &
                       s%xa(s%net_iso(ihe4), line_num), char(9), &
                       s%xa(s%net_iso(ic12), line_num), char(9), &
                       s%xa(s%net_iso(in14), line_num), char(9), &
                       s%xa(s%net_iso(io16), line_num), char(9), &
                       s%xa(s%net_iso(ine20), line_num), char(9), &
                       s%xa(s%net_iso(img24), line_num), char(9), &
                       s%xa(s%net_iso(isi28), line_num), char(9), &
                       s%xa(s%net_iso(is32), line_num), char(9), &
                       s%xa(s%net_iso(iar36), line_num), char(9), &
                       s%xa(s%net_iso(ica40), line_num), char(9), &
                       s%xa(s%net_iso(iti44), line_num), char(9), &
                       s%xa(s%net_iso(icr48), line_num), char(9), &
                       s%xa(s%net_iso(icr56), line_num), char(9), &
                       s%xa(s%net_iso(ife52), line_num), char(9), &
                       s%xa(s%net_iso(ife54), line_num), char(9), &
                       s%xa(s%net_iso(ife56), line_num), char(9), &
                       s%xa(s%net_iso(ico56), line_num), char(9), &
                       s%xa(s%net_iso(ini56), line_num)
               end do
               close(unit)
               call free_iounit(unit)
               write(photo_fname, '(a,i0.3)') 'pulse', pulse_num
               photo_fname = trim(photo_fname)
               call star_save_for_restart(id, 'photos/'//photo_fname, ierr)
	       if(ierr/=0) call mesa_error(__FILE__,__LINE__)
               call save_profile(id, id_extra, 55, ierr)

               write(*,*) "Relaxing model to lower mass!"
               time = s% time
               model_number = s% model_number

               write(*,*) "removing cells", k, s% m(k)/Msun
               species = s% species
               num_pts = s% nz - k + 1
               allocate(xq(s% nz - k + 1), xa(species, s% nz - k + 1), &
                  entropy(s% nz - k + 1), conv_vel(s% nz - k + 1))
               !xq(:) = (s% q(k) - s% q(k:s% nz))/s% q(k)
               !need to compute cell-centered q for remaining mass
               xq(1) = s% dq(k)/2/s% q(k)
               do k0 = 1, s% nz - k
                  xq(1+k0) = xq(1+k0-1) + (s% dq(k+k0) + s% dq(k+k0-1))/s% q(k)/2
                  !write(*,*) k0, s% dq(k+k0), xq(1+k0), (s% q(k) - s% q(k0+k))/s% q(k)
               end do

               xa(:,:) = s% xa(:,k:s% nz)
               entropy(:) = s% entropy(k:s% nz)*avo*kerg
               conv_vel(:) = s% conv_vel(k:s% nz)
               max_center_cell_dq = s% max_center_cell_dq
               s% max_center_cell_dq = s% dq(s% nz)
               dt = s% dt
               dt_next = s% dt_next

               call star_read_controls(id, 'inlist_hydro_off', ierr)
	       if(ierr/=0) call mesa_error(__FILE__,__LINE__)

               s% job% pgstar_flag = .false.
               call star_set_u_flag(id, .false., ierr)
               call star_set_RTI_flag(id, .false., ierr)
               s% have_previous_conv_vel = .false.
               s% initial_z = 0.02
               s% initial_y = 0.28
               s% initial_mass = s% m(k)/Msun
               ! avoid making photos
               s% photo_interval = 10000000
               ! write(*,*) "no predictive mixing for relaxation phase"
               ! s% predictive_mix(1) = .false.
               ! s% predictive_mix(2) = .false.
               ! s% predictive_mix(3) = .false.
               ! s% predictive_mix(4) = .false.

               call star_change_to_new_net(id, .true., 'basic.net', ierr)
	       if(ierr/=0) call mesa_error(__FILE__,__LINE__)
               call star_load_zams(id, ierr)
	       if(ierr/=0) call mesa_error(__FILE__,__LINE__)
               write(*,*) "*********************"
               write(*,*) "check ierr", ierr
               if (ierr /= 0) STOP "relaxing mass failed"
               write(*,*) "Done relaxing mass!!"
               write(*,*) "*********************"

               call star_change_to_new_net(id, .true., s% job% new_net_name, ierr)
	       if(ierr/=0) call mesa_error(__FILE__,__LINE__)
               call star_relax_composition( &
                  id, s% job% num_steps_to_relax_composition, num_pts, species, xa, xq, ierr)
	       if(ierr/=0) call mesa_error(__FILE__,__LINE__)

               write(*,*) "*********************"
               write(*,*) "check ierr", ierr
               if (ierr /= 0) STOP "relaxing composition failed"
               write(*,*) "Done relaxing composition!!"
               write(*,*) "*********************"

               call star_relax_entropy( &
                  id, s% job% max_steps_to_relax_entropy, num_pts, entropy, xq, ierr)
	       if(ierr/=0) call mesa_error(__FILE__,__LINE__)

               write(*,*) "*********************"
               write(*,*) "check ierr", ierr
               if (ierr /= 0) STOP "relaxing entropy failed"
               write(*,*) "Done relaxing entropy!!"
               write(*,*) "*********************"


               s% photo_interval = 50
               s% dt_next = 1d0!dt_next
               s% dt = 1d0!dt
               s% time = time
               s% model_number = model_number

               ! this avoid the history file from being rewritten
               s% doing_first_model_of_run = .false.
               ! putting xtra29 negative avoids the surface pressure timestep limit
               s% xtra29 = -1d0
               ! reloead relax counter
               s% ixtra2 = 0

               !take care of convective velocities
               allocate(xq_new(s% nz))
               xq_new(1) = s% dq(1)/2
               do k0 = 2, s% nz
                  xq_new(k0) = xq_new(k0-1) + (s% dq(k0) + s% dq(k0-1))/2
               end do
               !quick test for convective vels
               do k=1, s% nz
                  do k0=1, num_pts
                     if(xq_new(k) > xq(k0)) then
                        s% conv_vel(k) = conv_vel(k0)
                        s% conv_vel_old(k) = conv_vel(k0)
                        s% conv_vel_older(k) = conv_vel(k0)
                        exit
                     end if
                  end do
               end do
               s% have_previous_conv_vel = .true.

               s% lxtra1 = .false.
               s% max_center_cell_dq = max_center_cell_dq
               just_did_relax = .true.
               
               deallocate(xq, xa, entropy, conv_vel, xq_new)
               ! write(*,*) "re-turn on predictive mixing"
               ! s% predictive_mix(1) = .true.
               ! s% predictive_mix(2) = .true.
               ! s% predictive_mix(3) = .true.
               ! s% predictive_mix(4) = .true.

            end if
         end if

         ! turn on Riemman hydro if close to instability
         integral_norm = 0d0
         gamma1_integral = 0d0
         do k=1,s% nz
            integral_norm = integral_norm + s% P(k)*s% dm(k)/s% rho(k)
            gamma1_integral = gamma1_integral + &
               (s% gamma1(k)-4.d0/3.d0)*s% P(k)*s% dm(k)/s% rho(k)
         end do
         gamma1_integral = gamma1_integral/max(1d-99,integral_norm)
         write(*,*) "check gamma integral", gamma1_integral
         if (.not. just_did_relax .and. .not. s% lxtra1 .and. ((log10_cr(s% T(s% nz)) > s% x_ctrl(2) .and. &
            gamma1_integral < s% x_ctrl(3)) .or. log10_cr(s% T(s% nz)) > s% x_ctrl(10) )) then
            write(*,*) "Turning on Riemann hydro and RTI!"
            call star_read_controls(id, 'inlist_hydro_on', ierr)
	    if(ierr/=0) call mesa_error(__FILE__,__LINE__)
            call star_read_controls(id, 'inlist_after_first_pulse', ierr)
	    if(ierr/=0) call mesa_error(__FILE__,__LINE__)

            s% job% pgstar_flag = .false.
            call star_set_v_flag(id, .false., ierr)
	    if(ierr/=0) call mesa_error(__FILE__,__LINE__)
            call star_set_u_flag(id, .true., ierr)
	    if(ierr/=0) call mesa_error(__FILE__,__LINE__)
            call star_set_RTI_flag(id, .true., ierr)
	    if(ierr/=0) call mesa_error(__FILE__,__LINE__)
            !if(.not. s% doing_first_model_of_run) then
            !   s% dt_next = 0.001d0
            !end if
            s% lxtra1 = .true.
            s% lxtra2 = .true.
            s% xtra1 = -1
         end if

         ! check time to relax model
         ! the star is considered to be in a pulse if the velocity is high enough
         ! only check velocity in x_ctrl(7) of the mass of the star
         do k0 = 1, s% nz
            if (s% q(k0) < s% x_ctrl(7)) then
               exit
            end if
         end do
         write(*,*) s% lxtra1, k0, s% x_ctrl(7)
         if (s% lxtra1 .and. k0 < s% nz) then
            write(*,*) "check max vel", s% q(k0), maxval(abs(s% u(k0:s% nz)))/1e5, s% x_ctrl(4), s% xtra1
         end if
         if (s% lxtra1 .and. s% xtra1 < 0d0 .and. maxval(abs(s% u(k0:s% nz)))/1e5 > s% x_ctrl(4)) then
            core_mass = s% m(1)
            if(s% o_core_mass > 0d0) then
               core_mass = s% o_core_mass*Msun
            else if(s% c_core_mass > 0d0) then
               core_mass = s% c_core_mass*Msun
            else if (s% he_core_mass > 0d0) then
               core_mass = s% he_core_mass*Msun
            end if
            do k=1, s% nz
               if (s% m(k)/core_mass < s% x_ctrl(5)) then
                  exit
               end if
            end do
            tdyn = 1/sqrt(s% m(k)/(4d0/3d0*pi*s% r(k)**3)*standard_cgrav)
            s% xtra1 = s% time
            s% xtra2 = tdyn
            write(*,*) "reached high velocities, stopping in at least", s% x_ctrl(6)*tdyn, "seconds"
         end if

         extras_start_step = keep_going
      end function extras_start_step
      
      subroutine null_binary_controls(id, binary_id, ierr)
         integer, intent(in) :: id, binary_id
         integer, intent(out) :: ierr
         ierr = 0
      end subroutine null_binary_controls

      ! returns either keep_going or terminate.
      integer function extras_finish_step(id, id_extra)
         use run_star_support
         integer, intent(in) :: id, id_extra
         integer :: ierr,k
         type (star_info), pointer :: s
         real(dp) :: err    
         include 'formats'
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return

         if(abs(s% xtra25_old - s% center_h1) > 0.001) then
            s% dt_next = min(s% dt_next, s% dt * s% min_timestep_factor)
            write(*,*) "reducing dt due to large change in central hydrogen"
         else if(abs(s% xtra26_old - s% center_he4) > 0.001) then
            s% dt_next = min(s% dt_next, s% dt * s% min_timestep_factor)
            write(*,*) "reducing dt due to large change in central helium"
         else if(abs(s% xtra27_old - s% center_c12) > 0.001) then
            s% dt_next = min(s% dt_next, s% dt * s% min_timestep_factor)
            write(*,*) "reducing dt due to large change in central carbon"
         else if(abs(s% xtra28_old - s% center_o16) > 0.001) then
            s% dt_next = min(s% dt_next, s% dt * s% min_timestep_factor)
            write(*,*) "reducing dt due to large change in central oxygen"
         end if
         if(s% xtra29 > 0d0 .and. log10_cr(abs(s% P(1)/s% xtra29)) > 0.1) then
            s% dt_next = min(s% dt_next, s% dt * s% min_timestep_factor)
            write(*,*) "reducing dt due to large change in surface pressure"
         end if
         s% xtra25 = s% center_h1
         s% xtra26 = s% center_he4
         s% xtra27 = s% center_c12
         s% xtra28 = s% center_o16
         s% xtra29 = s% P(1)

         ! limit velocity
         if (s% u_flag) then
            do k=1, s% nz
               s% xh(s% i_u,k) = min(s% xh(s% i_u,k), 1d5*s% x_ctrl(11))!clight*s% x_ctrl(11))
            end do
         end if


         ! if the core is oxygen rich, there isn't silicon yet (so no oxygen burning has occurred)
         ! and the hydro is on (lxtra1 = .true.)
         ! and we haven't already done this (lxtra4 = .false.)
         ! then save a profile: we are getting closer to a pulse!
         if (s% xa(s%net_iso(io16),s%nz) >= 0.7) then
            if (s%xa(s%net_iso(isi28), s%nz) <= 0.1 .and. s% lxtra1 .eqv. .true.) then
               if (.not. s% lxtra4) then
                  s% lxtra4 = .true.
                  call star_write_model(id, 'pre-first-pulse.mod' , ierr)
                  s% need_to_update_history_now = .true.
                  s% need_to_save_profiles_now = .true.
                  s% save_profiles_model_priority = 54
               end if
            end if
         end if
         
         extras_finish_step = keep_going


         ! This doesn't work, because the infalling core accelerates much faster 
         ! than the shock moves out, therefore, MESA crashes because the input
         ! physics for the core loses sense.
         ! if (s% fe_core_infall < 1d8 .and. &
         !      s% fe_core_infall > 0.99d0*1d8) then
         !    write(*,1) 'nearing fe_core_infall limit', &
         !         s% fe_core_infall, s% fe_core_infall_limit
         ! end if
         
         ! if (s% fe_core_infall >= 1d8) then
         !    ! we have reached the core-collapse velocity threshold
         !    err = abs(s% fe_core_infall - 1d8)/ &
         !         (s% when_to_stop_atol + s% when_to_stop_rtol*max(abs(s% fe_core_infall),abs(1d8)))
         !    if (err > 1) then
         !       s% dt_next = 0.5*s%dt_next
         !       s% dt = 0.5*s%dt
         !       extras_finish_step = keep_going
         !       write(*,'(/,a,5e20.10)') &
         !             'redo with smaller timestep to get closer to stopping target', &
         !             s% fe_core_infall, 1d8
         !    else
         !       ! check if a pulse is going on:
         !       if (s% lxtra1 .eqv. .true.) then ! hydro is on
         !          ! if the maximum velocity is below x_ctrl(8)
         !          ! then it's fine (N.B. the infall velocity is negative, it will not be a max
         !          if (maxval(s% u(1:s% nz))/1e5 < s% x_ctrl(8)) then
         !             extras_finish_step = terminate                     
         !             print *, 'termination code: fe_core_infall_limit'
         !          end if
         !       end if
         !    end if
         ! end if


         if (s% u_flag .and. (s% total_linear_kinetic_energy - s% total_internal_energy) >= 0) then 
            if (minval(s% xh(s% i_u,1:s%nz)) >= 0d0) then 
               if (s% lxtra1 .eqv. .true.) then
                  print *, "minimum velocity is positive ..." 
                  if (maxval(s% T(1:s%nz)) <= 5d8) then
                     print *, '****************************'  
                     print *, ' star is going PISN!' 
                     print *, '****************************'
                     extras_finish_step = terminate  
                     print *, ' '  
                     print *, 'Terminated because:'
                     print *, '  Ekin - Etot > 1d5'
                     print *, ' min(u) > 0d0'  
                     print *, ' max(T) < 5d8' 
                  end if
               end if
            end if
         end if

      end function extras_finish_step
      
      
      ! routines for saving and restoring extra data so can do restarts
         
         ! put these defs at the top and delete from the following routines
         !integer, parameter :: extra_info_alloc = 1
         !integer, parameter :: extra_info_get = 2
         !integer, parameter :: extra_info_put = 3
      
      
      subroutine alloc_extra_info(s)
         integer, parameter :: extra_info_alloc = 1
         type (star_info), pointer :: s
         call move_extra_info(s,extra_info_alloc)
      end subroutine alloc_extra_info
      
      
      subroutine unpack_extra_info(s)
         integer, parameter :: extra_info_get = 2
         type (star_info), pointer :: s
         call move_extra_info(s,extra_info_get)
      end subroutine unpack_extra_info
      
      
      subroutine store_extra_info(s)
         integer, parameter :: extra_info_put = 3
         type (star_info), pointer :: s
         call move_extra_info(s,extra_info_put)
      end subroutine store_extra_info
      
      
      subroutine move_extra_info(s,op)
         integer, parameter :: extra_info_alloc = 1
         integer, parameter :: extra_info_get = 2
         integer, parameter :: extra_info_put = 3
         type (star_info), pointer :: s
         integer, intent(in) :: op
         
         integer :: i, j, num_ints, num_dbls, ierr
         
         i = 0
         ! call move_int or move_flg
         call move_int(redo_count)   
         num_ints = i
         
         i = 0
         num_dbls = i
         
         if (op /= extra_info_alloc) return
         if (num_ints == 0 .and. num_dbls == 0) return
         
         ierr = 0
         call star_alloc_extras(s% id, num_ints, num_dbls, ierr)
         if (ierr /= 0) then
            write(*,*) 'failed in star_alloc_extras'
            write(*,*) 'alloc_extras num_ints', num_ints
            write(*,*) 'alloc_extras num_dbls', num_dbls
            stop 1
         end if
         
         contains
         
         subroutine move_dbl(dbl)
            double precision :: dbl
            i = i+1
            select case (op)
            case (extra_info_get)
               dbl = s% extra_work(i)
            case (extra_info_put)
               s% extra_work(i) = dbl
            end select
         end subroutine move_dbl
         
         subroutine move_int(int)
            integer :: int
            i = i+1
            select case (op)
            case (extra_info_get)
               int = s% extra_iwork(i)
            case (extra_info_put)
               s% extra_iwork(i) = int
            end select
         end subroutine move_int
         
         subroutine move_flg(flg)
            logical :: flg
            i = i+1
            select case (op)
            case (extra_info_get)
               flg = (s% extra_iwork(i) /= 0)
            case (extra_info_put)
               if (flg) then
                  s% extra_iwork(i) = 1
               else
                  s% extra_iwork(i) = 0
               end if
            end select
         end subroutine move_flg
      
      end subroutine move_extra_info

      end module run_star_extras
      
