MODULE mod_tricks

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!  MOD_TRICKS : a bunch of CEPAM optional tricks to ease the calculations
!
!  Available tricks:
!       - repartition: set repartition constants to a given value
!
!  Authors: Mathieu Havel (mathieu.havel@oca.eu)
!
!  Language: Fortran 90
!
!  Changes: 15/07/2014 -- creation
!#######################################################################

      IMPLICIT NONE
      logical, save :: init=.true.

      ! Internal: repartition variables (with default values)
      character(len=15), save :: rep_mode
      double precision, save :: cte_switch=-99.d0
      double precision, save :: ctep, ctet, ctel=0.0, cter, ctem
      logical, save :: spl_p=.true., spl_t=.true., spl_l=.false.,       &
     & spl_r=.true., spl_m=.true., use_cte=.false.

      ! namelists
      namelist/nl_repartition/rep_mode, ctep, ctet, ctel, cter, ctem

      INTERFACE repartition
            MODULE PROCEDURE repartition_bp_S, repartition_xchims
      END INTERFACE

      PRIVATE
      PUBLIC :: repartition

      CONTAINS

! -------------------------
! -     NAMELIST FILE     -
! -------------------------
      SUBROUTINE read_options()

!   Read the .don file
        implicit none
        logical :: found=.false.

        include './cepam.inc'
        include './Communs/modelp.inc'

        open(unit=55,file='cepam_tricks.don', status='old',             &
     &      form='formatted', err=500)

! file exists
        if (silent.le.0) write(*,*)'Reading file "cepam_tricks.don" '
        read(55, nl_repartition)
        close(55)
        found = .true.
    
 500    if (found) then
! set the parameters
          if (rep_mode(:3).eq.'def') then
            ctel = 0.0
          elseif (rep_mode(:3).eq.'mix') then
            if (ctep.ne.cte_switch) spl_p=.false.
            if (ctet.ne.cte_switch) spl_t=.false.
            if (ctel.ne.cte_switch) spl_l=.false.
            if (cter.ne.cte_switch) spl_r=.false.
            if (ctem.ne.cte_switch) spl_m=.false.
          elseif (rep_mode(:3).eq.'spl') then
            spl_l=.true.
          else
            use_cte=.true.
            spl_p=.false.
            spl_t=.false.
            spl_l=.false.
            spl_r=.false.
            spl_m=.false.
          endif
        endif

      END SUBROUTINE read_options

! -------------------------
! -      REPARTITION      -
! -------------------------
      SUBROUTINE repartition_bp_S(bp_S, en_rayon, rcp, rct, rcl, rcr, rcm)

!   Version: 15-Jul-2014
!   Author: Mathieu Havel
!
!   DESCRIPTION
!       This subroutine computes the repartition constants used in the file Coeur/resoutp.f
!
!   CONFIGURATION FILE
!       In the file 'cepam_tricks.don', the namelist 'NL_REPARTITION' is used to control the behavior of this
!       subroutine. The following options are available:
!
!           rep_mode = 'default'
!           ctep = 1.0
!           ctet = 0.0
!           ctel = 0.0
!           cter = 0.0
!           ctem = 0.0
!
!   NOTES
!       If rep_mode = 'default', it will behave as CEPAM versions prior to this new addition to the code
!       If rep_mode = 'spline', then spline functions will be used to compute the repartition constants.
!       If rep_mode = 'mix', spline function are used only if the value of the constants is equal to -99.
!       Otherwise, the given constant values will always be used (eg. rep_mode = 'cte').
!
!   INPUTS
!       bp_S      -- spline solution coefficients
!       en_rayon  -- ...
!
!   OUTPUT
!       rcp       -- for the pressure
!       rct       -- for the temperature
!       rcl       -- for the luminosity
!       rcr       -- for the radius
!       rcm       -- for the mass

        implicit none

        double precision, intent(in), dimension(*) :: bp_S                   ! inputs
        logical, intent(in) :: en_rayon
        double precision, intent(out) :: rcp, rct, rcl, rcr, rcm             ! outputs

        include './cepam.inc'
        include './Communs/modelp.inc'      ! get: neqs, n, m_S, r_S

