MODULE lmdz_thermcell_dq
CONTAINS

  SUBROUTINE thermcell_dq(ngrid, nlay, impl, ptimestep, fm, entr, &
          masse, q, dq, qa, lev_out)
    USE lmdz_print_control, ONLY: prt_level
    USE lmdz_abort_physic, ONLY: abort_physic

    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

    ! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
    !  Introduction of an implicit computation of vertical advection in
    !  the environment of thermal plumes in thermcell_dq
    !  impl =     0 : explicit, 1 : implicit, -1 : old version

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

    ! arguments
    INTEGER, INTENT(IN) :: ngrid, nlay, impl
    REAL, INTENT(IN) :: ptimestep
    REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: masse
    REAL, INTENT(INOUT), DIMENSION(ngrid, nlay) :: entr, q
    REAL, INTENT(IN), DIMENSION(ngrid, nlay + 1) :: fm
    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: dq, qa
    INTEGER, INTENT(IN) :: lev_out                           ! niveau pour les print

    ! Local
    REAL, DIMENSION(ngrid, nlay) :: detr, qold
    REAL, DIMENSION(ngrid, nlay + 1) :: wqd, fqa
    REAL zzm
    INTEGER ig, k
    REAL cfl

    INTEGER niter, iter
    CHARACTER (LEN = 20) :: modname = 'thermcell_dq'
    CHARACTER (LEN = 80) :: abort_message


    ! Old explicite scheme
    IF (impl<=-1) THEN
      CALL thermcell_dq_o(ngrid, nlay, impl, ptimestep, fm, entr, &
              masse, q, dq, qa, lev_out)

    else


      ! Calcul du critere CFL pour l'advection dans la subsidence
      cfl = 0.
      DO k = 1, nlay
        DO ig = 1, ngrid
          zzm = masse(ig, k) / ptimestep
          cfl = max(cfl, fm(ig, k) / zzm)
          IF (entr(ig, k)>zzm) THEN
            PRINT*, 'entr*dt>m,1', k, entr(ig, k) * ptimestep, masse(ig, k)
            abort_message = 'entr dt > m, 1st'
            CALL abort_physic (modname, abort_message, 1)
          endif
        enddo
      enddo

      qold = q

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

      !   calcul du detrainement
      DO k = 1, nlay
        DO ig = 1, ngrid
          detr(ig, k) = fm(ig, k) - fm(ig, k + 1) + entr(ig, k)
          !           PRINT*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
          !test
          IF (detr(ig, k)<0.) THEN
            entr(ig, k) = entr(ig, k) - detr(ig, k)
            detr(ig, k) = 0.
            !               PRINT*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
            !     s         'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
          endif
          IF (fm(ig, k + 1)<0.) THEN
            !               PRINT*,'fm2<0!!!'
          endif
          IF (entr(ig, k)<0.) THEN
            !               PRINT*,'entr2<0!!!'
          endif
        enddo
      enddo

      ! Computation of tracer concentrations in the ascending plume
      DO ig = 1, ngrid
        qa(ig, 1) = q(ig, 1)
      enddo

      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

      ! Plume vertical flux
      DO k = 2, nlay - 1
        fqa(:, k) = fm(:, k) * qa(:, k - 1)
      enddo
      fqa(:, 1) = 0. ; fqa(:, nlay) = 0.


      ! Trace species evolution
      IF (impl==0) THEN
        DO k = 1, nlay - 1
          q(:, k) = q(:, k) + (fqa(:, k) - fqa(:, k + 1) - fm(:, k) * q(:, k) + fm(:, k + 1) * q(:, k + 1)) &
                  * ptimestep / masse(:, k)
        enddo
      else
        DO k = nlay - 1, 1, -1
          ! FH debut de modif : le calcul ci dessous modifiait numériquement
          ! la concentration quand le flux de masse etait nul car on divisait
          ! puis multipliait par masse/ptimestep.
          !        q(:,k)=(masse(:,k)*q(:,k)/ptimestep+fqa(:,k)-fqa(:,k+1)+fm(:,k+1)*q(:,k+1)) &
          !    &               /(fm(:,k)+masse(:,k)/ptimestep)
          q(:, k) = (q(:, k) + ptimestep / masse(:, k) * (fqa(:, k) - fqa(:, k + 1) + fm(:, k + 1) * q(:, k + 1))) &
                  / (1. + fm(:, k) * ptimestep / masse(:, k))
          ! FH fin de modif.
        enddo
      endif

      ! Tendencies
      DO k = 1, nlay
        DO ig = 1, ngrid
          dq(ig, k) = (q(ig, k) - qold(ig, k)) / ptimestep
          q(ig, k) = qold(ig, k)
        enddo
      enddo

    END IF ! impl=-1
    RETURN
  end


  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Obsolete version kept for convergence with Cmip5 NPv3.1 simulations
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  SUBROUTINE thermcell_dq_o(ngrid, nlay, impl, ptimestep, fm, entr, &
          masse, q, dq, qa, lev_out)
    USE lmdz_print_control, ONLY: prt_level
    USE lmdz_abort_physic, ONLY: abort_physic
    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 ngrid, nlay, impl

    REAL ptimestep
    REAL masse(ngrid, nlay), fm(ngrid, nlay + 1)
    REAL entr(ngrid, nlay)
    REAL q(ngrid, nlay)
    REAL dq(ngrid, nlay)
    INTEGER lev_out                           ! niveau pour les print

    REAL qa(ngrid, nlay), detr(ngrid, nlay), wqd(ngrid, nlay + 1)

    REAL zzm

    INTEGER ig, k
    REAL cfl

    REAL qold(ngrid, nlay)
    REAL ztimestep
    INTEGER niter, iter
    CHARACTER (LEN = 20) :: modname = 'thermcell_dq'
    CHARACTER (LEN = 80) :: abort_message



    ! Calcul du critere CFL pour l'advection dans la subsidence
    cfl = 0.
    DO k = 1, nlay
      DO ig = 1, ngrid
        zzm = masse(ig, k) / ptimestep
        cfl = max(cfl, fm(ig, k) / zzm)
        IF (entr(ig, k)>zzm) THEN
          PRINT*, 'entr*dt>m,2', k, entr(ig, k) * ptimestep, masse(ig, k)
          abort_message = 'entr dt > m, 2nd'
          CALL abort_physic (modname, abort_message, 1)
        endif
      enddo
    enddo

    !IM 090508     PRINT*,'CFL CFL CFL CFL ',cfl

