!
! $Id: calltherm.F90 1428 2010-09-13 08:43:37Z fairhead $
!
      subroutine calltherm_mars(dtime,zzlev,zzlay  &
     &      ,pplay,paprs,pphi  &
     &      ,u_seri,v_seri,t_seri,pq_therm,q2_therm  &
     &      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs,dq2_therm  &
     &      ,fm_therm,entr_therm,detr_therm,lmax,zmax,&
     &   zw2,fraca,zpopsk,ztla,heatFlux,heatFlux_down,&
     &     buoyancyOut,buoyancyEst,hfmax,wmax)

       USE ioipsl_getincom
      implicit none

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

      REAL dtime
      LOGICAL logexpr0, logexpr2(ngridmx,nlayermx), logexpr1(ngridmx)
      REAL fact
      INTEGER nbptspb

      REAL, INTENT(IN) :: zzlay(ngridmx,nlayermx)
      REAL, INTENT(IN) :: zzlev(ngridmx,nlayermx+1)

      REAL u_seri(ngridmx,nlayermx),v_seri(ngridmx,nlayermx)
      REAL t_seri(ngridmx,nlayermx),pq_therm(ngridmx,nlayermx,nqmx)
      REAL q2_therm(ngridmx,nlayermx)
      REAL paprs(ngridmx,nlayermx+1)
      REAL pplay(ngridmx,nlayermx)
      REAL pphi(ngridmx,nlayermx)
      real zlev(ngridmx,nlayermx+1) 
!test: on sort lentr et a* pour alimenter KE
      REAL zw2(ngridmx,nlayermx+1),fraca(ngridmx,nlayermx+1)
      REAL zzw2(ngridmx,nlayermx+1)

!FH Update Thermiques
      REAL d_t_ajs(ngridmx,nlayermx), d_q_ajs(ngridmx,nlayermx,nqmx)
      REAL d_u_ajs(ngridmx,nlayermx),d_v_ajs(ngridmx,nlayermx)
      REAL dq2_therm(ngridmx,nlayermx), dq2_the(ngridmx,nlayermx)
      real fm_therm(ngridmx,nlayermx+1)
      real entr_therm(ngridmx,nlayermx),detr_therm(ngridmx,nlayermx)

!********************************************************
!     declarations
      real zpopsk(ngridmx,nlayermx)
      real ztla(ngridmx,nlayermx)
      real wmax(ngridmx)
      real hfmax(ngridmx)
      integer lmax(ngridmx)
      real zmax(ngridmx)

!nouvelles variables pour la convection
!RC
      !on garde le zmax du pas de temps precedent
!********************************************************


! variables locales
      REAL d_t_the(ngridmx,nlayermx), d_q_the(ngridmx,nlayermx,nqmx)
      REAL d_u_the(ngridmx,nlayermx),d_v_the(ngridmx,nlayermx)
!
      integer isplit,nsplit_thermals
      real r_aspect_thermals

      real zfm_therm(ngridmx,nlayermx+1),zdt
      real zentr_therm(ngridmx,nlayermx),zdetr_therm(ngridmx,nlayermx)
      real heatFlux(ngridmx,nlayermx)
      real heatFlux_down(ngridmx,nlayermx)
      real buoyancyOut(ngridmx,nlayermx)
      real buoyancyEst(ngridmx,nlayermx)
      real zheatFlux(ngridmx,nlayermx)
      real zheatFlux_down(ngridmx,nlayermx)
      real zbuoyancyOut(ngridmx,nlayermx)
      real zbuoyancyEst(ngridmx,nlayermx)

      character (len=20) :: modname='calltherm'
      character (len=80) :: abort_message

      integer i,k
      logical, save :: first=.true.

      REAL tstart,tstop


!  Modele du thermique
!  ===================
 
         r_aspect_thermals=2.
         nsplit_thermals=40
         call getin("nsplit_thermals",nsplit_thermals)

         fm_therm(:,:)=0.
         detr_therm(:,:)=0.
         entr_therm(:,:)=0.

         heatFlux(:,:)=0.
         heatFlux_down(:,:)=0.
         buoyancyOut(:,:)=0.
         buoyancyEst(:,:)=0.

         zw2(:,:)=0.

         zdt=dtime/REAL(nsplit_thermals)

         do isplit=1,nsplit_thermals

