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

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

!--------------------------------------------------------------
! thermcell_updown_dq: calcul du transport d'un traceur en pr??sence 
! d'up/down drafts
!
!--------------------------------------------------------------

   ! Suite du travail :
   ! Calculer les tendances d'un traceur (ici theta) en tenant compte
   ! des up et down drafts et de la subsidence compensatoire.


   IMPLICIT NONE

! arguments

   integer,intent(in) :: ngrid,nlay
   real,intent(in) :: ptimestep
   real,intent(in), dimension(ngrid,nlay) :: eup,dup,edn,ddn,masse
   real,intent(inout), dimension(ngrid,nlay) :: theta
   integer, intent(in), dimension(ngrid) :: lmax

   
! Local

   real, dimension(ngrid,nlay+1) :: fup,fdn,fthu,fthd,fthe,fthtot
   real, dimension(ngrid,nlay) :: thetau,thetad,dtheta
   real :: www

   integer ig,ilay

   fdn(:,:)=0.
   fup(:,:)=0.
   fthu(:,:)=0.
   fthd(:,:)=0.
   fthe(:,:)=0.
   fthtot(:,:)=0.
   thetad(:,:)=0.
   thetau(:,:)=0.

   ! lmax : indice tel que fu(kmax+1)=0
   
   ! Dans ce cas, pas besoin d'initialiser thetad(lmax) ( =theta(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
               thetad(ig,ilay)=( fdn(ig,ilay+1)*thetad(ig,ilay+1) + edn(ig,ilay)*theta(ig,ilay) ) / (fdn(ig,ilay)+ddn(ig,ilay))
            else
               www=fdn(ig,ilay+1)/ (fdn(ig,ilay)+ddn(ig,ilay))
               thetad(ig,ilay)=www*thetad(ig,ilay+1) + (1.-www)*theta(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
               thetau(ig,ilay)=theta(ig,ilay)
            else
               !thetau(ig,ilay)=( fup(ig,ilay)*thetau(ig,ilay-1) + eup(ig,ilay)*theta(ig,ilay) ) / (fup(ig,ilay+1)+dup(ig,ilay))
               !eup(ig,ilay)=fup(ig,ilay+1)-fup(ig,ilay)+dup(ig,ilay)
               !thetau(ig,ilay)=( fup(ig,ilay)*thetau(ig,ilay-1) + (fup(ig,ilay+1)-fup(ig,ilay)+dup(ig,ilay))*theta(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))
               thetau(ig,ilay)=www*thetau(ig,ilay-1)+(1.-www)*theta(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)*thetau(ig,ilay-1)
       fthd(ig,ilay)=-fdn(ig,ilay)*thetad(ig,ilay)
       !!!!ATTENTION HYPOTHESE de FLUX COMPENSATOIRE DESCENDANT ET DONC comme schema amont on va chercher theta 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))*theta(ig,ilay)
       fthtot(ig,ilay)=fthu(ig,ilay)+fthd(ig,ilay)+fthe(ig,ilay)
     enddo
   enddo
   !Boucle pour calculer theta
   do ilay=1,nlay,1
     do ig=1,ngrid
       dtheta(ig,ilay)=(fthtot(ig,ilay)-fthtot(ig,ilay+1))*(ptimestep/masse(ig,ilay))