#undef CFL
#ifdef CFL
! On subdivise le calcul en niter pas de temps.
      niter=int(cfl)+1
#else
    niter = 1
#endif

    ztimestep = ptimestep / niter
    qold = q

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

      !   calcul du detrainement
      DO k = 1, nlay
        DO ig = 1, ngrid
          detr(ig, k) = fm(ig, k) - fm(ig, k + 1) + entr(ig, k)
          !           PRINT*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
          !test
          IF (detr(ig, k)<0.) THEN
            entr(ig, k) = entr(ig, k) - detr(ig, k)
            detr(ig, k) = 0.
            !               PRINT*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
            !     s         'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
          endif
          IF (fm(ig, k + 1)<0.) THEN
            !               PRINT*,'fm2<0!!!'
          endif
          IF (entr(ig, k)<0.) THEN
            !               PRINT*,'entr2<0!!!'
          endif
        enddo
      enddo

      !   calcul de la valeur dans les ascendances
      DO ig = 1, ngrid
        qa(ig, 1) = q(ig, 1)
      enddo

      DO k = 2, nlay
        DO ig = 1, ngrid
          IF ((fm(ig, k + 1) + detr(ig, k)) * ztimestep>  &
                  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
#undef centre
#ifdef centre
             wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
#else

#define CFL_plus_grand_que_un
#ifdef CFL_plus_grand_que_un
          ! Schema avec advection sur plus qu'une maille.
          zzm = masse(ig, k) / ztimestep
          IF (fm(ig, k)>zzm) THEN
            wqd(ig, k) = zzm * q(ig, k) + (fm(ig, k) - zzm) * q(ig, k + 1)
          else
            wqd(ig, k) = fm(ig, k) * q(ig, k)
          endif
#else
            wqd(ig,k)=fm(ig,k)*q(ig,k)
#endif
#endif

          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))  &
                  * ztimestep / masse(ig, k)
          !            if (dq(ig,k).lt.0.) THEN
          !               PRINT*,'dq<0!!!'
          !            endif
        enddo
      enddo

    END DO


    ! Calcul des tendances
    DO k = 1, nlay
      DO ig = 1, ngrid
        dq(ig, k) = (q(ig, k) - qold(ig, k)) / ptimestep
        q(ig, k) = qold(ig, k)
      enddo
    enddo

    RETURN
  END
END MODULE lmdz_thermcell_dq