! ***********************************************************************
!
!   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
      
      !include 'standard_run_star_extras.inc'
      
      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.inc)

         s % other_mixing => pavel_mixing_micro
         
         
      end subroutine extras_controls

      subroutine pavel_mixing_micro(id, ierr)

!         use chem_def, only: ih1, ihe4, in14, ic12, io16, ih2, ihe3
         use chem_def
         use chem_lib
         integer, intent(in) :: id
         integer, intent(out) :: ierr
         integer :: k, nz, h1, he4
         double precision :: x, dq00, dqm1, dqsum, T, kap, rho, cp, &
           molecular_diffusivity, mol_diff, zat, lnlam, &
           micro_diffusivity, y, z, mol_dif, amas 

         !integer, dimension(:), pointer :: Z_plus_N
         type (star_info), pointer :: s
         ierr = 0

         call get_star_ptr(id, s, ierr)
         if (ierr /= 0) return

         !! FOR DEBUGGING
         !print*,'NET_ISO',sum(s% net_iso(:))
         !print*,'****************************************'
         !print*,'CHEM_ID',s% chem_id
         !print*,'****************************************'
         !print*,'SPECIES',s% species
         !print*,'****************************************'

         nz = s% nz

       ! add microscopic diffusivity to D_mix
       ! clip_D_limit = 1 ! zero mixing diffusion coeffs that are smaller than this
       ! in inlist controls make clip_D_limt = 0.0001

         do k = 2, 4 ! start at 2 since no mixing across surface
           ! interpolate values at face

            !! FOR DEBUGGING
            !print*,'ABAR= ', s% abar(k)
            !print*,'ZBAR= ', s% zbar(k)
            !h1 = s% net_iso(ih1)
            !he4 = s% net_iso(ihe4)
            !x = s% xa(h1,k)
            !y = s% xa(he4,k)

            amas= s% abar(k)

            zat= s% zbar(k)

            !el_ch=1.60217653d-19

            !amu=1.660538921d-27

            !kb=1.3806488d-23
            
            dq00 = s% dq(k)
            dqm1 = s% dq(k-1)
            dqsum = dq00 + dqm1
            kap = (dqm1*s% opacity(k) + dq00*s% opacity(k-1))/dqsum
            rho = (dqm1*s% rho(k) + dq00*s% rho(k-1))/dqsum
            cp = (dqm1*s% cp(k) + dq00*s% cp(k-1))/dqsum
            T = (dqm1*s% T(k) + dq00*s% T(k-1))/dqsum

            lnlam=log((2./3.)*(1./(qe**3.))*(sqrt(((amas/zat)*amu*&
            ((boltzm*T)**3.))/(pi*rho*((zat)**5.)))))

            !!!radiative_diffusivity = 4*crad*T**4/(15*kap*clight*rho**2) !! THESE ARE THE OLD EXPRESSIONS USED TO ACCOUNT MICROSCOPIC DIFFUSION
            !!!molecular_diffusivity = 1.84d-17*(1+7*x)*(T**2.5/rho)

            mol_dif=2.21d-15*(amas**0.5/((zat**4)*lnlam)*(T**2.5/rho))

            !! FOR DEBUGGING !!
            !print*, '   '
            !print*, '********************'
            !print*, 'MOL_DIFF= ', mol_dif
            !print*, '********************'
            !print*, 'LN_LAM= ', lnlam
            !print*, '********************'
            !print*, '   '

            !s% D_mix(k) = s% D_mix(k) + mol_dif
            s% D_mix(k) = s% D_mix(k) + mol_dif + 100.0

         end do

      end subroutine pavel_mixing_micro

      end module run_star_extras