!
! $Id: calltherm.F90 1428 2010-09-13 08:43:37Z fairhead $
!
      subroutine calltherm_mars(ngrid,nlayer,dtime,nq,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,&
     &   zw2,fraca,zpopsk,ztla,heatFlux,heatFlux_down,&
     &     buoyancyOut,buoyancyEst,hfmax,wmax)

       USE thermcell, only : nsplit_thermals,r_aspect_thermals
       USE ioipsl_getincom
      implicit none

      INTEGER, INTENT(IN) :: ngrid,nlayer
      REAL dtime
      LOGICAL logexpr0, logexpr2(ngrid,nlayer), logexpr1(ngrid)
      REAL fact
      INTEGER nbptspb,nq

      REAL, INTENT(IN) :: zzlay(ngrid,nlayer)
      REAL, INTENT(IN) :: zzlev(ngrid,nlayer+1)

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

!FH Update Thermiques
      REAL d_t_ajs(ngrid,nlayer), d_q_ajs(ngrid,nlayer,nq)
      REAL d_u_ajs(ngrid,nlayer),d_v_ajs(ngrid,nlayer)
      REAL dq2_therm(ngrid,nlayer), dq2_the(ngrid,nlayer)
      real fm_therm(ngrid,nlayer+1)
      real entr_therm(ngrid,nlayer),detr_therm(ngrid,nlayer)

!********************************************************
!     declarations
      real zpopsk(ngrid,nlayer)
      real ztla(ngrid,nlayer)
      real wmax(ngrid)
      real hfmax(ngrid)
      integer lmax(ngrid)

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


! variables locales
      REAL d_t_the(ngrid,nlayer), d_q_the(ngrid,nlayer,nq)
      REAL d_u_the(ngrid,nlayer),d_v_the(ngrid,nlayer)
!
      integer isplit
      real zfm_therm(ngrid,nlayer+1),zdt
      real zentr_therm(ngrid,nlayer),zdetr_therm(ngrid,nlayer)
      real heatFlux(ngrid,nlayer)
      real heatFlux_down(ngrid,nlayer)
      real buoyancyOut(ngrid,nlayer)
      real buoyancyEst(ngrid,nlayer)
      real zheatFlux(ngrid,nlayer)
      real zheatFlux_down(ngrid,nlayer)
      real zbuoyancyOut(ngrid,nlayer)
      real zbuoyancyEst(ngrid,nlayer)

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

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

!  Modele du thermique
!  ===================

         nsplit_thermals=20
         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

! 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 (nq .ne. 0) then
            d_q_the(:,:,:)=0.
         endif

             CALL thermcell_main_mars(ngrid,nlayer,nq,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  &
     &      ,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 (nq .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(:,:)

         enddo ! isplit

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

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

      end