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

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

!-----------------------------------------------------------------
! 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(inout), dimension(ngrid,nlay) :: trac ! tracer 
   integer, intent(in), dimension(ngrid) :: lmax ! max level index at which downdraft are present

   
! Local

   real, dimension(ngrid,nlay+1) :: fup,fdn,fthu,fthd,fthe,fthtot
   real, dimension(ngrid,nlay) :: tracu,tracd,dtrac
   real :: www
   integer ig,ilay

   fdn(:,:)=0.
   fup(:,:)=0.
   fthu(:,:)=0.
   fthd(:,:)=0.
   fthe(:,:)=0.
   fthtot(:,:)=0.
   tracd(:,:)=0.
   tracu(:,:)=0.

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

   print*,'ON PASSE BIEN PAR LA CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC'
   ! Boucle pour le downdraft
   do ilay=nlay,1,-1
      do ig=1,ngrid
         if (ilay.le.lmax(ig) .and. lmax(ig)>1) then
            fdn(ig,ilay)=fdn(ig,ilay+1)+edn(ig,ilay)-ddn(ig,ilay)
            if ( 1 == 0 ) then
               tracd(ig,ilay)=( fdn(ig,ilay+1)*tracd(ig,ilay+1) + edn(ig,ilay)*trac(ig,ilay) ) / (fdn(ig,ilay)+ddn(ig,ilay))
            else
               www=fdn(ig,ilay+1)/ (fdn(ig,ilay)+ddn(ig,ilay))
               tracd(ig,ilay)=www*tracd(ig,ilay+1) + (1.-www)*trac(ig,ilay)
            endif
         endif
      enddo 
   enddo

   !Boucle pour l'updraft
   do ilay=1,nlay,1
      do ig=1,ngrid
         if (ilay.le.lmax(ig) .and. lmax(ig)>1) then
            fup(ig,ilay+1)=fup(ig,ilay)+eup(ig,ilay)-dup(ig,ilay)
            if (ilay == 1 ) then
               tracu(ig,ilay)=trac(ig,ilay)
            else
               !tracu(ig,ilay)=( fup(ig,ilay)*tracu(ig,ilay-1) + eup(ig,ilay)*trac(ig,ilay) ) / (fup(ig,ilay+1)+dup(ig,ilay))
               !eup(ig,ilay)=fup(ig,ilay+1)-fup(ig,ilay)+dup(ig,ilay)
               !tracu(ig,ilay)=( fup(ig,ilay)*tracu(ig,ilay-1) + (fup(ig,ilay+1)-fup(ig,ilay)+dup(ig,ilay))*trac(ig,ilay) ) / (fup(ig,ilay+1)+dup(ig,ilay))
               www=fup(ig,ilay)/(fup(ig,ilay+1)+dup(ig,ilay))
               !1-www=(fup(ig,ilay+1)+dup(ig,ilay)-fup(ig,ilay))/(fup(ig,ilay+1)+dup(ig,ilay))
               tracu(ig,ilay)=www*tracu(ig,ilay-1)+(1.-www)*trac(ig,ilay)
            endif
         endif
      enddo 
   enddo
   !Boucle pour calculer le flux up
   do ilay=2,nlay,1
     do ig=1,ngrid
       fthu(ig,ilay)=fup(ig,ilay)*tracu(ig,ilay-1)
       fthd(ig,ilay)=-fdn(ig,ilay)*tracd(ig,ilay)
       !!!!ATTENTION HYPOTHESE de FLUX COMPENSATOIRE DESCENDANT ET DONC comme schema amont on va chercher trac au dessus!!!!!
       !!!! si ce n'est pas le cas on stoppe le code !!!!
       if (fup(ig,ilay)-fdn(ig,ilay) .lt. 0.) then
          write(*,*) 'flux compensatoire montant, cas non traite par thermcell_updown_dq'
          stop          
       endif
       fthe(ig,ilay)=-(fup(ig,ilay)-fdn(ig,ilay))*trac(ig,ilay)
       fthtot(ig,ilay)=fthu(ig,ilay)+fthd(ig,ilay)+fthe(ig,ilay)
     enddo
   enddo
   !Boucle pour calculer trac
   do ilay=1,nlay,1
     do ig=1,ngrid
       dtrac(ig,ilay)=(fthtot(ig,ilay)-fthtot(ig,ilay+1))*(ptimestep/masse(ig,ilay))
!       trac(ig,ilay)=trac(ig,ilay) + (fthtot(ig,ilay)-fthtot(ig,ilay+1))*(ptimestep/masse(ig,ilay))
     enddo
   enddo
   if (1==0) then
    do ilay=1,nlay,1
     do ig=1,ngrid
       trac(ig,ilay)=trac(ig,ilay) + (fup(ig,ilay)*tracu(ig,ilay-1)-fup(ig,ilay+1)*tracu(ig,ilay) + &
       & (fup(ig,ilay+1)+fdn(ig,ilay+1))*trac(ig,ilay+1) - (fup(ig,ilay)+fdn(ig,ilay))*trac(ig,ilay) + &
       & fdn(ig,ilay+1)*tracd(ig,ilay+1)-fdn(ig,ilay)*tracd(ig,ilay))*(ptimestep/masse(ig,ilay))
     enddo
    enddo
   endif
! Il reste a coder : 
! d(rho trac)/dt = - d/dz(rho w'trac')
! d(trac)/dt = -1/rho * d/dz(rho w'trac')
! hydrostatique : dp - rho g dz
! dz = - dp /  (rho g)
! 1 / dz = - (rho g ) / dp

