source: trunk/LMDZ.MARS/libf/phymars/thermcell_dqup.F90 @ 621

      subroutine thermcell_dqup(ngrid,nlayer,ptimestep,fm0,entr0,  &
     &    masse0,q_therm,dq_therm)
      implicit none

! #include "iniprint.h"
!=======================================================================
!
!   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
!   Version modifiee pour prendre les downdrafts a la place de la 
!   subsidence compensatoire
!=======================================================================

#include "dimensions.h"
#include "dimphys.h"

! ============================ INPUTS ============================

      INTEGER, INTENT(IN) :: ngrid,nlayer
      REAL, INTENT(IN) :: ptimestep
      REAL, INTENT(IN) :: fm0(ngridmx,nlayermx+1)
      REAL, INTENT(IN) ::entr0(ngridmx,nlayermx)
      REAL, INTENT(IN) :: q_therm(ngridmx,nlayermx)
      REAL, INTENT(IN) :: masse0(ngridmx,nlayermx)

! ============================ OUTPUTS ===========================

      REAL, INTENT(OUT) :: dq_therm(ngridmx,nlayermx)  ! dq/dt -> derivative

! ============================ LOCAL =============================

      REAL q(ngridmx,nlayermx)
      REAL qa(ngridmx,nlayermx)
      INTEGER ig,k
      REAL detr(ngridmx,nlayermx),entr(ngridmx,nlayermx)
      REAL wqd(ngridmx,nlayermx+1)
      REAL zzm

! =========== Init ==============================================

      qa(:,:)=q_therm(:,:)
      q(:,:)=q_therm(:,:)
      entr(:,:)=entr0(:,:)
      detr(:,:)=0.

      do k=1,nlayermx
         do ig=1,ngridmx
            detr(ig,k)=fm0(ig,k)-fm0(ig,k+1)+entr(ig,k)
            if (detr(ig,k).lt.0.) then
               entr(ig,k)=entr(ig,k)-detr(ig,k)
               detr(ig,k)=0.
            endif
         enddo
      enddo

! =========== Updraft ============================================

      do ig=1,ngrid
         qa(ig,1)=q(ig,1)
      enddo

      do k=2,nlayermx
         do ig=1,ngridmx
            if ((fm0(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
     &         1.e-5*masse0(ig,k)) then
         qa(ig,k)=(fm0(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
     &     /(fm0(ig,k+1)+detr(ig,k))
            else
               qa(ig,k)=q(ig,k)
            endif
         enddo
      enddo

! =========== Subsidence ==========================================

      do k=2,nlayermx
         do ig=1,ngridmx

! Schema avec advection sur plus qu'une maille.
            zzm=masse0(ig,k)/ptimestep
            if (fm0(ig,k)>zzm) then
               wqd(ig,k)=(zzm*q(ig,k)+(fm0(ig,k)-zzm)*q(ig,k+1))
            else
               wqd(ig,k)=fm0(ig,k)*q(ig,k)
            endif
         enddo
      enddo
      do ig=1,ngrid
         wqd(ig,1)=0.
         wqd(ig,nlayermx+1)=0.
      enddo


! ====== dq ======================================================

      dq_therm(:,:)=0.

         do k=1,nlayermx-1
           q(:,k)=q(:,k)+         &
     & (detr(:,k)*qa(:,k)-entr(:,k)*q(:,k) &
     &       -wqd(:,k)+wqd(:,k+1))  &
     &               *ptimestep/masse0(:,k)
         enddo

         do k=1,nlayermx-1
          dq_therm(:,k)=(q(:,k)-q_therm(:,k))/ptimestep
         enddo

      return
      end