
! $Id: calltherm.F90 5101 2024-07-23 06:22:55Z abarral $

      subroutine calltherm(dtime  &
        ,pplay,paprs,pphi,weak_inversion  &
        ,u_seri_,v_seri_,t_seri_,q_seri_,t_env,q_env,zqsat,debut  &
        ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs  &
      ,fm_therm,entr_therm,detr_therm,zqasc,clwcon0,lmax,ratqscth,&
         ratqsdiff,zqsatth,ale_bl,alp_bl,lalim_conv,wght_th, &
         zmax0,f0,zw2,fraca,ztv,zpspsk,ztla,zthl &
!!! nrlmd le 10/04/2012
        ,pbl_tke,pctsrf,omega,airephy &
        ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
        ,n2,s2,strig,zcong,ale_bl_stat &
        ,therm_tke_max,env_tke_max &
        ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
        ,alp_bl_conv,alp_bl_stat &
!!! fin nrlmd le 10/04/2012
        ,zqla,ztva &
#ifdef ISO         
        ,xt_seri,d_xt_ajs &
#ifdef DIAGISO         
        ,q_the,xt_the &
#endif
#endif         
     )

      USE dimphy
      USE indice_sol_mod
      USE print_control_mod, ONLY: prt_level,lunout
      USE lmdz_thermcell_alp, ONLY: thermcell_alp
      USE lmdz_thermcell_main, ONLY: thermcell_main
      USE lmdz_thermcell_old, ONLY: thermcell, thermcell_2002, thermcell_eau, calcul_sec, thermcell_sec
#ifdef ISO
      use infotrac_phy, ONLY: ntiso
#ifdef ISOVERIF
      USE isotopes_mod, ONLY: iso_eau,iso_HDO
      USE isotopes_verif_mod, ONLY: iso_verif_aberrant_enc_vect2D, &
        iso_verif_egalite_vect2D
#endif   
#endif

      implicit none
      include "clesphys.h"
      include "thermcell_old.h"


!IM 140508
      INTEGER, SAVE ::  itap
!$OMP THREADPRIVATE(itap)
      REAL dtime
      LOGICAL debut
      LOGICAL logexpr0, logexpr2(klon,klev), logexpr1(klon)
      REAL fact(klon)
      INTEGER nbptspb

      REAL, DIMENSION(klon,klev), INTENT(IN) :: u_seri_,v_seri_
      REAL, DIMENSION(klon,klev), INTENT(IN) :: t_seri_,q_seri_
      REAL, DIMENSION(klon,klev), INTENT(IN) :: t_env,q_env
      REAL, DIMENSION(klon,klev) :: u_seri,v_seri
      REAL, DIMENSION(klon,klev) :: t_seri,q_seri
      REAL, DIMENSION(klon,klev) :: qmemoire
      REAL weak_inversion(klon)
      REAL paprs(klon,klev+1)
      REAL pplay(klon,klev)
      REAL pphi(klon,klev)
      real zlev(klon,klev+1) 
!test: on sort lentr et a* pour alimenter KE
      REAL wght_th(klon,klev)
      INTEGER lalim_conv(klon)
      REAL zw2(klon,klev+1),fraca(klon,klev+1)

!FH Update Thermiques
      REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
      REAL d_u_ajs(klon,klev),d_v_ajs(klon,klev)
      real fm_therm(klon,klev+1)
      real entr_therm(klon,klev),detr_therm(klon,klev)

!********************************************************
!     declarations
      LOGICAL flag_bidouille_stratocu
      real fmc_therm(klon,klev+1),zqasc(klon,klev)
      real zqla(klon,klev)
      real ztv(klon,klev),ztva(klon,klev)
      real zpspsk(klon,klev)
      real ztla(klon,klev)
      real zthl(klon,klev)
      real wmax_sec(klon)
      real zcong(klon)
      real zmax_sec(klon)
      real f_sec(klon)
      real detrc_therm(klon,klev)
! FH WARNING : il semble que ces save ne servent a rien
!     save fmc_therm, detrc_therm
      real clwcon0(klon,klev)
      real zqsat(klon,klev)
      real zw_sec(klon,klev+1)
      integer lmix_sec(klon)
      integer lmax(klon)
      real ratqscth(klon,klev)
      real ratqsdiff(klon,klev)
      real zqsatth(klon,klev)  
