MODULE lmdz_wake_vec_modulation
PUBLIC wake_vec_modulation
CONTAINS

SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon_loc, qb, d_qb, deltaqw, &
    d_deltaqw, sigmaw, d_sigmaw, alpha)
  ! ------------------------------------------------------
  ! Dtermination du coefficient alpha tel que les tendances
  ! corriges alpha*d_G, pour toutes les grandeurs G, correspondent
  ! a une humidite positive dans la zone (x) et dans la zone (w).
  ! ------------------------------------------------------
  IMPLICIT NONE

  ! Input
  INTEGER,                      INTENT(IN)               :: nl, nlon
  REAL, DIMENSION(nlon, nl),    INTENT(IN)               :: qb, d_qb
  REAL, DIMENSION(nlon, nl),    INTENT(IN)               :: deltaqw, d_deltaqw
  REAL, DIMENSION(nlon),        INTENT(IN)               :: sigmaw, d_sigmaw
  LOGICAL, DIMENSION(nlon),     INTENT(IN)               :: wk_adv
  ! Output
  REAL, DIMENSION(nlon),        INTENT(INOUT)            :: alpha
  ! Internal variables
  REAL zeta(nlon, nl)
  REAL alpha1(nlon)
  REAL x, a, b, c, discrim
  REAL epsilon_loc
  INTEGER i,k

  DO k = 1, nl
    DO i = 1, nlon
      IF (wk_adv(i)) THEN
        IF ((deltaqw(i,k)+d_deltaqw(i,k))>=0.) THEN
          zeta(i, k) = 0.
        ELSE
          zeta(i, k) = 1.
        END IF
      END IF
    END DO
    DO i = 1, nlon
      IF (wk_adv(i)) THEN
        x = qb(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qb(i, k) + &
          (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - d_sigmaw(i) * &
          (deltaqw(i,k)+d_deltaqw(i,k))
        a = -d_sigmaw(i)*d_deltaqw(i, k)
        b = d_qb(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - &
          deltaqw(i, k)*d_sigmaw(i)
        c = qb(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon_loc
        discrim = b*b - 4.*a*c
        ! print*, 'x, a, b, c, discrim', x, a, b, c, discrim
        IF (a+b>=0.) THEN !! Condition suffisante pour la positivite de ovap
          alpha1(i) = 1.
        ELSE
          IF (x>=0.) THEN
            alpha1(i) = 1.
          ELSE
            IF (a>0.) THEN
              alpha1(i) = 0.9*min( (2.*c)/(-b+sqrt(discrim)),  &
                                   (-b+sqrt(discrim))/(2.*a) )
            ELSE IF (a==0.) THEN
              alpha1(i) = 0.9*(-c/b)
            ELSE
              ! print*,'a,b,c discrim',a,b,c discrim
              alpha1(i) = 0.9*max( (2.*c)/(-b+sqrt(discrim)),  &
                                   (-b+sqrt(discrim))/(2.*a))
            END IF
          END IF
        END IF
        alpha(i) = min(alpha(i), alpha1(i))
      END IF
    END DO
  END DO

  RETURN
END SUBROUTINE wake_vec_modulation
END MODULE lmdz_wake_vec_modulation