! ------------------

        if (init) then
          init = .false.
          call read_options()
        endif

        if (use_cte) then
          rcp = ctep
          rct = ctet
          rcr = cter
          rcl = ctel
          rcm = ctem
          return
        endif

        if (spl_p) then
          rcp = 1./(bp_S(neqs*((n-1)*m_S+r_S-1)+1)-bp_S(1))
        else
          rcp = ctep
        endif
        if (spl_t) then
          rct = 1./(bp_S(neqs*((n-1)*m_S+r_S-1)+2)-bp_S(2))
        else
          rct = ctet
        endif
        if (spl_r) then
          if (.not.en_rayon) then
            rcr = -1./(bp_S(neqs*((n-1)*m_S+r_S-1)+3)-bp_S(3))
          else
            rcr = -1./(sqrt(abs(bp_S(neqs*((n-1)*m_S+r_S-1)+3)))-       &
     &       sqrt(abs(bp_S(3))))
          endif
        else
          rcr = cter
        endif
        if (spl_l) then
          rcl = 1./(bp_S(neqs*((n-1)*m_S+r_S-1)+4)-bp_S(4))
        else
          rcl = ctel
        endif
        if (spl_m) then
          rcm = 1./(bp_S(neqs*((n-1)*m_S+r_S-1)+5)-bp_S(5))
        else
          rcm = ctem
        endif

      END SUBROUTINE repartition_bp_S

!------------------------------------------------------

      SUBROUTINE repartition_xchims(xchims_Sc0, q0, n1, en_rayon, rcp, rct, rcl, rcr, rcm)

!   Version: 15-Jul-2014
!   Author: Mathieu Havel
!
!   DESCRIPTION
!       This subroutine computes the repartition constants used in the file Coeur/modelini.f
!
!   CONFIGURATION FILE
!       In the file 'cepam_tricks.don', the namelist 'NL_REPARTITION' is used to control the behavior of this
!       subroutine. The following options are available:
!
!           rep_mode = 'default'
!           ctep = 1.0
!           ctet = 0.0
!           ctel = 0.0
!           cter = 0.0
!           ctem = 0.0
!
!   NOTES
!       If rep_mode = 'default', it will behave as CEPAM versions prior to this new addition to the code
!       If rep_mode = 'spline', then spline functions will be used to compute the repartition constants.
!       If rep_mode = 'mix', spline function are used only if the value of the constants is equal to -99.
!       Otherwise, the given constant values will always be used (eg. rep_mode = 'cte').
!
!   INPUTS
!       xchims      -- spline solution coefficients
!       q0          -- ...
!       n1          -- ...
!       en_rayon    -- ...
!
!   OUTPUT
!       rcp       -- for the pressure
!       rct       -- for the temperature
!       rcl       -- for the luminosity
!       rcr       -- for the radius
!       rcm       -- for the mass

        implicit none

        double precision, intent(in), dimension(*) :: xchims_Sc0, q0         ! inputs
        logical, intent(in) :: en_rayon
        integer, intent(in) :: n1
        double precision, intent(out) :: rcp, rct, rcl, rcr, rcm             ! outputs

        include './cepam.inc'
        include './Communs/modelp.inc'      ! get: n1

! ------------------

        if (init) then
          init = .false.
          call read_options()
        endif

        if (use_cte) then
          rcp = ctep
          rct = ctet
          rcr = cter
          rcl = ctel
          rcm = ctem
          return
        endif

        if (spl_p) then
          rcp = 1./(xchims_Sc0(4*(n1-1)+1)-xchims_Sc0(1))
        else
          rcp = ctep
        endif
        if (spl_t) then
          rct = 1./(xchims_Sc0(4*(n1-1)+2)-xchims_Sc0(2))
        else
          rct = ctet
        endif
        if (spl_r) then
          if (.not.en_rayon) then
            rcr = -1./(xchims_Sc0(4*(n1-1)+3)-xchims_Sc0(3))
          else
            rcr = -1./(sqrt(abs(xchims_Sc0(4*(n1-1)+3)))-               &
     &        sqrt(abs(xchims_Sc0(3))))
          endif
        else
          rcr = cter
        endif
        if (spl_l) then
          rcl = 1./(xchims_Sc0(4*(n1-1)+4)-xchims_Sc0(4))
        else
          rcl = ctel
        endif
        if (spl_m) then
          rcm = 1./(q0(n1)-q0(1))
        else
          rcm = ctem
        endif

      END SUBROUTINE repartition_xchims

      END MODULE mod_tricks