! This file is part of GTC version 4.4 
! GTC version 4.4 is released under the 3-Clause BSD license:

! Copyright (c) 2002,2010,2016, GTC Team (team leader: Zhihong Lin, zhihongl@uci.edu)
! All rights reserved.

! Redistribution and use in source and binary forms, with or without 
! modification, are permitted provided that the following conditions are met:

! 1. Redistributions of source code must retain the above copyright notice, 
!    this list of conditions and the following disclaimer.

! 2. Redistributions in binary form must reproduce the above copyright notice, 
!    this list of conditions and the following disclaimer in the documentation 
!    and/or other materials provided with the distribution.

! 3. Neither the name of the GTC Team nor the names of its contributors may be 
!    used to endorse or promote products derived from this software without 
!    specific prior written permission.
! ==============================================================================
!=========================================================================
! subroutine ohm_law is used for solving non-adiabatic inductive parallel electric field (Eq. (31)) in J. Bao et al., Phys. Plasmas 2017
! subroutine ohm_law is only valid for conservative scheme
!=========================================================================
subroutine nonadiabatic_ohm_cs
  use system_env
  use precision,only: lk,machineEpsilon
  use global_parameters,only:  mpsi,mgrid,etemp0,eden0,r0,utime,ulength,rho0,deltabpara,&
                              mgridlow,mgridhigh,nonlinear,pi2,fielddir,nhybrid,eta,izonal,nfilter,gtcout,mype,istep,irk
  use particle_array,only: meshne,aelectron,qelectron,meshte,kapane,kapate,pressureepara,pressureeperp,betae
  use field_array,only: mtheta,deltat,zeta1,igrid,qmesh,qtinv,bmesh,jmesh,psimesh,gradbpara,gradpepara,gradphi,gradapara,gradue,phip00,d2apara,d4apara_para,d4apara_perp,sapara,dbmesh
  use equilibrium,only: lsp,lst,spdim,spdpsi_inv,spdtheta_inv,spdpsi,spdtheta,gpsi,cpsi,qpsi,bsp
  use petsc_array,only: luserb_ohm,luserx_ohm
  use conservative_scheme,only: d2phiind,gradphiind,gradapara_na,tapara_na,epara_field,gradphiind
  use spline_function,only: spb, spdbdt => dbdt, spdbdp=>dbdp, spdbdz=>dbdz   
#ifdef _PETSc
#include "petscversion.h"
#if PETSC_VERSION_GT(3,8,0)
  use petsc_array,only: petsc_ohm
#endif
#endif
  implicit none

  integer :: i,j,ij,isp,jst,n
  real(lk),dimension(mgrid) :: dentmp,para_grad_phiind
  real(lk) :: fullohm_na,perturb,cmratio_einv,temp,de_inv2,emfactor,b_inv,skin(0:mpsi),pdum,q_inv,dpx,dp2,&
              tdum,dtx,dt2,dbdp,dbdt,dbdz,&
              g,ri,dqdp,dgdp,dridp,gqi_inv,d2gdp2,d2ridp2,dsdp,&
              ne,ne_inv,te,dnedp,dtedp,jbb_inv,&
              dapara_nadp,dapara_nadt,dapara_nadz,&
              dphiinddp,dphiinddt,dphiinddz,&
              dbparadp,dbparadt,dbparadz,&
              dpressure_nadp,dpressure_nadt,dpressure_nadz,&
              dphidp,dphidt,dphidz,&
              dudp,dudt,dudz,&
              daparadp,daparadt,daparadz,&
              wpara_pressure,winertia,wdrive,wresistivity,whyper_para,whyper_perp,&
              wpressure_na_nl,wconvection_nl,wphiind_nl,wbpara,wphi_zonal_nl,&
              wconvection,wdrift,wmirror,tapara00_na(0:mpsi),&
              ei_collision_freq,ele_anomalous_para,ele_anomalous_perp,&
              d2aparatmp(0:1,mgrid)
  real(lk) :: tmp_field(mgrid)
  !$acc declare create(tmp_field, para_grad_phiind,skin, dentmp)

  real,external :: sprgpsi

  !$acc declare create(d2aparatmp, tapara00_na)

interface

