! ***********************************************************************
!
!   Copyright (C) 2011  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
!
! ***********************************************************************
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
!ly 8/9/15 Define global variables of the k-omega model  
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
      module define_convection_variables
      implicit none
!        cgma  = turbulent Prandtl number for k-omega model
!        cmu   = coefficient of eddy viscousity
!        cx3   = coefficient of diffusivity for convective mixing 
!        cta   = coefficient of temperature fluctuation model of turbulence
!        ct    = coefficient of buoyancy production of turbulence
!        ch    = coefficient of turbulent heat flux
!        afa   = coefficient of macro-length model for convective envelope
!        alf   = coefficient of macro-length model for convective core
      REAL*8,PARAMETER :: cgma=1.5d0, cmu=0.09d0
      real*8,parameter :: cta=0.50d0, cx3=1.00d0   
      real*8,parameter :: afa=0.70d0, ct=0.100d0       
      real*8,parameter :: alf=0.06d0, ch=2.344d0
      real*8,parameter :: ypn=0.00d0, ypm=1.00d0
!
      real*8,parameter :: ynl=1d0/(2d0+ypm)        
      real*8,parameter :: cl2=cmu*sqrt(cmu), cl=sqrt(cl2)        
      real*8,parameter :: cs=cmu/cl
      real*8,parameter :: cbe=1d0            !+cmu/(cl-cmu)
      real*8,parameter :: cc1=cbe*ct/ch
      real*8,parameter :: cc2=(cta+1d0)*ct
!        NG    = level of multi-grid method
!        JMAX  = maximum of mesh points in a stellar structure
!        ss    = variables of stellar structure used in the k-omega model
!        dxr   = dimension of the calculated region
      integer,parameter :: nv=2,NG=10,ngp=2**NG+1,JMAX=19999
      REAL*8 ss(5,ngp),dxr
!---
      real*8 ftlxyc(15,100000)
      integer klim,kcxh
!---
      end module define_convection_variables
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
      module run_star_extras
      use define_convection_variables
      use star_lib
      use star_def
      use const_def
      
      implicit none
      
      integer :: time0, time1, clock_rate

      
      ! these routines are called by the standard run_star check_model
      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_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  
!cxh
        s% other_D_mix =>null_other_D_mix
        !s% other_D_mix => my_other_D_mix
        !s% other_mlt => my_stc_mlt
        s% other_mlt => my_mlt
!cxh
      end subroutine extras_controls
      
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  


      
      subroutine my_mlt(  &
            id, k, cgrav, m, r, L, X, &            
            T_face, rho_face, P_face, &
            chiRho_face, chiT_face, &
            Cp_face, opacity_face, grada_face, &            
            alfa, beta, & ! f_face = alfa*f_00 + beta*f_m1
            T_00, T_m1, rho_00, rho_m1, P_00, P_m1, &
            chiRho_00, d_chiRho_00_dlnd, d_chiRho_00_dlnT, &
            chiRho_m1, d_chiRho_m1_dlnd, d_chiRho_m1_dlnT, &
            chiT_00, d_chiT_00_dlnd, d_chiT_00_dlnT, &
            chiT_m1, d_chiT_m1_dlnd, d_chiT_m1_dlnT, &
            Cp_00, d_Cp_00_dlnd, d_Cp_00_dlnT, &
            Cp_m1, d_Cp_m1_dlnd, d_Cp_m1_dlnT, &
            opacity_00, d_opacity_00_dlnd, d_opacity_00_dlnT, &
            opacity_m1, d_opacity_m1_dlnd, d_opacity_m1_dlnT, &
            grada_00, d_grada_00_dlnd, d_grada_00_dlnT, &
            grada_m1, d_grada_m1_dlnd, d_grada_m1_dlnT, &            
            gradr_factor, gradL_composition_term, &
            alpha_semiconvection, semiconvection_option, &
            thermohaline_coeff, thermohaline_option, &
            dominant_iso_for_thermohaline, &
            mixing_length_alpha, alt_scale_height, remove_small_D_limit, &
            MLT_option, Henyey_y_param, Henyey_nu_param, &
            gradT_smooth_low, gradT_smooth_high, smooth_gradT, &
            normal_mlt_gradT_factor, &
            prev_conv_vel, max_conv_vel, g_theta, dt, tau, just_gradr, &
            mixing_type, mlt_basics, mlt_partials1, ierr)
         use star_lib, only: star_mlt_eval
         integer, intent(in) :: id ! id for star         
         integer, intent(in) :: k ! cell number or 0 if not for a particular cell         
         real(dp), intent(in) :: &
            cgrav, m, r, L, X, &            
            T_face, rho_face, P_face, &
            chiRho_face, chiT_face, &
            Cp_face, opacity_face, grada_face, &            
            alfa, beta, &
            T_00, T_m1, rho_00, rho_m1, P_00, P_m1, &
            chiRho_00, d_chiRho_00_dlnd, d_chiRho_00_dlnT, &
            chiRho_m1, d_chiRho_m1_dlnd, d_chiRho_m1_dlnT, &
            chiT_00, d_chiT_00_dlnd, d_chiT_00_dlnT, &
            chiT_m1, d_chiT_m1_dlnd, d_chiT_m1_dlnT, &
            Cp_00, d_Cp_00_dlnd, d_Cp_00_dlnT, &
            Cp_m1, d_Cp_m1_dlnd, d_Cp_m1_dlnT, &
            opacity_00, d_opacity_00_dlnd, d_opacity_00_dlnT, &
            opacity_m1, d_opacity_m1_dlnd, d_opacity_m1_dlnT, &
            grada_00, d_grada_00_dlnd, d_grada_00_dlnT, &
            grada_m1, d_grada_m1_dlnd, d_grada_m1_dlnT, &
            gradr_factor, gradL_composition_term, &
            alpha_semiconvection, thermohaline_coeff, mixing_length_alpha, &
            Henyey_y_param, Henyey_nu_param, &
            prev_conv_vel, max_conv_vel, g_theta, dt, tau, remove_small_D_limit, &
            gradT_smooth_low, gradT_smooth_high, normal_mlt_gradT_factor
         logical, intent(in) :: alt_scale_height
         character (len=*), intent(in) :: thermohaline_option, MLT_option, semiconvection_option
         integer, intent(in) :: dominant_iso_for_thermohaline
         logical, intent(in) :: just_gradr, smooth_gradT
         integer, intent(out) :: mixing_type
         real(dp), intent(inout) :: mlt_basics(:) ! (num_mlt_results)
         real(dp), intent(inout), pointer :: mlt_partials1(:) ! =(num_mlt_partials, num_mlt_results)
         integer, intent(out) :: ierr

!---
         real(dp) :: q,Prcxh,grrad,grt,vt
         integer :: i, j
!-----

         ierr = 0
         call star_mlt_eval(  &
            id, k, cgrav, m, r, L, X, &            
            T_face, rho_face, P_face, &
            chiRho_face, chiT_face, &
            Cp_face, opacity_face, grada_face, &            
            alfa, beta, & ! f_face = alfa*f_00 + beta*f_m1
            T_00, T_m1, rho_00, rho_m1, P_00, P_m1, &
            chiRho_00, d_chiRho_00_dlnd, d_chiRho_00_dlnT, &
            chiRho_m1, d_chiRho_m1_dlnd, d_chiRho_m1_dlnT, &
            chiT_00, d_chiT_00_dlnd, d_chiT_00_dlnT, &
            chiT_m1, d_chiT_m1_dlnd, d_chiT_m1_dlnT, &
            Cp_00, d_Cp_00_dlnd, d_Cp_00_dlnT, &
            Cp_m1, d_Cp_m1_dlnd, d_Cp_m1_dlnT, &
            opacity_00, d_opacity_00_dlnd, d_opacity_00_dlnT, &
            opacity_m1, d_opacity_m1_dlnd, d_opacity_m1_dlnT, &
            grada_00, d_grada_00_dlnd, d_grada_00_dlnT, &
            grada_m1, d_grada_m1_dlnd, d_grada_m1_dlnT, &            
            gradr_factor, gradL_composition_term, &
            alpha_semiconvection, semiconvection_option, &
            thermohaline_coeff, thermohaline_option, &
            dominant_iso_for_thermohaline, &
            mixing_length_alpha, alt_scale_height, remove_small_D_limit, &
            MLT_option, Henyey_y_param, Henyey_nu_param, &
            gradT_smooth_low, gradT_smooth_high, smooth_gradT, &
            normal_mlt_gradT_factor, &
            prev_conv_vel, max_conv_vel, g_theta, dt, tau, just_gradr, &
            mixing_type, mlt_basics, mlt_partials1, ierr)