! d(trac)/dt = -1/rho * d/dp(rho w'trac') * (- rho g )
! d(trac)/dt =  d/dp(rho w'trac') * g 


! ----> calculer d/dz(w'trac')
! w'trac' = alpha_up*(w_up-w_bar)*(trac_up-trac_bar) 
!          + alpha_down*(w_down-w_bar)*(trac_down-trac_bar)
!          + (1-alpha_up-alpha_down)*(w_env-w_bar)*(trac_env-trac_bar)
! avec w_bar=0=alpha_up*w_up+alpha_down*w_dn+(1-alpha_up-alpha_down)w_env
! -> w_env= - (alpha_up*w_up+alpha_down*w_dn)/(1-alpha_up-alpha_down)
!  et on a : trac_bar=alpha_up*trac_up + alpha_down*trac_dn + (1-alpha_up-alpha_down)trac_env 
!
! w'trac' = alpha_up*w_up*(trac_up-trac_bar) 
!          + alpha_down*w_down*(trac_down-trac_bar)
!          + (1-alpha_up-alpha_down)*w_env*(trac_env-trac_bar)

! rho*w'trac' = rho*alpha_up*w_up*(trac_up-trac_bar) 
!              + rho*alpha_down*w_down*(trac_down-trac_bar)
!              + rho*(1-alpha_up-alpha_down)*w_env*(trac_env-trac_bar)
! 

! rho*w'trac' = fup*(trac_up-trac_bar) 
!              + fdn*(trac_down-trac_bar)
!              - rho*(alpha_up*w_up+alpha_down*w_dn)*(trac_env-trac_bar)

!              - rho*(alpha_up*w_up+alpha_down*w_dn)*((trac_bar-alpha_up*trac_up-alpha_down*trac_down) /(1-alpha_up-alpha_down)- trac_bar)
!              - rho*(alpha_up*w_up+alpha_down*w_dn)*(trac_bar*(1/(1-alpha_up-alpha_down)-1) - (alpha_up*trac_up+alpha_down*trac_down) /(1-alpha_up-alpha_down))

! rho*w'trac' = fup*(trac_up-trac_bar) 
!              + fdn*(trac_down-trac_bar)
!              - (fup+fdn)*(trac_env-trac_bar)


! rho*w'trac' = fup*(trac_up-trac_bar) 
!              + fdn*(trac_down-trac_bar)
!              - (fup+fdn)*( (trac_bar-alpha_up*trac_up-alpha_down*trac_down) /(1-alpha_up-alpha_down) - trac_bar)

!!! hypoth??se : alpha_up+alpha_down << 1 -> 1/(1-alpha_up-alpha_down) ~ 1

! rho*w'trac' = fup*(trac_up-trac_bar) 

!              + fdn*(trac_down-trac_bar)
!              + (fup+fdn)*(alpha_up*trac_up+alpha_down*trac_down)

! d(trac)/dt= -1/rho d/dz(rho*w'trac')
! choix de schema temporel (euler explicite)
! (trac(t,k)-trac(t-1,k))/dt = f(t-1,k)
! (trac(t,k)-trac(t-1,k))/dt = f(t,k) ---> euler implicite
! -----> on choisit explicite en temps

!dans le cas d'une ascendance sans downdraft : 
! d/dz(rho*w'trac') = fup(k)*qa(k-1) -fup(k+1)*qa(k) 

!(rho*w'trac'(k+1)-rho*w'trac'(k))/dz = fup(k)tracu(k-1)-fup(k+1)tracup(k) + (fup+fdn)(k+1)*trac_env(k+1) - (fup+fdn)(k)*trac_env(k) + fdn(k+1) tracdn(k+1)
!                   - fdn(k) tracdn(k)

! en continue on a  : d(trac)/dt = -1/rho * d/dz(rho w'trac')
! en discretis?? on a : 
! (trac(t,k)-trac(t-1,k))/dt = -1/rho * rho*w'trac'(k+1)-rho*w'trac'(k))/dz
! trac(t,k) = trac(t-1,k) - dt/rho * (rho*w'trac'(k+1)-rho*w'trac'(k))/dz 
! trac(t,k) = trac(t-1,k) - dt/rho * [fup(k)tracu(k-1)-fup(k+1)tracup(k) + (fup+fdn)(k+1)*trac_env(k+1) -
! (fup+fdn)(k)*trac_env(k) + fdn(k+1) tracdn(k+1) - fdn(k) tracdn(k)]



! d(trac)/dt =  d/dp(rho w'trac') * g
! -1/rho * d/dz(rho w'trac') = g*d/dp(rho w'trac') = g * [fup(k)tracu(k-1)-fup(k+1)tracup(k) + (fup+fdn)(k+1)*trac_env(k+1) -
! (fup+fdn)(k)*trac_env(k) + fdn(k+1) tracdn(k+1) - fdn(k) tracdn(k)]



! choix de sch??ma spatial 
! d/dz(rho*w'trac') = (rho*w'trac'(k+1)-rho*w'trac'(k))*dz
! d/dz(rho*w'trac') = (rho*w'trac'(k)-rho*w'trac'(k-1))*dz



! d/dz(rho*w'trac') = rho w'trac'(k+1)-rho w'trac'(k)

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


!
 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
   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)=0.5*dup(ig,ilay)
            ddn(ig,ilay)=0.5*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