subroutine laplacian_gpu(scalar, scalar_out, none0bound)
  use system_env
  use precision
  use global_parameters, only: mgrid,mype
  implicit none

  integer,optional, intent(in) :: none0bound
  real(lk),dimension(mgrid),intent(in) :: scalar
  real(lk),dimension(mgrid),intent(out) :: scalar_out
end subroutine laplacian_gpu

#ifdef _PETSc
#if PETSC_VERSION_GT(3,8,0)
     subroutine lapmat_pssolver(psdata,luserx_output,luserb_input)
       use precision, only:lk
       use petsc_array
       implicit none

       type (petscData),intent(in) :: psdata
       real(lk),dimension(0:),intent(inout) :: luserx_output
       real(lk),dimension(0:),intent(in) :: luserb_input
     end subroutine lapmat_pssolver
#endif
#endif

end interface

!=======================================================================
  if(nonlinear==1)then
     fullohm_na=0  ! Solve the reduced non-adiabatic ohm's law 
  elseif(nonlinear==0)then
     fullohm_na=1  ! Solve the full non-adiabatic ohm's law 
  else
     if(mype==0)write(gtcout,*)"wrong nonlinear parameter set for fullohm_na"
     stop
  endif
  if(mype==0 .and. istep==1 .and. irk==1)then
    write(gtcout,*)"equation form of the non_adiabatic ohm_law, fullohm_na=",fullohm_na

    if(fullohm_na==0)then
       write(gtcout,*)"non-adiabatic inductive electric field only contains nonlinear term {VI} in Eq. (31) of J. Bao POP 2017"
    elseif(fullohm_na==1)then
        write(gtcout,*)"non-adiabatic inductive electric field contains all terms in Eq. (31) of J.Bao et al., Phys. Plasmas 2017"
    endif

    flush(gtcout)
  endif
!=======================================================================

  perturb=real(nonlinear,lk)
  cmratio_einv=aelectron/qelectron
! Normalized square of debye length using axis equilibrium parameters
  temp=(etemp0/eden0)*7.43e2*7.43e2/(r0*r0)
! Normalized light speed
  emfactor=3.0*1e10*utime/ulength
! Normalized inverse of the square of electron skin depth using axis density
  de_inv2=rho0*rho0/temp/aelectron/emfactor/emfactor

! 10_lk is the approximate coulomb logarithm, ei_collision_freq is the normalized electron-to-ion collision frequency
  ei_collision_freq=0.0_lk*(2.91e-6_lk*eden0*10.0_lk/(etemp0*sqrt(etemp0)))*(pi2*utime)
! ele_anomalous_coe is the anomalous electron hyper viscosity,100.0 means 100cm^2/s=1m^2/s
  ele_anomalous_para=0.0_lk*(pi2*utime)/ulength/ulength
  ele_anomalous_perp=0.0_lk*(pi2*utime)/ulength/ulength


  if(mype==0 .and. istep==1 .and. irk==1)then
    write(gtcout,*)'eden0=',eden0,'etemp0=',etemp0
    write(gtcout,*)'utime=',utime,'ulength=',ulength
    write(gtcout,*)'electron-ion collision freq=',ei_collision_freq,'=',ei_collision_freq/(pi2*utime),'Hz'
    write(gtcout,*)'parallel anomalous coefficient ele_anomalous_para =',ele_anomalous_para*ulength*ulength/(pi2*utime),'cm^2/s'
    write(gtcout,*)'perpendicular anomalous coefficient ele_anomalous_perp =',ele_anomalous_perp*ulength*ulength/(pi2*utime),'cm^2/s'
    flush(gtcout)
  endif

! calculate 4-th order derivative of sapara in parallel direction for hyperviscosity in ohm's law
  !$acc kernels
  d4apara_para(:)=sapara(1,:)
  d2aparatmp=0.0_lk
  !$acc end kernels
  do n = 1, 2
    call grad2para(d4apara_para,d2aparatmp)
    !$acc kernels
    d4apara_para(:)=d2aparatmp(1,:)
    !$acc end kernels
  enddo


  !$acc kernels
  d4apara_perp(:)=sapara(1,:)
  d2aparatmp=0.0_lk
  !$acc end kernels
  do n = 1, 2
     call laplacian_gpu(d4apara_perp,tmp_field)
     !$acc kernels
     d2aparatmp(1,:)=tmp_field
     !$acc end kernels
     call smooth_gpu(d2aparatmp)
     !$acc kernels
     d4apara_perp(:)=d2aparatmp(1,:)
     !$acc end kernels
  enddo