!nouvelles variables pour la convection
      real ale_bl(klon)
      real alp_bl(klon)
      real ale(klon)
      real alp(klon)
!RC
      !on garde le zmax du pas de temps precedent
      real zmax0(klon), f0(klon)

!!! nrlmd le 10/04/2012
      real pbl_tke(klon,klev+1,nbsrf)
      real pctsrf(klon,nbsrf)
      real omega(klon,klev)
      real airephy(klon)
      real zlcl_th(klon),fraca0(klon),w0(klon),w_conv(klon)
      real therm_tke_max0(klon),env_tke_max0(klon)
      real n2(klon),s2(klon),strig(klon)
      real ale_bl_stat(klon)
      real therm_tke_max(klon,klev),env_tke_max(klon,klev)
      real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
!!! fin nrlmd le 10/04/2012

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

      real, dimension(klon) :: pcon
      real, dimension(klon,klev) :: rhobarz,wth3
      integer,dimension(klon) :: lalim
      real, dimension(klon,klev+1) :: fm
      real, dimension(klon,klev) :: alim_star
      real, dimension(klon) :: zmax




! variables locales
      REAL d_t_the(klon,klev), d_q_the(klon,klev)
      REAL d_u_the(klon,klev),d_v_the(klon,klev)

      real zfm_therm(klon,klev+1),zdt
      real zentr_therm(klon,klev),zdetr_therm(klon,klev)
! FH A VERIFIER : SAVE INUTILES
!      save zentr_therm,zfm_therm

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

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

#ifdef ISO
      REAL xt_seri(ntiso,klon,klev),xtmemoire(ntiso,klon,klev)
      REAL d_xt_ajs(ntiso,klon,klev)
      real d_xt_the(ntiso,klon,klev)
#ifdef DIAGISO
      real q_the(klon,klev)
      real xt_the(ntiso,klon,klev)
#endif
      real qprec(klon,klev)
      integer ixt
#endif


!$OMP THREADPRIVATE(first)
!********************************************************
      if (first) then
        itap=0
        first=.false.
      endif

      u_seri(:,:)=u_seri_(:,:)
      v_seri(:,:)=v_seri_(:,:)
      t_seri(:,:)=t_seri_(:,:)
      q_seri(:,:)=q_seri_(:,:)

! Incrementer le compteur de la physique
     itap   = itap + 1

!  Modele du thermique
!  ===================
!         print*,'thermiques: WARNING on passe t au lieu de t_seri'


! On prend comme valeur initiale des thermiques la valeur du pas
! de temps precedent
         zfm_therm(:,:)=fm_therm(:,:)
         zdetr_therm(:,:)=detr_therm(:,:)
         zentr_therm(:,:)=entr_therm(:,:)

! On reinitialise les flux de masse a zero pour le cumul en
! cas de splitting
         fm_therm(:,:)=0.
         entr_therm(:,:)=0.
         detr_therm(:,:)=0.

         ale_bl(:)=0.
         alp_bl(:)=0.
         if (prt_level>=10) then
          print*,'thermV4 nsplit: ',nsplit_thermals,' weak_inversion'
         endif

!   tests sur les valeurs negatives de l'eau
         logexpr0=prt_level>=10
         nbptspb=0
         do k=1,klev
            do i=1,klon
! Attention teste abderr 19-03-09
!               logexpr2(i,k)=.not.q_seri(i,k).ge.0.
                logexpr2(i,k)=.not.q_seri(i,k)>=1.e-15
               if (logexpr2(i,k)) then
#ifdef ISO
                qprec(i,k)=q_seri(i,k)
#endif
                q_seri(i,k)=1.e-15
                nbptspb=nbptspb+1
#ifdef ISO
                do ixt=1,ntiso
                  xt_seri(ixt,i,k)=1.e-15*(xt_seri(ixt,i,k)/qprec(i,k))
                  ! xt_seri(ixt,i,k)=1.e-15*(Rdefault(index_iso(ixt)))
                enddo
#endif
               endif
