MODULE lmdz_cv3_inip
  IMPLICIT NONE; PRIVATE
  PUBLIC cv3_inip
CONTAINS
  SUBROUTINE cv3_inip
    ! *******************************************************************
    ! *                                                                  *
    ! CV3_INIP Input = choice of mixing probability laws                 *
    !          Output = normalized coefficients of the probability laws. *
    ! *                                                                  *
    ! written by   : Jean-Yves Grandpeix, 06/06/2006, 19.39.27           *
    ! modified by :                                                      *
    ! *******************************************************************

    USE lmdz_print_control, ONLY: prt_level, lunout
    USE lmdz_abort_physic, ONLY: abort_physic
    ! IN: iflag_mix, gammas , alphas, betas, Fmax, scut, INOUT: qqa1 qqa2, OUT: Qcoef1max, Qcoef1max
    USE lmdz_yomcst2, ONLY: iflag_mix, gammas, alphas, betas, Fmax, scut, qqa1, qqa2, Qcoef1max, Qcoef2max

    IMPLICIT NONE

    !----------------------------------------------
    ! Local variables :
    CHARACTER (LEN = 20) :: modname = 'cv3_inip'
    CHARACTER (LEN = 80) :: abort_message

    REAL :: sumcoef
    REAL :: sigma, aire, pdf, mu, df
    REAL :: ff

    ! --   Mixing probability distribution functions

    REAL qcoef1, qcoef2, qff, qfff, qmix, rmix, qmix1, rmix1, qmix2, rmix2, f

    qcoef1(f) = tanh(f / gammas)
    qcoef2(f) = (tanh(f / gammas) + gammas * log(cosh((1. - f) / gammas) / cosh(f / gammas)))
    qff(f) = max(min(f, 1.), 0.)
    qfff(f) = min(qff(f), scut)
    qmix1(f) = (tanh((qff(f) - fmax) / gammas) + qcoef1max) / qcoef2max
    rmix1(f) = (gammas * log(cosh((qff(f) - fmax) / gammas)) + qff(f) * qcoef1max) / &
            qcoef2max
    qmix2(f) = -log(1. - qfff(f)) / scut
    rmix2(f) = (qfff(f) + (1. - qff(f)) * log(1. - qfff(f))) / scut
    qmix(f) = qqa1 * qmix1(f) + qqa2 * qmix2(f)
    rmix(f) = qqa1 * rmix1(f) + qqa2 * rmix2(f)

    IF (iflag_mix>0) THEN

      ! --      Normalize Pdf weights
      sumcoef = qqa1 + qqa2
      qqa1 = qqa1 / sumcoef
      qqa2 = qqa2 / sumcoef

      qcoef1max = qcoef1(fmax)
      qcoef2max = qcoef2(fmax)

      sigma = 0.
      aire = 0.0
      pdf = 0.0
      mu = 0.0
      df = 0.0001

      ! do ff = 0.0 + df, 1.0 - 2.*df, df
      ff = df
      DO WHILE (ff<=1.0 - 2. * df)
        pdf = (qmix(ff + df) - qmix(ff)) * (1. - ff) / df
        aire = aire + (qmix(ff + df) - qmix(ff)) * (1. - ff)
        mu = mu + pdf * ff * df
        IF (prt_level>9) WRITE (lunout, *) pdf, qmix(ff), aire, ff
        ff = ff + df
      END DO

      ! do ff=0.0+df,1.0 - 2.*df,df
      ff = df
      DO WHILE (ff<=1.0 - 2. * df)
        pdf = (qmix(ff + df) - qmix(ff)) * (1. - ff) / df
        sigma = sigma + pdf * (ff - mu) * (ff - mu) * df
        ff = ff + df
      END DO
      sigma = sqrt(sigma)

      IF (abs(aire - 1.0)>0.02) THEN
        WRITE (lunout, *) 'WARNING:: AREA OF MIXING PDF IS::', aire
        abort_message = ''
        CALL abort_physic(modname, abort_message, 1)
      ELSE
        PRINT *, 'Area, mean & std deviation are ::', aire, mu, sigma
      END IF
    END IF !  (iflag_mix .gt. 0)

  END SUBROUTINE cv3_inip
END MODULE lmdz_cv3_inip