calculs_skku_exp3 Subroutine

   subroutine calculs_skku_exp3(bounds, lg, ssk, sku)
   !================================================================================================
   !<@note Function to calculate the skewness and kurtosis of an **exponential** series.<br/>
   !< The principle is the same as [[calculs_skku_tan]], however it fits better some particular
   !< series quite binary (roughly two heights).
   !<
   !<@endnote
   !------------------------------------------------------------------------------------------------
   implicit none
   real   (kind=R8), intent(in), dimension(1:3) :: bounds  !! *interval limits* [-(1/bounds(1)-1), +(1/bounds(2)-1)]
   integer(kind=I4), intent(in)                 :: lg      !! *vec size*
   real   (kind=R8), intent(out)                :: ssk     !! *theoretical Ssk*
   real   (kind=R8), intent(out)                :: sku     !! *theoretical Sku*

      real(kind=R8) :: xa, xb, kk, kk2, kk3, kk4, mu, si, sk, ku, a, b, c
      real(kind=R8) :: mu2, mu3, mu4, si3, si4
      real(kind=R8) :: h, hh, b1, b2, alp, bet, gam
      real(kind=R8) :: exp1b, exp1a, exp2b, exp2a, exp3b, exp3a, exp4b, exp4a, tmp1a, tmp1b, tmp2a, tmp2b, tmp3a, tmp3b, tmp4a, tmp4b

      real(kind=R8) :: tmp1a2, tmp1b2, tmp1a3, tmp1b3,tmp1a4, tmp1b4, tmp1a5 , tmp1b5 , tmp1a6 , tmp1b6
      real(kind=R8) :: tmp1a7, tmp1b7, tmp1a8, tmp1b8,tmp1a9, tmp1b9, tmp1a10, tmp1b10, tmp1a13, tmp1b13

      integer(kind=I4) :: deb, fin

      deb = 1
      fin = lg

      a = bounds(1)
      b = bounds(2)
      c = bounds(3)

      hh = (-2._R8 + UN/a + UN/b) / (lg - 1)
      h  = hh

      xa = a
      xb = b

      b1 = -(UN - a)/a
      b2 = +(UN - b)/b

      alp = -(b2 - lg*b1)/(b2 - b1)
      bet =      (lg - 1)/(b2 - b1)

      kk  = c / (b2 - b1)**3
      gam = kk

      kk2 = kk**2
      kk3 = kk**3
      kk4 = kk**4

      tmp1a = -1*(UN - UN/xa) ; tmp1a = min(+0.9*HIG_E8, tmp1a) ; tmp1a = max(-0.9*HIG_E8, tmp1a)
      tmp1b = -1*(UN - UN/xb) ; tmp1b = min(+0.9*HIG_E8, tmp1b) ; tmp1b = max(-0.9*HIG_E8, tmp1b)

      tmp2a = -2*(UN - UN/xa) ; tmp2a = min(+0.9*HIG_E8, tmp2a) ; tmp2a = max(-0.9*HIG_E8, tmp2a)
      tmp2b = -2*(UN - UN/xb) ; tmp2b = min(+0.9*HIG_E8, tmp2b) ; tmp2b = max(-0.9*HIG_E8, tmp2b)

      tmp3a = -3*(UN - UN/xa) ; tmp3a = min(+0.9*HIG_E8, tmp3a) ; tmp3a = max(-0.9*HIG_E8, tmp3a)
      tmp3b = -3*(UN - UN/xb) ; tmp3b = min(+0.9*HIG_E8, tmp3b) ; tmp3b = max(-0.9*HIG_E8, tmp3b)

      tmp4a = -4*(UN - UN/xa) ; tmp4a = min(+0.9*HIG_E8, tmp4a) ; tmp4a = max(-0.9*HIG_E8, tmp4a)
      tmp4b = -4*(UN - UN/xb) ; tmp4b = min(+0.9*HIG_E8, tmp4b) ; tmp4b = max(-0.9*HIG_E8, tmp4b)

      tmp1a2  = tmp1a**2  ; tmp1b2  = tmp1b**2
      tmp1a3  = tmp1a**3  ; tmp1b3  = tmp1b**3
      tmp1a4  = tmp1a**4  ; tmp1b4  = tmp1b**4
      tmp1a5  = tmp1a**5  ; tmp1b5  = tmp1b**5
      tmp1a6  = tmp1a**6  ; tmp1b6  = tmp1b**6
      tmp1a7  = tmp1a**7  ; tmp1b7  = tmp1b**7
      tmp1a8  = tmp1a**8  ; tmp1b8  = tmp1b**8
      tmp1a9  = tmp1a**9  ; tmp1b9  = tmp1b**9
      tmp1a10 = tmp1a**10 ; tmp1b10 = tmp1b**10
      tmp1a13 = tmp1a**13 ; tmp1b13 = tmp1b**13

      exp1a = exp(-tmp1a) ! exp( 1 * (1 - 1/xa) )
      exp1b = exp(-tmp1b) ! exp( 1 * (1 - 1/xb) )
      exp2a = exp(-tmp2a) ! exp( 2 * (1 - 1/xa) )
      exp2b = exp(-tmp2b) ! exp( 2 * (1 - 1/xb) )
      exp3a = exp(-tmp3a) ! exp( 3 * (1 - 1/xa) )
      exp3b = exp(-tmp3b) ! exp( 3 * (1 - 1/xb) )
      exp4a = exp(-tmp4a) ! exp( 4 * (1 - 1/xa) )
      exp4b = exp(-tmp4b) ! exp( 4 * (1 - 1/xb) )

      mu = -exp1a+exp1b-(kk*tmp1a4)/4.0D0+(kk*tmp1b4)/4.0D0-1.0D0/xa+1.0D0/xb

      mu = (UN/h)*mu + add_expo3(1, deb, fin, alp, bet, gam, mu=0._R8, si=1._R8)

      mu = mu/lg

      mu2 = mu**2
      mu3 = mu**3
      mu4 = mu**4

      !. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

      si = kk*(-2.4D+1)-exp2a/2.0D0-exp2b/2.0D0+exp1a*(kk*1.2D+1+mu*2.0D0+2.0D0)+exp1b*(kk*1.2D+1-mu*2.0D0+2.0D0)+tmp1a*(mu*2.0D0+mu2+1.0D0)+tmp1b*(mu*(-2.0D0)+mu2+1.0D0)+(kk2*tmp1a7)/7.0D0+(kk2*tmp1b7)/7.0D0+(kk*tmp1a4*(mu+1.0D0))/2.0D0-(kk*tmp1b4*(mu-1.0D0))/2.0D0+kk*exp1a*tmp1a2*6.0D0+kk*exp1a*tmp1a3*2.0D0+kk*exp1b*tmp1b2*6.0D0+kk*exp1b*tmp1b3*2.0D0+kk*exp1a*tmp1a*1.2D+1+kk*exp1b*tmp1b*1.2D+1-3.0D0

      si = (UN/h)*si + add_expo3(2, deb, fin, alp, bet, gam, mu, si=1._R8)

      si = si/lg
      si = sqrt(si)

      si3 = si**3
      si4 = si**4

      !. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

      sk = -1.0D0/si3*(kk*(-2.79D+2/8.0D0)-mu*(9.0D0/2.0D0)+exp3a/3.0D0-kk*mu*3.6D+1-exp2a*(kk*(9.0D0/8.0D0)+mu*(3.0D0/2.0D0)+3.0D0/2.0D0)+tmp1a*(mu*3.0D0+mu2*3.0D0+mu3+1.0D0)+(kk3*tmp1a10)/1.0D+1-kk2*2.16D+3-mu2*3.0D0+exp1a*(kk*3.6D+1+mu*6.0D0+kk*mu*3.6D+1+kk2*2.16D+3+mu2*3.0D0+3.0D0)-kk*exp2a*tmp1a2*(9.0D0/4.0D0)-kk*exp2a*tmp1a3*(3.0D0/2.0D0)+kk*tmp1a4*(mu+1.0D0)**2*(3.0D0/4.0D0)+kk2*tmp1a7*(mu+1.0D0)*(3.0D0/7.0D0)+kk2*exp1a*tmp1a4*9.0D+1+kk2*exp1a*tmp1a5*1.8D+1+kk2*exp1a*tmp1a6*3.0D0-kk*exp2a*tmp1a*(9.0D0/4.0D0)+kk*exp1a*tmp1a*(kk*6.0D+1+mu+1.0D0)*3.6D+1+kk*exp1a*tmp1a2*(kk*6.0D+1+mu+1.0D0)*1.8D+1+kk*exp1a*tmp1a3*(kk*6.0D+1+mu+1.0D0)*6.0D0-1.1D+1/6.0D0)+1.0D0/si3*(kk*(-2.79D+2/8.0D0)+mu*(9.0D0/2.0D0)+exp3b/3.0D0+kk*mu*3.6D+1-exp2b*(kk*(9.0D0/8.0D0)-mu*(3.0D0/2.0D0)+3.0D0/2.0D0)-tmp1b*(mu*3.0D0-mu2*3.0D0+mu3-1.0D0)+(kk3*tmp1b10)/1.0D+1-kk2*2.16D+3-mu2*3.0D0+exp1b*(kk*3.6D+1-mu*6.0D0-kk*mu*3.6D+1+kk2*2.16D+3+mu2*3.0D0+3.0D0)-kk*exp2b*tmp1b2*(9.0D0/4.0D0)-kk*exp2b*tmp1b3*(3.0D0/2.0D0)+kk*tmp1b4*(mu-1.0D0)**2*(3.0D0/4.0D0)-kk2*tmp1b7*(mu-1.0D0)*(3.0D0/7.0D0)+kk2*exp1b*tmp1b4*9.0D+1+kk2*exp1b*tmp1b5*1.8D+1+kk2*exp1b*tmp1b6*3.0D0-kk*exp2b*tmp1b*(9.0D0/4.0D0)+kk*exp1b*tmp1b2*(kk*6.0D+1-mu+1.0D0)*1.8D+1+kk*exp1b*tmp1b3*(kk*6.0D+1-mu+1.0D0)*6.0D0+kk*exp1b*tmp1b*(kk*6.0D+1-mu+1.0D0)*3.6D+1-1.1D+1/6.0D0)

      sk = (UN/h)*sk + add_expo3(3, deb, fin, alp, bet, gam, mu, si)

      sk = sk/lg
      !. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

      ku =  1.0D0/si4*(kk*(-6.77962962962963D+1)-mu*(2.2D+1/3.0D0)-exp4a/4.0D0-kk*mu*(2.79D+2/2.0D0)+exp3a*(kk*(8.0D0/2.7D+1)+mu*(4.0D0/3.0D0)+4.0D0/3.0D0)+exp1a*(kk*7.2D+1+mu*1.2D+1+kk*mu*1.44D+2+kk*mu2*7.2D+1+kk2*mu*8.64D+3+kk2*8.64D+3+kk3*1.45152D+6+mu2*1.2D+1+mu3*4.0D0+4.0D0)-kk*mu2*7.2D+1-kk2*mu*8.64D+3+(kk4*tmp1a13)/1.3D+1-kk2*8.60625D+3-kk3*1.45152D+6-mu2*9.0D0-mu3*4.0D0-exp2a*(kk*(9.0D0/2.0D0)+mu*6.0D0+kk*mu*(9.0D0/2.0D0)+kk2*(1.35D+2/4.0D0)+mu2*3.0D0+3.0D0)+tmp1a*(mu*4.0D0+mu2*6.0D0+mu3*4.0D0+mu4+1.0D0)+kk*exp3a*tmp1a2*(4.0D0/3.0D0)+kk*exp3a*tmp1a3*(4.0D0/3.0D0)+kk*tmp1a4*(mu+1.0D0)**3+kk3*tmp1a10*(mu+1.0D0)*(2.0D0/5.0D0)-kk2*exp2a*tmp1a4*(4.5D+1/2.0D0)-kk2*exp2a*tmp1a5*9.0D0-kk2*exp2a*tmp1a6*3.0D0+kk3*exp1a*tmp1a7*2.88D+2+kk3*exp1a*tmp1a8*3.6D+1+kk3*exp1a*tmp1a9*4.0D0+kk2*tmp1a7*(mu+1.0D0)**2*(6.0D0/7.0D0)+kk*exp3a*tmp1a*(8.0D0/9.0D0)+kk*exp1a*tmp1a*(kk*1.2D+2+mu*2.0D0+kk*mu*1.2D+2+kk2*2.016D+4+mu2+1.0D0)*7.2D+1-kk*exp2a*tmp1a2*(kk*1.5D+1+mu*2.0D0+2.0D0)*(9.0D0/2.0D0)-kk*exp2a*tmp1a3*(kk*1.5D+1+mu*2.0D0+2.0D0)*3.0D0+kk2*exp1a*tmp1a4*(kk*1.68D+2+mu+1.0D0)*3.6D+2+kk2*exp1a*tmp1a5*(kk*1.68D+2+mu+1.0D0)*7.2D+1+kk2*exp1a*tmp1a6*(kk*1.68D+2+mu+1.0D0)*1.2D+1+kk*exp1a*tmp1a2*(kk*1.2D+2+mu*2.0D0+kk*mu*1.2D+2+kk2*2.016D+4+mu2+1.0D0)*3.6D+1+kk*exp1a*tmp1a3*(kk*1.2D+2+mu*2.0D0+kk*mu*1.2D+2+kk2*2.016D+4+mu2+1.0D0)*1.2D+1-kk*exp2a*tmp1a*(kk*1.5D+1+mu*2.0D0+2.0D0)*(9.0D0/2.0D0)-2.5D+1/1.2D+1)+1.0D0/si4*(kk*(-6.77962962962963D+1)+mu*(2.2D+1/3.0D0)-exp4b/4.0D0+kk*mu*(2.79D+2/2.0D0)+exp3b*(kk*(8.0D0/2.7D+1)-mu*(4.0D0/3.0D0)+4.0D0/3.0D0)+exp1b*(kk*7.2D+1-mu*1.2D+1-kk*mu*1.44D+2+kk*mu2*7.2D+1-kk2*mu*8.64D+3+kk2*8.64D+3+kk3*1.45152D+6+mu2*1.2D+1-mu3*4.0D0+4.0D0)-kk*mu2*7.2D+1+kk2*mu*8.64D+3+(kk4*tmp1b13)/1.3D+1-kk2*8.60625D+3-kk3*1.45152D+6-mu2*9.0D0+mu3*4.0D0-exp2b*(kk*(9.0D0/2.0D0)-mu*6.0D0-kk*mu*(9.0D0/2.0D0)+kk2*(1.35D+2/4.0D0)+mu2*3.0D0+3.0D0)+tmp1b*(mu*(-4.0D0)+mu2*6.0D0-mu3*4.0D0+mu4+1.0D0)+kk*exp3b*tmp1b2*(4.0D0/3.0D0)+kk*exp3b*tmp1b3*(4.0D0/3.0D0)-kk*tmp1b4*(mu-1.0D0)**3-kk3*tmp1b10*(mu-1.0D0)*(2.0D0/5.0D0)-kk2*exp2b*tmp1b4*(4.5D+1/2.0D0)-kk2*exp2b*tmp1b5*9.0D0-kk2*exp2b*tmp1b6*3.0D0+kk3*exp1b*tmp1b7*2.88D+2+kk3*exp1b*tmp1b8*3.6D+1+kk3*exp1b*tmp1b9*4.0D0+kk2*tmp1b7*(mu-1.0D0)**2*(6.0D0/7.0D0)+kk*exp3b*tmp1b*(8.0D0/9.0D0)+kk*exp1b*tmp1b*(kk*1.2D+2-mu*2.0D0-kk*mu*1.2D+2+kk2*2.016D+4+mu2+1.0D0)*7.2D+1-kk*exp2b*tmp1b2*(kk*1.5D+1-mu*2.0D0+2.0D0)*(9.0D0/2.0D0)-kk*exp2b*tmp1b3*(kk*1.5D+1-mu*2.0D0+2.0D0)*3.0D0+kk*exp1b*tmp1b2*(kk*1.2D+2-mu*2.0D0-kk*mu*1.2D+2+kk2*2.016D+4+mu2+1.0D0)*3.6D+1+kk*exp1b*tmp1b3*(kk*1.2D+2-mu*2.0D0-kk*mu*1.2D+2+kk2*2.016D+4+mu2+1.0D0)*1.2D+1+kk2*exp1b*tmp1b4*(kk*1.68D+2-mu+1.0D0)*3.6D+2+kk2*exp1b*tmp1b5*(kk*1.68D+2-mu+1.0D0)*7.2D+1+kk2*exp1b*tmp1b6*(kk*1.68D+2-mu+1.0D0)*1.2D+1-kk*exp2b*tmp1b*(kk*1.5D+1-mu*2.0D0+2.0D0)*(9.0D0/2.0D0)-2.5D+1/1.2D+1)

      ku = (UN/h)*ku + add_expo3(4, deb, fin, alp, bet, gam, mu, si)

      ku = ku/lg

      ssk = sk
      sku = ku

   contains

      real(kind=R8) function expo3(xi, n, alp, bet, gam, mu, si)
      !================================================================================================
      !<@note Profile function based on the exponential function
      !<
      !<@endnote
      !------------------------------------------------------------------------------------------------
      implicit none
      real   (kind=R8), intent(in) :: alp    !! *offset so that points are in [b1,b2]*
      real   (kind=R8), intent(in) :: bet    !! *reduction so that points are in [b1,b2]*
      real   (kind=R8), intent(in) :: gam    !! **
      real   (kind=R8), intent(in) :: mu     !! *numerical mean*
      real   (kind=R8), intent(in) :: si     !! *numerical standard deviation*
      real   (kind=R8), intent(in) :: xi     !! *abscissa*
      integer(kind=I4), intent(in) :: n      !! *statistical moment degree, n=3 for sk and n=4 for ku*

         real(kind=R8) :: tmp

         tmp  = (xi + alp)/bet

         expo3 = ( (sign(UN, tmp) * (UN - exp(-abs(tmp))) + gam * tmp**3 - mu)/si )**n

      return
      endfunction expo3

      real(kind=R8) function add_expo3(n, deb, fin, alp, bet, gam, mu, si)
      !================================================================================================
      !<@note Function that adds to the series mean the border integrals as explained in the modules
      !< presentation.
      !<
      !<@endnote
      !------------------------------------------------------------------------------------------------
      implicit none
      real   (kind=R8), intent(in) :: alp    !! *offset so that points are in [b1,b2]*
      real   (kind=R8), intent(in) :: bet    !! *reduction so that points are in [b1,b2]*
      real   (kind=R8), intent(in) :: gam    !! **
      real   (kind=R8), intent(in) :: mu     !! *numerical mean*
      real   (kind=R8), intent(in) :: si     !! *numerical standard deviation*
      integer(kind=I4), intent(in) :: n      !! *statistical moment degree, n=3 for sk and n=4 for ku*
      integer(kind=I4), intent(in) :: fin    !! *last integration point*
      integer(kind=I4), intent(in) :: deb    !! *first integration point*

         real(kind=R8) :: xdeb, xfin

         xdeb = deb
         xfin = fin
         add_expo3 = (UN/12)*( +9*(expo3(xdeb +0.0_R8, n, alp, bet, gam, mu, si)+expo3(xfin -0.0_R8, n, alp, bet, gam, mu, si)) &
                               +1*(expo3(xdeb +1.0_R8, n, alp, bet, gam, mu, si)+expo3(xfin -1.0_R8, n, alp, bet, gam, mu, si)) &
                               -4*(expo3(xdeb +0.5_R8, n, alp, bet, gam, mu, si)+expo3(xfin -0.5_R8, n, alp, bet, gam, mu, si)) )
      return
      endfunction add_expo3

   endsubroutine calculs_skku_exp3