SUBROUTINE calcratqs_inter(klon,klev,pdtphys, &
           ratqsbas, wake_deltaq, wake_s, q_seri,     &
           ratqs_inter)
USE ioipsl_getin_p_mod, ONLY : getin_p
implicit none

!========================================================================
! Subroutine écrite par L. d'Alençon le 25/02/2021
! Cette subroutine calcule une valeur de ratqsbas interactive dépendant de la présence de poches froides dans l'environnement. 
! Elle est appelée par la subroutine calcratqs lorsque iflag_ratqs = 10. 
!========================================================================

! Declarations


LOGICAL, SAVE :: first = .TRUE.
!$OMP THREADPRIVATE(first)
! Input
integer,intent(in) :: klon,klev
real,intent(in) :: pdtphys,ratqsbas
real, dimension(klon,klev),intent(in) :: wake_deltaq, q_seri
real, dimension(klon),intent(in) :: wake_s

! Output
real, dimension(klon,klev),intent(inout) :: ratqs_inter

! local
integer i,k
real, dimension(klon,klev) :: wake_dq
REAL, SAVE             :: tau_ratqs_wake
!$OMP THREADPRIVATE(tau_ratqs_wake) 
REAL, SAVE             :: a_ratqs_wake
!$OMP THREADPRIVATE(a_ratqs_wake)
real, dimension(klon) :: max_wake_dq, max_dqconv
!-------------------------------------------------------------------------
!  Caclul de ratqs_inter
!-------------------------------------------------------------------------

!  
      if (first) then
         tau_ratqs_wake = 3600. ! temps de relaxation de la variabilité
         a_ratqs_wake = 1.85    ! paramètre pilotant l'importance du terme dépendant des poches froides
         CALL getin_p('tau_ratqs_wake', tau_ratqs_wake)
         CALL getin_p('a_ratqs_wake', a_ratqs_wake)         
         first=.false.
      endif       
      max_wake_dq(:) = 0.
      max_dqconv (:) = 0
      
      do k=1,klev
        do i=1,klon
         max_wake_dq(i) = max(abs(wake_deltaq(i,k)),max_wake_dq(i))
        enddo
      enddo
      
      do k=1,klev
        do i=1,klon
         ratqs_inter(i,k)= ratqs_inter(i,k)*exp(-pdtphys/tau_ratqs_wake) +   &
         a_ratqs_wake*(max_wake_dq(i)*(wake_s(i)**0.5/(1.-wake_s(i))))*(1.-exp(-pdtphys/tau_ratqs_wake))/q_seri(i,1)   
         if (ratqs_inter(i,k)<ratqsbas) then
            ratqs_inter(i,k) = ratqsbas
         endif       
        enddo
      enddo



return
end