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