!               if (logexpr0) &
!    &             print*,'WARN eau<0 avant therm i=',i,'  k=',k  &
!    &         ,' dq,q',d_q_the(i,k),q_seri(i,k)
            enddo
         enddo
         if(nbptspb>0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb


         new_thermcell=iflag_thermals>=15.and.iflag_thermals<=18
#ifdef ISO
      if (.not.new_thermcell) then
           CALL abort_gcm('calltherm 234','isos pas prevus ici',1)
      endif
#ifdef ISOVERIF
      if (iso_eau.gt.0) then
       CALL iso_verif_egalite_vect2D( &
             xt_seri,q_seri, &
             'calltherm 174',ntiso,klon,klev)
      endif !if (iso_eau.gt.0) then
#endif    
#endif
         zdt=dtime/REAL(nsplit_thermals)


         do isplit=1,nsplit_thermals

          if (iflag_thermals>=1000) then
            CALL thermcell_2002(klon,klev,zdt,iflag_thermals   &
        ,pplay,paprs,pphi  &
        ,u_seri,v_seri,t_seri,q_seri  &
        ,d_u_the,d_v_the,d_t_the,d_q_the  &
        ,zfm_therm,zentr_therm,fraca,zw2  &
        ,r_aspect_thermals,30.,w2di_thermals  &
        ,tau_thermals)
          else if (iflag_thermals==2) then
            CALL thermcell_sec(klon,klev,zdt  &
        ,pplay,paprs,pphi,zlev  &
        ,u_seri,v_seri,t_seri,q_seri  &
        ,d_u_the,d_v_the,d_t_the,d_q_the  &
        ,zfm_therm,zentr_therm  &
        ,r_aspect_thermals,30.,w2di_thermals  &
        ,tau_thermals)
          else if (iflag_thermals==3) then
            CALL thermcell(klon,klev,zdt  &
        ,pplay,paprs,pphi  &
        ,u_seri,v_seri,t_seri,q_seri  &
        ,d_u_the,d_v_the,d_t_the,d_q_the  &
        ,zfm_therm,zentr_therm  &
        ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
        ,tau_thermals)
          else if (iflag_thermals==10) then
            CALL thermcell_eau(klon,klev,zdt  &
        ,pplay,paprs,pphi  &
        ,u_seri,v_seri,t_seri,q_seri  &
        ,d_u_the,d_v_the,d_t_the,d_q_the  &
        ,zfm_therm,zentr_therm  &
        ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
        ,tau_thermals)
          else if (iflag_thermals==11) then
              abort_message = 'cas non prevu dans calltherm'
              CALL abort_physic (modname,abort_message,1)
          else if (iflag_thermals==12) then
            CALL calcul_sec(klon,klev,zdt  &
        ,pplay,paprs,pphi,zlev  &
        ,u_seri,v_seri,t_seri,q_seri  &
        ,zmax_sec,wmax_sec,zw_sec,lmix_sec  &
        ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
        ,tau_thermals)
          else if (iflag_thermals==13.or.iflag_thermals==14) then
              abort_message = 'thermcellV0_main enleve svn>2084'
              CALL abort_physic (modname,abort_message,1)
          else if (new_thermcell) then
            CALL thermcell_main(itap,klon,klev,zdt  &
        ,pplay,paprs,pphi,debut  &
        ,u_seri,v_seri,t_seri,q_seri,t_env,q_env  &
        ,d_u_the,d_v_the,d_t_the,d_q_the  &
        ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax  &
        ,ratqscth,ratqsdiff,zqsatth  &
        ,zmax0,f0,zw2,fraca,ztv,zpspsk &
        ,ztla,zthl,ztva &
        ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax,zcong &
#ifdef ISO         
        ,xt_seri,d_xt_the &
#endif         
     )

            CALL thermcell_alp(klon,klev,zdt  &                      ! in
          ,pplay,paprs  &                                        ! in
          ,zfm_therm,zentr_therm,lmax  &                         ! in
          ,pbl_tke,pctsrf,omega,airephy &                        ! in
          ,zw2,fraca &                                           ! in
          ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &  ! in
          ,zcong,ale,alp,lalim_conv,wght_th &                          ! out
          ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &! out
          ,n2,s2,strig,ale_bl_stat &                                   ! out
          ,therm_tke_max,env_tke_max &                           ! out
          ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &          ! out
          ,alp_bl_conv,alp_bl_stat &                             ! out
          )

           if (prt_level>10) write(lunout,*)'Apres thermcell_main OK'
         else
           abort_message = 'Cas des thermiques non prevu'
           CALL abort_physic (modname,abort_message,1)
         endif

! Attention : les noms sont contre intuitif.
! flag_bidouille_stratocu est .true. si on ne fait pas de bidouille.
! Il aurait mieux valu avoir un nobidouille_stratocu
! Et pour simplifier :
! nobidouille_stratocu=.not.(iflag_thermals==13.or.iflag_thermals=15)
! Ce serait bien de changer, mai en prenant le temps de vérifier que ca
! fait bien ce qu'on croit.

       flag_bidouille_stratocu=iflag_thermals<=12.or.iflag_thermals==14.or.iflag_thermals==16.or.iflag_thermals==18

! Calcul a posteriori du niveau max des thermiques pour les schémas qui
! ne la sortent pas.
      if (iflag_thermals<=12.or.iflag_thermals>=1000) then
         lmax(:)=1
         do k=1,klev-1
            zdetr_therm(:,k)=zentr_therm(:,k)+zfm_therm(:,k)-zfm_therm(:,k+1)
         enddo
         do k=1,klev-1
            do i=1,klon
               if (zfm_therm(i,k+1)>0.) lmax(i)=k
            enddo
         enddo
      endif

      fact(:)=0.
      DO i=1,klon
       logexpr1(i)=flag_bidouille_stratocu.or.weak_inversion(i)>0.5
       IF(logexpr1(i)) fact(i)=1./REAL(nsplit_thermals)
      ENDDO

     DO k=1,klev
!  transformation de la derivee en tendance
            d_t_the(:,k)=d_t_the(:,k)*dtime*fact(:)
            d_u_the(:,k)=d_u_the(:,k)*dtime*fact(:)
            d_v_the(:,k)=d_v_the(:,k)*dtime*fact(:)
            d_q_the(:,k)=d_q_the(:,k)*dtime*fact(:)
            fm_therm(:,k)=fm_therm(:,k)  &
        +zfm_therm(:,k)*fact(:)
            entr_therm(:,k)=entr_therm(:,k)  &
         +zentr_therm(:,k)*fact(:)
            detr_therm(:,k)=detr_therm(:,k)  &
         +zdetr_therm(:,k)*fact(:)
#ifdef ISO
            do ixt=1,ntiso
              d_xt_the(ixt,:,k)=d_xt_the(ixt,:,k)*dtime*fact(:)
            enddo
#endif
      ENDDO
       fm_therm(:,klev+1)=0.



!  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(:,:)
#ifdef ISO
            d_xt_ajs(:,:,:)=d_xt_ajs(:,:,:)+d_xt_the(:,:,:)
#endif

!  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(:,:)
            qmemoire(:,:)=q_seri(:,:)
            q_seri(:,:) = q_seri(:,:) + d_q_the(:,:)
#ifdef ISO
            xtmemoire(:,:,:)=xt_seri(:,:,:)
            xt_seri(:,:,:) = xt_seri(:,:,:) + d_xt_the(:,:,:)
#ifdef ISOVERIF
!      write(*,*) 'calltherm 350 tmp: ajout d_xt_the'
      if (iso_HDO.gt.0) then
!      i=479
!      k=4
!      write(*,*) 'xt_seri(iso_hdo,i,k),q_seri(i,k)=', &
!     &   xt_seri(iso_hdo,i,k),q_seri(i,k)
!      write(*,*) 'd_xt_the(iso_hdo,i,k),d_q_the(i,k)=', &
!     &   d_xt_the(iso_hdo,i,k),d_q_the(i,k)
      CALL iso_verif_aberrant_enc_vect2D( &
          xt_seri,q_seri, &
          'calltherm 353, apres ajout d_xt_the',ntiso,klon,klev)
      endif      
#endif
#endif 
           if (prt_level>10) write(lunout,*)'Apres apres thermcell_main OK'

       DO i=1,klon
            fm_therm(i,klev+1)=0.
            ale_bl(i)=ale_bl(i)+ale(i)/REAL(nsplit_thermals)
!            write(22,*)'ALE CALLTHERM',ale_bl(i),ale(i)
            alp_bl(i)=alp_bl(i)+alp(i)/REAL(nsplit_thermals)
!            write(23,*)'ALP CALLTHERM',alp_bl(i),alp(i)
        if(prt_level>=10) print*,'calltherm i alp_bl alp ale_bl ale',i,alp_bl(i),alp(i),ale_bl(i),ale(i)
       ENDDO

!IM 060508 marche pas comme cela !!!        enddo ! isplit

!   tests sur les valeurs negatives de l'eau
         nbptspb=0
            DO k = 1, klev
            DO i = 1, klon
               logexpr2(i,k)=.not.q_seri(i,k)>=0.
               if (logexpr2(i,k)) then
                q_seri(i,k)=1.e-15
                nbptspb=nbptspb+1
#ifdef ISO
                do ixt=1,ntiso
                  xt_seri(ixt,i,k)=1.e-15*(xtmemoire(ixt,i,k)/qmemoire(i,k))
                enddo
#endif
!                if (prt_level.ge.10) then
!                  print*,'WARN eau<0 apres therm i=',i,'  k=',k  &
!    &         ,' dq,q',d_q_the(i,k),q_seri(i,k),  &
!    &         'fm=',zfm_therm(i,k),'entr=',entr_therm(i,k)
                 endif
            ENDDO
            ENDDO
#ifdef ISO
#ifdef ISOVERIF
      if (iso_HDO.gt.0) then
      CALL iso_verif_aberrant_enc_vect2D( &
          xt_seri,q_seri, &
          'calltherm 393, apres bidouille q<0',ntiso,klon,klev)
      endif      
#endif
#endif

        IF(nbptspb>0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
! tests sur les valeurs de la temperature
        nbptspb=0
            DO k = 1, klev
            DO i = 1, klon
               logexpr2(i,k)=t_seri(i,k)<50..or.t_seri(i,k)>370.
               if (logexpr2(i,k)) nbptspb=nbptspb+1
!              if ((t_seri(i,k).lt.50.) .or.  &
!    &              (t_seri(i,k).gt.370.)) then
!                 print*,'WARN temp apres therm i=',i,'  k=',k  &
!    &         ,' t_seri',t_seri(i,k)
!              CALL abort
!              endif
            ENDDO
            ENDDO
        IF(nbptspb>0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
         enddo ! isplit

!***************************************************************
!     calcul du flux ascencant conservatif
!            print*,'<<<<calcul flux ascendant conservatif'

      fmc_therm=0.
               do k=1,klev
            do i=1,klon
                  if (entr_therm(i,k)>0.) then
                     fmc_therm(i,k+1)=fmc_therm(i,k)+entr_therm(i,k)
                  else
                     fmc_therm(i,k+1)=fmc_therm(i,k)
                  endif
                  detrc_therm(i,k)=(fmc_therm(i,k+1)-fm_therm(i,k+1))  &
                   -(fmc_therm(i,k)-fm_therm(i,k))
               enddo
            enddo
      
     
!****************************************************************
!     calcul de l'humidite dans l'ascendance
!      print*,'<<<<calcul de lhumidite dans thermique'
!CR:on ne le calcule que pour le cas sec
      if (iflag_thermals<=11) then
      do i=1,klon
         zqasc(i,1)=q_seri(i,1)
         do k=2,klev
            if (fmc_therm(i,k+1)>1.e-6) then
               zqasc(i,k)=(fmc_therm(i,k)*zqasc(i,k-1)  &
                +entr_therm(i,k)*q_seri(i,k))/fmc_therm(i,k+1)
!CR:test on asseche le thermique
!               zqasc(i,k)=zqasc(i,k)/2.
!            else
!               zqasc(i,k)=q_seri(i,k)
            endif
         enddo
       enddo
      

!     calcul de l'eau condensee dans l'ascendance
!             print*,'<<<<calcul de leau condensee dans thermique'
             do i=1,klon
                do k=1,klev
                   clwcon0(i,k)=zqasc(i,k)-zqsat(i,k)
                   if (clwcon0(i,k)<0. .or.   &
               (fm_therm(i,k+1)+detrc_therm(i,k))<1.e-6) then
                      clwcon0(i,k)=0.
                   endif
                enddo
             enddo
       else
              do i=1,klon
                do k=1,klev
                   clwcon0(i,k)=zqla(i,k)  
                   if (clwcon0(i,k)<0. .or.   &
               (fm_therm(i,k+1)+detrc_therm(i,k))<1.e-6) then
                   clwcon0(i,k)=0. 
                   endif
                enddo
             enddo
       endif
!*******************************************************************    


!jyg  Protection contre les temperatures nulles
          do i=1,klon
             do k=1,klev
                if (ztla(i,k) < 1.e-10) fraca(i,k) =0.
             enddo
          enddo


      return

      end