!Calculation of the inverse of square of electron skin depth with local density
#ifdef _OPENACC
  !$acc parallel loop gang vector
#else
  !$omp parallel do private(i)
#endif
  do i=0,mpsi
     skin(i)=de_inv2*meshne(i)
  enddo
! Calulation of d2phiind: laplacian of the parallel gradient of phiind
  !$acc kernels
  para_grad_phiind=0.0_lk
  !$acc end kernels
#ifdef _OPENACC
  !$acc parallel loop gang
#else
  !$omp parallel do private(i,j,ij,b_inv)
#endif
  do i=0,mpsi
     !$acc loop vector
     do j=1,mtheta(i)
        ij=igrid(i)+j
        b_inv=1.0_lk/bmesh(ij)
        para_grad_phiind(ij)=qmesh(i)*gradphiind(3,1,ij)*b_inv*jmesh(ij)
     enddo
     para_grad_phiind(igrid(i))=para_grad_phiind(igrid(i)+mtheta(i))
  enddo
  call laplacian_gpu(para_grad_phiind,d2phiind)

  !$acc kernels
  dentmp=0.0_lk
  !$acc end kernels


#ifdef _OPENACC
  !$acc parallel loop gang
#else
  !$omp parallel do private(i,j,ij,pdum,q_inv,isp,dpx,dp2,&
  !$omp& tdum,jst,dtx,dt2,dbdp,dbdt,dbdz,&
  !$omp& g,ri,dqdp,dgdp,dridp,gqi_inv,d2gdp2,d2ridp2,dsdp,&
  !$omp& ne,ne_inv,te,dnedp,dtedp,b_inv,jbb_inv,&
  !$omp& dapara_nadp,dapara_nadt,dapara_nadz,&
  !$omp& daparadp,daparadt,daparadz,&
  !$omp& dpressure_nadp,dpressure_nadt,dpressure_nadz,&
  !$omp& dphiinddp,dphiinddt,dphiinddz,&
  !$omp& dphidp,dphidt,dphidz,&
  !$omp& dbparadp,dbparadt,dbparadz,&
  !$omp& dudp,dudt,dudz,&
  !$omp& wpara_pressure,winertia,wdrive,wresistivity,whyper_para,whyper_perp,&
  !$omp& wconvection,wdrift,wmirror,&
  !$omp& wpressure_na_nl,wbpara,wphiind_nl,wphi_zonal_nl,wconvection_nl)
#endif
  do i=0,mpsi
     pdum=psimesh(i)
     q_inv=1/qmesh(i)

     isp=max(1,min(lsp-1,ceiling(pdum*spdpsi_inv)))
     dpx=pdum-spdpsi*real(isp-1)
     dp2=dpx*dpx

     g=   gpsi(1,isp)   +gpsi(2,isp)*dpx   +gpsi(3,isp)*dp2
     ri=  cpsi(1,isp)   +cpsi(2,isp)*dpx   +cpsi(3,isp)*dp2

     dqdp=               qpsi(2,isp)       +2.0*qpsi(3,isp)*dpx
     dgdp=               gpsi(2,isp)       +2.0*gpsi(3,isp)*dpx
     dridp=              cpsi(2,isp)       +2.0*cpsi(3,isp)*dpx

     gqi_inv=1.0/(g*qmesh(i)+ri)