! add

!---
             q=chiT_face/chiRho_face
!            Prcxh=1d0/3d0*crad*T*T*T*T
!            grrad= getcxh_Paczynski_gradr(P,opacity,L,m,cgrav,Prcxh,tau,T,r,rho)

            
            grrad=mlt_basics(mlt_gradr)
            !print *,'call kettc!'

            call kettc(rho_face,T_face,r,m,P_face,Cp_face,q,opacity_face,grada_face,grrad,grt,vt,1d0)
            mlt_basics(mlt_gradT)=grt


!add

      end subroutine my_mlt



    
      subroutine my_stc_mlt( &
            id, k, cgrav, m, r, L, X, &            
            T_face, rho_face, P_face, &
            chiRho_face, chiT_face, &
            Cp_face, opacity_face, grada_face, &            
            alfa, beta, & ! f_face = alfa*f_00 + beta*f_m1
            T_00, T_m1, rho_00, rho_m1, P_00, P_m1, &
            chiRho_00, d_chiRho_00_dlnd, d_chiRho_00_dlnT, &
            chiRho_m1, d_chiRho_m1_dlnd, d_chiRho_m1_dlnT, &
            chiT_00, d_chiT_00_dlnd, d_chiT_00_dlnT, &
            chiT_m1, d_chiT_m1_dlnd, d_chiT_m1_dlnT, &
            Cp_00, d_Cp_00_dlnd, d_Cp_00_dlnT, &
            Cp_m1, d_Cp_m1_dlnd, d_Cp_m1_dlnT, &
            opacity_00, d_opacity_00_dlnd, d_opacity_00_dlnT, &
            opacity_m1, d_opacity_m1_dlnd, d_opacity_m1_dlnT, &
            grada_00, d_grada_00_dlnd, d_grada_00_dlnT, &
            grada_m1, d_grada_m1_dlnd, d_grada_m1_dlnT, &            
            gradr_factor, gradL_composition_term, &
            alpha_semiconvection, semiconvection_option, &
            thermohaline_coeff, thermohaline_option, &
            dominant_iso_for_thermohaline, &
            mixing_length_alpha, alt_scale_height, remove_small_D_limit, &
            MLT_option, Henyey_y_param, Henyey_nu_param, &
            gradT_smooth_low, gradT_smooth_high, smooth_gradT, &
            normal_mlt_gradT_factor, &
            prev_conv_vel, max_conv_vel, g_theta, dt, tau, just_gradr, &
            mixing_type, mlt_basics, mlt_partials1, ierr)
         !use star_lib, only: star_mlt_eval
!---
         use star_lib, only: star_ptr
         type(star_info), pointer :: s
!---
         integer, intent(in) :: id ! id for star         
         integer, intent(in) :: k ! cell number or 0 if not for a particular cell         
         real(dp), intent(in) :: &
            cgrav, m, r, L, X, &            
            T_face, rho_face, P_face, &
            chiRho_face, chiT_face, &
            Cp_face, opacity_face, grada_face, &            
            alfa, beta, &
            T_00, T_m1, rho_00, rho_m1, P_00, P_m1, &
            chiRho_00, d_chiRho_00_dlnd, d_chiRho_00_dlnT, &
            chiRho_m1, d_chiRho_m1_dlnd, d_chiRho_m1_dlnT, &
            chiT_00, d_chiT_00_dlnd, d_chiT_00_dlnT, &
            chiT_m1, d_chiT_m1_dlnd, d_chiT_m1_dlnT, &
            Cp_00, d_Cp_00_dlnd, d_Cp_00_dlnT, &
            Cp_m1, d_Cp_m1_dlnd, d_Cp_m1_dlnT, &
            opacity_00, d_opacity_00_dlnd, d_opacity_00_dlnT, &
            opacity_m1, d_opacity_m1_dlnd, d_opacity_m1_dlnT, &
            grada_00, d_grada_00_dlnd, d_grada_00_dlnT, &
            grada_m1, d_grada_m1_dlnd, d_grada_m1_dlnT, &
            gradr_factor, gradL_composition_term, &
            alpha_semiconvection, thermohaline_coeff, mixing_length_alpha, &
            Henyey_y_param, Henyey_nu_param, &
            prev_conv_vel, max_conv_vel, g_theta, dt, tau, remove_small_D_limit, &
            gradT_smooth_low, gradT_smooth_high, normal_mlt_gradT_factor
         logical, intent(in) :: alt_scale_height
         character (len=*), intent(in) :: thermohaline_option, MLT_option, semiconvection_option
         integer, intent(in) :: dominant_iso_for_thermohaline
         logical, intent(in) :: just_gradr, smooth_gradT
         integer, intent(out) :: mixing_type
         real(dp), intent(inout) :: mlt_basics(:) ! (num_mlt_results)
         real(dp), intent(inout), pointer :: mlt_partials1(:) ! =(num_mlt_partials, num_mlt_results)
         integer, intent(out) :: ierr
!---
         real(dp) :: q,Prcxh,grrad,grt,vt
         integer :: i, j

         ierr = 0
         call star_ptr(id, s, ierr)
         if(ierr /=0) return
!---
            q=chiT_face/chiRho_face
!            Prcxh=1d0/3d0*crad*T*T*T*T
!            grrad= getcxh_Paczynski_gradr(P,opacity,L,m,cgrav,Prcxh,tau,T,r,rho)
!
         !call star_mlt_eval(  &
         !   id, k, cgrav, m, r, L, X, &            
         !   T_face, rho_face, P_face, &
         !   chiRho_face, chiT_face, &
         !   Cp_face, opacity_face, grada_face, &            
         !   alfa, beta, & ! f_face = alfa*f_00 + beta*f_m1
         !   T_00, T_m1, rho_00, rho_m1, P_00, P_m1, &
         !   chiRho_00, d_chiRho_00_dlnd, d_chiRho_00_dlnT, &
         !   chiRho_m1, d_chiRho_m1_dlnd, d_chiRho_m1_dlnT, &
         !   chiT_00, d_chiT_00_dlnd, d_chiT_00_dlnT, &
         !   chiT_m1, d_chiT_m1_dlnd, d_chiT_m1_dlnT, &
         !   Cp_00, d_Cp_00_dlnd, d_Cp_00_dlnT, &
         !   Cp_m1, d_Cp_m1_dlnd, d_Cp_m1_dlnT, &
         !   opacity_00, d_opacity_00_dlnd, d_opacity_00_dlnT, &
         !   opacity_m1, d_opacity_m1_dlnd, d_opacity_m1_dlnT, &
         !   grada_00, d_grada_00_dlnd, d_grada_00_dlnT, &
         !   grada_m1, d_grada_m1_dlnd, d_grada_m1_dlnT, &            
         !   gradr_factor, gradL_composition_term, &
         !   alpha_semiconvection, semiconvection_option, &
         !   thermohaline_coeff, thermohaline_option, &
         !   dominant_iso_for_thermohaline, &
         !   mixing_length_alpha, alt_scale_height, remove_small_D_limit, &
         !   MLT_option, Henyey_y_param, Henyey_nu_param, &
         !   gradT_smooth_low, gradT_smooth_high, smooth_gradT, &
         !   prev_conv_vel, max_conv_vel, g_theta, dt, tau, just_gradr, &
         !   mixing_type, mlt_basics, mlt_partials1, ierr)