!       theta(ig,ilay)=theta(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
       theta(ig,ilay)=theta(ig,ilay) + (fup(ig,ilay)*thetau(ig,ilay-1)-fup(ig,ilay+1)*thetau(ig,ilay) + &
       & (fup(ig,ilay+1)+fdn(ig,ilay+1))*theta(ig,ilay+1) - (fup(ig,ilay)+fdn(ig,ilay))*theta(ig,ilay) + &
       & fdn(ig,ilay+1)*thetad(ig,ilay+1)-fdn(ig,ilay)*thetad(ig,ilay))*(ptimestep/masse(ig,ilay))
     enddo
    enddo
   endif
! Il reste a coder : 
! d(rho theta)/dt = - d/dz(rho w'theta')
! d(theta)/dt = -1/rho * d/dz(rho w'theta')
! hydrostatique : dp - rho g dz
! dz = - dp /  (rho g)
! 1 / dz = - (rho g ) / dp

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


! ----> calculer d/dz(w'theta')
! w'theta' = alpha_up*(w_up-w_bar)*(theta_up-theta_bar) 
!          + alpha_down*(w_down-w_bar)*(theta_down-theta_bar)
!          + (1-alpha_up-alpha_down)*(w_env-w_bar)*(theta_env-theta_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 : theta_bar=alpha_up*theta_up + alpha_down*theta_dn + (1-alpha_up-alpha_down)theta_env 
!
! w'theta' = alpha_up*w_up*(theta_up-theta_bar) 
!          + alpha_down*w_down*(theta_down-theta_bar)
!          + (1-alpha_up-alpha_down)*w_env*(theta_env-theta_bar)

! rho*w'theta' = rho*alpha_up*w_up*(theta_up-theta_bar) 
!              + rho*alpha_down*w_down*(theta_down-theta_bar)
!              + rho*(1-alpha_up-alpha_down)*w_env*(theta_env-theta_bar)
! 

! rho*w'theta' = fup*(theta_up-theta_bar) 
!              + fdn*(theta_down-theta_bar)
!              - rho*(alpha_up*w_up+alpha_down*w_dn)*(theta_env-theta_bar)

!              - rho*(alpha_up*w_up+alpha_down*w_dn)*((theta_bar-alpha_up*theta_up-alpha_down*theta_down) /(1-alpha_up-alpha_down)- theta_bar)
!              - rho*(alpha_up*w_up+alpha_down*w_dn)*(theta_bar*(1/(1-alpha_up-alpha_down)-1) - (alpha_up*theta_up+alpha_down*theta_down) /(1-alpha_up-alpha_down))

! rho*w'theta' = fup*(theta_up-theta_bar) 
!              + fdn*(theta_down-theta_bar)
!              - (fup+fdn)*(theta_env-theta_bar)


! rho*w'theta' = fup*(theta_up-theta_bar) 
!              + fdn*(theta_down-theta_bar)
!              - (fup+fdn)*( (theta_bar-alpha_up*theta_up-alpha_down*theta_down) /(1-alpha_up-alpha_down) - theta_bar)

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

! rho*w'theta' = fup*(theta_up-theta_bar) 

!              + fdn*(theta_down-theta_bar)
!              + (fup+fdn)*(alpha_up*theta_up+alpha_down*theta_down)

! d(theta)/dt= -1/rho d/dz(rho*w'theta')
! choix de schema temporel (euler explicite)
! (theta(t,k)-theta(t-1,k))/dt = f(t-1,k)
! (theta(t,k)-theta(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'theta') = fup(k)*qa(k-1) -fup(k+1)*qa(k) 

!(rho*w'theta'(k+1)-rho*w'theta'(k))/dz = fup(k)thetau(k-1)-fup(k+1)thetaup(k) + (fup+fdn)(k+1)*theta_env(k+1) - (fup+fdn)(k)*theta_env(k) + fdn(k+1) thetadn(k+1)
!                   - fdn(k) thetadn(k)

! en continue on a  : d(theta)/dt = -1/rho * d/dz(rho w'theta')
! en discretis?? on a : 
! (theta(t,k)-theta(t-1,k))/dt = -1/rho * rho*w'theta'(k+1)-rho*w'theta'(k))/dz
! theta(t,k) = theta(t-1,k) - dt/rho * (rho*w'theta'(k+1)-rho*w'theta'(k))/dz 
! theta(t,k) = theta(t-1,k) - dt/rho * [fup(k)thetau(k-1)-fup(k+1)thetaup(k) + (fup+fdn)(k+1)*theta_env(k+1) -
! (fup+fdn)(k)*theta_env(k) + fdn(k+1) thetadn(k+1) - fdn(k) thetadn(k)]



! d(theta)/dt =  d/dp(rho w'theta') * g
! -1/rho * d/dz(rho w'theta') = g*d/dp(rho w'theta') = g * [fup(k)thetau(k-1)-fup(k+1)thetaup(k) + (fup+fdn)(k+1)*theta_env(k+1) -
! (fup+fdn)(k)*theta_env(k) + fdn(k+1) thetadn(k+1) - fdn(k) thetadn(k)]



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



! d/dz(rho*w'theta') = rho w'theta'(k+1)-rho w'theta'(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