source: LMDZ6/trunk/libf/phylmd/thermcell_down.F90 @ 4570

SUBROUTINE thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,eup,dup,edn,ddn,masse,trac,dtrac)

USE thermcell_ini_mod, ONLY: iflag_thermals_down


!-----------------------------------------------------------------
! thermcell_updown_dq: computes the tendency of tracers associated
! with the presence of convective up/down drafts
! This routine that has been collectively written during the 
! "ateliers downdrafts" in 2022/2023
! Maelle, Frédéric, Catherine, Fleur, Florent, Etienne
!------------------------------------------------------------------


   IMPLICIT NONE

! declarations
!==============================================================

! input/output

   integer,intent(in)  :: ngrid ! number of horizontal grid points
   integer, intent(in) :: nlay  ! number of vertical layers
   real,intent(in) :: ptimestep ! time step of the physics [s]
   real,intent(in), dimension(ngrid,nlay) :: eup ! entrainment to updrafts * dz [same unit as flux]
   real,intent(in), dimension(ngrid,nlay) :: dup ! detrainment from updrafts * dz [same unit as flux]
   real,intent(in), dimension(ngrid,nlay) :: edn ! entrainment to downdrafts * dz [same unit as flux]
   real,intent(in), dimension(ngrid,nlay) :: ddn ! detrainment from downdrafts * dz [same unit as flux]
   real,intent(in), dimension(ngrid,nlay) :: masse ! mass of layers = rho dz 
   real,intent(in), dimension(ngrid,nlay) :: trac ! tracer 
   integer, intent(in), dimension(ngrid) :: lmax ! max level index at which downdraft are present
   real,intent(out),dimension(ngrid,nlay) ::dtrac ! tendance du traceur

   