!cxh
!            grrad=mlt_basics(mlt_gradr)
!cxh
            call kettc(rho_face,T_face,r,m,P_face,Cp_face,q,opacity_face,grada_face,grrad,grt,vt,1d0)
            mlt_basics(mlt_gradT)=grt
        end subroutine my_stc_mlt

        real(dp) function getcxh_Paczynski_gradr( &
            P,opacity,L,m,cgrav,Pr,tau,T,r,rho)
         use crlibm_lib, only: pow_cr
         real(dp), intent(in) :: P,opacity,L,m,cgrav,Pr,tau,T,r,rho
         real(dp) :: dilution_factor, s, f
         getcxh_Paczynski_gradr = P*opacity*L / (16*pi*clight*m*cgrav*Pr)
         if (tau < 2d0/3d0) then ! B. Paczynski, 1969, Acta Astr., vol. 19, 1., eqn 14.
            s = (2*crad*T*T*T*sqrt(r))/(3*cgrav*m*rho)*pow_cr(L/(8*pi*boltz_sigma), 0.25d0) ! eqn 15
            f = 1 - 1.5d0*tau ! Paczynski, 1969, eqn 8
            dilution_factor = (1 + f*s*(4*pi*cgrav*clight*m)/(opacity*L))/(1 + f*s)
            getcxh_Paczynski_gradr = getcxh_Paczynski_gradr*dilution_factor
         end if
      end function getcxh_Paczynski_gradr
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine null_other_D_mix(id, ierr)
         integer, intent(in) :: id
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         integer k,mz,m,i,j,ktop,mp2,mp17
         character stenv*30,stcor*30,stname*30
         integer,parameter:: fcxh01=2016,fcxh02=2017,fcxh03=2018
         real*8 fthe4
!         ftlxyc(13,:)=1d-99 
!         ftlxyc(14,:)=1d-99
!         ftlxyc(15,:)=1d-99
!---
        
         ierr = 0
         call star_ptr(id, s, ierr)
         if(ierr /=0) return
!---       
        m=0
        do k = s% nz,1,-1
        m=m+1
        ftlxyc(1,m)= s% q(k)*s% star_mass 
        ftlxyc(2,m)= s% r(k)      
        ftlxyc(3,m)= s% rho(k)
        ftlxyc(4,m)= s% t(k)
        ftlxyc(5,m)= s% p(k)
        ftlxyc(6,m)= s% L(k)
        ftlxyc(7,m)= s% chiRho(k) 
        ftlxyc(8,m)= s% chiT(k)
        ftlxyc(9,m)= s% gradr(k)
        ftlxyc(10,m)=s% grada(k)
        ftlxyc(11,m)=s% gradT(k)
        ftlxyc(12,m)=s% q(k)
        enddo
!---
      mp2=2
      if((ftlxyc(9,2)-ftlxyc(10,2)).lt.0d0)goto 2012
      do i=1, s% nz  
      if((ftlxyc(9,i)-ftlxyc(10,i)).lt.0d0)then
      mp2=i
      goto 2012
      endif    
      enddo
2012  continue 
!
      ktop=s% nz
      do i=s% nz-1,1,-1
      if((ftlxyc(9,i)-ftlxyc(10,i)).gt.0d0)then
      ktop=i
      goto 2013
      endif
      enddo
2013  continue
      
      mp17=s% nz
      do i=ktop-1,1,-1
      if((ftlxyc(9,i)-ftlxyc(10,i)).lt.0d0)then
      mp17=i
      goto 2014
      endif
      enddo
2014  continue
!Envelope overshoot
         fthe4=1d0-1.5d0*s% initial_z
         if(kcxh.eq.100)then
         !if(ftlxyc(4,mp17).gt.5d3)then
         if(ftlxyc(4,mp17).gt.6d5.and.s% center_h1.lt.1d-10)then
         call KOMEGA_env(s% nz)        
         do k= s% nz,1,-1
         if(k.ge.klim.and.ftlxyc(13,k).gt.1d0.and.ftlxyc(13,k).lt.1d20)s% D_mix(s% nz-k+1)=ftlxyc(13,k)
         enddo
!---writeSTenv
         !write (stenv,"('STE',i6.6,'.txt')") s% model_number
         !open(fcxh01,file=stenv,position='append')
         !do k= klim,s% nz
         !write(fcxh01,"(16(1pes27.16e3,1x))")s% q(k),s% p(k),s% D_mix(k)
         !write(fcxh01,"(16(1pes27.16e3,1x))")(ftlxyc(j,k),j=1,15),s% D_mix(k)
         !enddo
         !close(fcxh01)
         endif
          
!core overshoot
         if(mp2.gt.50.and.s% conv_bdy_q(1).gt.s% min_overshoot_q.and.s% star_mass.gt.1.2d0)then
         if(s% center_h1.gt.5d-4.or.(s% center_h1.lt.1d-10.and.s% center_he4.lt.fthe4.and.s% center_he4.gt.5d-4))then
         call KOMEGA_cor(s% nz)
         do k=1, s% nz
         if(k.le.klim)s% D_mix(s% nz-k+1)=ftlxyc(13,k)
         enddo
!---writeSTcor
         !write (stcor,"('STC',i6.6,'.txt')") s% model_number
         !open(fcxh02,file=stcor,position='append')
         !do k=1,klim
         !write(fcxh01,"(16(1pes27.16e3,1x))")s% q(k),s% p(k),s% D_mix(k)
         !write(fcxh02,"(16(1pes27.16e3,1x))")(ftlxyc(j,k),j=1,12),(st(j,k),j=1,3)
         !enddo
         !close(fcxh02)
         endif
         endif
!---
         !write (stname,"('ST',i6.6,'.txt')") s% model_number
         !open(fcxh01,file=stname,position='append')
         !do k= 1,s% nz
         !do k=1,klim
         !do k= klim,s% nz
         !write(fcxh01,"(16(1pes27.16e3,1x))")s% q(k),s% p(k),s% D_mix(k)
         !write(fcxh01,"(16(1pes27.16e3,1x))")(ftlxyc(j,k),j=1,12),(st(j,k),j=1,3)
         !enddo
         !close(fcxh01)
         endif
      end subroutine null_other_D_mix
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       integer function extras_startup(id, restart, ierr)        
	 integer, intent(in) :: id
         logical, intent(in) :: restart
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         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)
         else ! it is a restart
            call unpack_extra_info(s)
         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
         real(dp) :: 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
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         extras_check_model = keep_going         
!cxh
        kcxh=100
        if(s% center_h1.lt.1d-8.and.s% center_he4.lt.1d-4)s% varcontrol_target = 1d-2
        if(s% power_h_burn/s% photosphere_L.lt.s% Lnuc_div_L_zams_limit.and.s% center_h1.gt.0.5d0)then
        kcxh=10
        endif
        
        !s% profile_interval = 10000
        !if(s% model_number.gt.19950.and.s% model_number.lt.20100)then
        !s% profile_interval = 1
      !endif
        !s% xtra_coef_czb_full_on= 0.965d0
        !s% xtra_coef_czb_full_off= 0.975d0
        !s% xtra_coef_a_u_hb_czb= 0.01d0
        !s% xtra_coef_b_u_hb_czb= 0.01d0
        !s% xtra_dist_a_u_hb_czb= 0.35d0
        !s% xtra_dist_b_u_hb_czb= 0.05d0
        !s% xtra_coef_a_u_heb_czb= 0.01d0
        !s% xtra_coef_b_u_heb_czb= 0.01d0
        !s% xtra_dist_a_u_heb_czb= 0.35d0
        !s% xtra_dist_b_u_heb_czb= 0.05d0 
         if(s% center_h1.lt.2d-8)then
         !s% log_g_lower_limit = 3.0d0
         !s% Teff_lower_limit = 7000d0
         !s% xa_central_lower_limit_species(1)='h1'
         !s% xa_central_lower_limit(1)=1d-2
         !s% xa_central_lower_limit_species(1)='he4'
         !s% xa_central_lower_limit(1)=0.3d0
         endif
      end function extras_check_model


      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 = 0
      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
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
      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 = 2   !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

         if(n /= 2)stop 'data_for_extra_profile_columns wrong !'
         names(1)='kconv'
         names(2)='omega'
         do k=1, nz
            vals(nz-k+1,1)=log10(ftlxyc(14,k))
            vals(nz-k+1,2)=log10(ftlxyc(15,k))
         enddo
      end subroutine data_for_extra_profile_columns
      

      ! returns either keep_going or terminate.
      integer function extras_finish_step(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
         extras_finish_step = keep_going
         call store_extra_info(s)
      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    
         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)
            real(dp) :: 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
      
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
!ly 8/9/15 Interface of the k-omega model for convective core  
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
      subroutine KOMEGA_cor(nh)
      use define_convection_variables
      real,parameter:: sunm=1.9892d33,sunr=6.9598d10,sunl=3.8418d33
      real,parameter:: gconst=6.67428d-8, pi=3.141592653589793d0
      real,parameter::fln10=2.302585092994046d0
      save
