MODULE lmdz_thermcell_dtke
CONTAINS

      SUBROUTINE thermcell_dtke(ngrid,nlay,nsrf,ptimestep,fm0,entr0,  &
             rg,pplev,tke)
      USE lmdz_print_control, ONLY: prt_level
      IMPLICIT NONE

!=======================================================================

!   Calcul du transport verticale dans la couche limite en presence
!   de "thermiques" explicitement representes
!   calcul du dq/dt une fois qu'on connait les ascendances

!=======================================================================

      INTEGER, INTENT(IN) :: ngrid,nlay,nsrf
      REAL, INTENT(IN) :: ptimestep
      REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: entr0
      REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: fm0,pplev
      REAL, INTENT(IN) :: rg
      REAL, INTENT(INOUT) :: tke(ngrid,nlay+1,nsrf)

      REAL, DIMENSION(ngrid,nlay) :: masse0,detr0, masse,entr,detr
      REAL, DIMENSION(ngrid,nlay+1) :: fm,wqd,q,qa
      INTEGER lev_out                           ! niveau pour les print


      REAL :: zzm

      INTEGER ig,k
      INTEGER isrf


      lev_out=0

!PRINT*,'thermcell_dtke'

      IF (prt_level>=1) PRINT*,'Q2 THERMCEL_DQ 0'

!   calcul du detrainement
      do k=1,nlay
         detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
         masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
      enddo


! Decalage vertical des entrainements et detrainements.
      masse(:,1)=0.5*masse0(:,1)
      entr(:,1)=0.5*entr0(:,1)
      detr(:,1)=0.5*detr0(:,1)
      fm(:,1)=0.
      do k=1,nlay-1
         masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
         entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
         detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
         fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
      enddo
      fm(:,nlay+1)=0.



DO isrf=1,nsrf

   q(:,:)=tke(:,:,isrf)
!   calcul de la valeur dans les ascendances
      do ig=1,ngrid
         qa(ig,1)=q(ig,1)
      enddo


    IF (1==1) THEN
      do k=2,nlay
         do ig=1,ngrid
            IF ((fm(ig,k+1)+detr(ig,k))*ptimestep>  &
           1.e-5*masse(ig,k)) THEN
         qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
           /(fm(ig,k+1)+detr(ig,k))
            else
               qa(ig,k)=q(ig,k)
            endif
            IF (qa(ig,k)<0.) THEN
!               PRINT*,'qa<0!!!'
            endif
            IF (q(ig,k)<0.) THEN
!               PRINT*,'q<0!!!'
            endif
         enddo
      enddo

! Calcul du flux subsident

      do k=2,nlay
         do ig=1,ngrid
            wqd(ig,k)=fm(ig,k)*q(ig,k)
            IF (wqd(ig,k)<0.) THEN
!               PRINT*,'wqd<0!!!'
            endif
         enddo
      enddo
      do ig=1,ngrid
         wqd(ig,1)=0.
         wqd(ig,nlay+1)=0.
      enddo
     

! Calcul des tendances
      do k=1,nlay
         do ig=1,ngrid
            q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)  &
                 -wqd(ig,k)+wqd(ig,k+1))  &
                 *ptimestep/masse(ig,k)
         enddo
      enddo

 END IF

   tke(:,:,isrf)=q(:,:)

END DO

      RETURN
      END
END MODULE lmdz_thermcell_dtke