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

Last change on this file since 4393 was 4383, checked in by evignon, 2 years ago

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

File size: 11.0 KB
RevLine 
[4383]1   SUBROUTINE thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,eup,dup,edn,ddn,masse,trac)
[4351]2
[4383]3!-----------------------------------------------------------------
4! thermcell_updown_dq: computes the tendency of tracers associated
5! with the presence of convective up/down drafts
6! This routine that has been collectively written during the
7! "ateliers downdrafts" in 2022/2023
8! Maelle, Frédéric, Catherine, Fleur, Florent, Etienne
9!------------------------------------------------------------------
[4351]10
11
[4365]12   IMPLICIT NONE
13
[4383]14! declarations
15!==============================================================
[4365]16
[4383]17! input/output
[4365]18
[4383]19   integer,intent(in)  :: ngrid ! number of horizontal grid points
20   integer, intent(in) :: nlay  ! number of vertical layers
21   real,intent(in) :: ptimestep ! time step of the physics [s]
22   real,intent(in), dimension(ngrid,nlay) :: eup ! entrainment to updrafts * dz [same unit as flux]
23   real,intent(in), dimension(ngrid,nlay) :: dup ! detrainment from updrafts * dz [same unit as flux]
24   real,intent(in), dimension(ngrid,nlay) :: edn ! entrainment to downdrafts * dz [same unit as flux]
25   real,intent(in), dimension(ngrid,nlay) :: ddn ! detrainment from downdrafts * dz [same unit as flux]
26   real,intent(in), dimension(ngrid,nlay) :: masse ! mass of layers = rho dz
27   real,intent(inout), dimension(ngrid,nlay) :: trac ! tracer
28   integer, intent(in), dimension(ngrid) :: lmax ! max level index at which downdraft are present
29
[4365]30   
31! Local
32
[4377]33   real, dimension(ngrid,nlay+1) :: fup,fdn,fthu,fthd,fthe,fthtot
[4383]34   real, dimension(ngrid,nlay) :: tracu,tracd,dtrac
[4377]35   real :: www
[4365]36   integer ig,ilay
37
38   fdn(:,:)=0.
[4377]39   fup(:,:)=0.
40   fthu(:,:)=0.
41   fthd(:,:)=0.
42   fthe(:,:)=0.
43   fthtot(:,:)=0.
[4383]44   tracd(:,:)=0.
45   tracu(:,:)=0.
[4365]46
47   ! lmax : indice tel que fu(kmax+1)=0
48   
[4383]49   ! Dans ce cas, pas besoin d'initialiser tracd(lmax) ( =trac(lmax) )
[4365]50
[4377]51   print*,'ON PASSE BIEN PAR LA CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC'
52   ! Boucle pour le downdraft
[4365]53   do ilay=nlay,1,-1
54      do ig=1,ngrid
[4377]55         if (ilay.le.lmax(ig) .and. lmax(ig)>1) then
[4365]56            fdn(ig,ilay)=fdn(ig,ilay+1)+edn(ig,ilay)-ddn(ig,ilay)
[4377]57            if ( 1 == 0 ) then
[4383]58               tracd(ig,ilay)=( fdn(ig,ilay+1)*tracd(ig,ilay+1) + edn(ig,ilay)*trac(ig,ilay) ) / (fdn(ig,ilay)+ddn(ig,ilay))
[4377]59            else
60               www=fdn(ig,ilay+1)/ (fdn(ig,ilay)+ddn(ig,ilay))
[4383]61               tracd(ig,ilay)=www*tracd(ig,ilay+1) + (1.-www)*trac(ig,ilay)
[4377]62            endif
[4365]63         endif
64      enddo
65   enddo
66
[4377]67   !Boucle pour l'updraft
68   do ilay=1,nlay,1
69      do ig=1,ngrid
70         if (ilay.le.lmax(ig) .and. lmax(ig)>1) then
71            fup(ig,ilay+1)=fup(ig,ilay)+eup(ig,ilay)-dup(ig,ilay)
72            if (ilay == 1 ) then
[4383]73               tracu(ig,ilay)=trac(ig,ilay)
[4377]74            else
[4383]75               !tracu(ig,ilay)=( fup(ig,ilay)*tracu(ig,ilay-1) + eup(ig,ilay)*trac(ig,ilay) ) / (fup(ig,ilay+1)+dup(ig,ilay))
[4377]76               !eup(ig,ilay)=fup(ig,ilay+1)-fup(ig,ilay)+dup(ig,ilay)
[4383]77               !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))
[4377]78               www=fup(ig,ilay)/(fup(ig,ilay+1)+dup(ig,ilay))
79               !1-www=(fup(ig,ilay+1)+dup(ig,ilay)-fup(ig,ilay))/(fup(ig,ilay+1)+dup(ig,ilay))
[4383]80               tracu(ig,ilay)=www*tracu(ig,ilay-1)+(1.-www)*trac(ig,ilay)
[4377]81            endif
82         endif
83      enddo
84   enddo
85   !Boucle pour calculer le flux up
86   do ilay=2,nlay,1
87     do ig=1,ngrid
[4383]88       fthu(ig,ilay)=fup(ig,ilay)*tracu(ig,ilay-1)
89       fthd(ig,ilay)=-fdn(ig,ilay)*tracd(ig,ilay)
90       !!!!ATTENTION HYPOTHESE de FLUX COMPENSATOIRE DESCENDANT ET DONC comme schema amont on va chercher trac au dessus!!!!!
[4381]91       !!!! si ce n'est pas le cas on stoppe le code !!!!
92       if (fup(ig,ilay)-fdn(ig,ilay) .lt. 0.) then
93          write(*,*) 'flux compensatoire montant, cas non traite par thermcell_updown_dq'
94          stop         
95       endif
[4383]96       fthe(ig,ilay)=-(fup(ig,ilay)-fdn(ig,ilay))*trac(ig,ilay)
[4377]97       fthtot(ig,ilay)=fthu(ig,ilay)+fthd(ig,ilay)+fthe(ig,ilay)
98     enddo
99   enddo
[4383]100   !Boucle pour calculer trac
[4377]101   do ilay=1,nlay,1
102     do ig=1,ngrid
[4383]103       dtrac(ig,ilay)=(fthtot(ig,ilay)-fthtot(ig,ilay+1))*(ptimestep/masse(ig,ilay))
104!       trac(ig,ilay)=trac(ig,ilay) + (fthtot(ig,ilay)-fthtot(ig,ilay+1))*(ptimestep/masse(ig,ilay))
[4377]105     enddo
106   enddo
107   if (1==0) then
108    do ilay=1,nlay,1
109     do ig=1,ngrid
[4383]110       trac(ig,ilay)=trac(ig,ilay) + (fup(ig,ilay)*tracu(ig,ilay-1)-fup(ig,ilay+1)*tracu(ig,ilay) + &
111       & (fup(ig,ilay+1)+fdn(ig,ilay+1))*trac(ig,ilay+1) - (fup(ig,ilay)+fdn(ig,ilay))*trac(ig,ilay) + &
112       & fdn(ig,ilay+1)*tracd(ig,ilay+1)-fdn(ig,ilay)*tracd(ig,ilay))*(ptimestep/masse(ig,ilay))
[4377]113     enddo
114    enddo
115   endif
116! Il reste a coder :
[4383]117! d(rho trac)/dt = - d/dz(rho w'trac')
118! d(trac)/dt = -1/rho * d/dz(rho w'trac')
[4377]119! hydrostatique : dp - rho g dz
120! dz = - dp /  (rho g)
121! 1 / dz = - (rho g ) / dp
[4365]122
[4383]123! d(trac)/dt = -1/rho * d/dp(rho w'trac') * (- rho g )
124! d(trac)/dt =  d/dp(rho w'trac') * g
[4365]125
[4377]126
[4383]127! ----> calculer d/dz(w'trac')
128! w'trac' = alpha_up*(w_up-w_bar)*(trac_up-trac_bar)
129!          + alpha_down*(w_down-w_bar)*(trac_down-trac_bar)
130!          + (1-alpha_up-alpha_down)*(w_env-w_bar)*(trac_env-trac_bar)
[4377]131! avec w_bar=0=alpha_up*w_up+alpha_down*w_dn+(1-alpha_up-alpha_down)w_env
132! -> w_env= - (alpha_up*w_up+alpha_down*w_dn)/(1-alpha_up-alpha_down)
[4383]133!  et on a : trac_bar=alpha_up*trac_up + alpha_down*trac_dn + (1-alpha_up-alpha_down)trac_env
[4377]134!
[4383]135! w'trac' = alpha_up*w_up*(trac_up-trac_bar)
136!          + alpha_down*w_down*(trac_down-trac_bar)
137!          + (1-alpha_up-alpha_down)*w_env*(trac_env-trac_bar)
[4377]138
[4383]139! rho*w'trac' = rho*alpha_up*w_up*(trac_up-trac_bar)
140!              + rho*alpha_down*w_down*(trac_down-trac_bar)
141!              + rho*(1-alpha_up-alpha_down)*w_env*(trac_env-trac_bar)
[4377]142!
143
[4383]144! rho*w'trac' = fup*(trac_up-trac_bar)
145!              + fdn*(trac_down-trac_bar)
146!              - rho*(alpha_up*w_up+alpha_down*w_dn)*(trac_env-trac_bar)
[4377]147
[4383]148!              - 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)
149!              - 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))
[4377]150
[4383]151! rho*w'trac' = fup*(trac_up-trac_bar)
152!              + fdn*(trac_down-trac_bar)
153!              - (fup+fdn)*(trac_env-trac_bar)
[4377]154
155
[4383]156! rho*w'trac' = fup*(trac_up-trac_bar)
157!              + fdn*(trac_down-trac_bar)
158!              - (fup+fdn)*( (trac_bar-alpha_up*trac_up-alpha_down*trac_down) /(1-alpha_up-alpha_down) - trac_bar)
[4377]159
160!!! hypoth??se : alpha_up+alpha_down << 1 -> 1/(1-alpha_up-alpha_down) ~ 1
161
[4383]162! rho*w'trac' = fup*(trac_up-trac_bar)
[4377]163
[4383]164!              + fdn*(trac_down-trac_bar)
165!              + (fup+fdn)*(alpha_up*trac_up+alpha_down*trac_down)
[4377]166
[4383]167! d(trac)/dt= -1/rho d/dz(rho*w'trac')
[4377]168! choix de schema temporel (euler explicite)
[4383]169! (trac(t,k)-trac(t-1,k))/dt = f(t-1,k)
170! (trac(t,k)-trac(t-1,k))/dt = f(t,k) ---> euler implicite
[4377]171! -----> on choisit explicite en temps
172
173!dans le cas d'une ascendance sans downdraft :
[4383]174! d/dz(rho*w'trac') = fup(k)*qa(k-1) -fup(k+1)*qa(k)
[4377]175
[4383]176!(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)
177!                   - fdn(k) tracdn(k)
[4377]178
[4383]179! en continue on a  : d(trac)/dt = -1/rho * d/dz(rho w'trac')
[4377]180! en discretis?? on a :
[4383]181! (trac(t,k)-trac(t-1,k))/dt = -1/rho * rho*w'trac'(k+1)-rho*w'trac'(k))/dz
182! trac(t,k) = trac(t-1,k) - dt/rho * (rho*w'trac'(k+1)-rho*w'trac'(k))/dz
183! 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) -
184! (fup+fdn)(k)*trac_env(k) + fdn(k+1) tracdn(k+1) - fdn(k) tracdn(k)]
[4377]185
186
187
[4383]188! d(trac)/dt =  d/dp(rho w'trac') * g
189! -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) -
190! (fup+fdn)(k)*trac_env(k) + fdn(k+1) tracdn(k+1) - fdn(k) tracdn(k)]
[4377]191
192
193
194! choix de sch??ma spatial
[4383]195! d/dz(rho*w'trac') = (rho*w'trac'(k+1)-rho*w'trac'(k))*dz
196! d/dz(rho*w'trac') = (rho*w'trac'(k)-rho*w'trac'(k-1))*dz
[4377]197
198
199
[4383]200! d/dz(rho*w'trac') = rho w'trac'(k+1)-rho w'trac'(k)
[4377]201
[4365]202!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
203! Initialisations :
204!------------------
205
206
207!
208 RETURN
209   END
210
211!=========================================================================
212
213   SUBROUTINE thermcell_down(ngrid,nlay,po,pt,pu,pv,pplay,pplev,  &
214     &           lmax,fup,eup,dup,theta)
215
216!--------------------------------------------------------------
[4377]217!thermcell_down: calcul des propri??t??s du panache descendant.
[4365]218!--------------------------------------------------------------
219
220
[4351]221   USE thermcell_ini_mod, ONLY : prt_level,RLvCp,RKAPPA,RETV
222   IMPLICIT NONE
223
224! arguments
225
226   integer,intent(in) :: ngrid,nlay
[4365]227   real,intent(in), dimension(ngrid,nlay) :: po,pt,pu,pv,pplay,eup,dup
228   real,intent(in), dimension(ngrid,nlay) :: theta
[4351]229   real,intent(in), dimension(ngrid,nlay+1) :: pplev,fup
230   integer, intent(in), dimension(ngrid) :: lmax
231
232
233   
234! Local
235
[4365]236   real, dimension(ngrid,nlay) :: edn,ddn,thetad
[4351]237   real, dimension(ngrid,nlay+1) :: fdn
238
239   integer ig,ilay
240   real dqsat_dT
241   logical mask(ngrid,nlay)
242
243   edn(:,:)=0.
244   ddn(:,:)=0.
245   fdn(:,:)=0.
[4365]246   thetad(:,:)=0.
[4351]247
[4365]248   ! lmax : indice tel que fu(kmax+1)=0
[4351]249   
[4365]250   ! Dans ce cas, pas besoin d'initialiser thetad(lmax) ( =theta(lmax) )
[4351]251
[4365]252! FH MODIFS APRES REUNIONS POUR COMMISSIONS
253! quelques erreurs de declaration
[4377]254! probleme si lmax=1 ce qui a l'air d'??tre le cas en d??but de simu. Devrait ??tre 0 ?
[4365]255! Remarques :
[4377]256! on pourrait ??crire la formule de thetad
[4365]257!    www=fdn(ig,ilay+1)/ (fdn(ig,ilay)+ddn(ig,ilay))
258!    thetad(ig,ilay)= www * thetad(ig,ilay+1) + (1.-www) * theta(ig,ilay)
[4377]259! Elle a l'avantage de bien montr?? la conservation, l'id??e fondamentale dans le
[4365]260!   transport qu'on ne fait que sommer des "sources" au travers d'un "propagateur"
261!   (Green)
[4377]262! Elle montre aussi beaucoup plus clairement pourquoi on n'a pas ?? se souccier (trop)
263!   de la possible nulit?? du d??nominateur
[4365]264
265
266   do ilay=nlay,1,-1
[4351]267      do ig=1,ngrid
[4365]268         if (ilay.le.lmax(ig).and.lmax(ig)>1) then
269            edn(ig,ilay)=0.5*dup(ig,ilay)
270            ddn(ig,ilay)=0.5*eup(ig,ilay)
271            fdn(ig,ilay)=fdn(ig,ilay+1)+edn(ig,ilay)-ddn(ig,ilay)
272            thetad(ig,ilay)=( fdn(ig,ilay+1)*thetad(ig,ilay+1) + edn(ig,ilay)*theta(ig,ilay) ) / (fdn(ig,ilay)+ddn(ig,ilay))
[4351]273         endif
274      enddo
275   enddo
276
277   ! Suite du travail :
278   ! Ecrire la conservervation de theta_l dans le panache descendant
279   ! Eventuellement faire la transformation theta_l -> theta_v
[4377]280   ! Si l'air est sec (et qu'on oublie le c??t?? theta_v) on peut
[4351]281   ! se contenter de conserver theta.
282   !
[4377]283   ! Connaissant thetadn, on peut calculer la flotabilit??.
284   ! Connaissant la flotabilit??, on peut calculer w de proche en proche
285   ! On peut calculer le detrainement de facon ?? garder alpha*rho = cste
286   ! On en d??duit l'entrainement lat??ral
287   ! C'est le mod??le des mini-projets.
[4351]288
289!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
290! Initialisations :
291!------------------
292
293
294!
295 RETURN
296   END
Note: See TracBrowser for help on using the repository browser.