!        u     = unknowns in the k-omega model
!        nv    = number of unknows, i.e. k and omega
!        xr    = the independent variable of differential equations
!        cc    = variables of the stellar structure
       REAL*8  u(nv,ngp),cc(5,JMAX),xr(JMAX)
!
      INTEGER i,j,k,mp2,jl,nh,ii
      real*8  uj(2),sj(5),pl(2),pk(2)
      REAL*8  xrlim,fmr,fl,r1,r2,rho,T,P,gr,delta,hp 
      REAL*8  dxj,h,xrj,yj,zj,A1,hh,hh1,cx,cy,cz
!   find the boundary of the convective core, and add an overshooting region
      mp2=2
      do i=1, nh  
      if((ftlxyc(9,i)-ftlxyc(10,i)).lt.0d0)then
      mp2=i
      goto 2012
      endif    
      enddo
2012  continue
      xrlim=log(ftlxyc(5,mp2))-0.3d0
      Hcore=ftlxyc(2,mp2)
!   calculate coefficients of differential equations of the k-omega model
      maxr=0
      do k=1,nh
      if(log(ftlxyc(5,k)).lt.xrlim)then
         klim=k-1
         goto 1
      endif
      maxr=maxr+1
      fmr=ftlxyc(1,k)*sunm 
      fl =ftlxyc(6,k) 
      r2 =ftlxyc(2,k)**2d0
      rho=ftlxyc(3,k)
      T  =ftlxyc(4,k)
      P  =ftlxyc(5,k)
      gr =gconst*fmr/r2
      hp =P/rho/gr
      delta=ftlxyc(8,k)/ftlxyc(7,k)
      cc(1,maxr)=(4d0*pi*rho*r2)**2
      cc(2,maxr)=1d0/(1d0-ftlxyc(12,k))/fmr    
      cc(3,maxr)=alf*Hcore 
      cc(4,maxr)=gr*delta/hp*(ftlxyc(9,k)-ftlxyc(10,k))
      cc(5,maxr)=4d0*pi*gconst*rho*fmr*delta*ftlxyc(9,k)/ftlxyc(10,k)/fl
      xr(  maxr)=log(ftlxyc(12,k)/(1d0-ftlxyc(12,k)))
      enddo
    1 continue
!   Solve the k-omega model
      call getvar_cor(cc,xr,maxr)
      call iniguess_cor(u,ngp,nv)
      call nsnde(u,ngp,nv,99)
!   interpolate the turbulent diffusivity on the grid of stellar model
      h=dxr/(ngp-1d0)
      pk(:)=log(u(:,ngp-30)/u(:,ngp-10))/(20.*h)
      print*,' --- pk = ',pk(1),pk(2)
      do k=1,klim
!   Locate position xrj, and interpolate turbulent solution at point j.
         xrj=log(ftlxyc(12,k)/(1d0-ftlxyc(12,k)))
         dxj=(xrj-xr(1))/h
         jl =int(dxj)
         dxj=dxj-jl
         kl=jl+1
         if(kl.le.1024)then
         uj(:)= u(:,kl)+dxj*(u(:,kl+1)- u(:,kl))                  
         sj(:)=ss(:,kl)+dxj*(ss(:,kl+1)-ss(:,kl))
         else
         uj(:)= u(:,kl)!+dxj*(u(:,kl+1)- u(:,kl))                  
         sj(:)=ss(:,kl)!+dxj*(ss(:,kl+1)-ss(:,kl))
         endif
!   evaluate the turbulent diffusivity using solutions of k-omega model
         zj =sj(3)**2*sj(5)*uj(2)
         yj =uj(2)/uj(1)/sj(5)
         cq =0.5d0+cl2/yj-(cl2/yj)**2*(1d0-1d0/exp(1d0/cl2*yj))
         yj5=ch*cq/yj**ypm
         A1=cta*ct*sj(4)/uj(2)**2+yj5-yj-1d0
         hh=(A1+sqrt(A1**2+4d0*yj5*(1d0+yj)))/(1d0+yj)/2d0
         cx=1d0+yj-ct*cta*sj(4)/uj(2)**2/(1d0+hh)
         ftlxyc(13,k)=cx3*uj(1)/uj(2)/cx 
         ftlxyc(14,k)=uj(2)
         ftlxyc(15,k)=uj(1)
      enddo
      end subroutine KOMEGA_cor
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
!ly 8/9/15 Interpolate coefficients of differential equations   
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
      subroutine getvar_cor(cc,xr,maxr)
      use define_convection_variables
      implicit none
      INTEGER maxr
      REAL*8  cc(5,JMAX),xr(JMAX)
!
      INTEGER j,k
      REAL*8  h,xrj
      dxr=xr(maxr)-xr(1)
      h=dxr/(ngp-1d0)
      k=1
      do 12 j=1,ngp
         xrj=xr(1)+(j-1)*h
11       if(k.ge.maxr)stop ' -- k > maxr in getvar !'
         if(xrj-xr(k).lt.0d0)stop ' -- wrong position in getvar !'
         if(xrj-xr(k+1).gt.1d-14)then
            k=k+1
            goto 11
         endif
         ss(:,j)=cc(:,k)+(cc(:,k+1)-cc(:,k))*(xrj-xr(k))/(xr(k+1)-xr(k))
12    enddo
      return
      end subroutine getvar_cor
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
!ly 8/9/15 set up an initial solution of the k-omega model   
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
      SUBROUTINE iniguess_cor(u,n,np)
      use define_convection_variables
      implicit none
      INTEGER n,np
      REAL*8  u(np,n)
!
      INTEGER j,k,jlim1,jlim2
      REAL*8  h,xlim1,xlim2,sj(5),uj(2)
      h=dxr/(n-1d0)
      xlim1=0d0
      xlim2=0d0
      do j=1,n
         sj(:)=ss(:,j)
         if(sj(4).le.0d0)then
           jlim2=j
           goto 1
         endif
         CALL IMLT(uj,sj)
         u(:,j)=uj(:)
      enddo
1     jlim2=jlim2-2 !INT(0.0/h)
      do j=jlim2+1,n
         u(1,j)=1d0*u(1,jlim2)/exp(h*(j-jlim2)*56.)
         u(2,j)=1.0*u(2,jlim2)/exp(h*(j-jlim2)*28.)
      enddo
      return
      end SUBROUTINE iniguess_cor
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
!ly 8/9/15 calculate the local solutions of the k-omega model   
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
      subroutine IMLT(uj,cj)
      use define_convection_variables
      implicit none
      REAL*8  uj(2),cj(5)
!
      REAL*8  xx,xx2,hh,yy
!
      hh=cj(3)**2
      if(cj(4).gt.0d0)then
        xx=1d0/hh/cj(5)/sqrt(cj(4))
      else
        xx=1d0/hh/cj(5)/sqrt(1d-12)
      endif
      xx2=xx**2
      yy=1d0/((xx2/cc2)+(xx2/cc1)**ynl)
      uj(2)=yy/cj(5)/hh
      uj(1)=yy/cj(5)*uj(2)
      return
      end subroutine IMLT
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
!ly 8/9/15 calculate the diffusivities of the k-omega model   
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
      subroutine dfco(j,n,np,DY,u)
      use define_convection_variables
      implicit none
      INTEGER j,n,np
      REAL*8  DY(np,np+1),u(np,n)
!
      REAL*8  dxr2,BCmu
      DY=0
      dxr2=2d0*dxr*dxr
      BCmu=cmu*ss(1,j)*ss(2,j)/dxr2
