source: trunk/LMDZ.MARS/libf/phymars/calltherm_mars.F90 @ 339

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

       USE ioipsl_getincom
      implicit none

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

      REAL ptimestep
      LOGICAL logexpr0, logexpr2(ngridmx,nlayermx), logexpr1(ngridmx)
      REAL fact
      INTEGER nbptspb,iq,l

      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 pplev(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)
      REAL masse(ngridmx,nlayermx)

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

!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
      character (len=80) :: abort_message

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

      REAL tstart,tstop


!  Modele du thermique
!  ===================
 
!         r_aspect_thermals     ! ultimately conrols the amount of mass going through the thermals
                                ! decreasing it increases the thermals effect. Tests at gcm resolution
                                ! shows that too low values destabilize the model
                                ! when changing this value, one should check that the surface layer model
                                ! outputs the correct Cd*u and Ch*u through changing the gustiness coefficient beta


#ifdef MESOSCALE
         !! valid for timesteps < 200s 
         nsplit_thermals=2
         r_aspect_thermals=0.7
#else
         nsplit_thermals=35
         r_aspect_thermals=1.5
#endif

         call getin("nsplit_thermals",nsplit_thermals)
         call getin("r_aspect_thermals",r_aspect_thermals)

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

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

         zw2(:,:)=0.
         zmaxth(:)=0.
         lmax_real(:)=0.

         zdt=ptimestep/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.
         zmax(:)=0.
         lmax(:)=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  &
!             CALL thermcell_main_mars_coupled_v2(zdt  &
     &      ,pplay,pplev,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(:,:)*ptimestep*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

             zmaxth(:)=zmaxth(:)+zmax(:)*fact
             lmax_real(:)=lmax_real(:)+float(lmax(:))*fact
            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

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

! Now that we have computed total entrainment and detrainment, we can
! advect u, v, and q in thermals. (theta already advected). We can do 
! that separatly because u,v,and q are not used in thermcell_main for
! any thermals-related computation : they are purely passive.

!calcul de la masse
      do l=1,nlayermx
         masse(:,l)=(pplev(:,l)-pplev(:,l+1))/g
      enddo

!calcul de l'epaisseur des couches
      do l=1,nlayermx
         zdz(:,l)=zzlev(:,l+1)-zzlev(:,l)
      enddo


      modname='momentum'
      call thermcell_dqup(ngridmx,nlayermx,ptimestep                &
     &      ,fm_therm,entr_therm,detr_therm,  &
     &     masse,u_seri,d_u_ajs,modname,zdz)

      call thermcell_dqup(ngridmx,nlayermx,ptimestep    &
     &       ,fm_therm,entr_therm,detr_therm,  &
     &     masse,v_seri,d_v_ajs,modname,zdz)

      if (nqmx .ne. 0.) then
      modname='tracer'
      DO iq=1,nqmx
      call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
     &     ,fm_therm,entr_therm,detr_therm,  &
     &    masse,pq_therm(:,:,iq),d_q_ajs(:,:,iq),modname,zdz)

      ENDDO
      endif

      DO i=1,ngridmx
         hfmax(i)=MAXVAL(heatFlux(i,:)+heatFlux_down(i,:))
         wmax(i)=MAXVAL(zw2(i,:))
      ENDDO

      lmax(:)=nint(lmax_real(:))
         
      return

      end