subroutine chiv_ip(RLT,RLN,RLNe,q,kappa,
     .     shat,zth,nbeam,tau,eps,gnu,beta,
     .     g_perp,rmajor,rho_i,v_ti,rho_e,v_te,etg,
     .     RLTcrit,RLTcritz,chi_0,g,chi_i,chi_e,
     .     delGradT,RLTe,Ti_uprat,Te_uprat,ichi,kstep,rhohalf,
     .     chi_ihat,u_ihat,chi_eehat,chi_eihat,u_ehat,gamma)

c
c Calculate Turbulent Chi's and V's from the IFS-PPPL transport model 
c based on gyrofluid/gyrokinetic simulations.  Splits up the heat flux
c into effective convective and diffusive terms for numerical stability
c of transport codes.
c
c The first 4 lines of arguments are the same as for Dorland's original 
c ip_chi subroutine (see ip_chi.f for their full definitions):
c
c The additional arguments are:
c
c     delGradT  (input) relative change in RLT used in calculating 
c               d(heat flux)/d(RLT).  0.03 is a good default value.
c
c     RLTe  (input) R/L_Te, where R = rmajor and 1/L_Te = -dT_e/dr/T_e
c     Ti_uprat  (input)  Ti(k+1)/Ti(k) (Ti ratio across zone)
c     Te_uprat  (input)  Te(k+1)/Te(k) (Te ratio across zone)
c
c     chi_ihat  (output) effective ion heat diffusivity
c     u_ihat    (output) effective ion heat convection velocity
c     chi_ehat  (output) effective electron heat diffusivity
c     chi_eihat (output) effective electron-ion off-diagonal heat diffusivity
c     u_ehat    (output) effective electron heat convection velocity
c
c
c***********************************************************************
c
c Below are futher details on the definitions of the arguments and the
c logic behind this subroutine.
c
c The chi_hat's are in the same units as rho_i**2*v_ti/rmajor, and the
c u_hat convection velocities are in the same units as v_ti.
c
c The power balance equations are of the form:
c
c     (3/2) n_e d(T_e)/dt = -grad.(q_e) + P_e + ...
c     (3/2) n_i d(T_i)/dt = -grad.(q_i) + P_i + ...
c
c In Dorland's original ip_chi subroutine, the heat fluxes were:
c
c     q_e = -n_e*chi_e*grad(T_e)
c     q_i = -n_i*chi_i*grad(T_i)
c
c Here we split these up into effective pinch and diffusion terms:
c
c     q_e = -n_e*(chi_ehat*grad(T_e)+chi_eihat*grad(T_i)-u_ehat*T_e)
c     q_i = -n_i*(                  +chi_ihat *grad(T_i)-u_ihat*T_i)
c
c The two definitions of the heat fluxes are identical for the present
c value of grad(T)'s and T's.  I.e., this split into effective pinch
c and diffusion terms is exact in equilibrium, and is actually more 
c accurate for small perturabations around equilibrium.
c
c ?? Need to double check for factors of 3/2 in u_hat term?
c ?? Need to double check for n_e or n_i in power balance equations.
c
c In the present version of the code, there is no chi_ie term.  Though
c the variables chi_eihat and u_ehat exist, they are set to zero at present.
c (See the comments below for details, but in short, they didn't help much
c so I turned them off.)   This may change in future versions of the code 
c as we expand to a more complete transport matrix.
c
c For numerical stability and accuracy, it is important to be precise
c about some of the definitions.  The temperature profiles are assumed
c to be specified on a grid, and it is desired to find the heat flux
c half-way between two grid points k and k+1.  Then:
c
c     Ti_avg=(Ti(k+1)+Ti(k))/2
c     RLT = -Rmaj*(Ti(k+1)-Ti(k))/(r(k+1)-r(k))/Ti_avg
c
c and similarly for RLTe, RLN, and RLNe.  To get the best numerical
c stability, it appears important to use a centered definition of tau:
c
c     tau=Ti_avg/Te_avg
c
c It is recommended that the transport code be written with an upwind 
c differencing of the convective terms (otherwise, the temperatures can go
c negative if the diffusion terms aren't large enough).  The effective 
c convective terms u_ihat and u_ehat are calculated assuming this 
c upwind differencing (including small corrections based on the value of
c Ti_uprat and Te_uprat).
c
c***********************************************************************
c
c Finally, we describe the physics logic of this code.
c
c To avoid the numerical instability associated with critical-gradient
c chi models, we split the heat flux into an effective heat pinch term
c and an effective diffusivity term.  I.e., instead of the usual 
c convention of writing the heat flux as
c
c     q = - chi dT/dr
c
c We use:
c
c     q = q_0      + q' dT/dr 
c
c       = n hat{u} T - n hat{chi} dT/dr
c
c to define an effective heat pinch and diffusivity, where
c
c     q' = dq/d(dT/dr)
c
c is the local derivative of q in the vicinity of the present plasma
c parameters.
c
c In principle, one would include off-diagonal terms in this expansion
c of q.  For example, one could be interested in (d q_e)/(d T_i/dr),
c the change in the electron heat flux as the ion-temperature gradient
c varies.  We find that although such a term improves the convergence
c rate a little bit for small time steps, it is actually a bit less
c robust for very large time steps.  (Also, it is not clear to me that
c the off-diagonal terms can be guaranteed to lead to positive definite 
c temperatures, unlike the standard diagonal terms which can be proven
c to give positive definite temperatures no matter how large the time
c step.)  Thus we will ignore the off-diagonal terms for now, so that
c chi_eihat=u_ehat=0 is returned by this version of the code.
c
c To calculate q0 and q' numerically, we will calculate the heat flux
c at 2 different values of RLT:  the present RLT, and at (1+delGradT)*RLT.
c We find that delt=0.03 works fairly well in practice.  Convergence
c might be improved by using a larger delt (=0.2?) initially when things
c are changing rapidly, and a smaller delt for accuracy as the equilibrium
c is approached.  delt should be set to roughly the relative change
c in grad(T) expected in a single time step, as it represents the range
c in RLT over which the average chi and its derivative is sampled...
c
      use f_radial, only: f_amin, f_LTi
      implicit none
      real RLT,RLN,RLNe,shat,zth,tau,rmajor,chi_i,q,
     .     nbeam,rho_i,v_ti,eps,chi_e,g,g_perp,gnu,beta
      real rho_e,v_te
      real RLTcrit,RLTcritz,chi_0
      real Dmix, dDmix, qeqi, H_RLTcrit, H_RLTcritv
      real kappa,H_chi,aLT

      integer ichi,kstep,ifail
      logical etg


c more arguments:
      real delGradT,RLTe,chi_ihat,u_ihat,chi_eehat,chi_eihat,u_ehat
      real Ti_uprat,Te_uprat,rhohalf,gamma,gamma2

c more local variables:
      real RLT2,chi_i2,chi_e2,rmaj,RLT1
      integer ihcode_substeps

      call get_ihcode_substeps (ihcode_substeps)

      Rmaj=Rmajor

      RLT1=RLT
      if(delGradT.gt.0) then
         RLT2=RLT1+delGradT*abs(RLT1)
         if(abs(RLT1) .lt. 0.01) then
            RLT2=0.02
         endif
c
c Call Dorland's subroutine to calculate the IFS-PPPL transport model:
c
         if(ichi.eq.1 .and. ihcode_substeps.ge.0) then
            call ip_chi(RLT2,RLN,RLNe,q,kappa,
     .           shat,zth,nbeam,tau,eps,gnu,beta,
     .           g_perp,Rmaj,rho_i,v_ti,rho_e,v_te,etg,rlte,
     .           RLTcrit,RLTcritz,chi_0,g,chi_i2,chi_e2,gamma2)
            
            call ip_chi(RLT1,RLN,RLNe,q,kappa,
     .           shat,zth,nbeam,tau,eps,gnu,beta,
     .           g_perp,Rmaj,rho_i,v_ti,rho_e,v_te,etg,rlte,
     .           RLTcrit,RLTcritz,chi_0,g,chi_i,chi_e,gamma)

         elseif (ichi.eq.1) then
            call sub_ipchi(kstep,ihcode_substeps,
     >           RLT1,RLT2,RLN,RLNe,q,kappa,
     >           shat,zth,nbeam,tau,eps,gnu,beta,
     >           g_perp,Rmaj,rho_i,v_ti,rho_e,v_te,etg,rlte,
     >           RLTcrit,RLTcritz,chi_0,g,
     >           chi_i,chi_i2,chi_e,chi_e2,gamma,gamma2)
         else if(ichi.eq.3) then
            call setjr(kstep)
            aLT=f_amin(rhohalf)/f_LTi(rhohalf)
            if(aLT.ge.0.75 .and. q.gt.1) then
               call if_Dmix (kstep,rhohalf,Dmix,dDmix,qeqi,H_RLTcrit,
     .              H_RLTcritv,ifail)
               if(ifail.lt.0) then
                  write(50,*) 'H-code result may be flaky ',rhohalf
               endif
               chi_i2=H_chi(dDmix,RLT2,H_RLTcrit,H_RLTcritv,rhohalf)
     .              *rho_i**2*v_ti/f_amin(rhohalf)
               chi_i= H_chi(dDmix,RLT1,H_RLTcrit,H_RLTcritv,rhohalf)
     .              *rho_i**2*v_ti/f_amin(rhohalf)
               chi_e2=qeqi*chi_i2 
               chi_e =qeqi*chi_i
               RLTcrit = H_RLTcritv
               RLTcritz = H_RLTcritv
            else
               RLT1=0.
               chi_i=0.
               chi_i2=1.  ! m**2/sec (Ignore the center)
               chi_e=0.
               chi_e2=1.  ! m**2/sec (Ignore the center)
            endif
         else if (ichi.eq.-10) then
            ! bateman code
            call alt_chi (rhohalf,RLT1,chi_i,chi_e,RLTcrit)
            call alt_chi (rhohalf,RLT2,chi_i2,chi_e2,RLTcrit)
         else if(ichi.ge.4) then
            write(*,*) 'Set maxdel_RLT=0. if ichi>=4'
            stop
         endif
c     

c Ion Effective Thermal diffusivity and pinch velocity:
         chi_ihat=(chi_i2*RLT2-chi_i*RLT1)/(RLT2-RLT1)
         u_ihat=-(chi_ihat-chi_i)*RLT1/Rmaj

c Electron Effective Thermal diffusivities and pinch velocity:
         chi_eehat=chi_e
         
         chi_eihat=0.0
         u_ehat=0.0

c     write(*,*) RLT,RLT1,RLT2,RLTcrit
c BD 7/9/95 Another algorithm (not updated to full stability):
c (RLT1=1.01*RLT, RLT2=1.02*RLT)
c         chi_ihat=(chi_i2*RLT2-chi_i1*RLT1)/(RLT2-RLT1)
c         u_ihat=(chi_i1-chi_ihat)*RLT1/Rmaj
c         chi_eehat=max(0.,(chi_e2-chi_e1)/(RLT2-RLT1)
c     .        *(RLT-RLT1)+chi_e1)
c         chi_eihat=(chi_e2-chi_e1)/tau*RLTe/(RLT2-RLT1)
c         u_ehat=-chi_eihat*RLT/Rmaj*tau
c
c This is what I think the off-diagonal electron-ion terms should be,
c but in practice this is numerically unstable.  So instead rely on
c just the diagonal parts of the effective diffusivities and pinches,
c and on the smoothing of threshold function in ip_chi.f:
cc      chi_eihat=(chi_e2-chi_e)/tau*RLTe/(RLT2-RLT)
cc      u_ehat=-chi_eihat*RLT/Rmaj*tau
cc

c Debug:  Replace with dummy results:
c         chi_eehat=0.5+Te_r(k)*sqrt(Ti_r(k))*RLT/30.0
cc         u_ehat=-chi_eehat*rmin(k,2)/rmin(ngrid,2)**2
c         chi_eihat=0.5*chi_eehat/RLT*Te_r(k)/Ti_r(k)*RLTe
c         u_ehat=-chi_eihat*RLT/Rmaj*Ti_r(k)/Te_r(k)
         
         if(u_ihat .lt. 0.0) then
            u_ihat=u_ihat*(1.0+1.0/Ti_uprat)/2
         else
            u_ihat=u_ihat*(1.0+Ti_uprat)/2
         endif
         if(u_ehat .lt. 0.0) then
            u_ehat=u_ehat*(1.0+1.0/Te_uprat)/2
         else
            u_ehat=u_ehat*(1.0+Te_uprat)/2
         endif         
      else
         if(ichi.eq.1) then
            call ip_chi(RLT1,RLN,RLNe,q,kappa,
     .           shat,zth,nbeam,tau,eps,gnu,beta,
     .           g_perp,Rmaj,rho_i,v_ti,rho_e,v_te,etg,rlte,
     .           RLTcrit,RLTcritz,chi_0,g,chi_i,chi_e,gamma)
         else if(ichi.eq.3) then
c            call if_Dmix (kstep, rhohalf, Dmix, dDmix, qeqi, H_RLTcrit)
c            chi_i=H_chi(dDmix, RLT1, H_RLTcrit, rhohalf)
c            chi_e=qeqi*chi_i
            write(*,*) 'ichi=3 option not available'
         else if(ichi.eq.4) then
            chi_i=1.0   ! m**2/sec
            chi_e=1.0   ! m**2/sec
         else if(ichi.eq.5) then
            chi_i=1.0   ! m**2/sec
            chi_e=1.0   ! m**2/sec
         else if(ichi.eq.6) then
            chi_i=0.0   ! m**2/sec
            chi_e=1.0   ! m**2/sec
         else if(ichi.eq.-10) then
            call alt_chi (rhohalf,RLT1,chi_i,chi_e,RLTcrit)
         endif

         chi_ihat=chi_i
         chi_eehat=chi_e
         chi_eihat=0.
         u_ihat=0.
         u_ehat=0.
      endif


      return
      end

      real function H_chi(dDmix, RLT, RLTcrit, RLTcritv, rho)

      use f_radial, only: shat__r, f_RLT
      real dDmix, RLT, RLTcrit, RLTcritv, rho, dgrad
      real sfac

      sfac=(1+abs(shat__r(rho))**2.5)/(1+shat__r(rho))
      dgrad=max(1.e-10,f_RLT(rho)-RLTcrit)
      H_chi=6.6*min(1.0,1/sqrt(dgrad))
     >     *dDmix*(f_RLT(rho)-RLTcritv)
     >     *sfac
      
      return
      end

      subroutine get_ihcode_substeps (jhcode_substeps)
      use nt_data
      implicit none
cout
      integer jhcode_substeps
c

c
      jhcode_substeps = ihcode_substeps
      return
      end

      subroutine sub_ipchi(kstep,ihcode_substeps,
     >     RLT1,RLT2,RLN,RLNe,q,kappa,
     >     shat,zth,nbeam,tau,eps,gnu,beta,
     >     g_perp,Rmaj,rho_i,v_ti,rho_e,v_te,etg,rlte,
     >     RLTcrit,RLTcritz,chi_0,g,
     >     chi_i,chi_i2,chi_e,chi_e2,gamma,gamma2)
      implicit none
cin
      integer kstep,ihcode_substeps
      real RLT1,RLT2,RLN,RLNe,q,kappa,
     >     shat,zth,nbeam,tau,eps,gnu,beta,
     >     g_perp,Rmaj,rho_i,v_ti,rho_e,v_te,etg,rlte
cout
      real RLTcrit,RLTcritz,chi_0,g
      real chi_i,chi_i2,chi_e,chi_e2,gamma,gamma2
clocal
      include 'nt.par'
      integer istep
      data istep/0/
      save istep
      real rlt1_save(mxgrid)
      real rlt2_save(mxgrid)
      real chii1_save(mxgrid)
      real chii2_save(mxgrid)
      real chie1_save(mxgrid)
      real chie2_save(mxgrid)
      real gamma1_save(mxgrid)
      real gamma2_save(mxgrid)
      save rlt1_save, rlt2_save
      save chii1_save, chii2_save
      save chie1_save, chie2_save
      save gamma1_save, gamma2_save
      real rltc_save(mxgrid)
      real rltcz_save(mxgrid)
      real chi0_save(mxgrid)
      real g_save(mxgrid)
      real gperp_save(mxgrid)
      save rltc_save, rltcz_save, chi0_save, g_save, gperp_save
      real dmix1_old, dmix2_old
      integer i
c
c     Subcycle
c
      if (kstep .lt. 0) then
         istep = istep + 1
         if (istep .gt. ihcode_substeps) istep = 1
         return
      endif
      if (istep .eq. 1) then
         call ip_chi(RLT2,RLN,RLNe,q,kappa,
     .           shat,zth,nbeam,tau,eps,gnu,beta,
     .           g_perp,Rmaj,rho_i,v_ti,rho_e,v_te,etg,rlte,
     .           RLTcrit,RLTcritz,chi_0,g,chi_i2,chi_e2,gamma2)
            
         call ip_chi(RLT1,RLN,RLNe,q,kappa,
     .           shat,zth,nbeam,tau,eps,gnu,beta,
     .           g_perp,Rmaj,rho_i,v_ti,rho_e,v_te,etg,rlte,
     .           RLTcrit,RLTcritz,chi_0,g,chi_i,chi_e,gamma)

         rlt1_save(kstep) = rlt1
         rlt2_save(kstep) = rlt2
         chii1_save(kstep) = chi_i
         chii2_save(kstep) = chi_i2
         chie1_save(kstep) = chi_e
         chie2_save(kstep) = chi_e2
         gamma1_save(kstep) = gamma
         gamma2_save(kstep) = gamma2
         rltc_save(kstep) = rltcrit
         rltcz_save(kstep) = rltcritz
         chi0_save(kstep) = chi_0
         g_save(kstep) = g
         gperp_save(kstep) = g_perp
      endif

      rltcrit = rltc_save(kstep)
      rltcritz = rltcz_save(kstep)
      chi_0 = chi0_save(kstep)
      g = g_save(kstep)

      gamma = gamma1_save(kstep)
     >     + (rlt1-rlt1_save(kstep))
     >     *(gamma2_save(kstep)-gamma1_save(kstep))
     >     /(rlt2_save(kstep)  -rlt1_save(kstep))
      gamma2= gamma1_save(kstep)
     >     + (rlt2-rlt1_save(kstep))
     >     *(gamma2_save(kstep)-gamma1_save(kstep))
     >     /(rlt2_save(kstep)  -rlt1_save(kstep))
      
      dmix1_old = chii1_save(kstep)
     >     *sqrt(max(1.0,rlt1_save(kstep)-rltc_save(kstep)))
      dmix2_old = chii2_save(kstep)
     >     *sqrt(max(1.0,rlt2_save(kstep)-rltc_save(kstep)))

      chi_i = (dmix1_old
     >     + (rlt1-rlt1_save(kstep))
     >     *   (dmix1_old-dmix2_old)
     >     /   (rlt1_save(kstep)-rlt2_save(kstep)))
     >     /sqrt(max(rlt1 - rltc_save(kstep),1.0))
      chi_i2 = (dmix1_old
     >     + (rlt2-rlt1_save(kstep))
     >     *   (dmix1_old-dmix2_old)
     >     /   (rlt1_save(kstep)-rlt2_save(kstep)))
     >     /sqrt(max(rlt2 - rltc_save(kstep),1.0))
      chi_e = chie1_save(kstep)*chi_i/max(chii1_save(kstep),1e-10)
      chi_e2 = chie2_save(kstep)*chi_i2/max(chii2_save(kstep),1e-10)

      return
      end