!   Define turbulent diffusivities.
      DY(1,3)= BCmu*u(1,j)/u(2,j)
      DY(2,3)= BCmu*u(1,j)/u(2,j)/cgma
      DY(1,1)= DY(1,3)/u(1,j)
      DY(1,2)=-DY(1,3)/u(2,j)
      DY(2,1)= DY(2,3)/u(1,j)
      DY(2,2)=-DY(2,3)/u(2,j)
      return
      end subroutine dfco
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
!ly 8/9/15 calculate the coefficients of difference equations   
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
      subroutine dcff(j,n,np,DK,F,u)
      use define_convection_variables
      implicit none
      INTEGER j,n,np
      REAL*8  DK(np,np+1,-1:1),F(np,np+1),u(np,n)
!
      REAL*8  FX1,FX2,AA,AA1,AA2,BB,DD,yy,yy1,yy2,yy5,cq,cq1,cq2,qq
      REAL*8  CmD,CmE,DG2,CLD2,CmEE,AFM,AFM1,AFM2,AFN,AFN1,AFN2,qq1,qq2
!   Evaluate the coefficients of difference equations.
      call dfco(j-1,n,np,DK(:,:,-1),u)
      call dfco(j+1,n,np,DK(:,:,1),u)
      call dfco(j,n,np,DK(:,:,0),u)
      DK=DK*ss(2,j)
      F=0d0
!   Evaluate source functions
      CmE=ct*ss(4,j)                                ! ct*E
      CmD=cta*ct*ss(4,j)                            ! cta*ct*E
      DG2=ss(3,j)**2*ss(5,j)                        ! H**2*D
      CLD2=ss(5,j)/DG2                              ! 1/H**2
!
      FX1=u(1,j)					                      ! k
      FX2=u(2,j)                                    ! omega
      yy=FX2/FX1/ss(5,j)                            ! y**-1
      yy1=-yy/FX1
      yy2= yy/FX2
      yy5=cl2/yy
      if(yy5.gt.20.)then
      cq =1d0-1d0/(6d0*yy5+1.5d0)
      cq1=   -0.16666667d0/(yy5+0.25)**2*yy5/yy*yy1
      cq2=   -0.16666667d0/(yy5+0.25)**2*yy5/yy*yy2
      else
      cq =0.5d0+yy5-yy5**2*(1d0-exp(-1d0/yy5))
      cq1=(-yy5+2.*yy5**2*(1d0-exp(-1d0/yy5)))/yy*yy1-yy5**2/exp(1d0/yy5)/cl2*yy1
      cq2=(-yy5+2.*yy5**2*(1d0-exp(-1d0/yy5)))/yy*yy2-yy5**2/exp(1d0/yy5)/cl2*yy2
      endif
      yy5=ch*cq/yy**ypm
      qq=(cl2*DG2*FX2)**0.66666667   
      qq1=0.
      qq2=0.67*qq/FX2
!
      AFN=1d0+yy+yy5-CmD/FX2**2
      if(AFN.lt.0d0)AFN=0.01d99
      AFN1=yy1-ypm*yy5/yy*yy1                +ch/yy**ypm*cq1
      AFN2=yy2-ypm*yy5/yy*yy2+2d0*CmD/FX2**3 +ch/yy**ypm*cq2
      AA = 4.0d0*yy5*CmD/FX2**2/AFN**2
      AA1=-ypm*AA/yy*yy1-2.0*AA/AFN*AFN1
      AA2=-ypm*AA/yy*yy2-2.0*AA/AFN*AFN2-2.0*AA/FX2
      AA1=-1d0/(1d0+sqrt(1d0+AA))**2/sqrt(1d0+AA)*AA1
      AA2=-1d0/(1d0+sqrt(1d0+AA))**2/sqrt(1d0+AA)*AA2
      AA = 2d0/(1d0+sqrt(1d0+AA))
      AFM=(1d0+qq**2-cs*qq)/(1d0+qq**2)
      AFM1=(2.*qq-cs-AFM*2.*qq)*qq1/(1d0+qq**2)
      AFM2=(2.*qq-cs-AFM*2.*qq)*qq2/(1d0+qq**2)
      BB=AFM*FX1*FX2
      DD=AA*CmE/AFN*FX1/FX2
!   Define sorce function S_k.
      F(1,1)=BB/FX1+AFM1*FX1*FX2-AA1*CmE/AFN*FX1/FX2-DD/FX1+DD/AFN*AFN1
      F(1,2)=BB/FX2+AFM2*FX1*FX2-AA2*CmE/AFN*FX1/FX2+DD/FX2+DD/AFN*AFN2
      F(1,3)=BB-DD
!   Define sorce function S_omega.
      AFM=CLD2*FX1
      AFM1=AFM/FX1
      AFM2=0. 
      F(2,1)=-AFM1
      F(2,2)=-AFM2+2d0*FX2
      F(2,3)=-AFM +FX2**2
      return
      end subroutine dcff
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
!ly 8/9/15 interface to difference equations solver  
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
      SUBROUTINE nsnde(u,n,np,maxcyc)
      implicit none
      INTEGER n,np,maxcyc
      REAL*8  u(np,n)
!
      INTEGER i,j,jcycle,k,iumax(np),iemax(np)
      REAL*8  du(np,n),dumax(np),demax(np)
      REAL*8  dutol,detol,xx,fac
      PARAMETER (dutol=1.d-5,detol=1.d-5)
1     FORMAT(I5,4(I6,F17.10,I8,E14.4))
      DO 20 jcycle=1,maxcyc
         CALL nsolver(u,du,demax,iemax,n,np)
         DO 13 i=1,np
            dumax(i)=0d0
            DO 11 j=1,n
               xx=du(i,j)/u(i,j)
               IF (ABS(xx).GT.ABS(dumax(i))) THEN
                  dumax(i)=xx
                  iumax(i)=j
               ENDIF
!               if(xx.gt.9.9)du(i,j)=0.2*u(i,j)
               if(xx.gt.0.9)du(i,j)=0.9*u(i,j)
               if(xx.lt.-.5)du(i,j)=-.5*u(i,j)
               if(xx.lt.-.9)du(i,j)=-.1*u(i,j)
11          ENDDO
            fac=1d0
            DO 12 j=1,n
               u(i,j)=u(i,j)+du(i,j)*fac
12          ENDDO
13       ENDDO
         IF(ABS(dumax(1)).LT.dutol)GOTO 30
20    ENDDO
      print *, ' --- Newton iterations for k-omega solution diverge !!!'
30    RETURN
      END SUBROUTINE nsnde
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
!ly 8/9/15 Newton solver for nonlinear diffussion equations:
!          d/dx ( D du/dx ) - F(u) = rhs.
!   The current value of the solution u(1:np,1:n) is updated, using the
!   right-hand-side function rhs(1:np,1:n).
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
      SUBROUTINE nsolver(u,du,demax,iemax,n,np)
      implicit none
      INTEGER n,np,iemax(np)
      REAL*8  u(np,n),du(np,n),demax(np)
!
      INTEGER i,j,k,l,indx(np),np1
      REAL*8  DK(np,np+1,-1:1),F(np,np+1)
      REAL*8  G(np,np,-1:1),P(np,np),Q(np,np+1),S(np,np+1,n)
      REAL*8  h,h2i,dudx,du1dx1,dDdx,dD1dx1,par
      h=1d0/(n-1d0)
      h2i=1d0/h**2
      demax=0d0
      np1=np+1
!  Inner boundary conditions: X=X0
      S(:,:,1)=0d0
!  Inner boundary conditions: dX/dQ=0, when convection develops at the center
      if(u(1,1).gt.1d6)then
         DO i=1,np
            S(i,i,1)=1d0
            S(i,np1,1)=u(i,2)-u(i,1)
         ENDDO
      endif
!  Loop for the recursive relation
      DO 50 k=2,n-1
!  Coeffients for the differential equations
         CALL dcff(k,n,np,DK,F,u)
!  Coefficients of the correction equations
         DO 34 i=1,np
            dudx  =(u(i,k)-u(i,k-1))*h2i
            du1dx1=(u(i,k+1)-u(i,k))*h2i
            dDdx  =(DK(i,np1,0)+DK(i,np1,-1))*h2i
            dD1dx1=(DK(i,np1,0)+DK(i,np1, 1))*h2i
            DO 31 j=1,np
               G(i,j,-1)=-DK(i,j,-1)*dudx
               G(i,j, 1)= DK(i,j, 1)*du1dx1
               G(i,j, 0)= DK(i,j,0)*(du1dx1-dudx)-F(i,j)
