! ***********************************************************************
!
!   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 alert_lib
      use mlt_def
      
      implicit none
      
      ! these routines are called by the standard run_star check_model
      contains
      
      subroutine extras_controls(s, ierr)
         type (star_info), pointer :: s
         integer, intent(out) :: ierr
         ierr = 0
         
         ! this is the place to set any procedure pointers you want to change
         ! e.g., other_wind, other_mixing, other_energy  (see star_data.dek)

         s % other_mixing => RCB_mixing

      end subroutine extras_controls
      
      
      integer function extras_startup(s, id, restart, ierr)
         type (star_info), pointer :: s
         integer, intent(in) :: id
         logical, intent(in) :: restart
         integer, intent(out) :: ierr
         ierr = 0
         extras_startup = 0
         if (.not. restart) then
            call alloc_extra_info(s)
         else ! it is a restart
            call unpack_extra_info(s)
         end if
      end function extras_startup
      

      ! returns either keep_going, retry, backup, or terminate.
      integer function extras_check_model(s, id, id_extra)
         type (star_info), pointer :: s
         integer, intent(in) :: id, id_extra
         extras_check_model = keep_going         
         if (.false. .and. s% star_mass_h1 < 0.35d0) then
            ! stop when star hydrogen mass drops to specified level
            extras_check_model = terminate
            write(*, *) 'have reached desired hydrogen mass'
            return
         end if
      end function extras_check_model


      integer function how_many_extra_log_columns(s, id, id_extra)
         type (star_info), pointer :: s
         integer, intent(in) :: id, id_extra
         how_many_extra_log_columns = 0
      end function how_many_extra_log_columns
      
      
      subroutine data_for_extra_log_columns(s, id, id_extra, n, names, vals, ierr)
         type (star_info), pointer :: s
         integer, intent(in) :: id, id_extra, n
         character (len=maxlen_log_column_name) :: names(n)
         double precision :: vals(n)
         integer, intent(out) :: ierr
         ierr = 0
      end subroutine data_for_extra_log_columns

      
      integer function how_many_extra_profile_columns(s, id, id_extra)
         type (star_info), pointer :: s
         integer, intent(in) :: id, id_extra
         how_many_extra_profile_columns = 0
      end function how_many_extra_profile_columns
      
      
      subroutine data_for_extra_profile_columns(s, id, id_extra, n, nz, names, vals, ierr)
         type (star_info), pointer :: s
         integer, intent(in) :: id, id_extra, n, nz
         character (len=maxlen_profile_column_name) :: names(n)
         double precision :: vals(nz,n)
         integer, intent(out) :: ierr
         integer :: k
         ierr = 0
         ! here is an example for adding a profile column
         !if (n /= 1) stop 'data_for_extra_profile_columns'
         !names(1) = 'beta'
         !do k = 1, nz
         !   vals(k,1) = s% Pgas(k)/s% P(k)
         !end do
      end subroutine data_for_extra_profile_columns
      

      ! returns either keep_going or terminate.
      integer function extras_finish_step(s, id, id_extra)
         type (star_info), pointer :: s
         integer, intent(in) :: id, id_extra
         integer :: ierr
         extras_finish_step = keep_going
         call store_extra_info(s)

         ! to save a profile, 
            ! s% need_to_save_profiles_now = .true.
         ! to update the star log,
            ! s% need_to_update_logfile_now = .true.
            
      end function extras_finish_step
      
      
      subroutine extras_after_evolve(s, id, id_extra, ierr)
         type (star_info), pointer :: s
         integer, intent(in) :: id, id_extra
         integer, intent(out) :: ierr
         ierr = 0
      end subroutine extras_after_evolve
      
      
      ! 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    
         num_ints = i
         
         i = 0
         ! call move_dbl       
         
         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


      ! examples of other_mixing subroutines
 
      
      ! examples
                  
      
      subroutine RCB_mixing(id, ierr)

         use chem_def, only: io18

         integer, intent(in) :: id
         integer, intent(out) :: ierr         
         integer :: k, nz
         double precision :: lnT_mix, lnr_mix, ind_D     
         double precision :: T,m,o18,o18max,o18f,kmax,Tmax,m_kmax,kf1,m_kf1,kf2,m_kf2,r,res1,res2   
         double precision :: lgr_mix, Dr, lgT_mix, o18_mix, he_en, D_mlt, D_add, Dexp
         type (star_info), pointer :: s
         double precision :: gradr, grada, tri_alfa, rho, radius, D_mix, mlt_D, Tlim, dm_dr
         ierr = 0
         !double precision,dimension(size(gradr)):: D_add
         call get_star_ptr(id, s, ierr)
         if (ierr /= 0) return

         nz = s% nz
         o18 = s% net_iso(io18)
         
         !Dr = s% x_ctrl(1)
         !Tlim = s% x_ctrl(2)
         !lgT_mix = s% x_ctrl(4) 
         !lnT_mix = lgT_mix*ln10

        o18max=0d0
        o18f=0d0
        kmax=0d0
        kf1=0d0
        kf2=0d0
        m_kmax=0d0
        m_kf1=0d0
        m_kf2=0d0

        do k = 2, nz ! start at 2 since no mixing across surface
            dm_dr = 4*3.142*(s% r(k))**2*(s% rho(k))
            D_add = 1d51/(dm_dr)**2
            T = 10**(s% lnT(k)/ln10)
            m = s% m(k)            
            if (s% xa(o18,k) > o18max) then !find zone where o18 has max abundance
                kmax = k
                o18max = s% xa(o18,k)
                Tmax= T
                m_kmax = m
            end if
        end do  
       
        res2=abs(s% xa(o18,2) - 0.9*o18max) !find zone where o18 has 0.9*max abundance
        res1=0

        do r=kmax,nz
           m = s% m(r)            
           res1= abs(s% xa(o18,r) - 0.9*o18max) 
           if (res1 < res2) then
              res2 = res1
              kf1= r
              o18f= s% xa(o18,r)
              m_kf1 = m
           end if
        end do

        res2=abs(s% xa(o18,2) - 0.9*o18max) !find zone where o18 has 0.9*max abundance
        res1=0
        
        do r=2,kmax
           m = s% m(r)            
           res1= abs(s% xa(o18,r) - 0.9*o18max) 
           if (res1 < res2) then
              res2 = res1
              kf2= r
              o18f= s% xa(o18,r)
              m_kf2 = m
           end if
        end do
    
        do k = 2, nz ! Assign D
           D_add = 1d60/(dm_dr)**2
           T = 10**(s% lnT(k)/ln10)
           m = s% m(k)
           if (m > m_kf2) then
                 s% D_mix(k)= D_add
                 s% mixing_type(k) = semiconvective_mixing 
           else if (m<=m_kf2 .and. m>=m_kf1) then
                 s% D_mix(k)= (1d50*exp(-(T-Tmax)/(0.01*Tmax)))/(dm_dr)**2
                 s% mixing_type(k) = thermo_haline_mixing 
           else if (m < m_kf1) then
                s% D_mix(k)= 0.0
                s% mixing_type(k) =  semiconvective_mixing
           end if     
        end do     

      end subroutine RCB_mixing
                  

      end module run_star_extras