Index: LMDZ6/trunk/libf/phylmd/thermcell_down.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/thermcell_down.F90 (revision 4381) +++ LMDZ6/trunk/libf/phylmd/thermcell_down.F90 (revision 4383) @@ -1,24 +1,30 @@ - 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,7 +32,6 @@ 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,10 +42,10 @@ 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,8 +56,8 @@ 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,12 +71,12 @@ 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,7 +86,7 @@ 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,13 +94,13 @@ 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,95 +108,95 @@ 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) !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^