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

Last change on this file since 4484 was 4463, checked in by lguez, 20 months ago

Replace stop by call to abort_physiq

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