!         call cpu_time(tstart)


! On reinitialise les flux de masse a zero pour le cumul en
! cas de splitting

         zfm_therm(:,:)=0.
         zentr_therm(:,:)=0.
         zdetr_therm(:,:)=0.

         zheatFlux(:,:)=0.
         zheatFlux_down(:,:)=0.
         zbuoyancyOut(:,:)=0.
         zbuoyancyEst(:,:)=0.

         zzw2(:,:)=0.

         d_t_the(:,:)=0.
         d_u_the(:,:)=0.
         d_v_the(:,:)=0.
         dq2_the(:,:)=0.
         if (nqmx .ne. 0) then
            d_q_the(:,:,:)=0.
         endif

             CALL thermcell_main_mars(zdt  &
     &      ,pplay,paprs,pphi,zzlev,zzlay  &
     &      ,u_seri,v_seri,t_seri,pq_therm,q2_therm  &
     &      ,d_u_the,d_v_the,d_t_the,d_q_the,dq2_the  &
     &      ,zfm_therm,zentr_therm,zdetr_therm,lmax,zmax  &
     &      ,r_aspect_thermals &
     &      ,zzw2,fraca,zpopsk &
     &      ,ztla,zheatFlux,zheatFlux_down &
     &      ,zbuoyancyOut,zbuoyancyEst)

      fact=1./REAL(nsplit_thermals)
!  transformation de la derivee en tendance

            d_t_the(:,:)=d_t_the(:,:)*dtime*fact
            d_u_the(:,:)=d_u_the(:,:)*fact
            d_v_the(:,:)=d_v_the(:,:)*fact
            dq2_the(:,:)=dq2_the(:,:)*fact            

            if (nqmx .ne. 0) then
               d_q_the(:,:,:)=d_q_the(:,:,:)*fact
            endif

            fm_therm(:,:)=fm_therm(:,:)  &
     &      +zfm_therm(:,:)*fact
            entr_therm(:,:)=entr_therm(:,:)  &
     &       +zentr_therm(:,:)*fact
            detr_therm(:,:)=detr_therm(:,:)  &
     &       +zdetr_therm(:,:)*fact

            heatFlux(:,:)=heatFlux(:,:) &
     &       +zheatFlux(:,:)*fact
            heatFlux_down(:,:)=heatFlux_down(:,:) &
     &       +zheatFlux_down(:,:)*fact
            buoyancyOut(:,:)=buoyancyOut(:,:) &
     &       +zbuoyancyOut(:,:)*fact
            buoyancyEst(:,:)=buoyancyEst(:,:) &
     &       +zbuoyancyEst(:,:)*fact

            zw2(:,:)=zw2(:,:) + zzw2(:,:)*fact

!  accumulation de la tendance
      
            d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
            d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
            d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
            d_q_ajs(:,:,:)=d_q_ajs(:,:,:)+d_q_the(:,:,:)
            dq2_therm(:,:)=dq2_therm(:,:)+dq2_the(:,:)
!  incrementation des variables meteo
      
            t_seri(:,:) = t_seri(:,:) + d_t_the(:,:)
            u_seri(:,:) = u_seri(:,:) + d_u_the(:,:)
            v_seri(:,:) = v_seri(:,:) + d_v_the(:,:)
            pq_therm(:,:,:) = pq_therm(:,:,:) + d_q_the(:,:,:)
            q2_therm(:,:) = q2_therm(:,:) + dq2_therm(:,:)


!         call cpu_time(tstop)
!         print*,'elapsed time in thermals : ',tstop-tstart

         enddo ! isplit

     
!****************************************************************

!          do i=1,ngridmx
!             do k=1,nlayermx
!                if (ztla(i,k) .lt. 1.e-10) fraca(i,k) =0.
!               print*,'youpi je sers a quelque chose !'
!             enddo
!          enddo
        
          DO i=1,ngridmx
            hfmax(i)=MAXVAL(heatFlux(i,:)+heatFlux_down(i,:))
            wmax(i)=MAXVAL(zw2(i,:))
          ENDDO
 
      return

      end