Changeset 4383 for LMDZ6/trunk/libf

Timestamp:

Jan 13, 2023, 12:58:50 PM (2 years ago)

Author:

evignon

Message:

on renomme la variable theta de thermcell_updown_dq en trac et
on commente un peu l'entete du fichier

File:

• LMDZ6/trunk/libf/phylmd/thermcell_down.F90 (modified) (8 diffs)

LMDZ6/trunk/libf/phylmd/thermcell_down.F90

@@ -1 @@
-   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.
+   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
 
-! 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
+! 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
 
    
@@ -26 +32 @@
 
    real, dimension(ngrid,nlay+1) :: fup,fdn,fthu,fthd,fthe,fthtot
-   real, dimension(ngrid,nlay) :: thetau,thetad,dtheta
+   real, dimension(ngrid,nlay) :: tracu,tracd,dtrac
    real :: www
-
    integer ig,ilay
 
@@ -37 +42 @@
    fthe(:,:)=0.
    fthtot(:,:)=0.
-   thetad(:,:)=0.
-   thetau(:,:)=0.
+   tracd(:,:)=0.
+   tracu(:,:)=0.
 
    ! lmax : indice tel que fu(kmax+1)=0
    
-   ! Dans ce cas, pas besoin d'initialiser thetad(lmax) ( =theta(lmax) )
+   ! Dans ce cas, pas besoin d'initialiser tracd(lmax) ( =trac(lmax) )
 
    print*,'ON PASSE BIEN PAR LA CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC'
@@ -51 +56 @@
             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))
+               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))
-               thetad(ig,ilay)=www*thetad(ig,ilay+1) + (1.-www)*theta(ig,ilay)
+               tracd(ig,ilay)=www*tracd(ig,ilay+1) + (1.-www)*trac(ig,ilay)
             endif
          endif
@@ -66 +71 @@
             fup(ig,ilay+1)=fup(ig,ilay)+eup(ig,ilay)-dup(ig,ilay)
             if (ilay == 1 ) then
-               thetau(ig,ilay)=theta(ig,ilay)
+               tracu(ig,ilay)=trac(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))
+               !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)
-               !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))
+               !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))
-               thetau(ig,ilay)=www*thetau(ig,ilay-1)+(1.-www)*theta(ig,ilay)
+               tracu(ig,ilay)=www*tracu(ig,ilay-1)+(1.-www)*trac(ig,ilay)
             endif
          endif
@@ -81 +86 @@
    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!!!!!
+       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
@@ -89 +94 @@
           stop          
        endif
-       fthe(ig,ilay)=-(fup(ig,ilay)-fdn(ig,ilay))*theta(ig,ilay)
+       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 theta
+   !Boucle pour calculer trac
    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))
+       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
@@ -103 +108 @@
     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))
+       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 theta)/dt = - d/dz(rho w'theta')
-! d(theta)/dt = -1/rho * d/dz(rho w'theta')
+! 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(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)
+! 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 : theta_bar=alpha_up*theta_up + alpha_down*theta_dn + (1-alpha_up-alpha_down)theta_env 
+!  et on a : trac_bar=alpha_up*trac_up + alpha_down*trac_dn + (1-alpha_up-alpha_down)trac_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)
+! 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'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)
+! 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'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')
+! 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)
-! (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
+! (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'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')
+! 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 : 
-! (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)]
+! (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'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)
+! 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)
 
 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^