! To make the 2nd derivatives of I and g continuous,
! currently we use linear interpolations as a temporary solution.
! In the future, higher order spline functions of I and g need to be introduced
     d2gdp2=2.0*(dpx*spdpsi_inv*(gpsi(3,isp+1)-gpsi(3,isp))+gpsi(3,isp))
     d2ridp2=2.0*(dpx*spdpsi_inv*(cpsi(3,isp+1)-cpsi(3,isp))+cpsi(3,isp))

     dsdp=(g*d2ridp2-ri*d2gdp2)/(g*qmesh(i)+ri) &
        -(g*dridp-ri*dgdp)*(dgdp*qmesh(i)+g*dqdp+dridp)*((gqi_inv)**2)

     ne=meshne(i)
     ne_inv=1.0_lk/ne
     te=meshte(i)
     dnedp=-ne*kapane(i)
     dtedp=-te*kapate(i)
     !$acc loop vector
     do j=1,mtheta(i)
        ij=igrid(i)+j
        b_inv=1.0_lk/bmesh(ij)
        jbb_inv=jmesh(ij)*b_inv*b_inv


        if(fielddir==1 .or. fielddir==3)then
          tdum=modulo(deltat(i)*real(j)+(zeta1-pi2)*qtinv(i),pi2)
        else
          tdum=modulo(deltat(i)*real(j)+zeta1*qtinv(i),pi2)
        endif
        dbdt = dbmesh(2,ij)!spdbdt(pdum, tdum)
        dbdp = dbmesh(1,ij)!spdbdp(pdum, tdum)
        dbdz = dbmesh(3,ij)!spdbdz(pdum, tdum)

        dapara_nadp=gradapara_na(1,1,ij)
        dapara_nadt=gradapara_na(2,1,ij)
        dapara_nadz=gradapara_na(3,1,ij)-gradapara_na(2,1,ij)*q_inv

        daparadp=gradapara(1,1,ij)
        daparadt=gradapara(2,1,ij)
        daparadz=gradapara(3,1,ij)-gradapara(2,1,ij)*q_inv

        if(nhybrid==1)then
           dpressure_nadp=gradpepara(1,1,ij)
           dpressure_nadt=gradpepara(2,1,ij)
           dpressure_nadz=gradpepara(3,1,ij)-gradpepara(2,1,ij)*q_inv
        else
           dpressure_nadp=0.0_lk
           dpressure_nadt=0.0_lk
           dpressure_nadz=0.0_lk
        endif

        dphiinddp=gradphiind(1,1,ij)
        dphiinddt=gradphiind(2,1,ij)
        dphiinddz=gradphiind(3,1,ij)-gradphiind(2,1,ij)*q_inv
        if (deltabpara==1) then
        dbparadp=gradbpara(1,1,ij)
        dbparadt=gradbpara(2,1,ij)
        dbparadz=gradbpara(3,1,ij)-gradbpara(2,1,ij)*q_inv
        else
        dbparadp=0.0_lk
        dbparadt=0.0_lk
        dbparadz=0.0_lk
        endif
        dphidp=gradphi(1,1,ij)
        dphidt=gradphi(2,1,ij)
        dphidz=gradphi(3,1,ij)-gradphi(2,1,ij)*q_inv
 
        dudp=gradue(1,1,ij)
        dudt=gradue(2,1,ij)
        dudz=gradue(3,1,ij)-gradue(2,1,ij)*q_inv

!Following are the terms appeared in Eq. (31) of J. Bao et al., Phys. Plasmas2017

!Linear term{II}:electron inertia associated with phiind (adiabatic vector potential)
        winertia=d2phiind(ij)

!Linear term{III}: parallel non-adiabatic pressure force
        if(nhybrid==1)then
           wpara_pressure=skin(i)*jmesh(ij)*qmesh(i)*b_inv*gradpepara(3,1,ij)
        else
           wpara_pressure=0.0_lk
        endif

!Linear term{IV}: force of equilibrium pressure along the non-adiabatic magnetic field direction
        wdrive=skin(i)*jbb_inv*(g*dapara_nadt-ri*dapara_nadz)*(te*ne_inv*dnedp+dtedp)

!Linear convection in term{VIII}
        wconvection=-jmesh(ij)*b_inv*dsdp*(ri*dphidz-g*dphidt)

!Grad B drift dot grad delta upara in term{VIII}
        wdrift=4.0_lk*skin(i)*cmratio_einv*ne_inv*te*jbb_inv*(-ri*dbdp*dudz+g*dbdp*dudt+(ri*dbdz-g*dbdt)*dudp)

!Mirror force term associated with non-adiabatic pressure in term {X}
        if(nhybrid==1)then
           wmirror=-skin(i)*(pressureepara(1,ij)-pressureeperp(1,ij))*jbb_inv*(dbdt+qmesh(i)*dbdz)
        else
           wmirror=0.0_lk
        endif