31          ENDDO
            G(i,i,-1)=G(i,i,-1)+dDdx
            G(i,i, 1)=G(i,i, 1)+dD1dx1
            G(i,i, 0)=G(i,i, 0)-dDdx-dD1dx1
       du(i,k)=F(i,3)+dudx*(DK(i,np1,0)+DK(i,np1,-1))-du1dx1*(DK(i,np1,0)+DK(i,np1, 1))
            IF (ABS(du(i,k)).GT.demax(i)) THEN
               demax(i)=ABS(du(i,k))
               iemax(i)=k
            ENDIF
34       ENDDO
!  Solve for the recursive equations by LU decomposition
         DO 38 i=1,np
            DO 36 j=1,np
               Q(i,j)=-G(i,j,1)
               P(i,j)= G(i,j,0)
               DO 35 l=1,np
                  P(i,j)=P(i,j)+G(i,l,-1)*S(l,j,k-1)
35             ENDDO
36          ENDDO
            Q(i,np1)=du(i,k)
            DO 37 l=1,np
               Q(i,np1)=Q(i,np1)-G(i,l,-1)*S(l,np1,k-1)
37          ENDDO
38       ENDDO
            CALL ludcmp(P,np,np,indx,par)
         DO 39 j=1,np1
            CALL lubksb(P,np,np,indx,Q(1,j))
39       ENDDO
         DO 41 i=1,np
            DO 40 j=1,np1
               S(i,j,k)=Q(i,j)
40          ENDDO
41       ENDDO
50    ENDDO
!  Outer boundary conditions: X=X0
      du(:,n)=0d0
!  Outer boundary conditions: dX/dQ=0
!      DO i=1,np
!         DO j=1,np
!            P(i,j)=S(i,j,n-1)
!         ENDDO
!         P(i,i)=P(i,i)-1d0
!         du(i,n)=u(i,n)-u(i,n-1)-S(i,np+1,n-1)
!      ENDDO
!      CALL ludcmp(P,np,np,indx,par)
!      CALL lubksb(P,np,np,indx,du(1,n))
!  Find the corrections recursively
      DO 73 k=n-1,1,-1
         DO 72 i=1,np
            du(i,k)=S(i,np1,k)
            DO 71 l=1,np
               du(i,k)=du(i,k)+S(i,l,k)*du(l,k+1)
71          ENDDO
72       ENDDO
73    ENDDO
      RETURN
      END SUBROUTINE nsolver
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
      SUBROUTINE LUBKSB(A,N,NP,INDX,B)
!  .. Arguments ..
      INTEGER N,NP,INDX(NP)
      REAL*8 A(NP,NP),B(NP)
!  .. Local variables ..
      INTEGER I,II,J,LL
      REAL*8 SUM
!
      II = 0
      DO 12 I = 1,N
          LL = INDX(I)
          SUM = B(LL)
          B(LL) = B(I)
          IF (II.NE.0) THEN
              DO 11 J = II,I - 1
                  SUM = SUM - A(I,J)*B(J)
   11         CONTINUE
          ELSE IF (SUM.NE.0.) THEN
              II = I
          END IF
          B(I) = SUM
   12 CONTINUE
!
      DO 14 I = N,1,-1
          SUM = B(I)
          IF (I.LT.N) THEN
              DO 13 J = I + 1,N
                  SUM = SUM - A(I,J)*B(J)
   13         CONTINUE
          END IF
          B(I) = SUM/A(I,I)
   14 CONTINUE
      RETURN
      END SUBROUTINE LUBKSB
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
      SUBROUTINE LUDCMP(A,N,NP,INDX,D)
!  .. Arguments ..
      INTEGER N,NP,INDX(NP)
      REAL*8  A(NP,NP),D
!  .. Local variables ..
      INTEGER, parameter :: NMAX=100
      REAL*8,  parameter :: TINY=1d-20
      INTEGER  I,IMAX,J,K
      REAL*8   AAMAX,DUM,SUM,VV(NMAX)
!
      D = 1d0
      DO 12 I = 1,N
          AAMAX = 0d0
          DO 11 J = 1,N
              IF (ABS(A(I,J)).GT.AAMAX) AAMAX = ABS(A(I,J))
   11     CONTINUE
          IF (AAMAX.EQ.0d0) aamax=tiny    ! STOP 'Singular matrix in LUDCMP.'
          VV(I) = 1d0/AAMAX
   12 CONTINUE
!
      DO 19 J = 1,N
          IF (J.GT.1) THEN
              DO 14 I = 1,J - 1
                  SUM = A(I,J)
                  IF (I.GT.1) THEN
                      DO 13 K = 1,I - 1
                          SUM = SUM - A(I,K)*A(K,J)
   13                 CONTINUE
                      A(I,J) = SUM
                  END IF
   14         CONTINUE
          END IF
!
          AAMAX = 0d0
          DO 16 I = J,N
              SUM = A(I,J)
              IF (J.GT.1) THEN
                  DO 15 K = 1,J - 1
                      SUM = SUM - A(I,K)*A(K,J)
   15             CONTINUE
                  A(I,J) = SUM
              END IF
              DUM = VV(I)*ABS(SUM)
              IF (DUM.GE.AAMAX) THEN
                  IMAX = I
                  AAMAX = DUM
              END IF
   16     CONTINUE
!
          IF (J.NE.IMAX) THEN
              DO 17 K = 1,N
                  DUM = A(IMAX,K)
                  A(IMAX,K) = A(J,K)
                  A(J,K) = DUM
   17         CONTINUE
              D = -D
              VV(IMAX) = VV(J)
          END IF
!
          INDX(J) = IMAX
          IF (J.NE.N) THEN
              IF (A(J,J).EQ.0d0) A(J,J) = TINY
              DUM = 1d0/A(J,J)
              DO 18 I = J + 1,N
                  A(I,J) = A(I,J)*DUM
   18         CONTINUE
          END IF
   19 CONTINUE
      IF (A(N,N).EQ.0d0) A(N,N) = TINY
      RETURN
      END SUBROUTINE LUDCMP
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
!
      REAL*8 FUNCTION array_maximum(a,n)
      implicit none
      INTEGER n
      REAL*8  a(n)
!  Return the maximum element of matrix a(1:n).
      INTEGER i
      array_maximum=0d0
      DO i=1,n
         array_maximum=MAX(array_maximum,ABS(a(i)))
      ENDDO
      RETURN
      END
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
!ly 8/9/15 Interface of k-omega model for convective envelope  
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
      subroutine KOMEGA(nh)
      use define_convection_variables
      real,parameter:: sunm=1.9892d33,sunr=6.9598d10,sunl=3.8418d33
      real,parameter:: gconst=6.67428d-8, pi=3.141592653589793d0
      real,parameter::fln10=2.302585092994046d0
      save
!        u     = unknowns in k-omega model
!        nv    = number of unknows, i.e. k and omega
!        xr    = the independent variable of differential equations 
!        cc    = variables of the stellar structure
      REAL*8  u(nv,ngp),cc(5,JMAX),xr(JMAX)
!
      INTEGER i,j,k,mp17,jl,nh,kl,ktop
      real*8  uj(2),sj(5),pl(2),pk(2)
      REAL*8  xrlim,fmr,fl,r2,rho,T,P,gr,delta,hp 
      REAL*8  dxj,h,xrj,yj,zj,A1,hh,hh1,cx,cy,cz
!   find the boundary of the convective envelope, and add an overshooting region
      ktop=nh
      do i=nh-1,1,-1
      if((ftlxyc(9,i)-ftlxyc(10,i)).gt.0d0)then
      ktop=i
      goto 2012
      endif
      enddo
2012  continue
      
      mp17=nh
      do i=ktop-1,1,-1
      if((ftlxyc(9,i)-ftlxyc(10,i)).lt.0d0)then
      mp17=i
      goto 2013
      endif
      enddo