! Local

   real, dimension(ngrid,nlay+1) :: fup,fdn,fc,fthu,fthd,fthe,fthtot
   real, dimension(ngrid,nlay) :: tracu,tracd,traci,tracold
   real :: www, mstar_inv
   integer ig,ilay
   real, dimension(ngrid,nlay):: s1,s2,num !coefficients pour la resolution implicite
   integer :: iflag_impl=1 ! 0 pour explicite, 1 pour implicite "classique", 2 pour implicite avec entrainement et detrainement
   
   fdn(:,:)=0.
   fup(:,:)=0.
   fc(:,:)=0.
   fthu(:,:)=0.
   fthd(:,:)=0.
   fthe(:,:)=0.
   fthtot(:,:)=0.
   tracd(:,:)=0.
   tracu(:,:)=0.
   traci(:,:)=trac(:,:)
   tracold(:,:)=trac(:,:)
   s1(:,:)=0.
   s2(:,:)=0.
   num(:,:)=1.

   if ( iflag_thermals_down < 10 ) then
      call abort_physic("thermcell_updown_dq", &
           'thermcell_down_dq = 0 or >= 10', 1)
   else
        iflag_impl=iflag_thermals_down-10
   endif
      

   ! lmax : indice tel que fu(kmax+1)=0
   ! Dans ce cas, pas besoin d'initialiser tracd(lmax) ( =trac(lmax) )
   ! Boucle pour le downdraft
   do ilay=nlay,1,-1
      do ig=1,ngrid
         !if ( lmax(ig) > nlay - 2 ) stop "les thermiques montent trop haut"
         if (ilay.le.lmax(ig) .and. lmax(ig)>1 ) then
            fdn(ig,ilay)=fdn(ig,ilay+1)+edn(ig,ilay)-ddn(ig,ilay)
            if ( fdn(ig,ilay)+ddn(ig,ilay) > 0. ) then
               www=fdn(ig,ilay+1)/ (fdn(ig,ilay)+ddn(ig,ilay))
            else
               www=0.
            endif
            tracd(ig,ilay)=www*tracd(ig,ilay+1) + (1.-www)*trac(ig,ilay)
         endif
      enddo 
   enddo !Fin boucle sur l'updraft
   fdn(:,1)=0.

   !Boucle pour l'updraft
   do ilay=1,nlay,1
      do ig=1,ngrid
         if (ilay.lt.lmax(ig) .and. lmax(ig)>1) then
            fup(ig,ilay+1)=fup(ig,ilay)+eup(ig,ilay)-dup(ig,ilay)
            if (fup(ig,ilay+1)+dup(ig,ilay) > 0.) then
               www=fup(ig,ilay)/(fup(ig,ilay+1)+dup(ig,ilay))
            else
               www=0.
            endif
            if (ilay == 1 ) then
               tracu(ig,ilay)=trac(ig,ilay)
            else
               tracu(ig,ilay)=www*tracu(ig,ilay-1)+(1.-www)*trac(ig,ilay)
            endif
         endif
      enddo 
      enddo !fin boucle sur le downdraft

   ! Calcul des flux des traceurs dans les updraft et les downdrfat 
   ! et du flux de masse compensateur
   ! en ilay=1 et nlay+1, fthu=0 et fthd=0
   fthu(:,1)=0.
   fthu(:,nlay+1)=0.
   fthd(:,1)=0.
   fthd(:,nlay+1)=0.
   fc(:,1)=0.
   fc(:,nlay+1)=0.
   do ilay=2,nlay,1 !boucle sur les interfaces
     do ig=1,ngrid
       fthu(ig,ilay)=fup(ig,ilay)*tracu(ig,ilay-1)
       fthd(ig,ilay)=-fdn(ig,ilay)*tracd(ig,ilay)
       fc(ig,ilay)=fup(ig,ilay)-fdn(ig,ilay)
     enddo
   enddo
   

   !Boucle pour calculer le flux du traceur flux updraft, flux downdraft, flux compensatoire
   !Methode explicite : 
   if(iflag_impl==0) then
     do ilay=2,nlay,1
       do ig=1,ngrid
         !!!!ATTENTION HYPOTHESE de FLUX COMPENSATOIRE DESCENDANT ET DONC comme schema amont on va chercher trac au dessus!!!!!
         !!!! tentative de prise en compte d'un flux compensatoire montant  !!!!
         if (fup(ig,ilay)-fdn(ig,ilay) .lt. 0.) then
            call abort_physic("thermcell_updown_dq", 'flux compensatoire '&
                 // 'montant, cas non traite par thermcell_updown_dq', 1)
            !fthe(ig,ilay)=(fup(ig,ilay)-fdn(ig,ilay))*trac(ig,ilay-1)
         else
            fthe(ig,ilay)=-(fup(ig,ilay)-fdn(ig,ilay))*trac(ig,ilay)
         endif
         !! si on voulait le prendre en compte on
         !fthe(ig,ilay)=-(fup(ig,ilay)-fdn(ig,ilay))*trac(ig,ilay-1)
         fthtot(ig,ilay)=fthu(ig,ilay)+fthd(ig,ilay)+fthe(ig,ilay)
       enddo
     enddo
     !Boucle pour calculer trac
     do ilay=1,nlay
       do ig=1,ngrid
         dtrac(ig,ilay)=(fthtot(ig,ilay)-fthtot(ig,ilay+1))/masse(ig,ilay)
!         trac(ig,ilay)=trac(ig,ilay) + (fthtot(ig,ilay)-fthtot(ig,ilay+1))*(ptimestep/masse(ig,ilay))
       enddo
     enddo !fin du calculer de la tendance du traceur avec la methode explicite

   !!! Reecriture du schéma explicite avec les notations du schéma implicite
   else if(iflag_impl==-1) then
     write(*,*) 'nouveau schéma explicite !!!'
     !!! Calcul de s1
     do ilay=1,nlay
       do ig=1,ngrid
         s1(ig,ilay)=fthu(ig,ilay)-fthu(ig,ilay+1)+fthd(ig,ilay)-fthd(ig,ilay+1)
         s2(ig,ilay)=s1(ig,ilay)+fthe(ig,ilay)-fthe(ig,ilay+1)
       enddo
     enddo

     do ilay=2,nlay,1
       do ig=1,ngrid
         if (fup(ig,ilay)-fdn(ig,ilay) .lt. 0.) then
            call abort_physic("thermcell_updown_dq", 'flux compensatoire ' &
                 // 'montant, cas non traite par thermcell_updown_dq', 1)
         else
            fthe(ig,ilay)=-(fup(ig,ilay)-fdn(ig,ilay))*trac(ig,ilay)
         endif
         fthtot(ig,ilay)=fthu(ig,ilay)+fthd(ig,ilay)+fthe(ig,ilay)
       enddo
     enddo
     !Boucle pour calculer trac
     do ilay=1,nlay
       do ig=1,ngrid
         ! dtrac(ig,ilay)=(fthtot(ig,ilay)-fthtot(ig,ilay+1))/masse(ig,ilay)
         dtrac(ig,ilay)=(s1(ig,ilay)+fthe(ig,ilay)-fthe(ig,ilay+1))/masse(ig,ilay)
!         trac(ig,ilay)=trac(ig,ilay) + (fthtot(ig,ilay)-fthtot(ig,ilay+1))*(ptimestep/masse(ig,ilay))
!         trac(ig,ilay)=trac(ig,ilay) + (s1(ig,ilay)+fthe(ig,ilay)-fthe(ig,ilay+1))*(ptimestep/masse(ig,ilay))
       enddo
     enddo !fin du calculer de la tendance du traceur avec la methode explicite

   else if (iflag_impl==1) then
     do ilay=1,nlay
       do ig=1,ngrid
         s1(ig,ilay)=fthu(ig,ilay)-fthu(ig,ilay+1)+fthd(ig,ilay)-fthd(ig,ilay+1)
       enddo
     enddo
     
     !Boucle pour calculer traci = trac((t+dt)
     do ilay=nlay-1,1,-1
       do ig=1,ngrid
         if((fup(ig,ilay)-fdn(ig,ilay)) .lt. 0) then
            write(*,*) 'flux compensatoire montant, cas non traite par thermcell_updown_dq dans le cas d une resolution implicite, ilay : ', ilay
            call abort_physic("thermcell_updown_dq", "", 1)
         else
           mstar_inv=ptimestep/masse(ig,ilay)
           traci(ig,ilay)=((traci(ig,ilay+1)*fc(ig,ilay+1)+s1(ig,ilay))*mstar_inv+tracold(ig,ilay))/(1.+fc(ig,ilay)*mstar_inv)
         endif
       enddo
     enddo
     do ilay=1,nlay
       do ig=1,ngrid
         dtrac(ig,ilay)=(traci(ig,ilay)-tracold(ig,ilay))/ptimestep
       enddo
     enddo

   else
      call abort_physic("thermcell_updown_dq", &
           'valeur de iflag_impl non prevue', 1)
   endif

 RETURN
   END

!=========================================================================

   SUBROUTINE thermcell_down(ngrid,nlay,po,pt,pu,pv,pplay,pplev,  &
     &           lmax,fup,eup,dup,theta)

!--------------------------------------------------------------
!thermcell_down: calcul des propri??t??s du panache descendant.
!--------------------------------------------------------------


   USE thermcell_ini_mod, ONLY : prt_level,RLvCp,RKAPPA,RETV,fact_thermals_down
   IMPLICIT NONE

! arguments

   integer,intent(in) :: ngrid,nlay
   real,intent(in), dimension(ngrid,nlay) :: po,pt,pu,pv,pplay,eup,dup
   real,intent(in), dimension(ngrid,nlay) :: theta
   real,intent(in), dimension(ngrid,nlay+1) :: pplev,fup
   integer, intent(in), dimension(ngrid) :: lmax


   
! Local

   real, dimension(ngrid,nlay) :: edn,ddn,thetad
   real, dimension(ngrid,nlay+1) :: fdn

   integer ig,ilay
   real dqsat_dT
   logical mask(ngrid,nlay)

   edn(:,:)=0.
   ddn(:,:)=0.
   fdn(:,:)=0.
   thetad(:,:)=0.

   ! lmax : indice tel que fu(kmax+1)=0
   
   ! Dans ce cas, pas besoin d'initialiser thetad(lmax) ( =theta(lmax) )

! FH MODIFS APRES REUNIONS POUR COMMISSIONS
! quelques erreurs de declaration
! probleme si lmax=1 ce qui a l'air d'??tre le cas en d??but de simu. Devrait ??tre 0 ?
! Remarques :
! on pourrait ??crire la formule de thetad
!    www=fdn(ig,ilay+1)/ (fdn(ig,ilay)+ddn(ig,ilay))
!    thetad(ig,ilay)= www * thetad(ig,ilay+1) + (1.-www) * theta(ig,ilay) 
! Elle a l'avantage de bien montr?? la conservation, l'id??e fondamentale dans le 
!   transport qu'on ne fait que sommer des "sources" au travers d'un "propagateur"
!   (Green)
! Elle montre aussi beaucoup plus clairement pourquoi on n'a pas ?? se souccier (trop)
!   de la possible nulit?? du d??nominateur


   do ilay=nlay,1,-1
      do ig=1,ngrid
         if (ilay.le.lmax(ig).and.lmax(ig)>1) then
            edn(ig,ilay)=fact_thermals_down*dup(ig,ilay)
            ddn(ig,ilay)=fact_thermals_down*eup(ig,ilay)
            fdn(ig,ilay)=fdn(ig,ilay+1)+edn(ig,ilay)-ddn(ig,ilay)
            thetad(ig,ilay)=( fdn(ig,ilay+1)*thetad(ig,ilay+1) + edn(ig,ilay)*theta(ig,ilay) ) / (fdn(ig,ilay)+ddn(ig,ilay))
         endif
      enddo 
   enddo

   ! Suite du travail :
   ! Ecrire la conservervation de theta_l dans le panache descendant
   ! Eventuellement faire la transformation theta_l -> theta_v
   ! Si l'air est sec (et qu'on oublie le c??t?? theta_v) on peut
   ! se contenter de conserver theta.
   !
   ! Connaissant thetadn, on peut calculer la flotabilit??.
   ! Connaissant la flotabilit??, on peut calculer w de proche en proche
   ! On peut calculer le detrainement de facon ?? garder alpha*rho = cste
   ! On en d??duit l'entrainement lat??ral
   ! C'est le mod??le des mini-projets.

!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
! Initialisations :
!------------------


!
 RETURN
   END