! Resistivity     
        wresistivity=ei_collision_freq*d2apara(ij)

! Hyper-viscosity
        whyper_para=-ele_anomalous_para*d4apara_para(ij)
        whyper_perp=-ele_anomalous_perp*d4apara_perp(ij)

!Non-linear term{V}:non-adiabatic pressure force along the total perturbed magnetic field direction
        if(nhybrid==1)then
           wpressure_na_nl=skin(i)*jbb_inv*(dpressure_nadp*(g*daparadt-ri*daparadz)+daparadp*(ri*dpressure_nadz-g*dpressure_nadt))
        else 
           wpressure_na_nl=0.0_lk
        endif

!Non-linear term{VI}:adiabatic pressure force and electrostatic electric field force along total perturbed magnetic field direction
        wphiind_nl=skin(i)*jbb_inv*(dphiinddp*(g*daparadt-ri*daparadz)+daparadp*(ri*dphiinddz-g*dphiinddt))
        wbpara=skin(i)*jbb_inv*(dbparadp*(g*daparadt-ri*daparadz)+daparadp*(ri*dbparadz-g*dbparadt))

!Non-linear term{VII}:zonal electrostatic electric field force along total perturbed magentic field direction
        wphi_zonal_nl=-skin(i)*jbb_inv*phip00(i)*(g*daparadt-ri*daparadz)

!Non-linear term {IX}: Nonlinear convection (only important when kperp^2*de^2 is the order of 1)
        wconvection_nl=skin(i)*aelectron*ne_inv*jmesh(ij)*b_inv*(dphidp*(g*dudt-ri*dudz)+dudp*(ri*dphidz-g*dphidt))

        if(fullohm_na==1)then
           dentmp(ij)=wpara_pressure+winertia+wdrive+wresistivity+whyper_para+whyper_perp+&
                      wdrift+wconvection+wmirror+&          
                      perturb*(wpressure_na_nl+wphiind_nl+wbpara+wphi_zonal_nl+wconvection_nl)
        elseif(fullohm_na==0)then
           dentmp(ij)=perturb*(wphiind_nl+wbpara) ! remove the zonal phi component
        else
           if(mype==0)write(gtcout,*)"wrong set of fullohm_na, break down in ohm_law.F90"
           stop
        endif
     enddo
     dentmp(igrid(i))=dentmp(igrid(i)+mtheta(i))
  enddo

  if(fullohm_na==1)then

#ifdef _PETSc
#ifdef _OPENACC
     !$acc parallel loop gang vector
#else
     !$omp parallel do private(i)
#endif
     do i=mgridlow,mgridhigh
        luserb_ohm(i-mgridlow)=dentmp(i)
     enddo
     !$acc update host(luserb_ohm)
#if PETSC_VERSION_GT(3,8,0)
     call lapmat_pssolver(petsc_ohm,luserx_ohm,luserb_ohm)
#else
     call lapmat_ohmsolver
#endif
     !$acc update device(luserx_ohm)
    
#ifdef _OPENACC
     !$acc parallel loop gang vector
#else
     !$omp parallel do private(i)
#endif
     do i=mgridlow,mgridhigh
        tapara_na(1,i)=luserx_ohm(i-mgridlow)
     enddo
#else
     if (mype==0) then
        write(gtcout,*) 'hypre and amgx are not ready for lapmat_ohmsolver'
        flush(gtcout)
        call GTC_ABORT
     endif
#endif
  elseif(fullohm_na==0)then
#ifdef _OPENACC
     !$acc parallel loop gang
#else
     !$omp parallel do private(i,j,ij)
#endif
     do i=0,mpsi
        !$acc loop vector
        do j=1,mtheta(i)
           ij=igrid(i)+j
           tapara_na(1,ij)=dentmp(ij)/skin(i)
        enddo
        tapara_na(1,igrid(i))=tapara_na(1,igrid(i)+mtheta(i))
     enddo
     tapara_na(1,igrid(0):igrid(0)+mtheta(0))=0.0_lk
     tapara_na(1,igrid(mpsi):igrid(mpsi)+mtheta(mpsi))=0.0_lk

  endif

  call smooth_gpu(tapara_na)
  if(nfilter>0)then
     if(izonal==0)then
        call filter_gpu(tapara_na)
     elseif(izonal>0)then
        call filter_keep_zonal_gpu(tapara_na,tapara00_na)
     endif
  endif