2013  continue
      
      xrlim=log(ftlxyc(5,mp17))+3.0d0
      maxr=0
      do 1 k=1,nh
         if(log(ftlxyc(5,k)).gt.xrlim)then
         klim=k+1
         goto 1
         endif
!   calculate coefficients of the differential equations
      maxr=maxr+1
      fmr=ftlxyc(1,k)*sunm 
      fl =ftlxyc(6,k)
      r2 =ftlxyc(2,k)**2d0
      rho=ftlxyc(3,k)
      T  =ftlxyc(4,k)
      P  =ftlxyc(5,k)
      gr =gconst*fmr/r2
      hp =P/rho/gr
      delta=ftlxyc(8,k)/ftlxyc(7,k)
      cc(1,maxr)=(4d0*pi*rho*r2)**2
      cc(2,maxr)=1d0/hp/(4d0*pi*rho*r2)  
      cc(3,maxr)=afa*hp
      cc(4,maxr)=gr*delta/hp*(ftlxyc(9,k)-ftlxyc(10,k))
      cc(5,maxr)=4d0*pi*gconst*rho*fmr*delta*ftlxyc(9,k)/ftlxyc(10,k)/fl
      xr(  maxr)=log(P)
    1 continue

!      write(*,*)'Env:  ','mp2=',mp17,'klim=',klim
!   Solve the k-omega model
      call getvar(cc,xr,maxr)
      call iniguess(u,ngp,nv)
      call nsnde(u,ngp,nv,99)
!   interpolate the turbulent diffusivity on the grid of stellar model
      h=-dxr/(ngp-1d0)
      pl(:)=log(u(:,    60)/u(:,    20))/(40.*h)
      pk(:)=log(u(:,ngp-60)/u(:,ngp-20))/(40.*h)
      print*,' --- pl = ',pl(1),pl(2)
      print*,' --- pk = ',pk(1),pk(2)
!
      do k=klim,nh
         xrj=log(ftlxyc(5,k))
         dxj=(xr(1)-xrj)/h
!   Locate position xrj, and interpolate turbulent solution at point j.
         jl=int(dxj)
         dxj=dxj-jl
         kl=jl+1
         if(kl.le.1024)then
         uj(:)= u(:,kl)+dxj*(u(:,kl+1)- u(:,kl))                  
         sj(:)=ss(:,kl)+dxj*(ss(:,kl+1)-ss(:,kl))
         else
         uj(:)= u(:,kl)!+dxj*(u(:,kl+1)- u(:,kl))                  
         sj(:)=ss(:,kl)!+dxj*(ss(:,kl+1)-ss(:,kl))
         endif
!   evaluate the turbulent diffusivity using solutions of k-omega model
         zj =sj(3)**2*sj(5)*uj(2)
         yj =uj(2)/uj(1)/sj(5)
         cq =0.5d0+cl2/yj-(cl2/yj)**2*(1d0-1d0/exp(1d0/cl2*yj))
         yj5=ch*cq/yj**ypm
         A1=cta*ct*sj(4)/uj(2)**2+yj5-yj-1d0
         hh=(A1+sqrt(A1**2+4d0*yj5*(1d0+yj)))/(1d0+yj)/2d0
         cx=1d0+yj-ct*cta*sj(4)/uj(2)**2/(1d0+hh)
         ftlxyc(13,k)=cx3*uj(1)/uj(2)/cx 
         ftlxyc(14,k)=uj(2)
         ftlxyc(15,k)=uj(1)
      enddo
      end subroutine KOMEGA
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
!ly 8/9/15 Interpolate coefficients of differential equations   
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
      subroutine getvar(cc,xr,maxr)
      use define_convection_variables
      implicit none
      INTEGER maxr
      REAL*8  cc(5,JMAX),xr(JMAX)
!
      INTEGER j,k
      REAL*8  h,xrj
      dxr=xr(maxr)-xr(1)
      h=dxr/(ngp-1d0)
      k=1
      do 12 j=1,ngp
         xrj=xr(1)+(j-1)*h
11       if(k.ge.maxr)stop ' -- k > maxr in getvar !'
         if(xrj-xr(k).gt.1d-14)stop ' -- wrong position in getvar !'
         if(xrj-xr(k+1).lt.-1d-14)then
            k=k+1
            goto 11
         endif
         ss(:,j)=cc(:,k)+(cc(:,k+1)-cc(:,k))*(xrj-xr(k))/(xr(k+1)-xr(k))
12    enddo
      return
      end
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
!ly 8/9/15 set up an initial solution of the k-omega model   
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
      SUBROUTINE iniguess(u,n,np)
      use define_convection_variables
      implicit none
      INTEGER n,np
      REAL*8  u(np,n)
!
      INTEGER j,k,jlim1,jlim2,ibc
      REAL*8  h,xlim1,xlim2,sj(5),uj(2)
      h=dxr/(n-1d0)
      xlim1=0d0
      ibc=0
      do 2 j=1,n
         sj(:)=ss(:,j)
         if(sj(4).lt.0d0.and.ibc.eq.0)then
           jlim1=j
         else
           ibc=1
           CALL IMLT(uj,sj)
           u(:,j)=uj(:)
           if(uj(1).gt.xlim1)then
             xlim1=uj(1)
             xlim2=uj(2)
             jlim2=j
           endif
         endif
2     enddo
      jlim1=jlim1+3 !INT(0.1/h)
      do 13 j=jlim1-1,1,-1
         u(1,j)=1.0*u(1,jlim1)/exp(h*(j-jlim1)*9.0)
         u(2,j)=1.0*u(2,jlim1)/exp(h*(j-jlim1)*4.0)
13    enddo
!      jlim2=jlim2-INT(0.0/h)
      do 14 j=jlim2+1,n
         u(1,j)=1.0*u(1,jlim2)/exp(h*(jlim2-j)*3.5)
         u(2,j)=1.0*u(2,jlim2)/exp(h*(jlim2-j)*1.5)
14    enddo
      return
      end
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
      subroutine kettc(rh,t,r,fmr,p,cp,q,fkap,grad,grrad,grt,vt,tdc)
      use define_convection_variables
      implicit real*8 (a-h,o-z)
      real,parameter:: sunm=1.9892d33,sunr=6.9598d10,sunl=3.8418d33
      real,parameter:: gconst=6.67428d-8, pi=3.141592653589793d0
      real,parameter:: clight=2.99792458d10,arad=7.5659d-15
      real,parameter::fln10=2.302585092994046d0,sigcst=arad*clight/4d0
      g=gconst*fmr/r**2
      hp=p/rh/g
      cc3=afa*hp
      cc5=3d0/16d0/sigcst*cp*rh**2*fkap/t**3
      gamma=cc5*cc3**2
      if(grrad.gt.grad)then
         v=gamma*sqrt(g*q/hp*(grrad-grad))
      else
	 v=1d-39
      endif
      xx=max(1d0/v,1d-99)
      xx2=xx**2
      yy=1d0/((xx2/cc2)+(xx2/cc1)**ynl)
!   solve the local k-omega equation by Newton iteration method
      do j=1,10
         yy5=cl2*yy
         if(yy5.gt.20.)then
         cq =1d0-1d0/(6d0*yy5+1.5d0)
         cqy=0.16666667/(yy5+0.25d0)**2*cl2
         else
         cq =0.5d0+yy5-yy5**2*(1d0-exp(-1d0/yy5))
         cqy=cl2-2./yy*yy5**2*(1d0-exp(-1d0/yy5))+cl2/exp(1d0/yy5)        
         endif
         yy5=ch*cq*yy**ypm
	      qq =(cl2*yy)**0.66666667d0
         yy3=1d0+qq**2-cs*qq
         cb =(1d0+qq**2)/yy3
         cby=(1.3*qq**2-cb*(1.3*qq**2-0.67*cs*qq))/yy/yy3
         ff1=(1d0+1d0/yy+yy5*(1d0+cta/cb))*yy**2
         ffy=(-1./yy+ypm*yy5*(1d0+cta/cb))*yy+ch*yy**3*(1d0+cta/cb)*cqy&
                       -yy5*cta/cb**2*cby*yy**2+2d0*ff1/yy
         dF=ffy-ct/xx2*cby
         F0=ff1-ct/xx2*(cta+cb)
         yy=yy-F0/dF
         if(abs(F0/dF/yy).lt.1d-10) goto 1
         if(j.gt.5)print*, xx,cc3,cc5
      enddo
      stop ' --- no convergence in KETTC !'