! Diagnostic of epara_field and phiind in cs_model
  call diagnostic_cs_model

end subroutine nonadiabatic_ohm_cs
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!=====================================================================
!  Evaluate dnet for kinetic electron push in conservative scheme
!  Note that dnet is part of Eq. (29) in J. Bao et al., Phys. Plasmas 2017
!=====================================================================
subroutine dphieff_cs
  use system_env
  use global_parameters,only: mpsi
  use particle_array,only: dnet
  use field_array,only: mtheta,igrid,sdnedt
  implicit none

  integer :: i,j,ij
#ifdef _OPENACC
  !$acc parallel loop gang
#else
  !$omp parallel do private(i,j,ij)
#endif
  do i=0,mpsi
     !$acc loop vector
     do j=1,mtheta(i)
        ij=igrid(i)+j
        dnet(:,ij)=sdnedt(:,ij)
     enddo
     dnet(:,igrid(i))=dnet(:,igrid(i)+mtheta(i))
  enddo

end subroutine dphieff_cs
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine diagnostic_cs_model
  use precision
  use system_env
  use global_parameters
  use field_array
  use particle_array
  use conservative_scheme, only: phiind,epara_field,gradphiind,tapara_na
  implicit none

  integer i,j,ij
  real(lk) :: ne,ne_inv,te,dnedp,gradn_term,b_inv,phiind00(0:mpsi),tmp_field(mgrid)
  real(lk) :: tmp1, tmp2, tmp3,tmp4
  !$acc declare create(phiind00,tmp_field)

!=======================================================================
!     Solve phiind from density perturbation for conservative scheme    
!     Eq. (11) in J. Bao et al Phys. Plasmas 2017   
!=======================================================================

#ifdef _OPENACC
  !$acc parallel loop gang
#else
  !$omp parallel do private(i,j,ij,ne,ne_inv,te,dnedp,gradn_term)
#endif
  do i=0,mpsi
     ne=meshne(i)
     ne_inv=1.0/ne
     te=meshte(i)
     dnedp=-ne*kapane(i)
     !$acc loop vector
     do j=1,mtheta(i)
        ij=igrid(i)+j
        gradn_term=-te*ne_inv*dnedp*sdeltapsi(1,ij)

        phiind(1,ij)=-(phi(1,ij)-sdelfluidne(1,ij)*te*ne_inv)+gradn_term
     enddo
     phiind(1,igrid(i))=phiind(1,igrid(i)+mtheta(i))
  enddo

  CALL SMOOTH_GPU(phiind)
  if(nfilter>0)CALL FILTER_GPU(phiind)

! Safety guard: to guarantee phiind to be purely non-zonal
  !$acc kernels
  tmp_field=phiind(1,:)
  !$acc end kernels
  call zonal_field_gpu(tmp_field,phiind00)

#ifdef _OPENACC
  !$acc parallel loop gang
#else
  !$omp parallel do private(i,j,ij)
#endif
  do i=0,mpsi
     !$acc loop vector
     do j=1,mtheta(i)
        ij=igrid(i)+j
        phiind(0,ij)=phiind(0,ij)-phiind00(i)
        phiind(1,ij)=phiind(1,ij)-phiind00(i)
     enddo
     phiind(0,igrid(i))=phiind(0,igrid(i)+mtheta(i))
     phiind(1,igrid(i))=phiind(1,igrid(i)+mtheta(i))
  enddo

#ifdef _OPENACC
  !$acc parallel loop gang
#else
  !$omp parallel do private(i,j,ij,b_inv)
#endif
  do i=0,mpsi
     !$acc loop vector
     do j=1,mtheta(i)
        ij=igrid(i)+j
        b_inv=1.0_lk/bmesh(ij)
        epara_field(1,ij)=-jmesh(ij)*qmesh(i)*b_inv*(gradphi(3,1,ij)+gradphiind(3,1,ij))-tapara_na(1,ij)
     enddo
     epara_field(1,igrid(i))=epara_field(1,igrid(i)+mtheta(i))
  enddo

  call periodicity_gpu(epara_field)

end subroutine diagnostic_cs_model