1     continue
      hh=(1d0+cta/cb)*ch*cq*yy**ypm/(1d0+1d0/yy)
      grt=grad+(grrad-grad)/(1d0+hh)*tdc
      return
      end
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
      subroutine KOMEGA_env(nh)
      use define_convection_variables
      real,parameter:: sunm=1.9892d33,sunr=6.9598d10,sunl=3.8418d33
      real,parameter:: gconst=6.67428d-8, pi=3.141592653589793d0
      real,parameter::fln10=2.302585092994046d0
      save
!        u     = unknowns in k-omega model
!        nv    = number of unknows, i.e. k and omega
!        xr    = the independent variable of differential equations 
!        cc    = variables of the stellar structure
      REAL*8  u(nv,ngp),cc(5,JMAX),xr(JMAX)
!
      INTEGER i,j,k,mp17,jl,nh,kl,ktop
      real*8  uj(2),sj(5),pm(2),pn(2)
      REAL*8  xrlim,fmr,fl,r2,rho,T,P,gr,delta,hp 
      REAL*8  dxj,h,xrj,yj,zj,A1,hh,hh1,cx,cy,cz
!   find the boundary of the convective envelope, and add an overshooting region
      ktop=nh
      do i=nh-1,1,-1
      if((ftlxyc(9,i)-ftlxyc(10,i)).gt.0d0)then
      ktop=i
      goto 2012
      endif
      enddo
2012  continue
      
      mp17=nh
      do i=ktop-1,1,-1
         if((ftlxyc(9,i)-ftlxyc(10,i)).lt.0d0)then
         mp17=i
         if((ftlxyc(9,i-1)-ftlxyc(10,i-1)).lt.0d0)goto 2013
      endif
      enddo
2013  continue
      
      xrlim=log(ftlxyc(5,mp17))+3.0d0
      maxr=0
      !Henv1=ftlxyc(2,ktop)
      !Henv2=ftlxyc(2,mp17)
      !Henv=Henv2-Henv1
      do 1 k=1,nh
         if(log(ftlxyc(5,k)).gt.xrlim)then
         klim=k+1
         goto 1
         endif
!   calculate coefficients of the differential equations
      maxr=maxr+1
      fmr=ftlxyc(1,k)*sunm 
      fl =ftlxyc(6,k)
      r2 =ftlxyc(2,k)**2d0
      rho=ftlxyc(3,k)
      T  =ftlxyc(4,k)
      P  =ftlxyc(5,k)
      gr =gconst*fmr/r2
      hp =P/rho/gr
      delta=ftlxyc(8,k)/ftlxyc(7,k)
      cc(1,maxr)=(4d0*pi*rho*r2)**2
      cc(2,maxr)=1d0/hp/(4d0*pi*rho*r2)  
      cc(3,maxr)=afa*hp
      !cc(3,maxr)=afa*Henv
      cc(4,maxr)=gr*delta/hp*(ftlxyc(9,k)-ftlxyc(10,k))
      cc(5,maxr)=4d0*pi*gconst*rho*fmr*delta*ftlxyc(9,k)/ftlxyc(10,k)/fl
      xr(  maxr)=log(P)
    1 continue

!      write(*,*)'Env:  ','mp2=',mp17,'klim=',klim
!   Solve the k-omega model
      call getvar_env(cc,xr,maxr)
      call iniguess_env(u,ngp,nv)
      call nsnde(u,ngp,nv,99)
!   interpolate the turbulent diffusivity on the grid of stellar model
      h=-dxr/(ngp-1d0)
      pm(:)=log(u(:,    60)/u(:,    20))/(40.*h)
      pn(:)=log(u(:,ngp-60)/u(:,ngp-20))/(40.*h)
      print*,' --- pm = ',pm(1),pm(2)
      print*,' --- pn = ',pn(1),pn(2)
!
      do k=klim,nh
         xrj=log(ftlxyc(5,k))
         dxj=(xr(1)-xrj)/h
!   Locate position xrj, and interpolate turbulent solution at point j.
         jl=int(dxj)
         dxj=dxj-jl
         kl=jl+1
         if(kl.le.1024)then
         uj(:)= u(:,kl)+dxj*(u(:,kl+1)- u(:,kl))                  
         sj(:)=ss(:,kl)+dxj*(ss(:,kl+1)-ss(:,kl))
         else
         uj(:)= u(:,kl)!+dxj*(u(:,kl+1)- u(:,kl))                  
         sj(:)=ss(:,kl)!+dxj*(ss(:,kl+1)-ss(:,kl))
         endif
!   evaluate the turbulent diffusivity using solutions of k-omega model
         zj =sj(3)**2*sj(5)*uj(2)
         yj =uj(2)/uj(1)/sj(5)
         cq =0.5d0+cl2/yj-(cl2/yj)**2*(1d0-1d0/exp(1d0/cl2*yj))
         yj5=ch*cq/yj**ypm
         A1=cta*ct*sj(4)/uj(2)**2+yj5-yj-1d0
         hh=(A1+sqrt(A1**2+4d0*yj5*(1d0+yj)))/(1d0+yj)/2d0
         cx=1d0+yj-ct*cta*sj(4)/uj(2)**2/(1d0+hh)
         ftlxyc(13,k)=cx3*uj(1)/uj(2)/cx 
         ftlxyc(14,k)=uj(2)
         ftlxyc(15,k)=uj(1)
      enddo
      end subroutine KOMEGA_env
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
!ly 8/9/15 Interpolate coefficients of differential equations   
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
      subroutine getvar_env(cc,xr,maxr)
      use define_convection_variables
      implicit none
      INTEGER maxr
      REAL*8  cc(5,JMAX),xr(JMAX)
!
      INTEGER j,k
      REAL*8  h,xrj
      dxr=xr(maxr)-xr(1)
      h=dxr/(ngp-1d0)
      k=1
      do 12 j=1,ngp
         xrj=xr(1)+(j-1)*h
11       if(k.ge.maxr)stop ' -- k > maxr in getvar !'
         if(xrj-xr(k).gt.1d-14)stop ' -- wrong position in getvar !'
         if(xrj-xr(k+1).lt.-1d-14)then
            k=k+1
            goto 11
         endif
         ss(:,j)=cc(:,k)+(cc(:,k+1)-cc(:,k))*(xrj-xr(k))/(xr(k+1)-xr(k))
12    enddo
      return
      end
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
!ly 8/9/15 set up an initial solution of the k-omega model   
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
      SUBROUTINE iniguess_env(u,n,np)
      use define_convection_variables
      implicit none
      INTEGER n,np
      REAL*8  u(np,n)
!
      INTEGER j,k,jlim1,jlim2,ibc
      REAL*8  h,xlim1,xlim2,sj(5),uj(2)
      h=dxr/(n-1d0)
      xlim1=0d0
      ibc=0
      do 2 j=1,n
         sj(:)=ss(:,j)
         if(sj(4).lt.0d0.and.ibc.eq.0)then
           jlim1=j
         else
           ibc=1
           CALL IMLT(uj,sj)
           u(:,j)=uj(:)
           if(uj(1).gt.xlim1)then
             xlim1=uj(1)
             xlim2=uj(2)
             jlim2=j
           endif
         endif
2     enddo
      jlim1=jlim1+1 !INT(0.1/h)
      do 13 j=jlim1-1,1,-1
         u(1,j)=1.0*u(1,jlim1)/exp(h*(j-jlim1)*11.0)
         u(2,j)=1.0*u(2,jlim1)/exp(h*(j-jlim1)*5.0)
13    enddo
!      jlim2=jlim2-INT(0.0/h)
      do 14 j=jlim2+1,n
         u(1,j)=1.0*u(1,jlim2)/exp(h*(jlim2-j)*2.2)
         u(2,j)=1.0*u(2,jlim2)/exp(h*(jlim2-j)*0.5)
14    enddo
      return
      end
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++|

